樊華羽，詹浩，程詩信，米百剛，姚會勤

(1.西北工業(yè)大學 航空學院，西安 710072) (2.清華大學 航天航空學院，北京 100084)

0 引 言

飛行器外形多目標優(yōu)化設(shè)計是一個跨學科多目標的優(yōu)化任務(wù)，包含了總體、氣動、結(jié)構(gòu)、隱身等多個領(lǐng)域和專業(yè)的關(guān)鍵技術(shù)。這些因素之間相互交叉干擾，進一步增加了現(xiàn)代飛行器多目標優(yōu)化設(shè)計的復雜性。隨著電子計算機技術(shù)的蓬勃發(fā)展和計算流體力學(CFD)、電磁計算精度的提高，數(shù)值優(yōu)化算法結(jié)合先進的CFD與電磁計算，為氣動隱身一體化優(yōu)化設(shè)計提供了優(yōu)秀的平臺，進而針對各類型作戰(zhàn)飛機的氣動隱身一體化研究成為國內(nèi)外關(guān)注的熱點[1-4]。

相比于常規(guī)布局的飛行器氣動隱身優(yōu)化設(shè)計，飛翼布局自身的結(jié)構(gòu)特性使得其在隱身性能上優(yōu)于常規(guī)布局飛行器。但也因此在缺失尾翼、垂尾等結(jié)構(gòu)的情況下，其氣動及結(jié)構(gòu)載荷需要通過機翼來補足。國內(nèi)外針對飛翼布局無人機做了大量研究。K.Hyoungjin等[5]針對一體化翼身融合體進行了詳細地氣動分析；王豪杰等[6]通過風洞實驗對比，結(jié)合CFD完成對某型無尾飛翼布局無人機的氣動力設(shè)計；在飛翼構(gòu)型氣動優(yōu)化方面，王榮等[7]采用高精度粘性氣動計算模型建立了飛翼布局無人機外形精度下的氣動隱身多目標優(yōu)化；張樂等[8]針對雙發(fā)布局下的飛翼無人機大鼓包機身，提出了一種減小翼型前緣半徑的機身前緣類“鷹嘴”形飛翼布局優(yōu)化構(gòu)型；陳曦等[9]開展了飛翼布局無人機在考慮隱身迎角下的氣動隱身綜合優(yōu)化。

在針對飛翼布局無人機的多目標優(yōu)化中，由于需要頻繁調(diào)用目標函數(shù)，優(yōu)化過程需要投入較大的時間成本和計算機資源。如何在平衡計算精度與計算效率的情況下建立一套高效的多目標優(yōu)化算法，是目前亟待解決的問題。本文針對典型飛翼布局，結(jié)合基于EHVI(Expected Hyper-Volume Improvement)加點準則的高精度多目標優(yōu)化算法，建立一種新的能夠在較少調(diào)用目標函數(shù)的情況下完成優(yōu)化任務(wù)的高效多目標粒子群優(yōu)化算法，并對某型飛翼布局無人機進行隱身氣動多目標優(yōu)化設(shè)計研究。

1 高效多目標粒子群算法

1.1 多目標優(yōu)化問題

在n維搜索空間S∈Rn內(nèi)，定義一個目標數(shù)量為m的多目標優(yōu)化向量fi(x),i=1,…,m，那么帶約束的最小化多目標優(yōu)化問題的通用數(shù)學模型可以表示為

minf(x)=[f1(x),f2(x),…,fm(x)]

(1)

s.t.gj(x)≤0,j=1,…,l

hk(x)=0,k=1,…,p

式中：x=[x1,x2,…,xn]，為設(shè)計變量向量；gj(x)為第j個不等式約束；hk(x)為第k個等式約束；l和p分別為不等式約束和等式約束的數(shù)量。

1.2 多目標粒子群優(yōu)化算法

粒子群優(yōu)化(Particle Swarm Optimization，簡稱PSO)算法最初是為了模擬鳥群的飛行覓食行為而提出的一種群智能算法。其初始化為一群隨機粒子，通過迭代尋找最優(yōu)解。在每一次迭代中，粒子通過跟蹤兩個極值來更新自己：第一個極值是粒子個體所找到的歷史最優(yōu)解，稱為個體最優(yōu)值Pb；另一個極值是整個種群目前找到的最優(yōu)解，稱為全局最優(yōu)值Gb。在找到上述兩個最優(yōu)值時，粒子根據(jù)式(2)～式(3)來更新自己的速度和位置。

(2)

(3)

式中：c1和c2分別為認知和社會加速系數(shù)，取值范圍為[0,4]，一般取c1=c2=2；r1和r2為兩個在[0,1]內(nèi)服從均勻分布的隨機變量。

經(jīng)過眾多研究人員的深入研究，現(xiàn)在粒子群算法已可以方便地用于處理多目標優(yōu)化問題。

CDMOPSO(Crowding Distance Multi-objective PSO)算法是C.R.Raquel等在C.A.C.Coello等提出的多目標粒子群算法(MOPSO)[10]的基礎(chǔ)上，使用擁擠度算子代替超網(wǎng)格來維護外部檔案而改進(如圖1所示)的一種多目標優(yōu)化算法[11]。擁擠度算子[12]能夠通過選擇引導粒子運動的最優(yōu)解以及維護外部檔案來增加種群的多樣性，避免算法過早陷入局部最優(yōu)。這種算法的優(yōu)點是結(jié)構(gòu)簡單，易于實現(xiàn)，缺陷是基于均勻分布的變異操作使其易于陷入局部最優(yōu)，解決多模態(tài)優(yōu)化問題的能力較差。氣動隱身多目標優(yōu)化設(shè)計是一種高維強非線性問題，因此本文選擇α-stable變異函數(shù)來對其進行改進，通過α-stable分布[13-14]產(chǎn)生的隨機數(shù)，對PSO種群中的個體進行變異。在變異的過程中，通過動態(tài)調(diào)整函數(shù)的穩(wěn)定系數(shù)α來調(diào)整不同優(yōu)化階段變異的范圍和幅度。穩(wěn)定系數(shù)α描述了分布的尾跡大小，決定了隨機數(shù)的范圍。α的變化范圍是[1,2]。在開始階段使用較小的α值，增強全局搜索能力，隨著優(yōu)化循環(huán)的進程，α值越來越大，全局搜索能力減弱，局部搜索能力增強，為尋找高精度解提供幫助。其具體變異操作過程以及優(yōu)化結(jié)果對比詳見文獻[15]，流程圖如圖2所示。本文將其命名為ASMOPSO(Alpha-stable Multi-objective PSO)。

圖1 擁擠距離計算示意圖Fig.1 Schematic diagram of congestion distance calculation

圖2 ASMOPSO流程圖Fig.2 Flow chart of ASMOPSO

1.3 基于Kriging模型的EHVI加點準則

為了增加局部地區(qū)的代理模型精度，減少對初始樣本空間精度的依賴，采用基于自適應(yīng)加點的動態(tài)Kriging代理模型對粒子群中的未觀測點進行近似評估。在動態(tài)更新代理模型樣本點的過程中，使用EHVI準則選擇粒子群多目標優(yōu)化過程中若干個最可靠的未觀測點進行真實函數(shù)評估，更新Pareto解集，同時提高Kriging代理模型局部區(qū)域精度。

EHVI是M.Emmerich等[16-17]在超體積理論[18]的基礎(chǔ)上結(jié)合單目標EI(Expected Improvement)加點準則[19]提出的一種處理高耗時優(yōu)化問題的新型多目標加點準則。

給定一個新的解y，假設(shè)它不被Pareto解集P中的任意個體所支配，那么解集P的超體積改善為

H(y,P)=H(P∪{y})-H(P)

(4)

進而就可以定義超體積改善函數(shù)為

(5)

那么基于Kriging代理模型的多目標響應(yīng)可視為互相獨立且服從高斯分布的隨機變量Yj(x)，有：

(6)

式中：m為目標數(shù)量。

在此基礎(chǔ)上，可以得到多目標期望改善函數(shù)：

(7)

式中：Vnd為由Pareto解集確定的非支配區(qū)域，即如圖3所示的細實線與坐標軸封閉的區(qū)域。

圖3 超體積改善示意圖Fig.3 Diagram of hypervolume improvement

從式(6)可以看出：如果要精確計算EHVI值，需要在非支配區(qū)域進行多維積分。由于非支配區(qū)域的不規(guī)則性，將超體積區(qū)域分割成多個矩形單元是必不可少的步驟。當Pareto解較多或者目標維度較高時，精確地識別這些矩形單元并對其進行EI積分是一件非常耗時的工作。

為此，Cheng Shixin等[20]提出了一種動態(tài)的EHVI值計算方式，通過對多個標準函數(shù)的對比測試結(jié)果分析，相比其基本型CDMOPSO優(yōu)化算法，結(jié)合動態(tài)超體積期望改善的優(yōu)化算法能在大幅度減少調(diào)用真實函數(shù)次數(shù)的情況下保持計算精度，極大地提高了優(yōu)化效率。本文所使用的加點準則就是基于此動態(tài)計算方法的EHVI加點準則。

1.4 約束的處理

為了克服罰函數(shù)系數(shù)設(shè)定難的問題，優(yōu)化算法中采用不可行度(Infeasibility Degree，簡稱IFD)法來處理優(yōu)化問題中的約束[21]，當一個粒子處在不可行域內(nèi)時，IFD值可以表示為粒子與可行域邊界的接近程度，則定義第i個粒子的不可行度值為

(8)

式中：γ為等式約束的容限區(qū)間，γ≥0。

如果一個粒子位于可行域內(nèi)，則它的IFD值為0，這種粒子稱為可行粒子?；诓豢尚卸群蚉areto優(yōu)勝比較的粒子選擇機制是：①如果一個粒子可行，而另一個粒子不可行，那么選擇可行粒子；②如果兩個粒子都不可行，則選擇IFD值較小的那個粒子；③如果兩個粒子都可行，執(zhí)行Pareto優(yōu)勝比較，選擇非支配粒子。

基于式(8)將每個粒子的近似不可行度改寫為

(9)

1.5 基于EHVI的多目標粒子群算法

區(qū)別于傳統(tǒng)的尋找當前Kriging代理模型下的最大EHVI值的子優(yōu)化過程，本文計算當前種群中所有個體的EHVI值并對它們進行降序排列，選取前幾個個體進行真實函數(shù)評估并加入樣本集，用于更新外部非支配解檔案，引導ASMOPSO中個體的移動。在這個過程中Kriging代理模型只是作為提供未觀測點的均值和方差的工具，用于計算EHVI值，代理模型的精度不是首要考慮的因素。

綜上，結(jié)合多目標粒子群算法、Kriging代理模型、EHVI加點形成的高效多目標優(yōu)化算法(EHVIMOPSO)的流程如圖4所示。

圖4 EHVIMOPSO流程圖Fig.4 Flow chart of EHVIMOPSO

2 氣動隱身多目標優(yōu)化設(shè)計

2.1 優(yōu)化模型的建立

使用飛翼無人機作為初始外形，如圖5所示，以翼面為研究對象，以阻力和RCS(Radar Cross Section)為優(yōu)化目標，約束為升力系數(shù)不小于原始翼型和俯仰力矩系數(shù)絕對值不減小。設(shè)計點氣動計算狀態(tài)為：Ma=0.8,α=2.0°,Re=2.78×107。

圖5 飛翼布局無人機幾何外形Fig.5 Geometry of flywing UAV

由此可得優(yōu)化的數(shù)學模型為

(10)

式中：X為設(shè)計變量向量；CD為阻力系數(shù)；ARCS為RCS均值；CL為升力系數(shù)；CM為力矩系數(shù)；上標*表示初始翼型的計算值。

2.2 計算網(wǎng)格及仿真計算方法

流體計算方面，流場求解基于三維非定常雷諾平均N-S方程，湍流模型采用k-ωSST模型，物面邊界采用無滑移邊界條件，遠場采用壓力遠場邊界條件。

本文使用半模計算，采用由ANSYS ICEM CFD軟件生成的空間六面體網(wǎng)格對流動區(qū)域進行離散，網(wǎng)格數(shù)量為165萬，計算網(wǎng)格如圖6所示。在優(yōu)化設(shè)計循環(huán)中，使用該軟件自帶的腳本功能實現(xiàn)網(wǎng)格的運動[21]。

電磁隱身特性求解以RCS為衡量標準，它描述了物體因被電磁波照射而向各個方向散射，被雷達捕捉的雷達回波強度及其電磁特性。在計算目標RCS的過程中，需要平衡計算效率與計算精度。由于不同的計算方法對于不同尺寸的目標具有不同的適應(yīng)性，在選擇計算方法時需要考慮目標的電尺寸大小。

圖6 計算網(wǎng)格Fig.6 Computational grid

綜合考慮，本文選用大面元物理光學法LEPO(Large Element PO)配合一致性幾何繞射理論 (UTD)計算邊緣繞射場。計算狀態(tài)為單站，極化方式為水平極化，雷達頻率選擇10 GHz。采用半模計算，計算范圍(θ)為0°～180°，步長為1°。RCS計算示意圖如圖7所示。

圖7 RCS計算示意Fig.7 Diagram of RCS calculation

2.3 參數(shù)化方法

在優(yōu)化設(shè)計中，采用空間變形能力較強的FFD(Free Form Deformation)自由曲面變形法[22]作為參數(shù)化方法。該方法最早由T.W.Sederberg和S.R.Parry提出，是一種針對三維可變形物體的有效建模工具。其基本方法是通過構(gòu)建并操縱空間三維框架與被操縱面的映射關(guān)系，通過改變空間三維框架來改變被操縱面的形狀。由于此方法能夠適用于非常復雜的外形控制且模擬精度很高，在飛行器三維參數(shù)化中有較好的應(yīng)用。

飛翼優(yōu)化模型中的FFD參數(shù)化[23]空間控制點為翼面上下各25個控制點(如圖8(a)所示)，在坐標系中，x、y、z三個方向上布置的控制點數(shù)為(5,5,2)，由于控制變量明顯偏多，考慮在外形優(yōu)化中將整個變形框外框固定不變，選擇上下翼面中心位置的9個點(如圖8(b)所示)，一共18個點，坐標系中表示在x、y、z方向上的點數(shù)為(3,3,2)，變形方向限定為僅在z方向上。為了不讓變形過于劇烈，變形范圍設(shè)定在±0.05之內(nèi)。最后控制點一共為3×3×2=18個。

(b) FFD參數(shù)化翼面控制點示意圖圖8 翼面FFD參數(shù)化控制點Fig.8 Parameterized control point of airfoil FFD

在FFD參數(shù)化后，通過曲面差值生成新的飛翼翼面。參數(shù)化結(jié)果如表1所示，可以看出：FFD參數(shù)化后的代理模型與原始外形的翼面力學特性幾乎完全吻合，隱身性能存在0.3%的誤差，符合計算中的誤差范圍，故認為該參數(shù)化方法能夠精確表達原始模型的氣動外形。

表1 參數(shù)化結(jié)果比較Table 1 Comparison of parametric results

2.4 優(yōu)化過程及優(yōu)化結(jié)果分析

使用最優(yōu)拉丁超立方抽樣(Optimal Latin Hypercube Sampling，簡稱OLHS)[24]生成40個初始樣本點，建立初始代理模型并通過EHVIMOPSO多目標優(yōu)化算法搜索Pareto前沿。其中粒子群種群規(guī)模為200，外部檔案大小為50；慣性系數(shù)從0.75逐漸減小到0.25；認知加速系數(shù)和社會加速系數(shù)c1=c2=2.0；穩(wěn)定系數(shù)α的變化范圍為1.0～2.0。迭代步數(shù)為60步，每代使用EHVI值選擇3個最穩(wěn)定的個體加入到樣本空間中，并更新代理模型。優(yōu)化結(jié)果如圖9和表2所示。

圖9 Pareto前沿Fig.9 Forward position of Pareto 表2 Pareto前沿數(shù)值統(tǒng)計Table 2 Frontier point of Pareto

序號CDARCS/m210.008 5151.065 95120.008 5471.053 66030.008 5221.061 07740.008 5251.059 75450.008 5361.053 8496(A點)0.008 5321.055 57170.008 5271.057 25880.008 5301.056 98090.008 5421.053 797100.008 5451.053 757

從圖9和表2可以看出：Pareto解分布較為均勻，但是范圍還不夠?qū)拸V。

優(yōu)化前后數(shù)據(jù)比較如表3所示。由于此飛翼布局無人機為低可探測低阻外形，其原始外形的氣動及隱身性能非常優(yōu)秀(如表1所示)，因此在此優(yōu)化算法下的氣動優(yōu)化效果沒有特別明顯的提升。

表3 優(yōu)化前后數(shù)據(jù)比較Table 3 Data comparison before and after optimization

在平衡減阻與降低RCS兩者情況下選擇如圖9所示的ParetoA與原始構(gòu)型進行比較分析。

從表2可以看出：在ParetoA的構(gòu)型下阻力減少了0.34%，RCS縮減了4.41%。

在氣動計算結(jié)果對比方面，為了更好地說明，選取Y分別為1.0、2.5、4.5三個展向站位的翼型截面進行對比分析，如圖10所示。優(yōu)化后的翼型剖面形狀、壓力系數(shù)與原始外形對比結(jié)果如圖11～圖13所示。

圖10 翼面截取示意Fig.10 Sketch of wing interception

(a) 翼型對比

(b) 壓力系數(shù)對比圖11 對比結(jié)果(Y=1.0)Fig.11 Contrast result(Y=1.0)

(a) 翼型對比

(b) 壓力系數(shù)對比圖12 對比結(jié)果(Y=2.5)Fig.12 Contrast result(Y=2.5)

(a) 翼型對比

(b) 壓力系數(shù)對比圖13 對比結(jié)果(Y=4.5)Fig.13 Contrast result(Y=4.5)

從圖11～圖13可以看出：優(yōu)化構(gòu)型三個站位上的翼型剖面與原始翼型相比，在上翼面前部與下翼面后部都有小量縮減，翼面最大厚度向后移動且最大厚度值都減小、變薄，其中Y=2.5、Y=4.5兩個站位翼型變薄的情況非常明顯；翼型的前后緣位置壓力分布較陡峭，中段壓力系數(shù)曲線分布較為平緩，優(yōu)化后構(gòu)型仍能在此基礎(chǔ)上有一定程度的減緩翼面激波強度；Y=1.0站位，后緣附近有一個向下的加載力區(qū)域，使的整個優(yōu)化構(gòu)型的低頭力矩減小，提高了無人機的可操控性。

RCS優(yōu)化結(jié)果對比如圖14所示，可以看出：在入射角為57°和140°處有峰值存在，57°時的峰值為機翼前緣的電磁散射相干疊加而形成的，140°處的峰值為機翼端面造成的鏡面反射；RCS縮減部位主要集中在60°～90°范圍內(nèi)，這是由于本文選擇的翼面控制點在翼面中部位置，對于前緣和端面形成的電磁散射無法做到有效縮減。

圖14 優(yōu)化結(jié)果A與原始外形的RCS計算結(jié)果對比Fig.14 Comparison of RCS calculation results between Pareto A and original profile

3 結(jié) 論

本文使用計算流體力學和計算電磁學等數(shù)值模擬手段計算飛翼布局無人機的氣動和隱身性能，結(jié)合一種基于EHVI加點的高效多目標粒子群優(yōu)化算法使用FFD法進行參數(shù)化表達，對一種低阻、低可探測飛翼布局無人機進行氣動隱身多目標優(yōu)化設(shè)計。在200次調(diào)用目標函數(shù)的情況下，降低了無人機的阻力和雷達反射截面積，表明本文提出的高效優(yōu)化設(shè)計方法在解決類似飛翼式布局無人機氣動隱身多目標優(yōu)化設(shè)計等昂貴優(yōu)化問題時具有較大的應(yīng)用潛力。