基于共軛方程法的跨聲速機(jī)翼阻力優(yōu)化設(shè)計(jì)

張野平， 段卓毅， 張彥軍

（西安飛機(jī)設(shè)計(jì)研究所，陜西 西安710089）

0 引 言

在飛行器設(shè)計(jì)中，翼型、機(jī)翼、全機(jī)的氣動(dòng)設(shè)計(jì)技術(shù)歷來(lái)被視為各國(guó)航空領(lǐng)域的核心技術(shù)，傳統(tǒng)的設(shè)計(jì)方法主要有間接方法、余量修正方法幾種類型。間接設(shè)計(jì)法中典型的有虛擬氣體法和速度圖法。余量修正方法中具有代表性的方法是Takannashi在19世紀(jì)80年代中期提出的一種余量修正方法［1］。這兩類方法都計(jì)算量相對(duì)較小，但是自身存在不能直接控制氣動(dòng)特性或者解的適定性等諸多問(wèn)題。

20世紀(jì)70年代初至今，隨著計(jì)算機(jī)科學(xué)和數(shù)值方法的迅猛發(fā)展，計(jì)算流體力學(xué)（CFD）和數(shù)值優(yōu)化技術(shù)取得了長(zhǎng)足的進(jìn)步，基于CFD的氣動(dòng)外形優(yōu)化設(shè)計(jì)技術(shù)應(yīng)運(yùn)而生。根據(jù)目前的計(jì)算能力而言，以遺傳算法為代表的全局尋優(yōu)方法由于搜索面廣、計(jì)算量大，還只能用在比較簡(jiǎn)單的模型上；能夠應(yīng)用于工程實(shí)際的只有基于梯度計(jì)算的局部尋優(yōu)方法。在基于梯度的氣動(dòng)優(yōu)化設(shè)計(jì)方法中，對(duì)于有N個(gè)設(shè)計(jì)變量的設(shè)計(jì)問(wèn)題，往往進(jìn)行一次優(yōu)化設(shè)計(jì)就要進(jìn)行N＋1次流場(chǎng)求解，為獲得梯度所需進(jìn)行的流場(chǎng)計(jì)算次數(shù)與設(shè)計(jì)變量個(gè)數(shù)直接相關(guān)，當(dāng)設(shè)計(jì)變量眾多時(shí)計(jì)算代價(jià)難以接受。

20世紀(jì)80年代末期，A.Jameson提出的基于控制理論（共軛方程）的氣動(dòng)優(yōu)化設(shè)計(jì)技術(shù)［2－4］成功地解決了基于梯度的優(yōu)化方法中梯度獲得困難、計(jì)算量巨大這一難題。該方法以偏微分方程系統(tǒng)的控制理論為基礎(chǔ)，把物體邊界形狀作為控制函數(shù)，流場(chǎng)控制方程作為約束條件，在目標(biāo)函數(shù)中引入Lagrange因子，變條件約束問(wèn)題為無(wú)約束問(wèn)題，通過(guò)求解流場(chǎng)控制方程及其共軛方程進(jìn)行梯度求解，其計(jì)算量只相當(dāng)于兩倍的流場(chǎng)計(jì)算量，與設(shè)計(jì)變量數(shù)目無(wú)關(guān)，從而將設(shè)計(jì)問(wèn)題轉(zhuǎn)變成控制問(wèn)題。目前該設(shè)計(jì)理論的應(yīng)用領(lǐng)域不斷拓寬，對(duì)該方法的深入研究具有重要的意義和工程價(jià)值。

1 流動(dòng)控制方程

控制方程使用NS方程。在三維笛卡爾坐標(biāo)系（x1，x2，x3）中，有速度分量為（u1，u2，u3），守恒形式的NS方程形式為：

其中w是狀態(tài)矢量，fi是無(wú)粘矢通量，fvi是粘性矢通量。

在用結(jié)構(gòu)網(wǎng)格進(jìn)行求解時(shí)，需要將上述NS方程轉(zhuǎn)化為計(jì)算坐標(biāo)系（ξ1，ξ2，ξ3）下，其守恒形式轉(zhuǎn)化為：

其中F＝Sf，F(xiàn)＝Sf，S＝J，J為坐標(biāo)變換iijjviijvjij的Jacobin值。Sij的物理意義是包含控制體的邊界的外法向面矢量，J為控制體體積的倒數(shù)。

2 伴隨方程

引入目標(biāo)函數(shù)表示為：

各變量的含義如下：

w：流場(chǎng)參數(shù)

S：變換矩陣

B：物面邊界

D：計(jì)算體積空間

M：物面邊界變化對(duì)泛函的貢獻(xiàn)

P：計(jì)算區(qū)域體積空間變化對(duì)泛函的貢獻(xiàn)

很顯然，就空間劃分來(lái)說(shuō)，目標(biāo)函數(shù)受到的影響來(lái)自于物面邊界的變化和計(jì)算空間單元體積變化兩個(gè)部分。就物理上來(lái)說(shuō)，泛函的值取決于流場(chǎng)參數(shù)w和變換矩陣S。

物面邊界的變化最終將導(dǎo)致目標(biāo)函數(shù)的變化。則可以物面邊界作為控制函數(shù)，將流動(dòng)控制方程作為約束條件在目標(biāo)函數(shù)中引入，從而將約束問(wèn)題轉(zhuǎn)化為無(wú)約束問(wèn)題，設(shè)計(jì)問(wèn)題轉(zhuǎn)化為控制問(wèn)題。根據(jù)變分法則，設(shè)目標(biāo)函數(shù)對(duì)于邊界變形的變分為：

下標(biāo)Ⅰ和Ⅱ分別表示由流動(dòng)變量變化δw起的貢獻(xiàn)和由坐標(biāo)變換矩陣變化δS引起的貢獻(xiàn)。為了將流動(dòng)的控制方程作為約束條件引入到目標(biāo)函數(shù)里面，從而將有約束問(wèn)題轉(zhuǎn)化成無(wú)約束問(wèn)題，現(xiàn)在對(duì)定常情況下的流動(dòng)控制方程進(jìn)行變分處理：

在這里引入共軛變量Ψ，并對(duì)上式做全流場(chǎng)積分：＝0，假設(shè)“乘子”Ψ是可微分的，對(duì)上式做分部積分：

利用高斯定理，進(jìn)行體積分到面積分的轉(zhuǎn)化：

其中ni是面矢量。

將控制方程的積分形式的變分引入到目標(biāo)函數(shù)的變分當(dāng)中，用原來(lái)的變分減去現(xiàn)在的積分形式的控制方程的變分，再將δM、δP、δFi、δFvi帶入，整理可得：

選擇合適的共軛變量Ψ，使得δI和流動(dòng)變量δw不再顯性相關(guān)。目標(biāo)函數(shù)的計(jì)算可直接通過(guò)坐標(biāo)變換矩陣的變分得到，而不需要重新計(jì)算由擾動(dòng)每個(gè)設(shè)計(jì)變量產(chǎn)生的流動(dòng)變量變分δw。這樣就可以在全流場(chǎng)區(qū)域D的積分當(dāng)中首先得到全流場(chǎng)粘性伴隨方程：

然后在邊界積分區(qū)域B當(dāng)中得到粘性伴隨方程的物面邊界條件：

可得

對(duì)于本文涉及的正向設(shè)計(jì)問(wèn)題，目標(biāo)函數(shù)一般表達(dá)為：

式中CD為壓阻系數(shù)，CL為設(shè)計(jì)升力系數(shù)。

3 梯度求解

下面通過(guò)簡(jiǎn)單的變分計(jì)算問(wèn)題說(shuō)明改進(jìn)過(guò)的Sobolev梯度求解方法［5］，使用該方法能夠生成一系列連續(xù)的光滑外形，避免優(yōu)化所得的外形連續(xù)性降階。

選擇y（x）使得I（y，y′）dx最小，y（x）在兩端點(diǎn)的值分別為y（a）、y（b）。

使y（x）產(chǎn)生一個(gè)小的變量δy（x），可得：

定義梯度形式為：

內(nèi)積形式為

可以看出

假設(shè)：δy＝－λg，λ＞0

由此可得：

當(dāng)且僅當(dāng)g＝0時(shí)，等號(hào)成立，即取最小值。注意到g是y、y′、y″的函數(shù)，即g＝g（y，y′，y″）。

以著名的Brachistrone（最速降線）問(wèn)題為例：

并且：

每一步y(tǒng)n＋1＝y(tǒng)n－λngn都減小y的兩階光滑度。因此計(jì)算會(huì)導(dǎo)致越來(lái)越不光順，將導(dǎo)致計(jì)算的不穩(wěn)定性。

為了避免上述不利情況的發(fā)生，引入一種改進(jìn)過(guò)的Sobolev內(nèi)積：

這里ε是控制導(dǎo)數(shù)項(xiàng)權(quán)重的參數(shù)。如果定義一種梯度能夠使得：

可以得到：

其中

則有每一步y(tǒng)n＋1＝y(tǒng)n－λnn。

這樣就給出了一種改進(jìn)以后δI的形式：

這種情況下yn＋1和yn具有一致的連續(xù)性，可以得到穩(wěn)定的解。

使用這種求解梯度的方法，就可以只用網(wǎng)格點(diǎn)對(duì)模型進(jìn)行參數(shù)化，而不必使用其它復(fù)雜的模型參數(shù)化方法。一般來(lái)說(shuō)，使用每一個(gè)網(wǎng)格點(diǎn)來(lái)對(duì)模型進(jìn)行參數(shù)化，即把每一個(gè)網(wǎng)格點(diǎn)當(dāng)作設(shè)計(jì)變量，那么在三維機(jī)翼的情況下，設(shè)計(jì)變量將多達(dá)數(shù)千個(gè)，對(duì)于這么多的設(shè)計(jì)變量將產(chǎn)生一些問(wèn)題，使用這樣的設(shè)計(jì)空間，會(huì)導(dǎo)致所謂的病態(tài)問(wèn)題。通過(guò)使用上述梯度求解方法，可以光滑控制表面形狀，減弱了設(shè)計(jì)中的高頻震蕩，使病態(tài)問(wèn)題得到了改善。

4 數(shù)值方法

計(jì)算使用剪切拋物變換得到的C－H型網(wǎng)格；湍流模型采用BL模型；在流動(dòng)控制方程和共軛方程的數(shù)值求解中，采用格心格式的有限體積法進(jìn)行空間離散，時(shí)間推進(jìn)采用五步Runge－Kutta方法。

采用了目前常用的當(dāng)?shù)貢r(shí)間步長(zhǎng)、隱式殘值光順和多重網(wǎng)格法加速計(jì)算收斂，同時(shí)加入人工粘性抑制數(shù)值震蕩。

生成計(jì)算網(wǎng)格后直接使用網(wǎng)格點(diǎn)進(jìn)行幾何的參數(shù)化描述，每個(gè)網(wǎng)格點(diǎn)視作一個(gè)設(shè)計(jì)變量，通過(guò)文中所述的梯度求解方法可以避免優(yōu)化所得外形的光滑性降階；每完成一次優(yōu)化計(jì)算，重新生成網(wǎng)格進(jìn)行計(jì)算。

在求得梯度之后，使用最速下降法進(jìn)行優(yōu)化。

為保證設(shè)計(jì)所得的機(jī)翼在較大的范圍內(nèi)都具有良好的氣動(dòng)特性，將優(yōu)化方法和Pareto陣面的概念結(jié)合，進(jìn)行多目標(biāo)氣動(dòng)優(yōu)化［6］。

5 設(shè)計(jì)算例

優(yōu)化設(shè)計(jì)構(gòu)型為某飛機(jī)翼身組合體構(gòu)型，僅對(duì)機(jī)翼部分進(jìn)行優(yōu)化設(shè)計(jì)。

為保證設(shè)計(jì)所得的機(jī)翼在較大的范圍內(nèi)都具有良好的氣動(dòng)特性，提出四個(gè)設(shè)計(jì)狀態(tài)，如表1所示。

表1 設(shè)計(jì)狀態(tài)Table 1 Design points

第一個(gè)設(shè)計(jì)點(diǎn)用以對(duì)爬升狀態(tài)點(diǎn)進(jìn)行控制；第二、三個(gè)設(shè)計(jì)點(diǎn)用以對(duì)機(jī)翼巡航狀態(tài)進(jìn)行減阻設(shè)計(jì)；第四個(gè)設(shè)計(jì)點(diǎn)用以對(duì)較高馬赫數(shù)進(jìn)行減阻優(yōu)化，獲得較高的阻力發(fā)散馬赫數(shù)。

具體約束為機(jī)翼控制剖面的相對(duì)厚度、機(jī)翼的平面形狀不改變，考慮到機(jī)翼的低速特性與機(jī)翼控制剖面的頭部半徑關(guān)系密切，另外加入對(duì)機(jī)翼控制剖面的頭部半徑的約束。

選取一副經(jīng)過(guò)初步設(shè)計(jì)的機(jī)翼作為初始外形，采用大量實(shí)踐驗(yàn)證的幾何扭轉(zhuǎn)角分布進(jìn)行剖面扭轉(zhuǎn)以保證良好的幾何特性和翼尖失速特性，并且根據(jù)實(shí)際的剖面厚度要求對(duì)基本剖面形狀進(jìn)行厚度調(diào)整。初始機(jī)翼使用的控制剖面具有比較明顯的超臨界翼特征，這樣配置所得的初始機(jī)翼具有較好的跨聲速氣動(dòng)特性。

圖1為翼身組合體表面網(wǎng)格和空間剖面網(wǎng)格示意圖。

圖2為優(yōu)化前后機(jī)翼六個(gè)控制剖面外形對(duì)比。由于加入了剖面厚度、扭轉(zhuǎn)角以及頭部半徑的約束，可以看出，各剖面厚度未出現(xiàn)明顯變化，另外，優(yōu)化后剖面頭部半徑均與初始機(jī)翼持平或者略大于初始機(jī)翼。

圖3、圖4為兩個(gè)典型狀態(tài)下優(yōu)化前后剖面壓力分布對(duì)比，圖5、圖6為多個(gè)狀態(tài)下優(yōu)化前后機(jī)翼表面壓力云圖，從中可以看出，較大馬赫數(shù)下機(jī)翼上表面激波有明顯減弱，這樣可以大大降低波阻，提高阻力發(fā)散馬赫數(shù)。圖7為機(jī)翼優(yōu)化前后CL＝0.55下阻力隨馬赫數(shù)變化曲線，圖8為CL＝0.55下優(yōu)化前后巡航效率因子隨馬赫數(shù)變化曲線。可以看出初始機(jī)翼阻力發(fā)散馬赫數(shù)在0.76～0.77之間，經(jīng)過(guò)15次優(yōu)化后，阻力發(fā)散馬赫數(shù)提升至0.78～0.79之間。巡航效率因子也較初始機(jī)翼有明顯提高。

圖1 表面網(wǎng)格及機(jī)翼展向剖面網(wǎng)格Fig.1 Surface mesh and spanwise section mesh

圖2 機(jī)翼各展向位置剖面幾何外形變化Fig.2 Variation of section shape at different spanwise positions

圖3 Ma＝0.77，CL＝0.53剖面壓力分布對(duì)比Fig.3 Pressure comparison of section，Ma＝0.77，CL＝0.53

圖4 Ma＝0.79，CL＝0.5剖面壓力分布對(duì)比Fig.4 Pressure comparison of section，Ma＝0.79，CL＝0.5

圖5 Ma＝0.77，CL＝0.53機(jī)翼表面壓力云圖對(duì)比Fig.5 Comparison of pressure contours of wing，Ma＝0.77，CL＝0.53

圖6 Ma＝0.79，CL＝0.5機(jī)翼表面壓力云圖對(duì)比Fig.6 Comparison of pressure contours of wing，Ma＝0.79，CL＝0.5

圖7 CL＝0.55下阻力系數(shù)隨馬赫數(shù)變化曲線Fig.7 Variation of drag coefficient with Mach Number，CL＝0.55

圖8 CL＝0.55下巡航效率因子隨馬赫數(shù)變化曲線Fig.8 Variation of cruise efficient with Mach Number，CL＝0.55

表2為四個(gè)設(shè)計(jì)點(diǎn)優(yōu)化前后機(jī)翼主要?dú)鈩?dòng)力指標(biāo)對(duì)比，可以看出各個(gè)設(shè)計(jì)點(diǎn)的阻力均有所減小，升阻比相應(yīng)提高，大馬赫數(shù)下優(yōu)化效果更為明顯，由于優(yōu)化增大了機(jī)翼的后加載，因此機(jī)翼的俯仰力矩系數(shù)有相應(yīng)的增加，四個(gè)設(shè)計(jì)狀態(tài)下平均俯仰力矩系數(shù)增加0.02左右，在巡航狀態(tài)Ma＝0.75～0.77附近，俯仰力矩系數(shù)約為CM＝－0.145附近，量值對(duì)于超臨界機(jī)翼來(lái)說(shuō)屬于較為適中的范圍，工程上可以接受。

表2 優(yōu)化前后氣動(dòng)力對(duì)比Table 2 Comparison of aerodynamic performance before and after optimization

6 結(jié) 論

本文通過(guò)求解三維NS方程及其共軛方程對(duì)跨聲速機(jī)翼進(jìn)行氣動(dòng)優(yōu)化，文中針對(duì)某三維機(jī)翼氣動(dòng)外形優(yōu)化設(shè)計(jì)效果明顯，算例證明本文發(fā)展的基于粘性NS方程共軛方法的跨聲速機(jī)翼優(yōu)化設(shè)計(jì)具有良好的優(yōu)化效果，相比其它優(yōu)化方法，該方法具有優(yōu)化效率高、不受設(shè)計(jì)變量數(shù)目限制的明顯優(yōu)勢(shì)，優(yōu)化結(jié)果滿足設(shè)計(jì)要求，具有很高的工程實(shí)用價(jià)值。

［1］TAKANASHI S.An iterative procedure for three dimensional transonic wing design by the integral equation method［R］.AIAA Paper 84－2155，1985.

［2］JAMESON A.Aerodynamic design via control theory［R］.AGARD－CP－463，1989.

［3］JAMESON A，REUTHER J.Control theory based airfoil design using the euler equations［R］.AIAA Paper 94－4272－CP，1994.

［4］JAMESON A，PRINCE N A，MARTINELLI L.Optimum aerodynamic design using the navier－stokes equations［R］.AIAA Paper 97－0101，1997.

［5］JAMESON A，VASSBERG J C，SRIRAM SHANKARAN.Aerodynamic－structural design studies of lowsweep transonic wings［R］.AIAA Paper 2008－145，2008.

［6］TANG Z L.Pareto front capture using deterministic optimization methods in multi－criterion aerodynamic design［J］.TransactionofNanjingUniversityofAeronautics＆Astronautics，2006，23（2）：81－86.