呂志斌 萬志強 / LYU Zhibin WAN Zhiqiang

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

0 引言

飛機是高度綜合的現(xiàn)代科學(xué)技術(shù)的體現(xiàn)，飛機結(jié)構(gòu)設(shè)計又是飛機設(shè)計的主要階段[1]。但隨著現(xiàn)代飛機尤其是民用客機與運輸機追求高性能和低結(jié)構(gòu)重量等方面的要求，飛機結(jié)構(gòu)柔性較大，導(dǎo)致飛機結(jié)構(gòu)變形增大。加之近代工業(yè)在復(fù)合材料、智能材料等新型材料研究上的突破和應(yīng)用，使得氣動彈性問題在飛機設(shè)計中越來越凸顯出來[2]。傳統(tǒng)的“設(shè)計-校核-修改”的結(jié)構(gòu)設(shè)計模式越來越難以滿足現(xiàn)代飛機對于高性能和低重量的雙重要求。

基于此，隨著氣動彈性學(xué)科的發(fā)展，以克服氣動彈性效應(yīng)所帶來的不利影響的同時盡可能地降低飛機結(jié)構(gòu)重量為目標(biāo)的氣動彈性結(jié)構(gòu)優(yōu)化技術(shù)已經(jīng)逐漸成為一種主動設(shè)計的手段，應(yīng)用于飛機設(shè)計的各個階段[3]。

然而目前工程上在進行氣動彈性分析求解時氣動力數(shù)據(jù)來源仍然是采用低階面元法，該方法計算效率高但精度低；而CFD(Computational Fluid Dynamics，簡稱CFD)方法雖然計算精度高，但效率太低，而且難以實現(xiàn)優(yōu)化過程的迭代計算。當(dāng)前的氣動彈性優(yōu)化算法則主要集中于單一的優(yōu)化算法，如可行方向法、遺傳算法等。這些算法均有其自身的局限性，如敏度信息不易獲得、搜索精度低、收斂速度慢等問題[4]。同時考慮到氣動彈性結(jié)構(gòu)優(yōu)化本身往往具有計算量大、效率低、時間長、設(shè)備要求高等特點，難以滿足方案設(shè)計初期對時間的要求。

針對上述問題，本文提出了采用基于高階面元法的靜氣動彈性分析方法來求解約束響應(yīng)，使用遺傳/混合算法作為尋優(yōu)算法，并引入Kriging(克立格)代理模型來減少計算量和優(yōu)化時間。最后通過一個成熟的大展弦比機翼為算例，驗證本文提出的靜氣動彈性結(jié)構(gòu)優(yōu)化方法的可行性與可靠性，從而為飛機方案設(shè)計階段提供一種新的優(yōu)化方法。

1 基于高階面元法的靜氣動彈性分析

1.1 高階面元法

高階面元法與低階面元法一樣都是線性方法，都是通過在物面網(wǎng)格上進行布置偶極子或者源來模擬氣流流經(jīng)物體表面的流場，都是對線化后的小擾動位勢流方程進行求解[5]：

(1)

(2)

其中，M∞為來流馬赫數(shù)，?為擾動位函數(shù)。

而與低階面元法不同的是，低階面元法在平板氣動力網(wǎng)格上的壓力點需布置在每個氣動網(wǎng)格的四分之一弦線上，高階面元法則沒有類似限制。除可以按照幾何外形建立三維的氣動力計算模型外，高階面元法在氣動力計算模型的表面網(wǎng)格上布置連續(xù)的點源，在結(jié)點上布置連續(xù)的偶極子，線性分布奇點強度，如圖1所示[5]。

圖1 網(wǎng)格元素奇點分布示意圖

1.2 靜氣動彈性分析框架

與常用的基于低階面元法的Nastran軟件不同，本文使用的靜氣動彈性分析方法的氣動力數(shù)據(jù)來源于基于高階面元法的Panair開源程序，結(jié)構(gòu)計算基于柔度法的Nastran靜力分析模塊，剛性配平為求解基本的飛行力學(xué)平衡方程，載荷導(dǎo)數(shù)計算采用差分法。圖2所示即為基于高階面元法的靜氣動彈性分析方法框架。

圖2 基于高階面元法的靜氣動彈性分析方法框架

2 代理模型方法

代理模型方法主要包含兩方面的內(nèi)容：一是構(gòu)建模型的樣本點的取樣方法，又稱實驗設(shè)計方法；二是進行數(shù)據(jù)擬合和預(yù)測的模型構(gòu)建方法，又稱近似方法[6]。

2.1 拉丁超立方抽樣

拉丁超立方取樣方法于1979年由McKay提出，該方法最大的優(yōu)點是不需考慮樣本點的維度和取樣數(shù)目，通過控制取樣點的位置可以盡可能避免取樣點在小領(lǐng)域內(nèi)重合，從而實現(xiàn)取樣點在設(shè)計空間的均勻分布[7]。

假設(shè)需要在m維向量空間中抽取n個樣本，其具體步驟為：

1)將m維向量空間中的每一維均按策略(一般采用均勻分布的策略)分為互不重疊的n個區(qū)間，使得每個區(qū)間有相同的被選擇概率；

2)從向量空間每一維的每個區(qū)間中隨機抽取一個點作為該維該空間的值，按照劃分的區(qū)間共抽取m×n個值；

3)再從每一維中隨機抽取2種選擇的點組成m維向量，共抽取n個，至此樣本點抽取完成。

如圖3所示為用Matlab仿真出來的二維向量空間的拉丁超立方取樣的樣本點分布圖，取樣點數(shù)目為8個。

圖3 二維空間的拉丁超立方取樣

2.2 Kriging代理模型方法

Kriging模型是一種改進的線型回歸模型，其函數(shù)表達式中既包含了線性回歸部分，又包含了非參數(shù)部分[8]：

y=F(x)+z(x)

(3)

式中，F(xiàn)(x)為線性回歸部分，具有全局近似的特點；z(x)為非參數(shù)隨機部分，具有局部偏差的特點。

線性回歸部分可以表示為：

F(x)=f1(x)β1+…+fp(x)βp

(4)

按照f(x)形式的不同，可分為常數(shù)型、線性型和二次型。本文選擇如下所示的二次型，式中β可通過廣義最小二乘估計得到：

f1(x)=1

f2(x)=x1,…,fn+1(x)=xn

(5)

非參數(shù)部分具有如下的統(tǒng)計學(xué)特性：

E[z(x)]=0

Var(z(x)]=σ2

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

(6)

式中，G(θ,xj,xi)為i、j兩點的空間相關(guān)方程，該空間相關(guān)方程也有很多不同種類，對模型構(gòu)造影響很大。本文取Guass函數(shù)作為空間相關(guān)函數(shù)：

(7)

式中，ndv為樣本點的維度，θk是空間相關(guān)函數(shù)在樣本點第j個方向的參數(shù)，為待求解的量。

對于一組給定樣本點X=[x1…xn]T及其響應(yīng)值Y=[y1…yn]T，由樣本點得到待測點的響應(yīng)值為：

(8)

式中c是待求權(quán)系數(shù)向量。根據(jù)無偏條件及估計方差最小條件，結(jié)合拉格朗日乘子，可求得c：

(9)

具體推導(dǎo)過程見參考文獻[8]。

3 優(yōu)化算法

與其他優(yōu)化問題一樣，氣動彈性結(jié)構(gòu)優(yōu)化問題也可以用如下的公式進行描述，即在設(shè)計空間中尋找一組變量使目標(biāo)函數(shù)(一般為結(jié)構(gòu)重量)最小化[9]：

Min.F(v)

(10)

S.T.gj(v)≤0j=1,2,…,ncon

(11)

(12)

其中，式(10)為目標(biāo)函數(shù)，本文為結(jié)構(gòu)重量最??；式(11)為約束條件，包括靜氣彈約束(如發(fā)散速度約束、變形約束等)、動氣彈約束(如顫振速度約束等)、強度約束、操穩(wěn)約束以及氣動約束(如阻力約束等)等；式(12)為設(shè)計變量的上下限，指定設(shè)計空間的值域，又叫邊界約束。

3.1 遺傳算法

遺傳算法是一種借鑒自然選擇和遺傳進化機理而發(fā)明出來的自適應(yīng)概率搜索算法，具有全局尋優(yōu)的能力[10]。作為進化算法的一種，遺傳算法從初代個體開始，以適應(yīng)度函數(shù)為評價標(biāo)準(zhǔn)，通過遺傳算子實現(xiàn)對個體的自然選擇和遺傳進化，從而不斷改良個體直至滿足條件。

圖4所示為遺傳算法流程圖[9]，其基本步驟為：

圖4 遺傳算法流程圖

1)確定優(yōu)化策略，定義優(yōu)化參數(shù)，包括群體規(guī)模、交叉概率、變異概率、程序終止條件等遺傳參數(shù)，選擇編碼方式、選擇方式、交叉方式、變異方式等；

2)隨機產(chǎn)生指定規(guī)模的初代群體(即第一代父群體)；

3)計算群體中每個個體的目標(biāo)函數(shù)值以及約束響應(yīng)值，從而計算個體的適應(yīng)度值，從而對個體進行評估；

4)判斷是否滿足程序終止條件，如果滿足則計算結(jié)束，否則繼續(xù)進行計算；

5)通過基本遺傳操作(選擇、交叉和變異)產(chǎn)生新一代子群體，用子群體代替父群體形成新一代父群體，重復(fù)步驟3和步驟4。

3.2 分形算法

分形算法是一種借鑒樹的生長來實現(xiàn)搜索尋優(yōu)的啟發(fā)式算法，具有較強的局部搜索能力。圖5所示為分形算法流程圖，其中選點、分枝、剪枝為分形算法的三個主要操作運算[11]。分形算法的原理很簡單，通過選點來確定當(dāng)前區(qū)域的最優(yōu)解，為分枝和剪枝提供參考依據(jù)；通過分枝來產(chǎn)生新的樣本點，計算樣本點的適應(yīng)度(為了統(tǒng)一個體評價標(biāo)準(zhǔn)，此處依舊引入遺傳算法中的適應(yīng)度概念代替目標(biāo)函數(shù)值)以判斷是否更新當(dāng)前全局和局部最優(yōu)解；通過剪枝來縮小可行域范圍，使可行域逐漸逼近最優(yōu)解。在分形優(yōu)化算法中，分枝提高了算法的精度，剪枝加快了算法的收斂速度。

3.3 遺傳/分形算法

遺傳/分形混合算法結(jié)合了兩種算法的優(yōu)點，遺傳算法用于全局搜索以避免優(yōu)化算法陷入局部最優(yōu)解，分形算法用于最優(yōu)個體周圍的小范圍局部尋優(yōu)。圖6所示為遺傳/分形混合算法的流程圖，若該流程中不適用分形算法，則圖6為遺傳算法的流程圖。

確定遺傳策略和優(yōu)化參數(shù)后，首先使用遺傳算法對優(yōu)化問題進行一次基本遺傳操作。選擇遺傳操作后的最優(yōu)個體，以該個體為中心確定分形算法的初始可行域以及初始最優(yōu)解(即為該最優(yōu)個體適應(yīng)度)。然后對上述可行域使用分形算法進行尋優(yōu)，并用分形算法得到的最優(yōu)個體代替上一次遺傳算法得到的最優(yōu)個體。完成上述操作后重新評估該代群體，判斷是否滿足算法的終止條件。

圖6 遺傳/分形混合算法流程圖

3.4 靜氣動彈性結(jié)構(gòu)優(yōu)化框架

基于前文提到的基于高階面元法的靜氣動彈性分析方法、代理模型以及遺傳/分形算法，本文建立了一種高效高精度的靜氣動彈性結(jié)構(gòu)優(yōu)化框架，如圖7所示。在該優(yōu)化框架中，遺傳/分形算法為優(yōu)化算法(優(yōu)化時不考慮分形算法即為遺傳算法)，代理模型代替基于高階面元法的靜氣動彈性分析方法作為個體目標(biāo)函數(shù)和約束響應(yīng)的求解器。在實施本文提出的優(yōu)化方法時，首先應(yīng)該完成代理模型的構(gòu)建和精度校驗，如圖8所示。

圖7 靜氣動彈性結(jié)構(gòu)優(yōu)化框架

圖8 代理模型構(gòu)建與精度校驗流程圖

4 算例分析

4.1 算例描述

本文所用算例的氣動和結(jié)構(gòu)模型如圖9和圖10所示，機翼半展長16.27 m，翼尖弦長1.62 m。優(yōu)化時，本文僅對機翼進行優(yōu)化，而不改變機身的結(jié)構(gòu)重量。

圖9 算例氣動模型

圖10 算例結(jié)構(gòu)模型

表1給出了算例的優(yōu)化方案，優(yōu)化狀態(tài)為飛行高度在11 200 m、飛行馬赫數(shù)0.8時的定直平飛狀態(tài)，副翼和方向舵偏角均為0，依靠升降舵和攻角進行配平。

表1 算例的優(yōu)化方案

為了機翼的蒙皮滿足從翼根到翼尖、從后緣到前緣逐漸遞減的設(shè)計準(zhǔn)則，優(yōu)化時分別使用不同的線性函數(shù)來表示蒙皮的厚度變化。圖11和圖12為初始模型機翼上下表面蒙皮的厚度分布。

圖11 機翼上表面蒙皮厚度分度

圖12 機翼下表面蒙皮厚度分度

4.2 代理模型構(gòu)建與精度校驗

在優(yōu)化前，首先需要完成代理模型的構(gòu)建和精度校核。參照圖8的流程首先在機翼的設(shè)計變量空間中使用拉丁超立方取樣法選取4 000個樣本點作為樣本輸入用于模型構(gòu)建，使用隨機方法選取500個樣本點作為測試點用于精度校核，樣本點和測試點均采用基于高階面元法的靜氣動彈性分析方法進行目標(biāo)函數(shù)和約束響應(yīng)計算。

(13)

(14)

復(fù)相關(guān)系數(shù)R的值介于0和1之間，其值越大，說明預(yù)測值和真實值相關(guān)性越大即越接近。一般而言，0.5～1.0即為強相關(guān)。相對均方根誤差RMSE與響應(yīng)值無關(guān)，值介于0和1之間，其值越小，說明預(yù)測精度越高[12]。

圖13～圖16為校驗點的Kriging模型預(yù)測結(jié)果與高階面元法的靜氣動彈性程序計算結(jié)果的對比，表 2說明了Kriging模型的預(yù)測表現(xiàn)。從圖表的結(jié)果來看，Kriging模型的預(yù)測精度很高，可以完全滿足工程應(yīng)用和科學(xué)研究的精度要求。

圖13 結(jié)構(gòu)重量預(yù)測誤差百分比

圖14 翼尖位移預(yù)測誤差百分比

圖15 翼尖扭角預(yù)測誤差百分比

圖16 副翼效率預(yù)測誤差百分比

表2 Kriging模型預(yù)測表現(xiàn)

4.3 優(yōu)化結(jié)果對比與分析

本文使用了遺傳算法和遺傳/分形算法，結(jié)合Kriging代理模型對機翼進行了優(yōu)化，優(yōu)化共計進行了200代，圖17顯示了兩種算法的目標(biāo)函數(shù)隨迭代代數(shù)變化的曲線，表3為算例的優(yōu)化結(jié)果。由于算例使用的機翼是一個成熟的機翼，因此可優(yōu)化的程度并不大。

圖17 優(yōu)化目標(biāo)函數(shù)曲線

從目標(biāo)函數(shù)曲線來看，遺傳/分形混合優(yōu)化算法相比于單一的遺傳算法具有更快收斂速度和更好的全局尋優(yōu)能力，更容易收斂到最優(yōu)解，且最優(yōu)解也更接近全局最優(yōu)解。從算例的優(yōu)化結(jié)果來看，相比于初始模型，結(jié)構(gòu)重量減小，翼尖變形和扭角有一定程度的提高，副翼效率略微減小，但都在約束條件范圍內(nèi)(由于存在約束條件閾值，使得在優(yōu)化時約束響應(yīng)能略微超過約束條件)。

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

5 結(jié)論

本文建立了一種基于高階面元法和遺傳/分形算法的靜氣動彈性結(jié)構(gòu)優(yōu)化算法，并引入Kriging代理模型方法來加快優(yōu)化時間，降低優(yōu)化成本。其中，基于高階面元法的靜氣動彈性優(yōu)化方法用于約束響應(yīng)的求解，遺傳/分形算法用于目標(biāo)函數(shù)尋優(yōu)。最后以一個大展弦比客機機翼為例，通過與遺傳算法的對比，說明了本文提出方法的實用性和更強的尋優(yōu)效率。

值得一提的是，本文在對算例進行優(yōu)化時，大部分的優(yōu)化時間花費在了樣本數(shù)據(jù)的計算上，而代理模型的構(gòu)建和優(yōu)化計算總耗時不超過2h，也顯示了本文的方法在時間上的高效性。