趙承熙，葉正寅，華如豪

(西北工業(yè)大學(xué)翼型葉柵空氣動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，陜西西安 710072)

新型目標(biāo)壓力分布下的Licher雙翼反設(shè)計方法研究

趙承熙*，葉正寅，華如豪

(西北工業(yè)大學(xué)翼型葉柵空氣動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，陜西西安 710072)

為了提高超聲速氣動構(gòu)型的升阻比，減小激波阻力，使用反設(shè)計方法結(jié)合CFD技術(shù)優(yōu)化Licher雙翼，實(shí)際算例表明，來流馬赫數(shù)1.7，無粘情況下只需迭代15步即可得到優(yōu)化結(jié)果，優(yōu)化后翼型波阻可減小25.5%，升阻比提高30.5%。隨后依據(jù)翼型目標(biāo)壓力分布這一反設(shè)計關(guān)鍵點(diǎn)，分析了非零迎角下翼型各邊的受力情況，指出了原目標(biāo)壓力分布的不足，并提出了一種新的階梯形目標(biāo)壓力分布形式，該方法的優(yōu)化結(jié)果可使升阻比提高49.8%。此外基于NS方程的優(yōu)化結(jié)果表明，原目標(biāo)壓力分布的優(yōu)化效果被削弱，升阻比僅能提高17%，而新目標(biāo)壓力分布的優(yōu)化結(jié)果受到的影響較小，升阻比仍可提高49.2%，說明在考慮流動粘性特征時，階梯形目標(biāo)壓力分布形式更具實(shí)用價值。

Busemann消波原理;Licher雙翼;反設(shè)計方法;高升阻比;目標(biāo)壓力分布

0 引言

升阻比是衡量飛行器氣動效率的重要指標(biāo)之一。超聲速飛行器在飛行過程中的強(qiáng)波阻降低了飛行器的升阻比，影響了氣動效率，限制了民用和軍事領(lǐng)域超聲速飛行器的進(jìn)一步發(fā)展。因此，減小激波阻力、提高升阻比對超聲速飛行器的發(fā)展至關(guān)重要。

1935年，第五屆Volta空氣動力學(xué)國際大會上，Busemann首次提出Busemann雙翼的概念［1］，并于文獻(xiàn)［2］中詳細(xì)描述了消波原理。采用雙翼構(gòu)型，當(dāng)翼型間距離較近時，雙翼間的激波和膨脹波通過相互干涉可以明顯降低激波強(qiáng)度，如圖1。

圖1 Busemann消波原理Fig.1 Busemann concept

近幾年Busemann消波原理又引起廣泛關(guān)注，一些學(xué)者［3-7］利用該消波原理提出了一些新概念飛行器的氣動布局形式，Jameson等人使用伴隨方法優(yōu)化了該氣動構(gòu)型，解決了流動的壅塞問題［8-9］。

實(shí)際飛行中翼型具有一定迎角，在Busemann建立的理論基礎(chǔ)上，Licher提出了非對稱 Licher雙翼［10］。如圖2所示，根據(jù)疊加原理，將無厚度的平板雙翼和上下對稱的Busemann雙翼相結(jié)合，即保證了激波、膨脹波的相對位置，起到消波減阻的效果，又使該翼型具有一定的升力，得到了較大的升阻比。

圖2 Licher雙翼設(shè)計思路Fig.2 Licher biplane

但由于線化理論模型忽略了真實(shí)流動中的一些因素，Licher雙翼內(nèi)真實(shí)的波系分布不能達(dá)到理論預(yù)期的效果，因此，有必要對原翼型進(jìn)行優(yōu)化。目前梯度法、遺傳算法等優(yōu)化方法在工程上已經(jīng)大量推廣使用［11-13］，但這些方法需要大量的CFD數(shù)值模擬結(jié)果，因此計算的效率往往不能讓人滿意。另一種翼型優(yōu)化方法被稱為反設(shè)計方法，該方法一般通過給定的目標(biāo)壓力分布求解滿足這一壓力分布的氣動外形，計算效率較高，但也存在難以給定合理的目標(biāo)壓力分布的問題。

超聲速Busemann雙翼表面的壓力分布具有明顯的規(guī)律，故適合采用反設(shè)計方法對其進(jìn)行優(yōu)化。日本學(xué)者Kusunose教授團(tuán)隊(duì)成功的使用了反設(shè)計方法優(yōu)化了Licher翼型，進(jìn)一步減小了翼型的波阻［4，14-15］。

本文完善和發(fā)展了Licher雙翼的反設(shè)計方法，在保留原算法主要框架［4］的基礎(chǔ)上，提出翼型修正模塊解決優(yōu)化過程中的邊交叉、邊界不光滑，算法不收斂等問題;通過非零迎角下的翼型表面的受力分析，提出了一種新的目標(biāo)壓力分布方法，較文獻(xiàn)［4］進(jìn)一步提高了該方法的優(yōu)化效果;并引入基于NS方程的流場求解器，探究了粘性對該方法的影響。

1 Licher雙翼反設(shè)計方法原理

1.1 反設(shè)計法算法框架

圖3為反設(shè)計方法算法框架，該方法是迭代方法。首先確定初始翼型邊界函數(shù)f和目標(biāo)壓力系數(shù)Cp分布，采用CFD技術(shù)求解當(dāng)前翼型周圍流場，確定翼型表面Cp的分布，從而得到壓力系數(shù)的差值余量ΔCp，構(gòu)造ΔCp和翼型邊界函數(shù)變化量Δf之間的關(guān)系，確定新的翼型邊界函數(shù)f。反復(fù)重復(fù)上述步驟，直到ΔCp＜ε(ε為無量綱小量)

圖3 反設(shè)計優(yōu)化方法Fig.3 Inverse design method

但在實(shí)際優(yōu)化過程中，可能出現(xiàn)ΔCp過大造成的邊交叉、邊界不光滑或算法不收斂等現(xiàn)象(詳見圖4)，這里提出了一種復(fù)合光順方法對新生成的翼型進(jìn)行光順修形(見圖3下)，通過移動邊界解決邊交叉問題，采用多項(xiàng)式擬合不光滑邊界，具有自適應(yīng)調(diào)整多項(xiàng)式階數(shù)的功能，并可在局部使用直線消除擬合過程中曲線的Runge現(xiàn)象，大量數(shù)值實(shí)驗(yàn)表明，該光順方法有效提高了整個優(yōu)化系統(tǒng)的魯棒性。

圖4 無光順模塊的算法缺陷Fig.4 Setbacks of the method without smoothing airfoils

1.2 ΔCp和Δf間的關(guān)系

基于斜基波關(guān)系式［16］可以構(gòu)造壓力系數(shù)Cp和翼型表面邊界函數(shù)f之間的關(guān)系［4］，θ為壁面楔角(流動轉(zhuǎn)折角):

定義c1、c2為Busemann系數(shù)［16］，這里Ma∞為來流馬赫數(shù)，γ為比熱比。

可近似認(rèn)為超聲速翼型表面的曲線斜率與自由來流迎角的差等于θ(如圖5)，故可得翼型上下表面Cp與f之間的關(guān)系，這里 +代表上表面，-代表下表面，α為來流迎角。

1.3 目標(biāo)壓力分布的確定

翼型的目標(biāo)壓力分布是反設(shè)計方法的關(guān)鍵，圖6(a)展示了來流馬赫數(shù)為1.7情況下的Busemann雙翼表面的壓力系數(shù)云圖，(b)為壓力云圖對應(yīng)的Cp圖。分析圖6(b)可得，翼型表面Cp分布沿弦向可近似為直線，這是Busemann雙翼減小壓差阻力的原因。當(dāng)上下表面Cp完全變?yōu)橹本€，即Cp曲線圍成矩形時，翼型的壓差阻力為零。文獻(xiàn)［4］使用該原理，將優(yōu)化框架中的目標(biāo)壓力分布定義為矩形(后文簡稱矩形分布)，反復(fù)迭代優(yōu)化，達(dá)到減小波阻的效果。

圖6 Busemann雙翼壓力分布Fig.6 Pressure distribution of the Busemann airfoils

但是，當(dāng)翼型具有迎角時，矩形目標(biāo)壓力分布具有一定局限性。圖7為來流迎角為α的Licher雙翼受力分析圖，忽略粘性等復(fù)雜因素的影響，可分別列出上下翼的升力L和阻力D關(guān)系式(7)，這里p為翼型表面壓力，l為各邊長度，d為雙翼間距，下標(biāo)表示變量在翼型上的位置:

圖7 Licher雙翼受力分析圖Fig.7 Force analyses of Licher airfoils

根據(jù)壓力系數(shù)Cp=(p-p∞)/q∞，對式(7)進(jìn)行無量綱化得式(8)，其中無量綱長度珋l=l/lout，Cl、Cd分別表示升力系數(shù)和阻力系數(shù)。

采用文獻(xiàn)［4］做法時，目標(biāo)Cp分布定義為矩形，對應(yīng)的Cp1=Cp2;Cp3=Cp4，這顯然與式(10)結(jié)論相悖。故提出一種新的目標(biāo)壓力分布形式，使Cp1、Cp2、Cp3、Cp4之間的關(guān)系滿足式(10)，如圖8所示，Cp分布形似臺階，故簡稱該種目標(biāo)壓力分布為“階梯形分布”。

由于Cl_up＞0;Cl_low＜0，而該目標(biāo)壓力分布會使Cp2相對增大，Cp4相對減小，從而導(dǎo)致Cl_up和Cl_low同時增大。因此，階梯形分布不僅可讓Cd_up→0;Cd_low→0，還可以使雙翼升力增大，進(jìn)一步提高翼型升阻比。

圖8 上翼目標(biāo)Cp分布Fig.8 Target Cpdistribution of up airfoil

2 計算結(jié)果與分析

2.1 基于Euler方程的Licher翼型優(yōu)化算例

初始翼型選取Licher雙翼，相關(guān)參數(shù)見表1(參數(shù)具體含義參考圖7)，翼型表面控制點(diǎn)400個，各點(diǎn)流場狀態(tài)變量均滿足1.2節(jié)中的關(guān)系式。

表1 Licher翼型參數(shù)表Table 1 Parameters of Licher airfoils

圖9表明，優(yōu)化后升阻比顯著提高，且迭代15步算法基本收斂，說明該反設(shè)計方法具有極高的計算效率。表2統(tǒng)計了Licher雙翼優(yōu)化前后的氣動力參數(shù)，可以看出，采用矩形目標(biāo)壓力分布進(jìn)行優(yōu)化未能達(dá)到文獻(xiàn)［4］的優(yōu)化效果，對比發(fā)現(xiàn)本文優(yōu)化結(jié)果的波阻偏大。而階梯形分布的優(yōu)化結(jié)果具有高升力特性，顯著提高了翼型的升阻比，但阻力相對于矩形分布優(yōu)化結(jié)果有明顯增加。

圖9 基于Euler方程的優(yōu)化方法收斂性Fig.9 Convergence of the optimization method based on Euler equations

表2 基于Euler方程的優(yōu)化結(jié)果Table 2 Optimization results based on Euler Equations

為進(jìn)一步探究兩種目標(biāo)壓力分布下不同優(yōu)化效果的產(chǎn)生機(jī)理，分析優(yōu)化前后翼型周圍流場。圖10為優(yōu)化前流場，該圖表明:線化理論設(shè)計的Licher雙翼間的波系與理論估計有所差異，翼型前緣形成的激波打在了轉(zhuǎn)折點(diǎn)前方，形成了Cp圖上的第一個壓力峰值;第一個激波的前移造成了反射激波與膨脹波的錯位，反射的激波打在翼型后緣，形成了Cp圖的第二個壓力峰值;翼型后緣產(chǎn)生的膨脹波使得Cp值急劇減小。因此Licher雙翼未能達(dá)到預(yù)期的減阻效果。

圖10 Licher雙翼壓力分布Fig.10 Pressure distribution of the Licher airfoils

圖11 表明以矩形Cp分布為目標(biāo)優(yōu)化后，上下翼型被削薄，翼型內(nèi)側(cè)流動轉(zhuǎn)折點(diǎn)前移，前緣形成的激波打在轉(zhuǎn)折點(diǎn)處，消去了第一個Cp峰值;通過控制后緣翼型內(nèi)側(cè)形狀，使反射激波與膨脹波重合，消去第二個Cp峰值;同時改變翼型后緣曲率，降低了Cp曲線的衰減程度，從而達(dá)到減小阻力的作用。與矩形目標(biāo)壓力分布類似，階梯形目標(biāo)壓力分布也削平了Licher雙翼內(nèi)側(cè)的Cp峰值，減小了阻力，但翼型彎度較矩形優(yōu)化結(jié)果有所增加。

對比圖12(a)與圖12(b)，圖(a)Cp1、Cp3和Cp2、Cp4區(qū)域顏色基本相同，圖12(b)相對于圖12(a)，Cp2區(qū)云圖顏色變深，上翼內(nèi)側(cè)壓力系數(shù)增大，Cp4區(qū)云圖顏色變淺，下翼內(nèi)側(cè)壓力系數(shù)減小，該分布有效增加了雙翼升力，使升阻比進(jìn)一步提高，這與1.3節(jié)的理論分析的結(jié)果一致。

為探究階梯形目標(biāo)優(yōu)化結(jié)果阻力偏大的原因，分別將上下翼阻力系數(shù)列于表3中，可得，在矩形目標(biāo)分布下，下方翼型水平方向氣動力呈推力狀態(tài)，有效地減少了阻力。因此，矩形分布的優(yōu)化結(jié)果側(cè)重于低波阻，而階梯分布的優(yōu)化結(jié)果側(cè)重于高升力。

優(yōu)化結(jié)果僅給出設(shè)計迎角α=1°下的翼型升阻特性，保持來流馬赫數(shù)不變，分別計算三種翼型在0°～4°迎角下的升阻特性。由圖13可得:優(yōu)化設(shè)計的雙翼具有較好的穩(wěn)健性，在迎角0°～4°的變化范圍內(nèi)沒有出現(xiàn)流動壅塞現(xiàn)象［5，8-9］;同時，設(shè)計狀態(tài)下具有較高升阻比的翼型，在其它迎角下升阻比仍占據(jù)優(yōu)勢。

圖11 基于Euler方程的翼型優(yōu)化結(jié)果Fig.11 Optimization results based on Euler Equations

圖12 基于Euler方程優(yōu)化結(jié)果的壓力云圖Fig.12 Pressure contours of optimization airfoils based one Euler equations

表3 上下翼型阻力系數(shù)Table 3 Drag coefficient of airfoils

圖13 不同迎角下的升阻比Fig.13 Lift-drag ratio of different angles

2.2 基于NS方程的Licher雙翼優(yōu)化算例

為了探究粘性對雙翼性能的影響，使用基于NS方程的CFD求解器代替優(yōu)化框架中的Euler方程求解器，雷諾數(shù)Re=3.8×107，選擇SST湍流模型，翼型數(shù)據(jù)，來流馬赫數(shù)，來流迎角等初始條件均保持不變。

對比表2和表4的升阻比變化率可得，考慮粘性后，階梯分布優(yōu)化結(jié)果的優(yōu)勢更加明顯，相對矩形分布優(yōu)化的結(jié)果，其升阻比提高率達(dá)到了32.2%，遠(yuǎn)高于無粘情況下19.3%。為分析原因，將阻力分解成摩擦阻力和激波阻力(如圖14)。圖14表明:與菱形翼對比，雙翼構(gòu)型使激波阻力減小了約一個量級，但摩擦阻力卻增加了兩倍，階梯分布優(yōu)化結(jié)果的特點(diǎn)是高升力，而矩形分布優(yōu)化結(jié)果的特點(diǎn)是低波阻，考慮粘性后，摩擦阻力在阻力中占據(jù)主導(dǎo)的地位，減小波阻提高升阻比的效果受到抑制，因此矩形分布的優(yōu)化效果變差。但階梯形分布具有高升力特性，仍可保持良好的優(yōu)化效果，故以階梯形分布為目標(biāo)的反設(shè)計方法更具實(shí)用價值。

表4 基于NS方程的優(yōu)化結(jié)果Table 4 Optimization results based on NS Equations

圖14 阻力分解Fig.14 Drag analyses

圖15 為基于NS方程的翼型優(yōu)化結(jié)果，對比圖11中的無粘優(yōu)化結(jié)果可以看出，Cp分布規(guī)律基本保持不變，且階梯形目標(biāo)壓力分布優(yōu)化的翼型彎度仍大于矩形分布的優(yōu)化結(jié)果。但由于邊界層效應(yīng)，Licher雙翼內(nèi)側(cè)的兩個壓力峰明顯前移，導(dǎo)致優(yōu)化翼型的轉(zhuǎn)折點(diǎn)前移，此外，激波(膨脹波)邊界層干擾，邊界層流動轉(zhuǎn)捩等復(fù)雜流動現(xiàn)象使得計算出的翼型表面Cp值出現(xiàn)波動，影響了優(yōu)化穩(wěn)定性，削弱了減阻效果。

圖15 基于NS方程的翼型優(yōu)化結(jié)果Fig.15 Optimization results based on NS Equations

3 結(jié)論

采用反設(shè)計方法，結(jié)合CFD技術(shù)，對Licher雙翼進(jìn)行了優(yōu)化設(shè)計，得到以下結(jié)論:

1)反設(shè)計方法有效修正了線化理論和實(shí)際流動之間的差異，提高了Licher雙翼的升阻比;且該方法的計算效率極高，一般迭代不超過20步，即可得到收斂的結(jié)果。

2)目標(biāo)壓力分布是反設(shè)計方法中的關(guān)鍵因素，矩形目標(biāo)壓力分布可以有效減小激波阻力，來流馬赫數(shù)1.7、無粘情況下使升阻比提高30.5%。但該分布在翼型有迎角的情況下存在不足，本文提出階梯形目標(biāo)壓力分布，在保留消波減阻特性的基礎(chǔ)上，增加了雙翼升力，無粘假設(shè)下升阻比可提高49.8%。

3)當(dāng)上述優(yōu)化方法推廣到粘性情況下進(jìn)行設(shè)計時，由于摩擦阻力在阻力中占主導(dǎo)地位，矩形目標(biāo)壓力分布優(yōu)化效果受到抑制，來流馬赫數(shù)1.7時升阻比僅能提高17%，而階梯形分布的優(yōu)化效果受到影響較小，升阻比仍可提高49.2%，因此階梯形目標(biāo)壓力分布更具實(shí)用價值。

［1］Busemann A.Aerodynamic lift at supersonic speeds［C］.The 5th Volta Aerodynamic Conference，Volta，Italy，October 3，1935.

［2］Busemann A.The relation between minimized drag and noise at supersonic speed［C］.The High-speed Aeronautics Conference.Brooklyn，USA，January 20-22，1955.

［3］Xu Y Z，Xu Z Q，Li S G，et al.A hypersonic lift mechanism with decoupled lift and drag surfaces［J］.Science China:Physics，Mechanics＆Astronomy，2013，56(5):981-988.

［4］Kusunose K.A fundamental study for the development of boomless supersonic transport aircraft［C］.The 44th AIAA Aerospace Sciences Meeting and Exhibit，Reno，USA，January 9-12，2006.

［5］Maruyama D，Matsuzawa T.Consideration at off-design conditions of supersonic flows around biplane airfoils［C］.The 45th AIAA Aerospace Sciences Meeting and Exhibit，Reno，USA，January 8-11，2007.

［6］Yonezawa M，Yamashita H，Obayashi S.Investigation of supersonic wing shape using Busemann biplane airfoil［C］.The 45th AIAA Aerospace Sciences Meeting and Exhibit，Reno，USA，January 8-11，2007.

［7］Hua Ruhao，Ye Zhengyin.Drag reduction method for supersonic missile based on Busemann biplane concept［J］.Chinese Journal of Applied Mechanics，2012，29(5):536-540.(in Chinese)華如豪，葉正寅.基于Busemann雙翼構(gòu)型的超音速導(dǎo)彈減阻技術(shù)研究［J］.應(yīng)用力學(xué)學(xué)報，2012，29(5):536-540.

［8］Hu R，Jameson A.Wang Q Q.Adjoint-based aerodynamic optimization of supersonic biplane airfoils［J］.Journal of Aircraft，2012，49(3):802-814.

［9］Hu R.Supersonic biplane design via adjoint method［D］.Stanford: Institution of aeronautics and astronautics，2009.

［10］Licher R M.Optimum two-dimensional multiplanes in supersonic flow［C］.Douglass Aircraft Conference，USA，1955.

［11］Ye Youda，Luo Tianning，Lu Shen，et al.Application of optimization method based on dynamic boundary control in configuration optimization of high lift-to-drag vehicle［J］.Acta Aerodynamica Sinica，2009，27:29-33.(in Chinese)葉友達(dá)，羅天寧，盧笙，等.基于動態(tài)邊界控制的優(yōu)化方法在高升阻比氣動外形優(yōu)化中的初步應(yīng)用［J］.空氣動力學(xué)學(xué)報，2009，27:29-33.

［12］Chang Yanxing，Gao Zhenghong.Applicability analysis of genetic algorithms using in aerodynamic optimization［J］.Acta Aerodynamica Sinica，2009，27(2):167-172.(in Chinese)常彥鑫，高正紅.遺傳算法在氣動優(yōu)化問題中的適應(yīng)性分析［J］.空氣動力學(xué)學(xué)報，2009，27(2):167-172.

［13］Li Jing，Gao Zhenghong，Huang Jiangtao，et al.Airfoil optimization based on distributed panicle swarm algorithm［J］.Acta Aerodynamica Sinica，2011，29(4):464-469.(in Chinese)李靜，高正紅，黃江濤，等.基于分布式粒子群算法的翼型優(yōu)化設(shè)計［J］.空氣動力學(xué)學(xué)報，2011，29(4):464-469.

［14］Matsushima K，Takanashi S.An inverse design method for wings using integral equations and its recent progress［C］.Proceedings of Lectures and Workshop International Recent Advances in Multidisciplinary Technology and Modeling，Japan，2008.

［15］Kusunose K，Matsushima K，Maruyama D.Supersonic biplane—a review［J］.Progress in Aerospace Sciences，2011，47:53-87.

［16］Liepmann H W，Roshko A.Elements of gas dynamics［M］.New York:John Wiley＆Sons，Inc，1957.

［17］Yamashita H，Yonezawa M，Obayashi S，et al.A study of Busemann-type biplane for avoiding choked flow［C］.The 45th AIAA Aerospace Sciences Meeting and Exhibit，Reno，USA，2007.

Inverse design method for the Licher biplane with a new target pressure distribution

Zhao Chengxi*，Ye Zhengyin，Hua Ruhao
(National Key Laboratory of Science and Technology on Aerodynamic Design and Research，School of Aeronautics，Northwestern Polytechnical University，Xi’an 710072，China)

An inverse design method has been applied for the Licher biplane with modern computational fluid dynamics(CFD)technique to increase the lift-drag ratio and reduce the shock wave drag during its the supersonic flight regime.An example of numerical simulation demonstrates that the wave drag can be reduced by 25.5%and the lift-drag ratio rises by 30.5%simultaneously using this optimization algorithm with inviscid flow assumption as the free stream Mach number is 1.7.In addition，it takes about 15 iterations to get convergent results，which demonstrates that the inverse design method in this paper is efficient.A new stair target pressure distribution，as the key of the inverse design method，has been developed to overcome setbacks of the original one by some analyses of aerodynamic forces.With the new target，the lift-drag ratio can increase up to 49.8%.Furthermore，the Euler equations in the optimization procedure have been replaced by the Navier-Stokes equations to investigate the influence of viscous fluid.The results indicate that the effect of the optimization with the original target pressure distribution is beyond neglect，the lift-drag ratio increases only by 17.0%.However，the results of the new stair target pressure distribution，with more potential practical value，stay almost the same level.

Busemann concept;licher biplane;inverse design method;high lift-drag ratio;target pressure distribution

V211.3

:Adoi:10.7638/kqdlxxb-2014.0040

0258-1825(2015)05-0610-07

2014-05-19;

:2014-06-21

國家自然科學(xué)基金(11272262);航天支撐基金(12115001)

趙承熙*(1990-)，男，安徽六安人，碩士研究生，研究方向:計算流體力學(xué).E-mail:zhaochengxi1990@163.com

趙承熙，葉正寅，華如豪.新型目標(biāo)壓力分布下的Licher雙翼反設(shè)計方法研究［J］.空氣動力學(xué)學(xué)報，2015，33(5):610-616，630.

10.7638/kqdlxxb-2014.0040 Zhao C X，Ye Z Y，Hua R H.Inverse design method for the Licher biplane with a new target pressure distribution［J］.Acta Aerodynamica Sinica，2015，33(5):610-616，630.