基于替代模型的三維后體尾噴管優(yōu)化設(shè)計(jì)

苗 萌，曾 鵬，閻 超

（1.北京航空航天大學(xué) 國家計(jì)算流體力學(xué)實(shí)驗(yàn)室，北京 100191；2.北京臨近空間飛行器系統(tǒng)工程研究所，北京 100076；3.中國航天電子技術(shù)研究院無人機(jī)系統(tǒng)工程研究所，北京 100094）

0 引 言

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

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

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

1 優(yōu)化方法

基于替代模型的氣動(dòng)優(yōu)化設(shè)計(jì)需要有四大技術(shù)作為支撐：高精度的CFD求解器、高效可靠的替代模型、高效且方便的現(xiàn)代優(yōu)化算法、對(duì)優(yōu)化外形的參數(shù)化建模及CFD計(jì)算網(wǎng)格的自動(dòng)生成。

1.1 CFD求解器

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

式中，Q為守恒變量，F(xiàn)、G、H分別為三個(gè)方向上的無粘通量，F(xiàn)v、Gv、Hv分別為三個(gè)方向上的粘性通量。

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

1.2 Kriging模型

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

式中：f（x）是對(duì)全部設(shè)計(jì)空間的全局模擬，可看作一個(gè)常數(shù)β，β值可由已知響應(yīng)值進(jìn)行估計(jì)，z（x）是期望為0、方差為σ2的高斯隨機(jī)函數(shù)，表示全局模擬的插值。z（x）的協(xié)方差矩陣可表示為：

式中：R為相關(guān)矩陣；R（xi，xj）為相關(guān)函數(shù)；i，j＝1，2，…，n。n為已知的樣本中數(shù)據(jù)的個(gè)數(shù)，其表達(dá)式為：

式中：nk為設(shè)計(jì)變量的個(gè)數(shù)；θk為未知的相關(guān)參數(shù)矢量。一般地可用標(biāo)量θ來代替θk。

根據(jù)Kriging理論，未知點(diǎn)x處的響應(yīng)值y的估計(jì)值可通過下式給出：

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

相關(guān)系數(shù)θ可以由極大似然估計(jì)給出。

1.3 NSGAⅡ多目標(biāo)遺傳算法

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

2 三維后體尾噴管優(yōu)化設(shè)計(jì)

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

2.1 三維外形參數(shù)化及CAD模型自動(dòng)生成

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

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

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

圖2 對(duì)稱面外形及參數(shù)化Fig.2 Symmetry plane geometry map of 3Dnozzle

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

在圖3中，以‘o’點(diǎn)為坐標(biāo)原點(diǎn)的超橢圓函數(shù)可以表示為：

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

其中L和H表示尾噴管型面在某個(gè)流向截面上的展向?qū)挾群透叨?，L表示尾噴管側(cè)緣距基準(zhǔn)面的展向距離，H表示側(cè)緣距基準(zhǔn)面的法向距離，參數(shù)n可以控制超橢圓的彎曲程度，而φ＝L／H表示超橢圓的寬高比。同過以上的方程三維尾噴管展向截面的剖面線可以通過L，H，n這三個(gè)參數(shù)得到完全描述。L在尾噴管的內(nèi)膨脹段和發(fā)動(dòng)機(jī)入口處寬度相等，為150mm，在外膨脹段，噴管末端寬度定義為LL已知，為220mm。內(nèi)膨脹段與外膨脹段通過一段圓弧過渡，圓弧半徑為r，圓弧角為θ，圓弧與一段拐點(diǎn)配置在尾噴管末端且外凸的三次曲線相連。其中r和θ為設(shè)計(jì)變量。圖4為展向?qū)挾茸兓?guī)律示意圖。高度方面，在內(nèi)膨脹段尾噴管的高度由下檔板的位置確定，在外膨脹段沿流向的高度由以下公式確定：

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

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

圖4 三維尾噴管展向?qū)挾茸兓?guī)律Fig.4 Side wall width changing curve

在三維建模中，先根據(jù)具體外形在流向選定疏密適當(dāng)?shù)?00個(gè)站位，而后根據(jù)讀入的外形參數(shù)可計(jì)算得到該100個(gè)站位上展向切面上的控制型線，在每條控制型線生成21個(gè)控制點(diǎn)并連接為樣條線，則半個(gè)上表面曲面即可采用這100條樣條線結(jié)合為多截面曲面得到。在VB程序中，得到九個(gè)控制外形的參數(shù)后，即可計(jì)算得到個(gè)控制線上點(diǎn)的坐標(biāo)，并調(diào)用CATIA接口將點(diǎn)連成線并最終生成曲面，實(shí)現(xiàn)外形生成的自動(dòng)化。自動(dòng)生成的半模外形見圖5。

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

2.2 計(jì)算網(wǎng)格及其自動(dòng)生成

計(jì)算網(wǎng)格為分區(qū)對(duì)接結(jié)構(gòu)網(wǎng)格，采用Gridgen軟件生成，錄制宏使網(wǎng)格生成自動(dòng)化。由于考慮粘性流動(dòng)，網(wǎng)格采用內(nèi)O型拓?fù)湟员阌谶吔鐚拥牟贾?，網(wǎng)格見圖6，共有8個(gè)分區(qū)，網(wǎng)格總數(shù)為217280。

圖6 Gridgen軟件自動(dòng)生成的計(jì)算網(wǎng)格Fig.6 Grid generated automatically by Gridgen

2.3 試驗(yàn)設(shè)計(jì)及Kriging模型構(gòu)建

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

2.4 多目標(biāo)優(yōu)化結(jié)果及分析

遺傳算法參數(shù)選取為：實(shí)數(shù)編碼，種群大小為100，遺傳算法共計(jì)算100代；交叉操作采用單點(diǎn)交叉，交叉概率0.9，變異操作采用單點(diǎn)變異，變異概率為0.1。

圖7 CFD與Kriging模型計(jì)算推力結(jié)果Fig.7 Thrust computed by CFD and Kriging model

圖8 CFD與Kriging模型計(jì)算升力結(jié)果Fig.8 Lift computed by CFD and Kriging model

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

圖9 多目標(biāo)優(yōu)化結(jié)果和Pareto前緣Fig.9 Multi－objects optimization result for 3Dnozzle

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

表1 原始構(gòu)型與優(yōu)化構(gòu)型氣動(dòng)性能比較Table 1 Comparison between original and optimized 3Dafterbody nozzle

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

圖11 原始尾噴管與優(yōu)化尾噴管對(duì)稱面流動(dòng)等壓力線圖Fig.11 Pressure contour comparison of symmetry plane for optimization and original nozzle

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

3 結(jié) 論

本文通過試驗(yàn)設(shè)計(jì)方法抽取樣本點(diǎn)進(jìn)行CFD計(jì)算以構(gòu)建替代模型來代替優(yōu)化過程中的大規(guī)模CFD計(jì)算，在較短的時(shí)間內(nèi)完成了吸氣式高超聲速飛行器三維后體尾噴管優(yōu)化設(shè)計(jì)。在工程實(shí)際中，三維構(gòu)型氣動(dòng)優(yōu)化方法難以應(yīng)用的瓶頸是復(fù)雜外形參數(shù)化及自動(dòng)生成困難以及過大的計(jì)算量，本文的工作表明：

（1）采用CATIA二次開發(fā)能方便快捷地實(shí)現(xiàn)三維后體尾噴管參數(shù)化外形的自動(dòng)生成，結(jié)合Gridgen軟件的宏功能實(shí)現(xiàn)計(jì)算網(wǎng)格的自動(dòng)生成，能高效地實(shí)現(xiàn)前處理的自動(dòng)化，便于在工程中應(yīng)用。

（2）使用優(yōu)化拉丁方方法抽樣來構(gòu)建Kriging模型，能準(zhǔn)確地模擬三維后體尾噴管在不同外形參數(shù)下產(chǎn)生的推力及升力，對(duì)非線性系統(tǒng)具有很好的模擬能力。

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

［1］李曉宇.高超聲速飛行器一體化布局外形設(shè)計(jì)［D］.長沙：國防科學(xué)技術(shù)大學(xué)，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］羅世彬.高超聲速飛行器機(jī)體發(fā)動(dòng)機(jī)一體化及總體多學(xué)科設(shè)計(jì)優(yōu)化方法研究［D］.長沙：國防科學(xué)技術(shù)大學(xué)，2004.

［5］陳兵.空間推進(jìn)算法及超燃沖壓發(fā)動(dòng)機(jī)部件優(yōu)化設(shè)計(jì)研究［D］.北京：北京航空航天大學(xué)，2005.

［6］甘文彪，閻超.高超聲速飛行器后體／尾噴管優(yōu)化設(shè)計(jì)［J］.北京航空航天大學(xué)學(xué)報(bào)，2011，37（11）：1440－1445.

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

［8］胡挺，吳立軍.CATIA二次開發(fā)技術(shù)基礎(chǔ)［M］.北京：電子工業(yè)出版社，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.