基于替代模型的三維后體尾噴管優化設計

苗 萌，曾 鵬，閻 超

（1.北京航空航天大學 國家計算流體力學實驗室，北京 100191；2.北京臨近空間飛行器系統工程研究所，北京 100076；3.中國航天電子技術研究院無人機系統工程研究所，北京 100094）

0 引 言

以吸氣式沖壓發動機為動力的飛行器將成為本世紀航空航天技術的發展方向和里程碑。自20世紀50年代以來，美國與前蘇聯就開展了這些方面的基礎研究，80年代以后，歐洲、日本也逐漸加入到研究行列。X－43A、X－51A等試飛成功，更是預示著國外在高超聲速飛行器領域取得了相當大的成果。為了確保我國在該領域達到國際先進水平，如何解決高聲速飛行器氣動設計問題顯得尤為重要。與常規的亞聲速、低超聲速飛行器的設計不同，吸氣式高超聲速飛行器的機體與推進系統必須進行一體化優化設計［1］，而其后體尾噴管是飛行器產生推力、升力的關鍵部件，在一體化設計中占重要地位。

國內外對于吸氣式高超聲速飛行器后體尾噴管優化設計都做過相應的研究。對于二維平面或軸對稱外形的尾噴管常采用特征線法進行設計，但應用于三維外形則過為復雜，且沒有考慮粘性影響。20世紀70年代，Edwards C L Q等［2］將后體簡化為斜平面，以噴管推力系數尾性能目標，完成了尾噴管簡化構型的優化設計；Baysal O，Burgreen G W 等［3］將CFD技術引入尾噴管的優化設計領域，以軸向推力為目標，對二維噴管構型進行優化；羅世彬等［4］人基于二次和三次型線利用響應面和多目標優化方法，實現了二維構型的多目標優化；陳兵等［5］采用空間推進方法求解層流PNS方程，對二維尾噴管進行了優化設計；甘文彪等［6］結合試驗設計方法、替代模型等技術通過求解RANS方程對二維后體尾噴管進行了多點多目標優化設計。賀旭照等［7］采用采用空間推進解算器求解超聲速無粘流動，在二維尾噴管優化設計的基礎上，確定了對稱面的外形，并依此進行三維拓展，實現了三維后體尾噴管優化設計。

為了更加充分地挖掘后體尾噴管流動的三維膨脹性能，本文直接將尾噴管三維外形參數化并通過CATIA二次開發［8］使參數化外形生成自動化，通過求解層流NS方程模擬考慮粘性的真實情況，為降低計算量利用了試驗設計方法對樣本點在設計空間內抽樣并進行CFD計算，依此構建Kriging模型［9］來代替優化過程中海量的CFD計算，優化算法采用NSGAⅡ［10－11］多目標遺傳算法，對后體尾噴管三維構型進行了優化設計，并對優化出的結果進行了CFD計算驗證。

1 優化方法

基于替代模型的氣動優化設計需要有四大技術作為支撐：高精度的CFD求解器、高效可靠的替代模型、高效且方便的現代優化算法、對優化外形的參數化建模及CFD計算網格的自動生成。

1.1 CFD求解器

本文采用課題組編寫的CFD求解器進行CFD計算，該求解器已經在多項科研及工程項目中得到成功應用。通過求解三維可壓全N－S方程來獲得尾噴管的推力和升力：

式中，Q為守恒變量，F、G、H分別為三個方向上的無粘通量，Fv、Gv、Hv分別為三個方向上的粘性通量。

空間格式采用高分辨率的Roe的FDS格式，時間格式為高可靠性、高效率的LU－SGS格式。

1.2 Kriging模型

Kriging模型是一種全局替代模型。通常Kriging模型包括兩部分，具體模型為：

式中：f（x）是對全部設計空間的全局模擬，可看作一個常數β，β值可由已知響應值進行估計，z（x）是期望為0、方差為σ2的高斯隨機函數，表示全局模擬的插值。z（x）的協方差矩陣可表示為：

式中：R為相關矩陣；R（xi，xj）為相關函數；i，j＝1，2，…，n。n為已知的樣本中數據的個數，其表達式為：

式中：nk為設計變量的個數；θk為未知的相關參數矢量。一般地可用標量θ來代替θk。

根據Kriging理論，未知點x處的響應值y的估計值可通過下式給出：

式中：y為樣本的響應值；f為長度為n的單位列向量；rT（x）為未知向量x與樣本輸入數據之間的相關向量［x1，x2，…xn］，表達式為：

相關系數θ可以由極大似然估計給出。

1.3 NSGAⅡ多目標遺傳算法

NSGAⅡ在2002年提出。與基本遺傳算法相比，NSGAⅡ在適應度計算方面主要有非支配排序及擁擠度排序這兩個特性。非支配排序使處于前緣的個體具有更大的概率被選中執行交叉變異等操作進入下一代；擁擠度排序使個體分散更加均勻，并能更好地保持解的多樣性。NSGA II中，對每個解來說，需要確定多少解支配它和它支配的解集。NSGA II需要估計圍繞著種群中一個特定解的解密度，即沿著問題的每個目標計算兩個解之間的平均距離，這個值被稱為密集距離（crowding measure）。在選擇期間，NSGA II密集比較算子既考慮種群中個體的非劣解秩，也考慮密集距離。也就是說，優先選擇非劣解；但當兩個解具有相同的非劣解秩時，則選擇那個不處于擁擠距離區域內的解。與其它多目標遺傳算法相比，NSGA II的精英保留策略使用（μ＋λ）選擇，包含了最好的父代和子代個體。正是這種機制使新一代種群比前一代種群更有效，效果更好。本文采用ISIGHT軟件包中的NSGAⅡ程序進行優化計算。

2 三維后體尾噴管優化設計

吸氣式高超聲速飛行器后體尾噴管要求有盡可能大的推力及升力，本文即選取推力及升力作為優化目標。首先對三維噴管外形進行參數化，并編制VB程序調用CATIA接口進行二次開發使外形生成自動化；計算網格采用Gridgen軟件生成并通過其宏錄制功能使網格生成自動化，結合高精度CFD求解器及NSGAⅡ多目標遺傳算法對三維尾噴管進行優化設計。流動條件與文獻［7］一致，噴流入口條件為：Ma＝1.2，靜溫T＝2170K，靜壓p＝0.195MPa；外部來流條件為：Ma＝6，靜壓p＝2978Pa，靜溫T＝235.45K，攻角4.5°。

2.1 三維外形參數化及CAD模型自動生成

尾噴管整體外形如圖1所示，由內膨脹段（有下壁板）和外膨脹段（無下壁板）構成。對稱面形狀如圖2所示，在對稱面外形基礎上，沿流向在對稱面型線上定義一系列展向剖面線，就刻畫出整個三維尾噴管上表面型線。

尾噴管總長L和噴管入口高度h確定，分別為850mm和64mm。對稱面采用三次曲線的建模方式，三次曲線多變的外形及良好的適應性能較好地描述多樣的噴管外形，燃燒室出口和噴管上表面用圓弧過渡，圓弧半徑為R，即圖中的圓弧AB。BC段為三次曲線，為保證曲線上凸，把C點配置為拐點。BC段與AB段在B點處以初始膨脹角θ1相切，三次曲線BC的出口角為θ2；噴管下壁板OE長度為Lb，OE與水平線夾角為θc。其中圓弧半徑R，初始膨脹角與噴管出口角θ1、θ2，下壁板長度Lb、角度θc為變量，在確定變量之后，即可確定對稱面的整體外形。

圖1 三維尾噴管建模示意圖Fig.1 CAD model of 3Dnozzle

圖2 對稱面外形及參數化Fig.2 Symmetry plane geometry map of 3Dnozzle

尾噴管內膨脹段（OE段）上壁面型線，通過對基準面型線展向延拓得到。在每一個展向剖面上，采用超橢圓定義展向剖面型線。圖3為三維尾噴管展向剖面的建模示意圖。

在圖3中，以‘o’點為坐標原點的超橢圓函數可以表示為：

圖3 三維尾噴管展向剖面Fig.3 Section plane curve of 3Dnozzle

其中L和H表示尾噴管型面在某個流向截面上的展向寬度和高度，L表示尾噴管側緣距基準面的展向距離，H表示側緣距基準面的法向距離，參數n可以控制超橢圓的彎曲程度，而φ＝L／H表示超橢圓的寬高比。同過以上的方程三維尾噴管展向截面的剖面線可以通過L，H，n這三個參數得到完全描述。L在尾噴管的內膨脹段和發動機入口處寬度相等，為150mm，在外膨脹段，噴管末端寬度定義為LL已知，為220mm。內膨脹段與外膨脹段通過一段圓弧過渡，圓弧半徑為r，圓弧角為θ，圓弧與一段拐點配置在尾噴管末端且外凸的三次曲線相連。其中r和θ為設計變量。圖4為展向寬度變化規律示意圖。高度方面，在內膨脹段尾噴管的高度由下檔板的位置確定，在外膨脹段沿流向的高度由以下公式確定：

式中He為尾噴管內膨脹段結束處上下型面之間的高度，ah為一設計變量，決定側壁高度的變化，ah越大，則側壁高度隨x增加而減小越快，當ah為0時，側壁高度保持He不變。沿流向，噴管截面超橢圓曲線指數通過一下方程確定：

式中M為一正整數，本文取為200，c為一較小整數，本文取為3。an為設計變量，確定展向截面的形狀，an越大則展向截面越接近于矩形。通過上面的型面建模過程，三維尾噴管的型面可以由完全由θ1、θ2、θc、R、Lb、r、θ、ah、an九個參數描述，其中前五個為二維對稱面的參數，后四個為向三維延拓后所增加的參數。九個參數的取值下限為［1°，1°，－10°，50mm，80mm，50mm，1°，0，0.01］，取值上限為［40°，40°，10°，500mm，250mm，200mm，10°，0.5，0.2］。原始外形的九個參數分別為［20°，20°，－3°，200mm，100mm，100mm，5°，0.2，0.1］。

圖4 三維尾噴管展向寬度變化規律Fig.4 Side wall width changing curve

在三維建模中，先根據具體外形在流向選定疏密適當的100個站位，而后根據讀入的外形參數可計算得到該100個站位上展向切面上的控制型線，在每條控制型線生成21個控制點并連接為樣條線，則半個上表面曲面即可采用這100條樣條線結合為多截面曲面得到。在VB程序中，得到九個控制外形的參數后，即可計算得到個控制線上點的坐標，并調用CATIA接口將點連成線并最終生成曲面，實現外形生成的自動化。自動生成的半模外形見圖5。

圖5 采用CATIA二次開發自動生成的半模外形Fig.5 Half of the model automatically created by CATIA secondary development

2.2 計算網格及其自動生成

計算網格為分區對接結構網格，采用Gridgen軟件生成，錄制宏使網格生成自動化。由于考慮粘性流動，網格采用內O型拓撲以便于邊界層的布置，網格見圖6，共有8個分區，網格總數為217280。

圖6 Gridgen軟件自動生成的計算網格Fig.6 Grid generated automatically by Gridgen

2.3 試驗設計及Kriging模型構建

采用優化拉丁方方法［12］在整個設計區間內抽取100個樣本點進行CFD計算得出其推力和升力，以此構建Kriging替代模型。優化拉丁方方法在隨機拉丁方方法的基礎上做了改進，使樣本點散布更加均勻。為驗證替代模型精度，在設計空間內再次隨機抽取20個點分別采用CFD和替代模型計算以檢驗替代模型精度，結果如圖7和圖8所示，可見替代模型對升力FL的模擬效果好于對推力FT的模擬，替代模型的相對誤差多在5%以內，最大不超過10%，替代模型精度滿足使用要求。

2.4 多目標優化結果及分析

遺傳算法參數選取為：實數編碼，種群大小為100，遺傳算法共計算100代；交叉操作采用單點交叉，交叉概率0.9，變異操作采用單點變異，變異概率為0.1。

圖7 CFD與Kriging模型計算推力結果Fig.7 Thrust computed by CFD and Kriging model

圖8 CFD與Kriging模型計算升力結果Fig.8 Lift computed by CFD and Kriging model

圖9為優化得到的Pareto前沿，并標明了原始噴管和選優噴管在Pareto前沿中對應的位置，可見Pareto前緣較為飽滿，得到的優選噴管性能較原始噴管有了較大的提高。優選噴管的性能應保證推力和升力較原始噴管都有所提升，并綜合考慮使噴管性能最優。由圖9可見，原始構型的推力已足夠大，若使推力有大幅度提升則會大大損失升力，因此所選擇的優選構型推力僅有小幅度提高，但升力較原始構型提高了123.3%。優選噴管的九個參數為［17.65°，1.14°，－5.69°，211.68mm，199.92mm，70.03mm，1.36°，0.16，0.18］。由于這是替代模型得出的結果，再次對選優噴管進行CFD計算，結果見表1，可見兩者相差不大，基于替代模型的優化設計達到了所期望的結果。

圖9 多目標優化結果和Pareto前緣Fig.9 Multi－objects optimization result for 3Dnozzle

圖10為原始外形與優選外形壁面等壓力線圖，可見優選外型較原始構型沿流向壓力下降更緩慢，沿展向分布更加均勻。圖11為原始外形和優化外形對稱面等壓力線圖，可見優化外形內膨脹段較長，流動在下壁板出口處產生的膨脹波較原始噴管弱很多，因而使上壁面壓力沒有很快的下降，從而使升力大大增加。

表1 原始構型與優化構型氣動性能比較Table 1 Comparison between original and optimized 3Dafterbody nozzle

圖10 原始尾噴管與優化尾噴管壁面壓力比較Fig.10 Pressure contour comparison for optimization and original nozzle

圖11 原始尾噴管與優化尾噴管對稱面流動等壓力線圖Fig.11 Pressure contour comparison of symmetry plane for optimization and original nozzle

本次優化工作在辦公機上進行，使用單核CPU 2.4GHz進行計算。共進行CFD計算122次（構建替代模型100次、驗證替代模型20次、原始及優化外形計算2次），每次CFD計算耗時30min，加上遺傳算法運行時間10min，共耗時3670min。若不采用替代模型，優化過程共需進行CFD計算10000次，需耗時300000min，則優化工作根本無法進行。替代模型的引入大大縮短了優化過程所需時間；在精度上，由優化結果的CFD驗證可以看出，替代模型的相對誤差很小，精度滿足要求。

3 結 論

本文通過試驗設計方法抽取樣本點進行CFD計算以構建替代模型來代替優化過程中的大規模CFD計算，在較短的時間內完成了吸氣式高超聲速飛行器三維后體尾噴管優化設計。在工程實際中，三維構型氣動優化方法難以應用的瓶頸是復雜外形參數化及自動生成困難以及過大的計算量，本文的工作表明：

（1）采用CATIA二次開發能方便快捷地實現三維后體尾噴管參數化外形的自動生成，結合Gridgen軟件的宏功能實現計算網格的自動生成，能高效地實現前處理的自動化，便于在工程中應用。

（2）使用優化拉丁方方法抽樣來構建Kriging模型，能準確地模擬三維后體尾噴管在不同外形參數下產生的推力及升力，對非線性系統具有很好的模擬能力。

（3）基于替代模型的優化設計方法能大大縮短優化時間，高效地完成了三維后體尾噴管的優化設計，便于三維構型的優化設計工作廣泛應用于工程實際。由于后體尾噴管為吸氣式高超聲速飛行器的一個部件，設計時必須考慮飛行器整體力矩平衡并兼顧各種性能，進一步的研究將進行三維機體／推進系統的一體化優化設計，即在本文方法的基礎上，發展出一套對三維整機進行參數化的方法，并對其進行氣動性能優化設計。

［1］李曉宇.高超聲速飛行器一體化布局外形設計［D］.長沙：國防科學技術大學，2007.

［2］EDWARDS C L Q，SMALL W J，WEIDNER J P.Studies of scramjet／airframe integration techniques for hypersonic aircraft［R］.AIAA 75－58，1975.

［3］BAYSAL O，ELESHAKY M，BURGREEN G.Aerodynamic shape optimization using sensitivity analysis on third－order Euler equations［J］.JournalofAircraft，1993，30（6）：953－961.

［4］羅世彬.高超聲速飛行器機體發動機一體化及總體多學科設計優化方法研究［D］.長沙：國防科學技術大學，2004.

［5］陳兵.空間推進算法及超燃沖壓發動機部件優化設計研究［D］.北京：北京航空航天大學，2005.

［6］甘文彪，閻超.高超聲速飛行器后體／尾噴管優化設計［J］.北京航空航天大學學報，2011，37（11）：1440－1445.

［7］賀旭照，倪鴻禮，周正，等.吸氣式高超聲速飛行器三維后體尾噴管優化設計［J］.推進技術，2009，30（6）：687－690.

［8］胡挺，吳立軍.CATIA二次開發技術基礎［M］.北京：電子工業出版社，2006.

［9］ALEXANDER I J fORRESTER，ANDY J KEANE.Recent advances in surrogate－based optimization［J］.ProgressinAerospaceSciences，2009，45（3）：50－79.

［10］SRINIVAS N，DEB K.Multiobjective optimization using nondominated sorting in genetic algorithms［J］.EvolutionaryComputation，1994，2（3）：221－248.

［11］DEB K，PRATAP A，AGARWAL S，et al.A fast and elitist multiobjective genetic algorithm：NSGAⅡ［J］.IEEETransactionsonEvolutionaryComputation，2002，6（2）：182－197.

［12］ANTHONY A G.Use of data sampling，surrogate models，and numerical simulation optimization in engineering design［R］.AIAA 2002－0538，2002.