徐 嘉，蔡晉生，劉秋洪

（西北工業(yè)大學 翼型葉柵空氣動力學國防科技重點實驗室，陜西 西安 710072）

0 引 言

減阻技術，一直是飛機設計研究中的熱點問題?，F(xiàn)代噴氣戰(zhàn)斗機的后體阻力大約占到全機阻力的38%～50%［1］，所以有必要對這類阻力進行研究。研究后體阻力之前，首先需要對飛機后體流場進行準確的預測。早期科研人員通過采用風洞實驗的方法得到后體噴流流場信息［2－3］。隨著計算機與數(shù)值方法的發(fā)展，科研人員采用數(shù)值方法研究飛機后體噴流與外流相互干擾問題［4－6］。

戰(zhàn)斗機后體外形設計方法通常是將表面外形設計成流線型，避免后體流場對噴管的不利干擾，使得噴管表面得到很好的再壓縮［1］。Catt等［7－8］對 F－16戰(zhàn)斗機后體外形進行減阻設計，通過延長后體長度，改變后體外形曲線等方法降低后體阻力。這種設計的思路就是保證戰(zhàn)斗機后體外形表面光滑流線，不發(fā)生流動分離，從而達到良好的減阻效果。由于方案設計取決于實際經驗，因此具有一定的局限性。隨著CFD理論的不斷深入，科研人員開始探索新的減阻技術。Kentfield［9］提出了后體多臺階被動減阻方法，與原有的減阻思路不同，外形表面并不是設成光滑流線，而是用多個階梯代替。當氣流流過階梯時，階梯內產生環(huán)形渦，使得流向膨脹相對較小，表面壓力得到很好的恢復，從而降低了后體阻力。對于超聲速和高超聲速流動，這種方法能夠降低后體阻力，但是對于不同的飛行狀態(tài)，需要設計相應的臺階數(shù)目和位置。除此之外，還有多種主動減阻措施，如噴流、微吹吸氣、振動壁面等，這樣飛機必須添加額外機構或系統(tǒng)才能實現(xiàn)，往往會增加自重。被動減阻措施，如鼓包修形［10］，凹陷修形［1］等，在亞聲速飛行狀態(tài)，不但不能減阻反而還會使阻力增大，只有在特定的飛行狀態(tài)下的被動減阻措施才有明顯減阻效果。

本文通過計算二維軸對稱N－S方程、k－ω SST湍流模型和組分方程，對戰(zhàn)斗機后體繞流與尾噴流進行數(shù)值模擬，研究不同噴流介質和噴壓比對后體阻力的影響。為了降低后體阻力，提高戰(zhàn)斗機的氣動性能，采用梯度法［11－12］對后體外形進行減阻優(yōu)化設計。與其他優(yōu)化方法相比，梯度法簡單、實用，但是計算效率不高。為了提高梯度法的優(yōu)化效率，本文提出一種優(yōu)化設計點加速算法。

1 后體流場數(shù)值模擬

1.1 流場計算方法

流動控制方程為二維軸對稱Navier－Stokes方程，采用有限體積法對方程進行離散求解，方程中的對流項采用二階精度的Roe－FDS格式［13］進行離散，粘性項采用中心差分格式，通過隱式LU－SGS［14］進行時間推進。湍流模型選用k－ω SST模型，以k－ω為基礎的剪切應力輸運（SST）模型計算湍流剪切應力的輸運項，準確模擬在逆壓梯度下邊界層湍流剪切應力的影響［15］。

假設當前發(fā)動機內燃燒完全，流動過程中不考慮化學反應，燃氣噴流可視為多組分的凍結流?；旌衔镏忻恳环N組分氣體都滿足氣體狀態(tài)方程，假設相同控制體內每種組分氣體溫度相同，對濃度相對較大的水蒸氣、氧氣和氮氣進行數(shù)值模擬。不同氣體質量組分方程采用以下形式：

其中Yi為第i種物質質量分數(shù)，Ji為物質i的擴散通量。

1.2 后體幾何外型

本文采用的后體幾何外型與文獻［16］進行實驗的幾何外型相同，圖1給出了幾何外型尺寸示意圖。坐標軸x，y分別沿對稱軸及徑向，原點位于尾噴管喉道o點處，箭頭表示氣流流動方向。機身直徑d＝0.1524m，后體出口高度e＝0.03885m，內噴管喉道處高度t＝0.02935m，收縮段長度s＝0.1084m，擴張段長度f＝0.1506m。

1.3 計算網格和邊界條件

根據圖1所示后體幾何外型，生成如圖2所示的計算網格。后體外流場計算網格為H型結構網格，網格邊界到物面距離為15倍后體直徑，根據來流雷諾數(shù)Re∞確定物面邊界層網格的尺寸為2×10－5左右。噴管內流場網格為H型結構，在噴管內流場與外流場剪切層區(qū)域分布較密網格。

圖1 計算模型示意圖Fig.1 Computational model

物面采用絕熱、無滑移邊界條件，側向與前緣遠場邊界取無反射邊界條件，下游遠場邊界采用壓力出口邊界條件。在內流場的噴管進口處，指定噴流馬赫數(shù)Mj與噴流溫度Tj。當Mj＜1時，噴管進口靜壓值pj由計算域外推得到；當Mj＞1時，pj取指定值。

1.4 計算結果與分析

1.4.1 后體流場特性

為了驗證計算方法與網格的可行性，采用文獻［16］的實驗來流條件，馬赫數(shù)Ma∞為1.2，迎角α為0°，基于機身直徑的雷諾數(shù)Re∞為14.06×106，噴管出口馬赫數(shù)Mae為2.0，出口總溫Ttotal為1013K，噴管壓比pt，j／p∞（Nozzle Pressure Rates，簡稱 NPR）為8，噴流介質組分為57.7%水蒸氣和42.3%氧氣，水蒸氣和氧氣不再發(fā)生化學反應，噴流氣體的比熱比γ為1.265，氣體常數(shù)R 為376.19J／kg·K，外場氣體由77%氮氣和23%氧氣組成。

圖3為后體表面壓力分布與流線圖，圖4為后體流場馬赫數(shù)云圖。從圖3中可以看到，本文的計算結果與文獻的實驗結果幾乎吻合，由于外流在x／d為0.0的位置以后經歷了一個膨脹加速的過程，氣流加速使得表面壓力降低，所以在0.0＜x／d＜0.56區(qū)間內后體壓力值逐漸降低，在x／d為0.56左右的位置達到最低點之后迅速升高，這是由于跨聲速氣流在流過后體表面膨脹加速后在x／d為0.56左右的位置產生激波。由于后體表面激波導致邊界層流動分離而產生分離渦，如圖3的流線圖所示。

圖5為噴流馬赫數(shù)分布云圖，從圖中可以看到，噴流與外流相互作用形成剪切層，噴口處的膨脹波在中心線處相互作用后延伸至噴流邊界反射形成壓縮波，多個壓縮波疊加在一起形成一道截斷膨脹波的攔截激波以及馬赫盤。

圖5 噴流馬赫數(shù)云圖Fig.5 Contour of Mach number on jet

1.4.2 噴流介質對后體阻力的影響

噴流氣體介質的比熱比是影響噴流和外流相互作用的一個重要參數(shù)，多組分氣體的比熱比隨著流場氣體質量變化而變化，理想氣體的比熱比為常數(shù)，現(xiàn)在研究不同比熱比的噴流氣體對后體阻力的影響。來流邊界條件與1.4.1節(jié)相同，兩種噴流氣體介質和出口參數(shù)如表1所示。

表1 噴流介質出口狀態(tài)參數(shù)Table 1 The parameters for jet flow

圖6為理想氣體噴流介質實驗與計算，以及組分氣體噴流介質計算后體表面壓力分布比較圖。從圖中可以看到，與組分氣體噴流介質后體壓力分布相比，理想氣體噴流介質的后體表面壓力分布的激波位置略靠后。采用理想氣體作為噴流介質的后體外形阻力系數(shù)為0.199，組分氣體噴流介質的后體阻力系數(shù)為0.184，理想氣體噴流介質的阻力系數(shù)高于組分氣體噴流介質的阻力系數(shù)7.5%。由于采用組分氣體做噴流介質需要求解組分方程，造成數(shù)值計算中計算量的增加。在優(yōu)化設計的過程中，在不影響計算精度的情況下，為了降低計算量在優(yōu)化設計中可采用理想氣體作為噴流介質。

圖6 后體壓力分布比較Fig.6 Comparison of pressure distribution

1.4.3 噴管壓比對后體阻力系數(shù)的影響

噴管壓比對后體表面壓力分布有較大的影響，通過模擬不同噴管壓比的后體流場，研究噴管壓比對后體阻力系數(shù)的影響。采用1.4.2節(jié)相同的流場條件，表2給出不同噴管壓比下后體阻力系數(shù)的變化?？梢钥吹?，后體阻力系數(shù)隨噴管壓比的增大而減小。因為隨著噴管壓比的增大，噴流膨脹主要對后體尾部產生較大的壓力分布，從而降低后體壓差阻力。

表2 不同噴管壓比下的阻力系數(shù)Table 2 Drag coefficients with different pressure rates

2 后體減阻優(yōu)化設計

2.1 優(yōu)化設計方法

應用梯度法對后體外形進行減阻優(yōu)化設計。在優(yōu)化之前，需要構造設計變量與目標函數(shù)，如圖7所示，從后體外形上選取n個設計控制點，第i設計控制點的坐標值為（xi，yi），設計控制點軸向坐標xi值取后體弦長的i個平均值，設計控制點徑向坐標yi值為設計變量，采用三次樣條插值函數(shù)來表示后體外形，確定后體起始點O處的一階導數(shù)值和N點的自然邊界插值，以保證后體外形與機身平直段的光滑過渡。

圖7 后體外形及初始優(yōu)化控制點分布示意圖Fig.7 Design pointes and optimal model of afterbody

定義后體阻力系數(shù)為：

其中，U∞為來流速度，ρ∞為來流密度，Daft為后體阻力。以徑向坐標y值為設計變量，最小阻力系數(shù)Cd，aft為設計目標，定義目標函數(shù)為：

梯度法優(yōu)化后體外形的具體過程如下：

（1）給定初始點y0，初始搜索步長λ0和誤差限0＜ε≤1迭代步數(shù)i＝1

（2）計算搜索方向si＝ －▽f（yi），代入到迭代公式y(tǒng)i＋1＝y(tǒng)i＋λi·si；

（3）根據徑向坐標yi＋1值和軸向坐標x值得到后體外形優(yōu)化曲線，應用CFD計算出后體阻力系數(shù)Cd，aft，即得到f（yi＋1）；

（4）若‖f（yi＋1）－f（yi）‖≤ε，則停止計算，yi＋1所確定的后體外形為最終優(yōu)化外形；否則，求最優(yōu)步長λi＋1，使得不等式f（yi＋λi＋1·di）＜f（yi）成立。

（5）令i＝i＋1，轉（2）。

2.2 優(yōu)化設計加速算法

梯度法產生每一優(yōu)化步的搜索方向都要用到目標函數(shù)的梯度信息，求解流場的次數(shù)取決于設計控制點個數(shù)，當設計控制點個數(shù)比較多時，求解流場的次數(shù)會很多。選取少量的設計控制點雖然能夠降低優(yōu)化時間，但是不可能得到最優(yōu)的設計外形。所以，本文提出一種梯度法優(yōu)化設計的加速算法，即基于后體外形曲率值，逐步增加設計控制點數(shù)目并合理分布設計控制點的位置，從而提高優(yōu)化設計的效率和精度。

對于給定的后體外形曲線，其曲率分布為k（x），在曲線上布置n個優(yōu)化控制點（后體外形曲線的起始和終止位置相對固定，因此不作為設計控制點），可以得到設計控制點之間的區(qū)域數(shù)為m＝n＋1，定義區(qū)域曲率與長度積分的平均值為：

可以得到設計控制點的軸向坐標x1，…，xi，…，xn，使下式成立，

式（4）表示優(yōu)化控制點之間區(qū)域曲率相對長度積分值相同，使得設計控制點在外形曲率值較大的區(qū)域分布比較密集，而曲率值較小的區(qū)域分布比較稀少，這樣就可以采用較少的控制點來精確地描述較為復雜的外形曲線。

采用梯度法進行優(yōu)化設計，先選取少量的設計控制點對初始外形進行優(yōu)化設計，使目標函數(shù)基本達到最優(yōu)。然后增加設計控制點個數(shù)，通過計算外形曲線的曲率對所有設計控制點重新進行分布，對后體外形進行重新描述，采用梯度法進行新一輪優(yōu)化設計。重復上述逐步增加設計控制點過程，最終達到優(yōu)化設計目標。如圖8所示為基于曲率逐步增加設計控制點方法與梯度法相結合進行優(yōu)化設計的流程示意圖，其中Cd為新一輪優(yōu)化得到后體阻力值，Cd＿t為目標值，一般取沒有結合加速算法直接梯度法優(yōu)化得到的后體阻力值。

圖8 優(yōu)化設計流程示意圖Fig.8 The flow chart of the optimization

2.3 優(yōu)化設計結果與分析

2.3.1 后體減阻優(yōu)化結果

以1.4.1節(jié)算例中的計算外形為初始外形，采用相同的來流條件進行后體減阻優(yōu)化設計。流場數(shù)值計算結果發(fā)現(xiàn)噴流介質氣體組分對后體阻力影響較小，為了減少優(yōu)化設計中流場數(shù)值計算的時間，不考慮噴流介質中氣體組分的影響，以理想氣體為噴流介質。噴管出口馬赫數(shù)Mae為2.0，出口總溫Ttotal為1013K，噴管壓比NPR為8。取11個設計控制點進行后體外形減阻優(yōu)化設計，設計點軸向坐標xi值取后體弦長的i個平均值，即在初始外形上均勻分布設計點。

圖9對比優(yōu)化前后的后體外形，實線是優(yōu)化后的外形，虛線是優(yōu)化前的初始外形，實線上的點為均勻分布的設計控制點。與初始外形相比，優(yōu)化外形出現(xiàn)了一個類似“鼓包”形狀，外形其他部分表面光滑平直。

圖9 初始與優(yōu)化外形Fig.9 Initial and optimal shape

圖10為優(yōu)化前后外形表面壓力分布比較，實線（optimized）是優(yōu)化后的外形表面壓力分布，虛線（initial）是初始外形表面壓力分布。比起初始外形的壓力分布，在0＜x／d＜0.5區(qū)間內優(yōu)化外形的壓力不降反升，壓力最低點也要高于初始外形。優(yōu)化前后的后體阻力系數(shù)的計算結果分別為0.199和0.178，優(yōu)化后的外形使得后體阻力系數(shù)降低了10.6%。在初始外形與優(yōu)化外形的前段，外流都進行了一個膨脹加速的過程，比起初始外形，優(yōu)化外形表面表過渡的曲率變化小，后體表面壓力得到恢復，從而降低了后體阻力。在激波位置附近，優(yōu)化外形出現(xiàn)了一個類似“鼓包”形狀，弱化了激波強度，減小了壓力損失。

圖10 優(yōu)化前后的后體表面壓力分布Fig.10 Pressure distribution

2.3.2 優(yōu)化設計加速算法結果分析

采用本文提出的優(yōu)化設計加速算法法與梯度法相結合進行后體減阻優(yōu)化設計，可以提高優(yōu)化設計的效率和精度。先在初始外形上均勻布置2個設計控制點，得到2點優(yōu)化外形。然后增加3個優(yōu)化控制點，采用基于外形曲率的設計控制點分布方法重新分布5個優(yōu)化控制點位置，進一步得到5點優(yōu)化外形。為了更好地描述后體外形局部形狀，特別是曲率變化較大區(qū)域，在5個優(yōu)化控制點優(yōu)化外形的基礎上再增加6個優(yōu)化控制點，采用同樣的方法重新分布11個優(yōu)化控制點，通過梯度法進行優(yōu)化計算得到最終優(yōu)化外形。

圖11至圖13分別給出采用2、5與11個優(yōu)化控制點進行優(yōu)化設計得到的外形與初始外形的對比。圖中實線是優(yōu)化后外形，虛線是初始外形，實心點為設計點。從圖13中可以看到采用11個優(yōu)化控制點進行優(yōu)化設計，設計控制點主要集中分布在激波區(qū)域的“鼓包”附近以及后體起始位置附近。采用2、5和11個設計控制點優(yōu)化得到的后體阻力系數(shù)分別為0.184、0.178和0.175，與初始外形相比，最終得到的優(yōu)化外形使得后體阻力系數(shù)降低了13.0%。與直接均勻分布11個優(yōu)化控制點采用梯度法進行優(yōu)化設計的方法相比，本文提出的逐步增加設計控制點個數(shù)并根據外形曲率合理分布設計控制點位置與梯度法相結合的優(yōu)化設計加速算法能夠得到更好的優(yōu)化結果。

圖11 初始外形與2點優(yōu)化外形Fig.11 Initial and 2－point optimal shape

圖12 初始外形與5點優(yōu)化外形Fig.12 Initial and 5－point optimal shape

圖13 初始外形與11點優(yōu)化外形Fig.13 Initial and 11－point optimal shape

圖14為優(yōu)化過程中的收斂曲線圖，實心點劃線（11design）為直接11設計控制點梯度法迭代曲線，空心點劃線為逐步增加設計控制點方法迭代曲線。直接11設計控制點梯度法優(yōu)化迭代了26步，平均迭代1步需要3.07小時，完成11點優(yōu)化外形最終花費79.8小時。逐步增加設計控制點方法共計優(yōu)化迭代了32步，其中2設計控制點迭代18步，每步需要0.87小時，5設計控制點迭代7步，每步需要1.45小時，完成優(yōu)化計算最終花費47.3小時。與直接11設計控制點梯度法優(yōu)化方法相比，這種優(yōu)化加速算法使得優(yōu)化計算時間降低了40%左右。

圖14 收斂曲線Fig.14 Convergence curve

3 結 論

本文通過數(shù)值計算二維軸對稱Navier－Stokes方程、k－ω SST湍流模型和氣體組分方程，分析飛機后體流場特性，并進行后體減阻優(yōu)化設計。通過本文的數(shù)值分析，可以得到以下的幾點結論：

（1）由于后體表面激波導致邊界層流動分離而產生分離渦。噴流與外流相互作用形成剪切層，噴口處形成一系列膨脹波和壓縮波以及馬赫盤。在跨聲速來流下，理想氣體噴流介質的后體阻力系數(shù)略高于組分氣體噴流介質的后體阻力系數(shù)，而阻力系數(shù)隨著噴管壓比的增大而減小。

（2）通過優(yōu)化設計得到的后體外形曲率變化緩慢，使得表面壓力得到恢復，后體外形在激波位置附近出現(xiàn)“鼓包”形狀，弱化了激波強度，與原始外形相比后體阻力值降低13%。

（3）優(yōu)化設計加速算法逐步增加設計控制點個數(shù)，并根據外形曲率合理分布設計控制點的位置，更好地描述外形曲線。優(yōu)化結果表明，該算法不僅能夠有效的降低優(yōu)化設計時間，而且還能夠提高優(yōu)化設計精度。

［1］方寶瑞.飛機氣動布局設計［M］.北京：航空工業(yè)出版社，1997.

［2］CHAMBERLAIN D.Measurement of drag from interaction of jet exhaust and surrounding airframe［R］.AIAA Paper 68－395，1968.

［3］PRESZ W M，PITKIN E T.Flow separation over axisymmetric afterbody models［R］.AIAA Paper 74－17，1974.

［4］DEIWERT G S，AANDREWS A E，NAKAHASHI K.Theroretical analysis of aircraft afterbody flow［J］.J.Spacecraft ＆ Rockets，24（6）：496－503.

［5］蔡晉生，劉秋洪.超聲速流場中側向射流的數(shù)值研究［J］.空氣動力學學報，2010，28（5）：553－558.

［6］張正科.用N－S方程有限體積法計算彈體繞流／底噴流［J］.空氣動力學學報，2003，21（3）：288－294.

［7］CATT J A，WELTERLE T J.Decreasing F－16nozzle drag using computation fluid dynamics R .AIAA Paper 93－2572，1993.

［8］WELTERLEN T J，CATT J A C.Evaluating F－16nozzle drag using computational fluid dynamics［R］.AIAA Paper，94－0022，1994.

［9］KENTFIELD J A C.Drag reduction of controlled separated flows［R］.AIAA Paper，85－1800，1985.

［10］QIN N，WONG W S，LE Moigne A.Three－dimensional contour bumps for transonic wing dody drag reduction［A］.Proceedings of the Institution of Mechanical Engineers，Part G：Journal of Aerospace Engineering［C］.2008.

［11］YUAN YA－XIANG.A new stepsize for the steepest descent method［J］.Journal of Computational Mathe－matics 2006242149－156.

［12］JONATHAN BARZILAI，JONATHAN M Borwein.Two－point step size gradient methods［J］.IMA Journal of Numerical Analysis，1988，8：141148.

［13］閻超.計算流體力學方法及應用［M］.北京：北京航空航天大學出版社，2006.

［14］VOS J B，LEYLAND P，KEMEMADE V V，et al.NSMB handbook version 5.0［R］.2003：36－38.

［15］MENTER F R.Zonal two－equation k－ω turbulence model for aerodynamic flows［R］.AIAA Paper 93－2906，1993.

［16］WILLIAM B Comptom Ⅲ.Effects of jet exhaust gas properties on exhaust simulation and afterbody drag［R］.NASA TR R－444，1975.