基于伴隨算子的大飛機氣動布局精細優化設計

吳文華，范召林，陳德華，覃 寧，孟德虹

（1.中國空氣動力研究與發展中心空氣動力學國家重點實驗室，四川 綿陽 621000；2.中國空氣動力研究與發展中心高速空氣動力研究所，四川 綿陽 621000；3.謝菲爾德大學，英國 謝菲爾德S3 7JJ）

0 引言

在飛行器氣動布局設計的后期，布局的主要特征參數和外形都已經確定，比如機身的長度、圓柱段直徑、機翼的展弦比、前緣后掠角、根梢比、截面最大厚度，機翼面積，尾翼位置及面積等。在這個時期，這些參數都已經成為其他學科設計的依據，主要氣動性能也不宜有大的變動，以免全部設計工作推倒重來，因為這將帶來巨大的經濟損失并大大拖延進度。這時候的優化設計，可以對布局進行一些局部的、細節的調整，在不對其他學科設計造成較大的影響的前提下，進一步提高布局氣動性能，比如減阻，降低巡航力矩等。調整布局的局部曲面形狀，可以滿足這個要求。對布局的局部曲面形狀進行微調必然是多參數的以達到設計目的，只有設計參數達到足夠的數量，才能夠對曲面進行足夠精細的調整。這種布局曲面的細致調整也將使得曲面微調引發的布局性能變化很小，比如一次調整導致的阻力變化有可能在1阻力單位以內。然而1個阻力單位，對于大飛機等巡航距離很遠的飛行器仍然具有很大的意義。在優化過程中分辨這么小的阻力變化量，是對優化軟件中流場解算數的極高要求，必須在網格密度、網格拓撲結構、差分格式等方面進行精心的設計、試驗和調整才能做到。同時，由于參數多，使得曲面形狀的變化多種多樣，也會導致目標函數與設計參數之間形成復雜的函數關系，這種關系極有可能是非凸的、多極值的，這就使得傳統的牛頓法、二次規劃尋優算法效果變差，甚至無法使用。太多的參數又使得遺傳算法、響應曲面算法、粒子群算法等全局尋優算法無法直接使用，需要探索新的，能夠適應氣動布局設計后期多參數精細優化設計需要的新型尋優算法。

由于數值計算非常耗時，在數值優化技術發展的最初階段，設計參數的數量通常很少，這是因為不論采取哪種優化算法，優化設計的計算量都隨著設計參數的增多而急劇增長，使得參數太多的優化問題根本無法完成。這一狀況直到基于伴隨算子的敏感導數解算方法出現才得到改觀。敏感導數是指優化目標對設計參數的導數，可用于指導設計參數往哪個方向調整，以迅速獲得最優結果。這種方法可以一次性求出所有敏感導數，計算精度高，計算時間隨著設計參數的增加不會明顯增多。伴隨算子求導技術結合基于敏感導數的尋優算法，就可以完成極多參數的氣動布局優化設計，有可能快速并準確地獲得限定條件內的最優布局，而且該方法的計算量隨設計參數數量的增多變化不大。

伴隨算子求導技術在偏微分方程敏感導數求解中的應用至今為止已經超過30年。伴隨算子的應用最早出現在控制優化中，隨后在結構有限元優化中也有應用。在空氣動力學領域中最先應用伴隨算子技術的是Pironneau［1］，之后美國的Jameson［2－4］將其應用到機翼的優化設計中。隨后美國和西歐對基于伴隨算子求導算法的氣動布局優化技術進行了大量的研究［2－8］，將伴隨算子技術在氣動優化設計中的應用大大拓展，比如飛行器全機氣動布局的優化設計，最新的波音787、A380等飛機的研制中都應用了這項技術。

國內的氣動布局優化設計研究，目前主要還集中在布局設計前期的少量主參數優化設計上，用于后期多參數高精度優化設計的研究很少。氣動力的計算常用近似模型，無法滿足精細優化設計的需要。

中國空氣動力研究與發展中心總體技術部［9－10］、南京航空航天大學［11－12］和西工大［13－17］在基于伴隨算子的多參數優化設計方面開展了一些研究，但是研究應用水平與西方發達國家相比還有差距，開展的多是翼型或者單獨機翼的優化設計研究［11－17］。單獨的機翼優化設計沒能夠考慮發動機或者機身對布局的影響，在布局設計后期的應用會受到限制。

目前的大型飛機布局大都采用機翼加柱形機身再加垂尾、平尾的布局方式。波音公司、空客公司的客機、運輸機，從最早的波音737，到最新的A380、波音787，都采用了這種布局。經過了幾十年的研究，這種布局的性能潛力已經得到了充分的挖掘，經過傳統的設計手段的優化，其氣動性能在給定設計條件下已經達到了相當高的水準。即使采用新的設計手段，受到布局形式的限制，其氣動性能的提高也將很有限。國內外基于伴隨算子的多參數優化設計算例，優化對象的初始氣動特性大都較差，優化容易獲得較大的性能提高。本文的優化對象已經經過傳統設計手段的多輪優化，初始外形就具有很高的氣動性能，對這種布局的優化難度大大提高。

本文以伴隨算子優化方法為基礎，發展了一套多參數、高精度優化設計軟件，用于布局設計后期精細優化設計。該軟件由雷諾平均N－S方程解算器、伴隨方程解算器、動網格程序、二次規劃尋優程序、外形參數化程序等構成。該軟件系統在一種大型飛機布局全機構型的高精度優化設計中得到應用，優化過程中計入了短艙和機身等對機翼氣動特性的影響，優化獲得了良好的效果。

1 數值方法

常規的氣動布局氣動特性計算要求計算結果具有很高的精準度，雖然氣動布局優化設計也要求計算結果準確，但是更強調優化過程中，不同布局外形氣動特性——優化目標之間差量的準確性。如果計算過程中，解的振蕩性太大，那么在解滿足收斂條件時，有可能某個解在波峰，而另外一個解處在波谷，那么它們的差量就包含了較大的數值誤差，因此本研究需要的解算器不僅要速度快，還要穩定性高，能很好地抑制解的振蕩，以便數值誤差互相抵消，獲得高精度的差量。

用于優化設計的解算器還要具有很高的健壯性，即要求在求解過程中，最大限度地避免出現解不收斂等導致計算過程中斷的情況出現，因為優化是完全自動進行的，期間要進行時間很長的迭代并且完成大量網格不同的數值計算，如果在其中的某一步出現問題，那么整個優化過程就會中斷，或者導致錯誤的優化結果產生。

1.1 流場解算方法

根據優化設計平臺對流場解算器的技術要求，在借鑒和改進國內外研究的基礎上，采用基于黎曼近似基本解的Osher矢通量分裂格式和有限體積法方法來構造解算器［7］，空間離散采用MUSCL 格式，時間離散采用隱式差分格式。這種方法具有較高的精度和良好的穩定性，同時能夠有效抑制解的振蕩，能夠滿足優化計算對數值解算器的幾個要求。

任意控制體 上N－S方程組的積分形式為：

其中?Ω表示控制體邊界，F為N－S方程組除了時間項以外各計算項的矩陣形式［7］：

采用有限體積法將上述方程離散，意味著整個計算區域被分割成很多小的控制體，第i個控制體Vi上N－S方程組的形式為：

Ui為第i個控制體上狀態變量的平均值，Ri為流經控制體表面的通量總和的殘差矢量，

本計算程序是利用時間推進法求解穩態問題，對于第n個時間層，有：

采用當地時間步長，以加速收斂。

對流項的處理采用近似黎曼解的Osher格式。Osher格式對間斷問題的處理表現出色。由于采用有限體積法時，整個流場被劃分為許多小控制體，每個控制體和周圍相鄰的某個控制體共享一個界面，通過計算每個界面上的通量來得到相鄰控制體之間的影響，一旦每個控制體的變化得到，整個流場的變化也就知道了。

湍流模擬采用k－ωSST 湍流模型。

1.2 敏感導數解算方法

敏感導數是指優化目標對設計參數的導數，可用于指導設計參數往哪個方向調整，以獲得更好的結果。本文采用基于敏感導數的尋優算法，因此如何快速準確地獲得敏感導數，成為本研究的關鍵所在。

最直接的敏感導數計算方法是有限差分法，可以分為單邊差分和中心差分兩種方式，單邊差分格式具有一階精度，中心差分格式具有二階精度。有限差分法簡單，但是計算量太大，如果有n個設計參數，采用單邊差分格式需要進行n＋1次數值計算，中心差分格式需要進行2n次數值計算，才能得到全部敏感導數，如果有1000個設計參數，那么優化的第一步就至少需要進行1001次數值計算，因此這種方法只是在設計參數很少時才能應用。

伴隨算子敏感導數求解方法可以一次性求出所有敏感導數，計算精度高，計算時間隨著設計參數的增加不會明顯增多。該方法是通過求解伴隨方程，獲得伴隨算子，然后求得網格參數對設計參數的導數，目標函數對設計參數的偏導數、目標函數對網格的偏導數和流場殘差對設計參數的偏導數等，再使用這些計算結果進行一些代數運算獲得敏感導數。這種算法的好處在于可一次性獲得所有的敏感導數，而且精度很高，有效避免了有限差分法的缺點。由于伴隨算子只有在氣動布局外形改變之后才需要重新計算，所以每一個優化步，只需要求解一次伴隨方程、一次NS方程，使得計算效率大大提高，對于優化變量特別多的優化問題尤其有效，比如對于設計參數為1000的優化計算，只需要進行一次N－S方程求解，一次伴隨方程求解，伴隨方程求解的復雜程度和耗時與N－S方程基本類似，然后再進行一些代數計算，就可以獲得全部的敏感導數。

基于粘性方程的敏感導數計算公式推導如下。

目標函數可以表示為：

其中各項意義如下：f＝f（Q＊（β），X＊（β），β）為流場變量，上標＊表示為收斂的流場變量；X為網格變量組成的矢量；β為設計變量矢量。

上式的差分表達形式為：

采用伴隨方程求解時，加入伴隨向量λ，表達式變為：

為了避免求解dQ＊／dβk時需要多次求解流場，我們令：

所以在λ求解得到伴隨矢量以后，敏感性導數可以由以下公式計算得到：

伴隨方程的求解分為連續法和離散法，兩種方法的結果一致，本項研究中伴隨方法的求解采用離散法。首先需要求解伴隨方程，求得伴隨向量λ以后可以計算敏感導數。

2 參數化外形設計及網格變形方法

本文采用了Bezier－Bernstein外形參數化方法。該方法的優點是能夠以比較少的參數，比較精確地表示布局外形，同時又能夠對外形進行有效、精確而又細致地調整，還能保持布局表面的光滑性。這一點對于大型飛機氣動布局的參數化具有特別重要的意義。主要設計思想是將在每條Bezier－Bernstein曲線上設置幾個到幾十個控制點，這些點的位置（坐標）就成為設計參數。曲線的數量和每條曲線點的數量取決于布局外形曲面的復雜程度以及優化所要求的精細程度，一般光滑的表面外形，設置十幾個控制點就夠了。

對于二維的曲線，Bezier可以采用下式表示：

其中，S2（u）＝x（u）／y（u），PK＝Px／Py是Bezier曲線的控制點，Bernstein多項式BK，N（u）＝uK（1－u）N－KN！／［K！（N－K）！］中，u表示曲線的參考弧長，N表示控制點個數，也就數參數數量，Px、Py表示控制點的縱橫向坐標。以機翼為例，由于機翼通常在優化過程中，只需要改變控制點的Y坐標，因此設計變量是控制點的Y坐標Py。

由于優化設計中超臨界翼型的外形變化不會太大，而且本研究主要采用結構網格，因此采用代數法進行網格變形設計。代數法網格變形技術能夠基本滿足優化的需求，而且能夠比較方便地計算網格點對設計變量的導數，這些導數在使用伴隨算子計算敏感導數時需要用到。該方法先移動布局表面的網格點到新位置，然后將這個變動逐漸傳遞到外圍，在傳遞過程中，根據點的位置按比例調整節點位移量以保證計算域的外邊界保持不變。這樣還可以保證網格的拓撲結構和附面層網格等非常相似，從而抑制網格變化產生的數值誤差，提高優化結果的精度。

3 尋優算法

采用了二次規劃尋優算法。優化問題可以描述為：

最小化目標函數：F（β）；約束條件：gi（β）≤0，i＝1，l；hj（β）＝0，j＝1，m；βlk≤βk≤βnk，k＝1，NDV。其中，β＝（β1，β2，β3，β4……）T表示設計變量。在求解過程中，二次規劃尋優算法通過目標函數值及其敏感導數，用二次曲面在優化起始點附近擬合目標函數。給出能使得目標函數最優的設計參數變化方向和大小。計算獲得新的目標函數值，并判斷新的設計點是否滿足約束條件，是否優于老的設計點。重復這個循環，直到滿足優化收斂條件。這種尋優算法只能得到局部最優點，本文針對大飛機機翼的特點，對這種算法進行了一些改進，提高了優化效果。

4 W9大飛機布局的多參數優化設計

4.1 基本外形氣動特性

W9布局基本外形如圖1所示，機翼表面的壓力分布如圖2所示。

原始外形全機巡航升力系數CL＝0.5，升阻比K＝17.53，阻力系數CD＝0.02853。

圖1 W9布局外形及表面網格分布Fig.1 The original shape and surface grids of the W9

4.2 設計條件

以給定的初始大型飛機布局為基礎，在全機狀態下對機翼進行多參數的減阻優化，這樣所得到的優化結果就計入了短艙、機身、掛架等部件對機翼的影響。機翼初始外形及網格如圖1。設計目標與設計要求如下：

設計要求：CL≥0.5，Ma＝0.785，Re＝2.4×107，保持機翼最大相對厚度不減小；設計目標：減阻；

設計參數：設計參數168個，短艙掛架左右的機翼各使用5 個截面控制，每個截面上下表面各8 個Bezier控制參數，前后緣點不動，短艙上方使用8個設計參數控制，總計168個設計參數，控制面位置及原形壓力分布見圖2。

圖2 機翼展向控制面的位置Fig.2 The position distribution of wing span control section

4.3 優化結果與分析

優化前：CL＝0.5，CD＝0.02853，K＝17.53；優化后：CL＝0.5，CD＝0.02785，K＝17.95。優化總共進行了41步，阻力變化過程如圖3所示。

整個優化設計共計進行了41步，最大的一步阻力減小也不到2個阻力單位，其中很多優化步的阻力減小量小于0.2個阻力單位。那么，這么小的阻力優化量，其計算結果是否可靠呢？優化過程本身就能說明這個問題，雖然每一步的優化量不大，但是整個優化過程進行的很順利，每一步都較前一步阻力有所減小，雖然有些優化步的減小量很小。如果這些阻力的變化僅僅是數值誤差的話，那么受到數值誤差干擾的優化過程就無法順利的進行下去，中間必然會出現阻力增大的情況。由于網格的數量和質量較高，同時阻力收斂較好，優化過程中又保證了網格的相似性，因此阻力變化量的計算精度很高，能夠滿足大型飛機精細優化設計的計算精度要求。

圖3 優化過程中阻力的變化Fig.3 Variation of drag coefficient with the optimization process

相比優化前，阻力減小了6.8個阻力單位，占全機總阻力的2.38%。這個比重雖然比較小，但是考慮到這種優化只能減小機翼的壓差阻力，并不能減小摩擦阻力，而機翼的總阻力為0.0144，因此阻力減小量約為機翼阻力的4.72%，占壓差阻力的比重就更大，接近10%。為了驗證優化結果的可靠性，采用了相同的網格和湍流模型，使用空氣動力學國家重點實驗室的亞跨超平臺Trip 2.0進行了驗算。驗算得到的升力值和阻力值均與原值有一定的差別，阻力系數相差約4個阻力單位，但是阻力優化量卻變化很小，驗算得到的阻力優化量約為6.6個阻力單位，減小了約0.2個阻力單位。從驗算結果可以看出，本次優化的精度是很高的，結果是可信的。

圖4給出了優化前后機翼的表面壓力分布對比，從圖中可以看出，優化后機翼上表面的激波明顯減弱，消除了原型機翼上外翼部位存在的二次壓縮與膨脹。圖5到圖8給出了機翼展向30%到60%共計2個截面的壓力分布和形狀對比圖。從圖中可以看出，優化后的機翼壓力分布更加接近標準的超臨界壓力分布，激波更弱，部分截面激波消除。原形壓力分布中的凹坑被抹平。下表面壓力分布基本未變，整個機翼的截面外形改變較小，截面最大厚度基本保持不變。其他截面的壓力分布變化與這兩個截面類似，由于篇幅所限，不再一一列出。

圖4 優化前后壓力分布對比圖（168設計變量）Fig.4 The comparison of surface pressure distribution of original shape and optimized one（168design variables）

圖5 機翼展向30%截面處優化前后壓力分布對比Fig.5 Pressure distribution comparison between original and optimized at 30% wing span section for W9

圖6 機翼展向30%處翼型優化前后外形對比圖Fig.6 The comparison of airfoil shape of original and optimized at 30% wing span section for W9

圖7 機翼展向60%截面處優化前后壓力分布對比Fig.7 Pressure distribution comparison between original and optimized at 60% wing span section for W9

圖8 機翼展向60%處翼型優化前后外形對比圖Fig.8 The comparison of airfoil shape of original and optimized at 60% wing span section for W9

5 結論

從優化結果可以看出本文采用的優化方法是有效的，所采用的數值方法、伴隨算子解算方法、尋優方法能夠滿足大型飛機全機構型多參數高精度優化設計的要求，優化結果與理論分析相符。雖然原始外形已經經過常規設計手段的多輪優化，其性能已經很高，但是多參數精細優化仍然取得了明顯的優化效果，在保持升力不變，最大厚度不減小的前提下阻力降低了0.00068，減小了2.38%，其中機翼阻力減小了4.72%。從截面壓力分布來看，優化后的截面壓力分布明顯改善。這與總阻力的減小能夠互相印證。

從優化過程和驗證結果來看，本次優化設計對阻力變化量的計算精度較高，能夠滿足大飛機氣動布局優化精細設計的要求。

［1］IRONNEAU O.On optimum design in fluid mechanics［J］.JournalofFluidMechanics，1974，64（1）：97－110.

［2］JAMESON A.Aerodynamic design via control theory［J］.JournalofScientificComputing，1989，3（3）：233－260.

［3］JAMESON A.Optimum aerodynamic design using CFD and control theory［R］.AIAA 95－1729，1995.

［4］JAMESON A.Automatic design of transonic airfoils to reduce the shock induced pressure drag［A］／／Proceedings of the 31st Israel Annual Confer－ence on Aviation and Aeronautics［C］.Tel.Aviv.，1990.

［5］MOIGNE A L，QIN N.Variable－fidelity aerodynamic optimization for turbulent flows using a discrete adjoint formulation［J］.AIAAJournal，2004，42（7）：1281－1292.

［6］QIN N，WONG W S，MOIGNE A L.Three－dimensional contour bumps for transonic wing drag reduction［J］.ProceedingsoftheInstitutionofMechanicalEngineers，PartG：JournalofAerospaceEngineering，2008，222（5）：619－629.

［7］MOIGNE A L.A discrete Navier－Stokes adjoint method for aerodynamic optimization of blended wing－body configurations［D］.［PhD Thesis］.Cranfield，Cranfield University，2002.

［8］GILES M B，PIERCE N A.An introduction to the adjoint approach to design［J］.Flow，Turbulenceand Combustion，2000，65（3）：393－415.

［9］黃勇，陳作斌，劉剛.基于伴隨方程的翼型數值優化設計方法研究［J］.空氣動力學學報，1999，17（4）：413－422.

［10］周鑄，陳作斌.基于N－S方程的翼型氣動優化設計方法研究［J］.空氣動力學學報，2002，20（2）：141－149.

［11］唐智禮，黃明恪.基于控制理論的Euler方程翼型減阻優化設計［J］.空氣動力學學報，2001，19（3）：262－270.

［12］唐智禮，黃明恪.約束最優控制理論及其在氣動優化中的應用［J］.力學學報，2007，39（2）：273－277.

［13］喬志德，楊旭東，朱兵.亞、跨音速三維機翼氣動外形反設計的控制理論方法［J］.2003，21（1）：11－19.

［14］楊旭東，喬志德.基于共軛方程法的跨音速機翼氣動力優化設計［J］.航空學報，2003，24（1）：1－5.

［15］李少峰，喬志德，楊旭東.基于控制理論方法的翼型最大升力優化設計［J］.航空計算技術，2005，35（4）：98－102.

［16］楊旭東，喬志德，朱兵.基于控制理論和NS 方程的氣動優化設計方法研究［J］.空氣動力學學報，2005，23（1）：46－51.