基于伴隨方法的飛翼布局多目標(biāo)氣動優(yōu)化設(shè)計

劉曉冬, 張沛良, 何光洪, 王永恩, 楊旭東

(1.沈陽飛機設(shè)計研究所, 遼寧 沈陽 110035; 2.西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072)

目前結(jié)合數(shù)值模擬技術(shù)的氣動外形優(yōu)化設(shè)計方法主要分為兩類：隨機類方法和梯度類方法[1-2]。隨機類方法追蹤目標(biāo)函數(shù)值相關(guān)信息，帶有“隨機”特性，具有全局性好、不要求設(shè)計變量連續(xù)分布以及導(dǎo)數(shù)存在等假設(shè)的優(yōu)點，如遺傳算法、代理模型方法等[3]，但其缺點是在采用N-S方程進行多設(shè)計變量的多目標(biāo)氣動設(shè)計時，CFD計算量巨大。盡管在過去幾十年里計算機能力大大提高，但執(zhí)行大量工程實際設(shè)計仍然不切實際，目前遺傳算法應(yīng)用仍局限于設(shè)計變量較少的設(shè)計問題；代理模型計算量相對遺傳算法雖然有所減少，但其適用性受到設(shè)計變量取值范圍和空間樣本點分布的影響，不合理的設(shè)計變量及樣本點分布將難以獲得高精度的代理模型，從而導(dǎo)致設(shè)計結(jié)果偏離設(shè)計要求[4]。上述缺陷限制了隨機類方法在多設(shè)計變量設(shè)計問題中的應(yīng)用。

在梯度類方法中，基于伴隨理論的氣動優(yōu)化設(shè)計方法(控制理論、共軛方法)在兼顧梯度的精確、快速求解和計算量方面取得較大的進步[1,5]。其以偏微分方程系統(tǒng)的控制理論為基礎(chǔ)，把物面邊界作為控制函數(shù)，基于拉格朗日的觀點將線化流動方程作為約束引入到目標(biāo)函數(shù)與設(shè)計變量的表達式中，將設(shè)計問題轉(zhuǎn)化為控制問題，計算量只相當(dāng)于兩倍流場計算量，與設(shè)計變量數(shù)目無關(guān)。國外采用伴隨方法在翼型、機翼和翼身組合體乃至全機的氣動優(yōu)化設(shè)計中，取得了一系列的研究成果[5-11]。在國內(nèi)，比較具代表性的是西北工業(yè)大學(xué)以喬志德教授為核心的團隊，先后對翼型、機翼和翼身組合體開展了基于連續(xù)/離散伴隨方法的Euler和N-S方程的優(yōu)化設(shè)計，獲得了比較滿意的結(jié)果[12-14]；西安交通大學(xué)的豐鎮(zhèn)平教授將伴隨方法應(yīng)用到內(nèi)流領(lǐng)域的二維、三維透平葉柵氣動優(yōu)化設(shè)計中[1,15]；西北工業(yè)大學(xué)白俊強團隊近幾年采用伴隨方法在大型客機氣動減阻和聲爆特性優(yōu)化方面做了研究[16-17]。總的來說，目前國內(nèi)基于伴隨方法在氣動外形優(yōu)化設(shè)計中取得了一些較為成熟的結(jié)果，提高了伴隨方法在氣動外形設(shè)計方面的工程應(yīng)用前景，今后還將不斷發(fā)展完善。

飛翼布局由于良好的氣動效率及隱身特性，在軍用飛機得到了廣泛的應(yīng)用，同時飛翼布局又存在操縱效能低，配平損失等典型問題[18]。所以無尾飛翼布局的氣動外形設(shè)計是涉及總體、氣動、控制、結(jié)構(gòu)、隱身等專業(yè)約束的綜合設(shè)計問題，即多目標(biāo)多約束的設(shè)計問題。一般來說，要求飛機在巡航狀態(tài)時具有較高升阻比，同時具有較小的低頭力矩，從而不會引起較大的配平損失，而這兩種要求往往是矛盾的，尤其在飛翼布局上表現(xiàn)的較為突出，即提高巡航升阻比的同時會帶來較大的低頭力矩。為了解決飛翼布局設(shè)計中的這種問題，本文基于伴隨方法的基本原理，發(fā)展了一種考慮氣動、結(jié)構(gòu)約束的飛翼布局多目標(biāo)多約束氣動優(yōu)化方法，并進行了典型算例驗證。

1 流動控制方程

采用結(jié)構(gòu)化貼體網(wǎng)格進行空間數(shù)值離散，并按照求和約定，在計算域中，N-S方程可表示為：

(1)

2 伴隨方程及邊界條件

2.1 伴隨方法優(yōu)化設(shè)計原理

氣動優(yōu)化問題是以外形變化對氣動特性的影響為基礎(chǔ)而進行的,氣動特性的目標(biāo)函數(shù)可表述為:

(2)

dBξ,dDξ分別為計算空間中的表面與空間積分單元，M與P取決于流動變量w及計算空間的矩陣S。

在滿足流場控制方程約束條件下,氣動外形的變化將導(dǎo)致流動變量變分δw與矩陣變分δS,目標(biāo)函數(shù)的變化可表示為

(3)

式中:δM=[Mw]Ⅰδw+δMⅡ,δP=[Pw]Ⅰδw+δPⅡ；下標(biāo)Ⅰ表示由流動變量變化δw引起的貢獻；下標(biāo)Ⅱ表示由矩陣變化δS引起的貢獻。

定常狀態(tài)下,氣動外形變化的約束方程可表述為

(4)

引入伴隨矢量ψ=(ψ1,φ1,φ2,φ3,θ)T,與(4)式在整個計算空間求積,則有

(5)

假定ψ可微,(5)式分部積分及高斯定理可進一步寫成

(6)

整合(3)、(5)、(6)式,目標(biāo)函數(shù)變分可寫為

(7)

令(7)式中空間積分項流動物理量變分δw系數(shù)項組合為0,則可得伴隨方程

(8)

令(7)式邊界積分中流動物理量變分δw的系數(shù)項組合在一起為0,則得對應(yīng)的伴隨方程邊界條件

(9)

剩余項即為目標(biāo)函數(shù)的梯度求解公式

(10)

2.2 伴隨方程

鑒于伴隨方程理論推導(dǎo)的復(fù)雜性,詳細過程可參考文獻[12],此處僅給出最終的數(shù)學(xué)表達式。

(11)

式中：矢量Y中a為聲速；Pr為普朗特數(shù)。

2.3 減阻優(yōu)化多目標(biāo)函數(shù)

為了實現(xiàn)飛翼布局設(shè)計點有效減阻優(yōu)化,綜合考慮以下設(shè)計目標(biāo)。

通過加權(quán)組合方法定義如(12)式的統(tǒng)一目標(biāo)函數(shù)

(12)

2.4 伴隨方程邊界條件

定義自由來流馬赫數(shù)M∞,壓力P∞,迎角α以及參考面積Sref,力矩參考點(xref,yref),易知

則升力系數(shù)、阻力系數(shù)、俯仰力矩系數(shù)中流動變量的變分為

由(12)式得,目標(biāo)函數(shù)中流動變量的變分為

(19)

根據(jù)(9)式在參考文獻[12]中伴隨方程邊界條件的推導(dǎo)可知

(20)

由(20)式得方程組,并求解伴隨邊界條件為

(21)

2.5 梯度方程

同理,根據(jù)(10)式目標(biāo)函數(shù)及流動控制方程中對矩陣的變分項,即可獲得梯度求解式為

(22)

3 優(yōu)化設(shè)計流程與數(shù)值求解

圖1給出了優(yōu)化設(shè)計流程圖。

圖1 優(yōu)化設(shè)計流程圖

網(wǎng)格生成:采用無限插值法生成結(jié)構(gòu)化計算網(wǎng)格,同時采用正交控制、加權(quán)平均光順和法向量控制等措施確保網(wǎng)格質(zhì)量。

流場數(shù)值求解:采用Jameson的中心格式有限體積法進行空間離散,五步Runge-Kutta顯示格式時間推進,同時加入人工黏性抑制振蕩,采用當(dāng)?shù)貢r間步長、隱式殘值光順、多重網(wǎng)格等加速收斂措施;湍流模型采用B-L湍流模型。

伴隨方程數(shù)值求解:采用與N-S方程類似的數(shù)值解法。

設(shè)計變量及優(yōu)化算法:采用Hicks-Henne形狀函數(shù)描述設(shè)計變量對物體外形變化的影響,同時采用最速下降法進行梯度搜索。

4 算例分析

4.1 算例1:某小展弦比飛翼減阻設(shè)計

設(shè)計點:Ma=0.85,α=3°,控制剖面取機翼3個剖面,每個剖面26個設(shè)計變量,共78個設(shè)計變量。計算網(wǎng)格采用C-H網(wǎng)格,如圖2所示。先后開展了2種不同約束下的優(yōu)化設(shè)計。

圖2 C-H計算網(wǎng)格示意圖

設(shè)計1:升力、面積約束下的減阻優(yōu)化設(shè)計,根據(jù)設(shè)計狀態(tài),限定升力系數(shù)、各剖面面積相對約束值變化不超過5%。選取目標(biāo)函數(shù)中各部分的權(quán)值分別為:Ω1=120,Ω2=1,Ω3=0,Ω4=1。

表1給出了初始與設(shè)計外形的氣動特性與幾何特性具體數(shù)值對比。從中看到,初始阻力系數(shù)為0.013 1,升力系數(shù)為0.213 7,優(yōu)化設(shè)計25步后,阻力系數(shù)變?yōu)?.010 6,升力系數(shù)變?yōu)?.214 3,阻力系數(shù)下降19.5%,而升力系數(shù)、各控制剖面面積滿足約束條件。

表1 某小展弦比飛翼不同約束條件優(yōu)化前后特性對比

注意到相比初始外形,優(yōu)化后帶來了較大的低頭力矩,這樣就會帶來較大的配平損失,從而導(dǎo)致實際使用升阻比降低。

設(shè)計2:為改善設(shè)計1帶來的問題,進行了升力、俯仰力矩、面積共同約束下的減阻優(yōu)化設(shè)計,對應(yīng)目標(biāo)函數(shù)中各部分的權(quán)值分別取:Ω1=120,Ω2=1,Ω3=0.01,Ω4=1,保證|Cm|≤0.004。表1為優(yōu)化后具體的數(shù)值變化。迭代優(yōu)化12步后,阻力系數(shù)減小為0.010 8(與設(shè)計1保持相同量級),下降約18.0%;升力系數(shù)變?yōu)?.206 7,減小3.3%,變化相對較大,但也滿足約束條件;俯仰力矩系數(shù)由初始的0.001 8變?yōu)?0.001 6,滿足約束條件;各控制剖面面積變化同樣滿足約束指標(biāo)。

圖4給出了初始外形與2種約束下優(yōu)化設(shè)計外形機翼展向不同站位剖面壓力分布對比。可以清楚地看到2種設(shè)計外形機翼上表面的激波都被很大程度削弱,不同之處在于有力矩約束得到的外形壓力分布在50%弦長前負壓值較大,50%弦長后負壓值則較小,則對應(yīng)上表面的壓心靠前,因此不會帶來很大的低頭力矩,與氣動力計算結(jié)果一致。

圖3 初始外形與無/有力矩約束優(yōu)化設(shè)計外形表面壓力對比

通過算例1的2種優(yōu)化設(shè)計結(jié)果對比,驗證了所定義的目標(biāo)函數(shù),推導(dǎo)的伴隨方程邊界條件及梯度求解公式是正確并有效的。

4.2 算例2:某大展弦比飛翼減阻設(shè)計

采用本文所發(fā)展的方法,對某大展弦比飛翼進行跨聲速狀態(tài)升力、俯仰力矩和面積約束下的減阻優(yōu)化設(shè)計。設(shè)計點:Ma=0.75,α=4°,目標(biāo)函數(shù)中各部分的權(quán)值分別取:Ω1=50,Ω2=2,Ω3=0.001,Ω4=0.5,限定升力系數(shù)、各剖面面積相對約束值變化不超過5%,俯仰力矩-0.004≤Cm≤0.008;控制剖面取機翼的4個剖面,每個剖面26個設(shè)計變量,加上4個剖面扭轉(zhuǎn)角,共108個設(shè)計變量。

表2給出了優(yōu)化前后氣動力系數(shù)及控制剖面面積的具體數(shù)值變化。優(yōu)化迭代8步,阻力系數(shù)由初始的0.016 65減小為0.015 06,下降約9.55%;升力系數(shù)由0.361變?yōu)?.355,減小1.66%,滿足約束條件;俯仰力矩系數(shù)由初始的0.006 2變?yōu)?.003 6,滿足約束條件;各控制剖面面積變化也滿足約束指標(biāo)。

表2 某大展弦比飛翼優(yōu)化前后氣動、幾何特性對比

圖4給出了機翼展向不同站位剖面壓力分布及外形對比?？梢钥吹皆O(shè)計外形上表面壓力負壓峰值區(qū)域減小,逆壓梯度變小,展向不同位置的激波強度都有不同程度的減弱,尤其在展向60%～70%范圍激波削弱明顯。對應(yīng)剖面外形的主要變化趨勢是最大厚度略有減小,且弦向位置略有后移;扭轉(zhuǎn)角主要是靠近對稱面的剖面有了一個較小的正扭轉(zhuǎn)角,其余剖面變化不大。

圖4 初始外形與優(yōu)化外形展向不同剖面外形及壓力對比

5 結(jié) 論

本文針對飛翼布局氣動優(yōu)化設(shè)計中的多目標(biāo)多約束問題，基于伴隨方法和N-S方程發(fā)展了一種氣動優(yōu)化設(shè)計方法，并先后進行了2種不同展弦比飛翼布局的跨聲速減阻優(yōu)化設(shè)計，結(jié)果表明：

1) 通過構(gòu)建合理的統(tǒng)一目標(biāo)函數(shù)形式來解決氣動優(yōu)化設(shè)計中的多目標(biāo)多約束問題是合適的，根據(jù)伴隨方法基本原理推導(dǎo)的伴隨方程物面邊界條件以及梯度求解方程是正確、有效的；

2) 所發(fā)展的設(shè)計方法在飛翼布局的多約束氣動優(yōu)化設(shè)計問題上，具有良好的設(shè)計效果和優(yōu)化效率，因此在工程上具有廣闊的應(yīng)用前景。