宋 倩 萬志強 / SONG Qian WAN Zhiqiang

(北京航空航天大學(xué)，北京 100191)

0 引言

在現(xiàn)代客機的機翼設(shè)計中，需要對氣動、結(jié)構(gòu)和內(nèi)部容積等要求進(jìn)行綜合考慮，以氣動和結(jié)構(gòu)的耦合關(guān)系為基礎(chǔ)，找到一個最佳的設(shè)計方案。對于飛翼這種采用新型氣動布局的設(shè)計方案，相較傳統(tǒng)方案，則更需注重考慮氣動與結(jié)構(gòu)的耦合關(guān)系。

目前，對飛翼布局飛機的單項學(xué)科研究較多，研究的重點領(lǐng)域包含氣動、隱身、結(jié)構(gòu)、操穩(wěn)等學(xué)科；而對于飛翼布局飛機多學(xué)科綜合優(yōu)化研究的相對較少，且在對其進(jìn)行多學(xué)科優(yōu)化的過程中往往對于氣動與結(jié)構(gòu)的耦合問題做出了簡化，譬如只將氣動載荷單向傳遞給結(jié)構(gòu)模型，這樣獲得的結(jié)果往往并非準(zhǔn)確[1]。針對目前對飛翼客機機翼優(yōu)化過程中對氣動彈性效應(yīng)考慮不充分、氣動/結(jié)構(gòu)耦合不充分的問題，本文開展了針對此類機翼的氣動/結(jié)構(gòu)綜合優(yōu)化設(shè)計研究，同時獲得機翼最優(yōu)的氣動外形和結(jié)構(gòu)構(gòu)型。

1 理論基礎(chǔ)

1.1 靜氣動彈性響應(yīng)方程

靜氣動彈性響應(yīng)分析方程在f-set的位移向量集下一般可以寫為[2]：

(1)

1.2 氣動阻力計算方法

氣動外形優(yōu)化設(shè)計需要將氣動阻力約束考慮在內(nèi)，但是基于小擾動線化勢流理論的面元法無法計算誘導(dǎo)阻力、摩擦阻力以及激波這幾類阻力項。若采用基于N-S方程的CFD(Computational Fluid Dynamics，簡稱CFD)方法雖可以較為精確地得到高精度氣動力，但是計算量會劇增，時間耗費巨大，無法滿足效率的要求。本文考慮選擇基于Euler方程的無黏流方法，結(jié)合邊界層黏性修正，便可以在滿足壓力分布、激波阻力、誘導(dǎo)阻力及摩擦阻力的計算要求的同時，有效減少計算時間的耗費。

1.2.1 笛卡爾網(wǎng)格及MGAERO軟件

在氣動彈性優(yōu)化中進(jìn)行氣動與結(jié)構(gòu)的耦合分析時，要求結(jié)構(gòu)的變形必須可以準(zhǔn)確地反映在CFD網(wǎng)格上。笛卡爾網(wǎng)格是一種特殊的空間離散網(wǎng)格，它網(wǎng)格的生成不依賴于幾何外形，免除了貼體網(wǎng)格生成的各種問題。因此，它在需要反復(fù)修改網(wǎng)格的氣動彈性優(yōu)化問題中十分適用。

MGAERO是一個以求解歐拉方程為基礎(chǔ)的CFD軟件。相較其他的CFD軟件，MGAERO突出的特點便是采用了多重笛卡爾網(wǎng)格嵌套技術(shù)。另外，它所具備的的局部網(wǎng)格加密技術(shù)也可以提高效率，減少計算量，加快收斂。MGAERO軟件的計算流程如圖1所示。

圖1 MGAERO軟件的計算流程圖

1.2.2 黏性修正

在計算摩擦阻力之前需要計算無黏流理論的表面流線，再通過積分該表面流線進(jìn)行邊界層的黏性修正。得到邊界層的厚度后，修改原無黏流的邊界條件，得到考慮黏性后的壓力分布，同時表面流線的改變又會影響邊界層厚度，由此即可進(jìn)行有黏/無黏的迭代修正。

其中，層流邊界層計算采用Cohen-Reshotko方法，邊界層轉(zhuǎn)捩計算采用Granville方法，湍流邊界層計算采用Green方法[3]。

邊界層修正是一個迭代的過程，大概需要3～5次黏性/無黏循環(huán)，本文的黏性/無黏迭代的流程如圖2所示。

圖2 有黏/無黏迭代流程圖

1.3 曲面樣條插值

在使用MGAERO軟件計算機翼巡航阻力前，需要將NASTRAN對模型機翼進(jìn)行氣動彈性分析后得到的位移插值到MGAERO所使用的三維氣動模型之上，以求解機翼的巡航阻力。

本文使用曲面樣條函數(shù)插值方法，其主要思路是：利用已知條件(即已知的獨立結(jié)點及這些結(jié)點所對應(yīng)的位移)構(gòu)成一個逼近函數(shù)去代替真實函數(shù)，從而計算出任何插值點處的位移以及這些插值點處的導(dǎo)數(shù)[4]。逼近函數(shù)的形式見式(2)。

式(2)中，ε為給定的常數(shù)，一般平坦曲面ε=10-2～1，有奇性的曲面取ε=10-5～10-6；c1,c2,L,cN+1+n為待定系數(shù)；xpi為第i個結(jié)點的第p維坐標(biāo)。

逼近函數(shù)確定后，便可以求出出任一向量Xk對應(yīng)的函數(shù)近似值，此即為插值。

1.4 氣動彈性優(yōu)化

氣動彈性優(yōu)化是一個標(biāo)準(zhǔn)的帶約束優(yōu)化問題，即在n維空間中搜索一組設(shè)計變量是目標(biāo)函數(shù)F(v)最小[5]。

MinF(v)=F(x1,x2,K,xn)

ST.gj(v)≤0,j=1,2,K,m

vimin≤vi≤vimax,i=1,2,k,n

其中，v為設(shè)計變量，式(3)為目標(biāo)函數(shù)；m為約束的個數(shù)，通過式(4)定義約束；為設(shè)計變量的個數(shù)，式(5)定義了設(shè)計變量的上下界。

本文使用遺傳算法進(jìn)行優(yōu)化，相比敏度算法，遺傳算法具備較強的全局搜索能力，能有效避免陷入局部最優(yōu)解。

1.5 Kriging模型

代理模型可有效解決優(yōu)化問題中計算耗費大、效率低的問題。響應(yīng)面方法是一類比較常見的代理模型，其中Kriging模型構(gòu)建方便，精度高，在實際工程問題中應(yīng)用較為廣泛[6]。

它是一種改進(jìn)的線性回歸方法，由線性回歸和非參數(shù)兩部分構(gòu)成[6-7]。非參數(shù)部分通常采用隨機方法，可表示為：

y(x)=F(x)+z(x)

(6)

式(6)中，F(xiàn)(x)是線性回歸項，是對全部設(shè)計空間的模擬，通過已知的響應(yīng)進(jìn)行數(shù)據(jù)估計；z(x)為非參數(shù)隨機項，是對全部設(shè)計空間的背離，它的協(xié)方差矩陣可表示為：

Cov[z(xi),z(xj)]=σ2G(xi,xj)

(7)

式(7)中，G(xi,xj)為任意兩點的空間相關(guān)方程，可根據(jù)需求定義，對模型的精度起決定性作用，應(yīng)用最廣泛的是高斯相關(guān)方程。

(8)

式(8)中，θk為相關(guān)性參數(shù)，ndv為每組采樣點中的維數(shù)。

相關(guān)性參數(shù)θk可由極大似然法求得，最終可轉(zhuǎn)化為ndv維的無約束非線性最優(yōu)問題。求解后便可完全建立一個Kriging模型。

1.6 參數(shù)化建模方法

機翼外形最基本的參數(shù)是展弦比、根梢比、前緣后掠角、上反角等，這些參數(shù)可以唯一確定翼面的形狀[8]。機翼的外形參數(shù)改變后，雖然所有節(jié)點的位置都發(fā)生了改變，但機翼的基本結(jié)構(gòu)形式保持不變。所以變參只需改變點的位置即可，在NASTRAN中便體現(xiàn)為對于所有GRID卡片坐標(biāo)值的修改。變參新建模型的關(guān)鍵在于兩點：一是建立初始模型關(guān)鍵點的位置與外形參數(shù)的函數(shù)關(guān)系；二是找到初始模型的所有點和新模型點之間的位置函數(shù)關(guān)系，并以此關(guān)系生成新的GRID卡片。

在整個變參過程中，所涉及到的參數(shù)包括兩類：外形參數(shù)和初始關(guān)鍵點參數(shù)。

外形參數(shù)有內(nèi)翼根梢比η1、內(nèi)翼展弦比λ1、內(nèi)翼上反角τ1、外翼根梢比η2、外翼展弦比λ2、外翼上反角τ2、后掠角χ，共計七個參數(shù)，這七個參數(shù)也是對機翼氣動彈性性能產(chǎn)生影響的主要參數(shù)。

初始關(guān)鍵點參數(shù)(與初始模型相關(guān))有翼根前后緣點A1、A2、KINK位置(機翼后緣轉(zhuǎn)折處)前后緣點A3、A4、翼稍前后緣點A5、A6，總計6個點的xyz軸位置參數(shù)，如圖3所示。

需要指出的是，為計算方便，本文的展弦比中的弦長均為翼根處的弦長。同時，在計算中將內(nèi)翼展弦比修正為內(nèi)翼的絕對展長/根部弦長，外翼展弦比同理，外翼根梢比修正為KINK位置弦長/翼稍弦長。

圖3 機翼平面投影簡化圖與參數(shù)說明示意圖

圖4 機翼關(guān)鍵點及比例系數(shù)計算示意圖

如圖3所示，機翼可從KINK位置分為內(nèi)翼和外翼，這兩部分可看做兩個連接的梯形。就一個梯形區(qū)域而言，如圖4所示，對機翼進(jìn)行變參計算的流程如圖5所示。其中，各比例系數(shù)計算式見式(9)～式(11)。

k1=(yG0-yA1)/l1

(9)

k2=(xG0-xP0)/(xQ0-xP0)

(10)

K3=(zG0-zP0)/(xQ0-xP0)

(11)

其中，k1為展向比例系數(shù)，k2為弦向比例系數(shù)，k3為厚度方向比例系數(shù)。

圖5 變參計算流程圖

2 機翼氣動/結(jié)構(gòu)綜合優(yōu)化框架

2.1 常規(guī)優(yōu)化流程

本文的綜合優(yōu)化設(shè)計方法考慮了各個子學(xué)科間的耦合作用，同時得到最佳氣動外形、結(jié)構(gòu)布局以及結(jié)構(gòu)參數(shù)分布，避免反復(fù)設(shè)計，從而為機翼設(shè)計提供參考。在滿足精度及效率的前提下，氣動子學(xué)科選用亞聲速偶極子格網(wǎng)法以及基于Euler方程的CFD方法；結(jié)構(gòu)子學(xué)科選用有限元分析方法，優(yōu)化算法采用遺傳算法。

整體優(yōu)化流程如圖6所示。

圖6 常規(guī)綜合優(yōu)化流程圖

常規(guī)優(yōu)化流程中由于需要對每一個個體進(jìn)行CFD氣動阻力計算，時間耗費巨大，故引入Kriging代理模型進(jìn)行機翼氣動/結(jié)構(gòu)綜合優(yōu)化。此常規(guī)優(yōu)化流程則可將氣動阻力約束去掉進(jìn)行優(yōu)化(即不執(zhí)行圖3中紅色虛線框內(nèi)的操作)，可將此結(jié)果作為后者結(jié)果的對比。

2.2 基于Kriging模型的優(yōu)化流程

基于代理模型的優(yōu)化流程包括代理模型的建立和遺傳優(yōu)化兩個部分，如圖7所示。

具體步驟如下：

1)準(zhǔn)備初始模型(NASTRAN：機翼三維板桿模型和平板氣動模型；MGAERO：三維氣動模型)。

2)準(zhǔn)備輸入文件，輸入文件包括控制文件和模型文件。

3)通過拉丁超立方抽樣在設(shè)計變量空間內(nèi)選取2 000個樣本點。

4)根據(jù)平板氣動模型和結(jié)構(gòu)模型進(jìn)行氣動彈性分析，將氣動彈性分析得到的變形插值到CFD模型上，對插值處理后模型的巡航外形進(jìn)行流場分析。

5)對樣本中的每一個個體執(zhí)行步驟4)～5)的操作，得到樣本的響應(yīng)值。

6)構(gòu)建代理模型。

構(gòu)建完畢后，便可進(jìn)入到遺傳優(yōu)化中進(jìn)行計算，在滿足精度要求的前提之下，代理模型則類似于一個“黑匣子”，在輸入不同的個體后，便可以近似地計算出各個個體的響應(yīng)，極大地提高優(yōu)化設(shè)計效率。

在優(yōu)化滿足結(jié)束條件，輸出最優(yōu)個體后，再對該個體進(jìn)行真實響應(yīng)值的計算，若不滿足約束要求，則更新個體，重新生成初代群體，重新開始優(yōu)化；若滿足約束要求，則優(yōu)化結(jié)束。

3 算例模型

3.1 三維有限元模型

本文研究所用的飛機機翼半翼展約29 m，采用雙梁式結(jié)構(gòu)。機翼上下蒙皮、前后梁腹板采用板單元；前后梁凸緣、上下桁條采用桿單元。材料為鋁合金，彈性模量72 GPa，剪切模量27 GPa。有限元模型及機翼內(nèi)部結(jié)構(gòu)如圖8～圖9所示。

圖8 有限元模型

圖9 機翼內(nèi)部結(jié)構(gòu)示意圖

模型機身部分采用很硬的殼單元，未建立內(nèi)部結(jié)構(gòu)。機翼部分為模型的主要部分，建模則較為細(xì)致。

3.2 平板氣動模型

本文將模型的右半部分機翼及機身劃分成為具有不同網(wǎng)格密度的7塊平板升力面，如圖10所示。

圖10 平板氣動模型

3.3 CFD三維氣動模型

本文在MGAERO中建立的三維氣動模型及其空間網(wǎng)格如圖11所示。

圖11 三維氣動模型

4 優(yōu)化問題描述

本文對飛翼布局大型客機機翼的氣動/結(jié)構(gòu)綜合優(yōu)化共包含201個設(shè)計變量，其中包括7個機翼外形參數(shù)和194個結(jié)構(gòu)參數(shù)(上下蒙皮厚度及前后梁上下凸緣面積)。

優(yōu)化以結(jié)構(gòu)質(zhì)量最小為目標(biāo)函數(shù)，約束條件包括：翼尖位移、翼尖扭角、顫振速度及氣動阻力。

本文對飛翼布局大型客機機翼氣動/結(jié)構(gòu)進(jìn)行綜合優(yōu)化給定的巡航條件為11 000m的飛行高度，馬赫數(shù)為0.8，配平重量為185t。

5 代理模型的構(gòu)建及精度校驗

圖12 目標(biāo)函數(shù)相對誤差

圖13 翼尖扭角相對誤差

圖14 翼尖位移相對誤差

圖15 氣動阻力相對誤差

本文使用2 000個樣本用于模型構(gòu)建，另外50個樣本用于精度校驗，各響應(yīng)的相對誤差如圖12～圖15所示。目標(biāo)函數(shù)的平均相對誤差為0.43%，翼尖位移的平均相對誤差為4.71%，翼尖扭角的平均相對誤差為3.73%，氣動阻力的平均相對誤差為4.86%，顫振的相對誤差由于其特殊性恒為0%，滿足優(yōu)化設(shè)計的需求。

6 計算結(jié)果

本文對不考慮氣動阻力約束的常規(guī)優(yōu)化(以下簡稱常規(guī)優(yōu)化)以及基于代理模型考慮氣動阻力約束的優(yōu)化(以下簡稱代理優(yōu)化)，二者優(yōu)化的結(jié)果見表1。

表1 優(yōu)化結(jié)果

代理優(yōu)化結(jié)束條件設(shè)定為達(dá)到指定的最大迭代代數(shù)，本文共迭代200代，每代個體數(shù)為1 000。代理優(yōu)化每代最優(yōu)個體的目標(biāo)函數(shù)響應(yīng)及每代個體的適應(yīng)度情況如圖16、圖17所示。

圖16 目標(biāo)函數(shù)響應(yīng)隨迭代次數(shù)的變化趨勢

圖17 每代最優(yōu)個體適應(yīng)度和平均適應(yīng)度變化趨勢

采用這兩類方法對機翼進(jìn)行綜合優(yōu)化設(shè)計后的機翼外形參數(shù)見表2。

表2 外形參數(shù)對比

使用代理模型進(jìn)行氣動/阻力綜合優(yōu)化總耗時為3 360 s，不到一個小時，可以在滿足精度的前提下，極大地提高計算效率。

從設(shè)計結(jié)果的響應(yīng)值可以看出，不考慮氣動阻力約束的常規(guī)優(yōu)化雖有效地減小了目標(biāo)函數(shù)，但氣動阻力系數(shù)卻有所增加，不滿足設(shè)計要求，如此便造成了反復(fù)設(shè)計。在引入了Kriging代理模型進(jìn)行氣動/結(jié)構(gòu)綜合優(yōu)化，將氣動阻力約束考慮在后，可在精度允許范圍內(nèi)，顯著減少機翼的結(jié)構(gòu)質(zhì)量以及氣動阻力。

可以看出，常規(guī)優(yōu)化后的機翼外形整體變小，這也是機翼整體減重的重要原因之一；而代理優(yōu)化后的外形變化最為明顯的則是機翼外翼部分變長、后掠角增大。

初始機翼及在兩類優(yōu)化方法下，最優(yōu)個體的上下翼面壓強分布對比如表3所示。從表中各圖可以看出，由于機翼外形的變化，機翼表面靠近前緣部分的壓力分布變化較大，而其他區(qū)域變化較小。代理優(yōu)化后，飛機上下表面的壓力變化更為光順，階躍突變減少，這有利于在保證一定升力的同時，減小氣動阻力。

表3 各模型上下翼面Cp分布對比

在0.8馬赫數(shù)的巡航條件下，大型客機機翼上表面會出現(xiàn)局部激波，相比于黏性阻力，激波阻力和誘導(dǎo)阻力是飛行阻力中非常重要的。前緣后掠角的增加可有效減少上翼面的超聲速區(qū)域，從而降低局部激波阻力。表4給出了三個階段設(shè)計結(jié)果機翼與初始模型上翼面馬赫數(shù)的對比，可以看出，與初始模型及常規(guī)優(yōu)化結(jié)果比較，考慮阻力約束的代理優(yōu)化結(jié)果其超聲速區(qū)域可大大降低。

表4 各模型上下翼面Ma分布對比

7 結(jié)論

對于飛翼這種采用新型氣動布局的設(shè)計方案，相較于傳統(tǒng)方案，更需注重考慮氣動與結(jié)構(gòu)的耦合關(guān)系。由本文對飛翼布局客機機翼氣動/結(jié)構(gòu)綜合優(yōu)化所得結(jié)果的分析和對比中可以看出：一方面，將阻力約束考慮在內(nèi)的氣動/結(jié)構(gòu)綜合優(yōu)化，可以在減重的同時，顯著地減少氣動阻力，同時獲得機翼的最優(yōu)氣動外形和結(jié)構(gòu)布局；另一方面，基于代理模型的氣動/結(jié)構(gòu)綜合優(yōu)化極大地減少了時間的耗費，且其結(jié)果的精度滿足要求，具備有較好的工程適用性。