劉曉冬, 張沛良, 何光洪, 王永恩, 楊旭東
(1.沈陽飛機設計研究所, 遼寧 沈陽 110035; 2.西北工業大學 航空學院, 陜西 西安 710072)
目前結合數值模擬技術的氣動外形優化設計方法主要分為兩類:隨機類方法和梯度類方法[1-2]。隨機類方法追蹤目標函數值相關信息,帶有“隨機”特性,具有全局性好、不要求設計變量連續分布以及導數存在等假設的優點,如遺傳算法、代理模型方法等[3],但其缺點是在采用N-S方程進行多設計變量的多目標氣動設計時,CFD計算量巨大。盡管在過去幾十年里計算機能力大大提高,但執行大量工程實際設計仍然不切實際,目前遺傳算法應用仍局限于設計變量較少的設計問題;代理模型計算量相對遺傳算法雖然有所減少,但其適用性受到設計變量取值范圍和空間樣本點分布的影響,不合理的設計變量及樣本點分布將難以獲得高精度的代理模型,從而導致設計結果偏離設計要求[4]。上述缺陷限制了隨機類方法在多設計變量設計問題中的應用。
在梯度類方法中,基于伴隨理論的氣動優化設計方法(控制理論、共軛方法)在兼顧梯度的精確、快速求解和計算量方面取得較大的進步[1,5]。其以偏微分方程系統的控制理論為基礎,把物面邊界作為控制函數,基于拉格朗日的觀點將線化流動方程作為約束引入到目標函數與設計變量的表達式中,將設計問題轉化為控制問題,計算量只相當于兩倍流場計算量,與設計變量數目無關。國外采用伴隨方法在翼型、機翼和翼身組合體乃至全機的氣動優化設計中,取得了一系列的研究成果[5-11]。在國內,比較具代表性的是西北工業大學以喬志德教授為核心的團隊,先后對翼型、機翼和翼身組合體開展了基于連續/離散伴隨方法的Euler和N-S方程的優化設計,獲得了比較滿意的結果[12-14];西安交通大學的豐鎮平教授將伴隨方法應用到內流領域的二維、三維透平葉柵氣動優化設計中[1,15];西北工業大學白俊強團隊近幾年采用伴隨方法在大型客機氣動減阻和聲爆特性優化方面做了研究[16-17]。總的來說,目前國內基于伴隨方法在氣動外形優化設計中取得了一些較為成熟的結果,提高了伴隨方法在氣動外形設計方面的工程應用前景,今后還將不斷發展完善。
飛翼布局由于良好的氣動效率及隱身特性,在軍用飛機得到了廣泛的應用,同時飛翼布局又存在操縱效能低,配平損失等典型問題[18]。所以無尾飛翼布局的氣動外形設計是涉及總體、氣動、控制、結構、隱身等專業約束的綜合設計問題,即多目標多約束的設計問題。一般來說,要求飛機在巡航狀態時具有較高升阻比,同時具有較小的低頭力矩,從而不會引起較大的配平損失,而這兩種要求往往是矛盾的,尤其在飛翼布局上表現的較為突出,即提高巡航升阻比的同時會帶來較大的低頭力矩。為了解決飛翼布局設計中的這種問題,本文基于伴隨方法的基本原理,發展了一種考慮氣動、結構約束的飛翼布局多目標多約束氣動優化方法,并進行了典型算例驗證。
采用結構化貼體網格進行空間數值離散,并按照求和約定,在計算域中,N-S方程可表示為:

(1)

氣動優化問題是以外形變化對氣動特性的影響為基礎而進行的,氣動特性的目標函數可表述為:

(2)
dBξ,dDξ分別為計算空間中的表面與空間積分單元,M與P取決于流動變量w及計算空間的矩陣S。
在滿足流場控制方程約束條件下,氣動外形的變化將導致流動變量變分δw與矩陣變分δS,目標函數的變化可表示為

(3)
式中:δM=[Mw]Ⅰδw+δMⅡ,δP=[Pw]Ⅰδw+δPⅡ;下標Ⅰ表示由流動變量變化δw引起的貢獻;下標Ⅱ表示由矩陣變化δS引起的貢獻。
定常狀態下,氣動外形變化的約束方程可表述為

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

(5)
假定ψ可微,(5)式分部積分及高斯定理可進一步寫成
(6)
整合(3)、(5)、(6)式,目標函數變分可寫為
(7)
令(7)式中空間積分項流動物理量變分δw系數項組合為0,則可得伴隨方程

(8)
令(7)式邊界積分中流動物理量變分δw的系數項組合在一起為0,則得對應的伴隨方程邊界條件
(9)
剩余項即為目標函數的梯度求解公式
(10)
鑒于伴隨方程理論推導的復雜性,詳細過程可參考文獻[12],此處僅給出最終的數學表達式。

(11)

式中:矢量Y中a為聲速;Pr為普朗特數。
為了實現飛翼布局設計點有效減阻優化,綜合考慮以下設計目標。
通過加權組合方法定義如(12)式的統一目標函數
(12)

定義自由來流馬赫數M∞,壓力P∞,迎角α以及參考面積Sref,力矩參考點(xref,yref),易知
則升力系數、阻力系數、俯仰力矩系數中流動變量的變分為
由(12)式得,目標函數中流動變量的變分為
(19)
根據(9)式在參考文獻[12]中伴隨方程邊界條件的推導可知

(20)

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

(21)
同理,根據(10)式目標函數及流動控制方程中對矩陣的變分項,即可獲得梯度求解式為
(22)
圖1給出了優化設計流程圖。

圖1 優化設計流程圖
網格生成:采用無限插值法生成結構化計算網格,同時采用正交控制、加權平均光順和法向量控制等措施確保網格質量。
流場數值求解:采用Jameson的中心格式有限體積法進行空間離散,五步Runge-Kutta顯示格式時間推進,同時加入人工黏性抑制振蕩,采用當地時間步長、隱式殘值光順、多重網格等加速收斂措施;湍流模型采用B-L湍流模型。
伴隨方程數值求解:采用與N-S方程類似的數值解法。
設計變量及優化算法:采用Hicks-Henne形狀函數描述設計變量對物體外形變化的影響,同時采用最速下降法進行梯度搜索。
設計點:Ma=0.85,α=3°,控制剖面取機翼3個剖面,每個剖面26個設計變量,共78個設計變量。計算網格采用C-H網格,如圖2所示。先后開展了2種不同約束下的優化設計。

圖2 C-H計算網格示意圖
設計1:升力、面積約束下的減阻優化設計,根據設計狀態,限定升力系數、各剖面面積相對約束值變化不超過5%。選取目標函數中各部分的權值分別為:Ω1=120,Ω2=1,Ω3=0,Ω4=1。
表1給出了初始與設計外形的氣動特性與幾何特性具體數值對比。從中看到,初始阻力系數為0.013 1,升力系數為0.213 7,優化設計25步后,阻力系數變為0.010 6,升力系數變為0.214 3,阻力系數下降19.5%,而升力系數、各控制剖面面積滿足約束條件。

表1 某小展弦比飛翼不同約束條件優化前后特性對比
注意到相比初始外形,優化后帶來了較大的低頭力矩,這樣就會帶來較大的配平損失,從而導致實際使用升阻比降低。
設計2:為改善設計1帶來的問題,進行了升力、俯仰力矩、面積共同約束下的減阻優化設計,對應目標函數中各部分的權值分別取:Ω1=120,Ω2=1,Ω3=0.01,Ω4=1,保證|Cm|≤0.004。表1為優化后具體的數值變化。迭代優化12步后,阻力系數減小為0.010 8(與設計1保持相同量級),下降約18.0%;升力系數變為0.206 7,減小3.3%,變化相對較大,但也滿足約束條件;俯仰力矩系數由初始的0.001 8變為-0.001 6,滿足約束條件;各控制剖面面積變化同樣滿足約束指標。
圖4給出了初始外形與2種約束下優化設計外形機翼展向不同站位剖面壓力分布對比。可以清楚地看到2種設計外形機翼上表面的激波都被很大程度削弱,不同之處在于有力矩約束得到的外形壓力分布在50%弦長前負壓值較大,50%弦長后負壓值則較小,則對應上表面的壓心靠前,因此不會帶來很大的低頭力矩,與氣動力計算結果一致。

圖3 初始外形與無/有力矩約束優化設計外形表面壓力對比
通過算例1的2種優化設計結果對比,驗證了所定義的目標函數,推導的伴隨方程邊界條件及梯度求解公式是正確并有效的。
采用本文所發展的方法,對某大展弦比飛翼進行跨聲速狀態升力、俯仰力矩和面積約束下的減阻優化設計。設計點:Ma=0.75,α=4°,目標函數中各部分的權值分別取:Ω1=50,Ω2=2,Ω3=0.001,Ω4=0.5,限定升力系數、各剖面面積相對約束值變化不超過5%,俯仰力矩-0.004≤Cm≤0.008;控制剖面取機翼的4個剖面,每個剖面26個設計變量,加上4個剖面扭轉角,共108個設計變量。
表2給出了優化前后氣動力系數及控制剖面面積的具體數值變化。優化迭代8步,阻力系數由初始的0.016 65減小為0.015 06,下降約9.55%;升力系數由0.361變為0.355,減小1.66%,滿足約束條件;俯仰力矩系數由初始的0.006 2變為0.003 6,滿足約束條件;各控制剖面面積變化也滿足約束指標。

表2 某大展弦比飛翼優化前后氣動、幾何特性對比
圖4給出了機翼展向不同站位剖面壓力分布及外形對比。可以看到設計外形上表面壓力負壓峰值區域減小,逆壓梯度變小,展向不同位置的激波強度都有不同程度的減弱,尤其在展向60%~70%范圍激波削弱明顯。對應剖面外形的主要變化趨勢是最大厚度略有減小,且弦向位置略有后移;扭轉角主要是靠近對稱面的剖面有了一個較小的正扭轉角,其余剖面變化不大。

圖4 初始外形與優化外形展向不同剖面外形及壓力對比
本文針對飛翼布局氣動優化設計中的多目標多約束問題,基于伴隨方法和N-S方程發展了一種氣動優化設計方法,并先后進行了2種不同展弦比飛翼布局的跨聲速減阻優化設計,結果表明:
1) 通過構建合理的統一目標函數形式來解決氣動優化設計中的多目標多約束問題是合適的,根據伴隨方法基本原理推導的伴隨方程物面邊界條件以及梯度求解方程是正確、有效的;
2) 所發展的設計方法在飛翼布局的多約束氣動優化設計問題上,具有良好的設計效果和優化效率,因此在工程上具有廣闊的應用前景。