跨聲速層流翼型的混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)方法

陳靜，宋文萍，朱震，許朕銘，韓忠華,*

1. 西北工業(yè)大學(xué) 航空學(xué)院 翼型葉柵空氣動(dòng)力學(xué)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室，西安 710072 2. 武漢第二船舶設(shè)計(jì)研究所，武漢 430064

對(duì)于現(xiàn)代大型客機(jī)，接近50%氣動(dòng)阻力來(lái)自流動(dòng)邊界層內(nèi)黏性效應(yīng)，而層流邊界層內(nèi)黏性阻力遠(yuǎn)小于湍流邊界層。因此，旨在保持大范圍層流邊界層的層流技術(shù)具有非常巨大的減阻潛力，被美國(guó)《航空周刊》評(píng)選成為未來(lái)20～40年內(nèi)最有前景的前沿方向之一。

層流技術(shù)涵蓋自然層流技術(shù)、層流控制技術(shù)和混合層流控制技術(shù)3種。自然層流技術(shù)，作為其他層流技術(shù)的基礎(chǔ)，研究歷史已有80多年。2005年Honda Jet噴氣公務(wù)機(jī)[1-2]的問(wèn)世，標(biāo)志著自然層流技術(shù)在小型商業(yè)飛機(jī)上首次成功應(yīng)用。而波音737的層流吊艙[3]以及波音787 Max的層流翼尖小翼[4]則開(kāi)啟了自然層流技術(shù)在現(xiàn)代民用客機(jī)上局部應(yīng)用的先河。

針對(duì)自然層流翼型/機(jī)翼設(shè)計(jì)，基于CFD流動(dòng)的數(shù)值優(yōu)化方法被廣泛應(yīng)用于氣動(dòng)設(shè)計(jì)中，包括反設(shè)計(jì)和直接優(yōu)化設(shè)計(jì)兩種。Green[5]和Gopalarathnam[6]等均采用反設(shè)計(jì)方法開(kāi)展了自然層流翼型設(shè)計(jì)工作。Campbell和Lynde[7]發(fā)展了考慮流動(dòng)和幾何約束的自然層流機(jī)翼反設(shè)計(jì)方法，并初步開(kāi)展了適用于大型客機(jī)的層流機(jī)翼設(shè)計(jì)工作。德國(guó)宇航院在LamAiR項(xiàng)目[8]中以壓力分布反設(shè)計(jì)為主要方法，設(shè)計(jì)了跨聲速自然層流前掠機(jī)翼，其上下表面均可實(shí)現(xiàn)大范圍層流流動(dòng)。Lee和Jameson[9]進(jìn)行了自然層流翼型目標(biāo)壓力分布反設(shè)計(jì)，并以設(shè)計(jì)翼型為剖面配置三維機(jī)翼，進(jìn)行直接減阻優(yōu)化，但由于其梯度求解沒(méi)有考慮轉(zhuǎn)捩的影響，減阻效果主要來(lái)自激波的消除。Amoignon等[10]在跨聲速狀態(tài)下，將邊界層內(nèi)擾動(dòng)動(dòng)能作為直接優(yōu)化的最小化目標(biāo)函數(shù)，推遲了翼型層流-湍流轉(zhuǎn)捩，但翼面激波仍較強(qiáng)。Zingg團(tuán)隊(duì)[11-13]開(kāi)展了系列單設(shè)計(jì)點(diǎn)和多設(shè)計(jì)點(diǎn)自然層流翼型直接優(yōu)化設(shè)計(jì)工作。Cameron等[14]開(kāi)展了基于代理模型的自然層流翼型多目標(biāo)優(yōu)化設(shè)計(jì)。國(guó)內(nèi)對(duì)層流翼型設(shè)計(jì)也開(kāi)展了大量研究，華俊等[15]采用“正反迭代、余量修正”原理進(jìn)行層流翼型反設(shè)計(jì)工作；喬志德[16]采用人-機(jī)對(duì)話和壓力分布反設(shè)計(jì)的方法設(shè)計(jì)了自然層流超臨界翼型，首次提出了利用弱激波形成足夠順壓梯度的自然層流超臨界翼型設(shè)計(jì)思想。王迅等[17]利用本征正交分解(POD)代理模型實(shí)現(xiàn)快速轉(zhuǎn)捩判斷，并進(jìn)行翼型減阻優(yōu)化設(shè)計(jì)；韓忠華等[18-20]基于代理優(yōu)化算法分別開(kāi)展了跨聲速層流翼型和機(jī)翼減阻直接優(yōu)化設(shè)計(jì)。馬曉永等[21]采用序列二次規(guī)劃優(yōu)化算法開(kāi)展了某民用通用飛機(jī)自然層流機(jī)翼減阻設(shè)計(jì)。同時(shí)，為了提高非設(shè)計(jì)點(diǎn)氣動(dòng)性能，張宇飛等[22]研究了不同順壓梯度約束下優(yōu)化翼型氣動(dòng)特性對(duì)雷諾數(shù)等參數(shù)的敏感性，驗(yàn)證了優(yōu)化數(shù)學(xué)模型的選擇對(duì)于設(shè)計(jì)結(jié)果性能魯棒的重要性。高正紅等[23-25]開(kāi)展了大量層流翼型減阻設(shè)計(jì)和穩(wěn)健設(shè)計(jì)工作。國(guó)內(nèi)還開(kāi)展了大量層流翼型多目標(biāo)優(yōu)化設(shè)計(jì)研究[26-27]。

可見(jiàn)，在跨聲速層流翼型/機(jī)翼設(shè)計(jì)方面，國(guó)內(nèi)外基于反設(shè)計(jì)或直接優(yōu)化設(shè)計(jì)方法開(kāi)展了大量研究。在反設(shè)計(jì)方法中，需預(yù)設(shè)目標(biāo)流場(chǎng)特征(如壓力分布、擾動(dòng)增長(zhǎng)曲線等)，只有在設(shè)計(jì)領(lǐng)域積累了豐富經(jīng)驗(yàn)的專家才能根據(jù)設(shè)計(jì)指標(biāo)給出合適的目標(biāo)流場(chǎng)分布，因此，反設(shè)計(jì)方法只能局限于具有相當(dāng)豐富積累的領(lǐng)域，而當(dāng)新技術(shù)更新加速時(shí)，反設(shè)計(jì)方法很難適應(yīng)于工程設(shè)計(jì)。在直接優(yōu)化設(shè)計(jì)方法中，設(shè)定優(yōu)化目標(biāo)(如最小化阻力、最大化升阻比等)和約束，并采用數(shù)值優(yōu)化算法(如遺傳算法或梯度算法等)得到優(yōu)化解，該方法依賴于優(yōu)化數(shù)學(xué)模型的選擇，可能會(huì)過(guò)度追求優(yōu)化目標(biāo)最優(yōu)，而不能全面兼顧其他性能，針對(duì)該不足發(fā)展起來(lái)的多目標(biāo)優(yōu)化和魯棒優(yōu)化設(shè)計(jì)等方法計(jì)算量巨大，很難在工程設(shè)計(jì)中實(shí)現(xiàn)。因此，發(fā)展能夠綜合反設(shè)計(jì)和直接優(yōu)化設(shè)計(jì)方法優(yōu)點(diǎn)而避開(kāi)各自缺點(diǎn)的氣動(dòng)優(yōu)化設(shè)計(jì)新方法具有重要意義。

本文將局部反設(shè)計(jì)目標(biāo)，即基于經(jīng)驗(yàn)/物理機(jī)理的局部目標(biāo)流動(dòng)分布特性，與直接優(yōu)化設(shè)計(jì)目標(biāo)加權(quán)組合，實(shí)現(xiàn)基于代理優(yōu)化算法的混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)方法。針對(duì)適用于現(xiàn)代民機(jī)的跨聲速層流機(jī)翼，采用本文發(fā)展方法開(kāi)展基于代理優(yōu)化算法和耦合轉(zhuǎn)捩自動(dòng)判定的雷諾平均Navier-Stokes(RANS)方程求解器的跨聲速層流翼型氣動(dòng)設(shè)計(jì)研究，驗(yàn)證本文方法的有效性。進(jìn)一步將設(shè)計(jì)翼型配置到三維機(jī)翼上，檢驗(yàn)其氣動(dòng)性能改進(jìn)效果，并進(jìn)行設(shè)計(jì)結(jié)果減阻機(jī)理分析。

1 基于混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)的氣動(dòng)優(yōu)化設(shè)計(jì)方法

1.1 數(shù)學(xué)模型

混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)方法旨在兼顧局部流場(chǎng)特性和總體性能指標(biāo)，形成實(shí)用的設(shè)計(jì)結(jié)果。在單次優(yōu)化設(shè)計(jì)過(guò)程中，既引入部分設(shè)計(jì)經(jīng)驗(yàn)或設(shè)計(jì)準(zhǔn)則以指導(dǎo)優(yōu)化設(shè)計(jì)過(guò)程，確保設(shè)計(jì)結(jié)果的實(shí)用性和魯棒性，又通過(guò)優(yōu)化算法搜索設(shè)計(jì)空間，得到設(shè)計(jì)指標(biāo)最優(yōu)的氣動(dòng)外形。本文中，反設(shè)計(jì)問(wèn)題和直接優(yōu)化問(wèn)題被分別表達(dá)為單目標(biāo)優(yōu)化問(wèn)題，并通過(guò)權(quán)重系數(shù)組合成多目標(biāo)優(yōu)化問(wèn)題，其數(shù)學(xué)模型表示為

min obj=w1objinv+w2objdir

w.r.t.x∈(xl,xu)

(1)

式中：objinv表示反設(shè)計(jì)目標(biāo)函數(shù)，objdir表示直接優(yōu)化設(shè)計(jì)目標(biāo)函數(shù)，通過(guò)權(quán)重系數(shù)w1、w2組合成混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)總目標(biāo)obj；hp為氣動(dòng)或幾何相關(guān)等式約束；Np為等式約束數(shù)目；gm為氣動(dòng)或幾何相關(guān)不等式約束；Nm為不等式約束數(shù)目；設(shè)計(jì)變量為x序列，xl和xu分別為設(shè)計(jì)變量的變化上下邊界。

當(dāng)w1=1、w2=0時(shí)，該模型為直接優(yōu)化形式下的反設(shè)計(jì)問(wèn)題；當(dāng)w1=0、w2=1時(shí)，該模型為直接優(yōu)化設(shè)計(jì)問(wèn)題；當(dāng)w1、w2在0～1之間時(shí)，該數(shù)學(xué)模型為綜合直接優(yōu)化設(shè)計(jì)和反設(shè)計(jì)的混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)問(wèn)題。針對(duì)不同的氣動(dòng)優(yōu)化問(wèn)題，反設(shè)計(jì)目標(biāo)objinv是根據(jù)設(shè)計(jì)經(jīng)驗(yàn)或物理背景等確定的流場(chǎng)特性分布函數(shù)，如局部壓力分布、載荷分布或某種氣動(dòng)設(shè)計(jì)準(zhǔn)則等，而直接優(yōu)化設(shè)計(jì)目標(biāo)objdir通常為升阻特性或其他有待改善的性能指標(biāo)。通過(guò)將混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)問(wèn)題構(gòu)造為上述多目標(biāo)優(yōu)化問(wèn)題，確保優(yōu)化外形具備目標(biāo)流場(chǎng)特征，且性能指標(biāo)得到優(yōu)化。

1.2 幾何參數(shù)化及設(shè)計(jì)空間選擇

采用基于型函數(shù)/類函數(shù)變換的參數(shù)化方法(CST)描述翼型幾何形狀。該方法最早由Kulfan[28]提出，表達(dá)式為

(2)

(3)

采用擾動(dòng)CST方法得到描述翼型的設(shè)計(jì)參數(shù)后，還需確定各參數(shù)在優(yōu)化設(shè)計(jì)中的取值范圍，從而確定設(shè)計(jì)空間大小。在優(yōu)化前并不知道最優(yōu)外形所在設(shè)計(jì)空間，因此，先以基準(zhǔn)翼型為設(shè)計(jì)空間中心，給定設(shè)計(jì)參數(shù)范圍，如果優(yōu)化翼型設(shè)計(jì)參數(shù)達(dá)到上下邊界，則將優(yōu)化翼型作為新的設(shè)計(jì)空間中心，開(kāi)展新一輪優(yōu)化設(shè)計(jì)，直到優(yōu)化翼型設(shè)計(jì)參數(shù)全部被包含在設(shè)計(jì)空間內(nèi)。

1.3 代理優(yōu)化設(shè)計(jì)

代理優(yōu)化算法作為一種新型的高效優(yōu)化設(shè)計(jì)方法，在航空航天界受到了重視，并快速發(fā)展起來(lái)。該方法通過(guò)建立一種計(jì)算量小的近似數(shù)學(xué)模型“代替”原有的高精度數(shù)值模擬，大幅降低計(jì)算量和計(jì)算時(shí)間，實(shí)現(xiàn)高效的數(shù)值優(yōu)化。本文在代理優(yōu)化器“SurroOpt”[30-32]基礎(chǔ)上，實(shí)現(xiàn)了基于混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)方法的氣動(dòng)優(yōu)化設(shè)計(jì)框架，如圖1所示，具體步驟如下：

圖1 代理優(yōu)化設(shè)計(jì)方法流程圖Fig.1 Flow chart of surrogate optimization design method

步驟1采用試驗(yàn)設(shè)計(jì)方法(DoE，如拉丁超立方等)在給定設(shè)計(jì)空間內(nèi)隨機(jī)抽樣產(chǎn)生初始樣本點(diǎn)。

步驟2調(diào)用CFD分析程序，計(jì)算混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)目標(biāo)函數(shù)值及約束函數(shù)值。

步驟3基于已有樣本數(shù)據(jù)分別建立目標(biāo)函數(shù)和約束函數(shù)的代理模型(如Kriging等)。

步驟4采用不同的自適應(yīng)加點(diǎn)技術(shù)(如Expected Improvement(EI)、Minimizing Surrogate Prediction(MSP)、Probability of Improvement(PI)、Lower Confidence Bound(LCB)等)在代理模型上進(jìn)行子優(yōu)化，增加新樣本點(diǎn)。

步驟5調(diào)用CFD分析程序，計(jì)算新增樣本點(diǎn)混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)目標(biāo)函數(shù)值及約束值，更新代理模型，重復(fù)步驟4和步驟5，直到優(yōu)化過(guò)程終止于優(yōu)化收斂準(zhǔn)則。

2 流動(dòng)數(shù)值模擬方法及驗(yàn)證

2.1 耦合轉(zhuǎn)捩自動(dòng)判定的RANS流場(chǎng)求解器

采用耦合轉(zhuǎn)捩自動(dòng)判定的RANS流場(chǎng)求解器(PMNS2D/3D)[33-34]開(kāi)展翼型/機(jī)翼繞流數(shù)值模擬計(jì)算。RANS求解器與轉(zhuǎn)捩預(yù)測(cè)模塊的耦合框架如圖2所示。RANS主流場(chǎng)求解模塊采用S-A(Spalart-Allmaras)湍流模型,中心格式空間離散，LU-SGS(Lower-Upper Symmeffic Gauss-Seidel)隱式時(shí)間推進(jìn)；轉(zhuǎn)捩預(yù)測(cè)模塊包括可提供高精度邊界層信息的層流邊界層求解和基于線性穩(wěn)定性理論的eN方法兩部分。在eN方法中，擾動(dòng)放大因子增長(zhǎng)到預(yù)設(shè)閾值時(shí)，判定層流-湍流轉(zhuǎn)捩發(fā)生。

圖2 RANS求解器與轉(zhuǎn)捩模塊耦合策略示意圖Fig.2 Coupling scheme of RANS solver and transition module

2.2 流場(chǎng)求解器驗(yàn)證

2.2.1 NLF(1)-0416翼型繞流自由轉(zhuǎn)捩計(jì)算

對(duì)NLF(1)-0416翼型開(kāi)展考慮轉(zhuǎn)捩影響的數(shù)值模擬計(jì)算。在計(jì)算狀態(tài)馬赫數(shù)Ma=0.3,雷諾數(shù)Re=4×106,迎角α=2.03°下，層流邊界層內(nèi)位移厚度計(jì)算值與文獻(xiàn)[9]吻合較好，如圖3所示，驗(yàn)證了求解器中邊界層求解程序的可靠性。在計(jì)算狀態(tài)Ma=0.1,Re=4×106下，采用eN方法對(duì)NLF(1)-0416翼型邊界層進(jìn)行轉(zhuǎn)捩判斷。對(duì)于半經(jīng)驗(yàn)的eN方法，航空界經(jīng)過(guò)大量研究，綜合考慮影響轉(zhuǎn)捩的各種因素，將轉(zhuǎn)捩放大因子閥值(NTS)tr取為9，作為設(shè)計(jì)參考值。然而在具體風(fēng)洞試驗(yàn)中，試驗(yàn)條件、模型加工精度等差異會(huì)導(dǎo)致不同試驗(yàn)標(biāo)定得到的轉(zhuǎn)捩放大因子取值存在差異。本驗(yàn)證算例中(NTS)tr=11，參考德國(guó)宇航院對(duì)該試驗(yàn)數(shù)據(jù)的標(biāo)定結(jié)果[35]。轉(zhuǎn)捩放大因子(NTS)tr=9，11時(shí)，NLF(1)-0416翼型0°迎角下壓力分布與轉(zhuǎn)捩位置Xtr對(duì)比見(jiàn)圖4，可見(jiàn)兩個(gè)不同轉(zhuǎn)捩放大因子值計(jì)算得到的壓力系數(shù)Cp分布和所預(yù)測(cè)的轉(zhuǎn)捩位置差異不大。

圖5是上表面Tollmien-Schlichting (TS)波擾動(dòng)NTS增長(zhǎng)曲線，可以看出TS波擾動(dòng)因子(NTS)tr增長(zhǎng)劇烈，曲線陡峭。針對(duì)此特點(diǎn)，(NTS)tr在合理范圍內(nèi)取值時(shí)，不會(huì)對(duì)TS波誘導(dǎo)轉(zhuǎn)捩位置有較大影響。不同迎角下，翼型上表面轉(zhuǎn)捩位置預(yù)測(cè)值與試驗(yàn)結(jié)果的對(duì)比如圖6所示?？梢?jiàn)(NTS)tr=9，11時(shí)，所預(yù)測(cè)的轉(zhuǎn)捩位置差異不大，均基本在試驗(yàn)值合理范圍內(nèi)。根據(jù)參考文獻(xiàn)標(biāo)定取11時(shí)，計(jì)算結(jié)果更接近試驗(yàn)值，表明本文轉(zhuǎn)捩判定方法能準(zhǔn)確預(yù)測(cè)翼型表面層流-湍流轉(zhuǎn)捩。

圖3 NLF(1)-0416翼型上表面位移厚度δ*比較Fig.3 Comparison of displacement thickness δ* for NLF(1)-0416 airfoil upper surface

圖4 NLF(1)-0416翼型壓力分布和轉(zhuǎn)捩位置比較Fig.4 Comparison of pressure distribution and transition locations on NLF(1)-0416 airfoil

圖5 NLF(1)-0416翼型上表面TS波擾動(dòng)增長(zhǎng)曲線Fig.5 Disturbance growth curves of TS wave on NLF(1)-0416 airfoil upper surface

圖6 NLF(1)-0416翼型上表面轉(zhuǎn)捩位置與試驗(yàn)結(jié)果比較Fig.6 Comparison of transition locations on NLF(1)-0416 airfoil upper surface with experiment results

2.2.2 NPU-LSC-72613翼型跨聲速繞流自由轉(zhuǎn)捩計(jì)算

本節(jié)對(duì)NPU-LSC-72613[36]翼型開(kāi)展考慮轉(zhuǎn)捩影響的數(shù)值模擬計(jì)算，驗(yàn)證跨聲速流動(dòng)狀態(tài)下本文使用的轉(zhuǎn)捩預(yù)測(cè)技術(shù)的可靠性。由于現(xiàn)有公開(kāi)資料中，沒(méi)有帶自然轉(zhuǎn)捩的跨聲速翼型試驗(yàn)數(shù)據(jù),因此本文通過(guò)與MSES軟件結(jié)果對(duì)比來(lái)驗(yàn)證跨聲速狀態(tài)下轉(zhuǎn)捩判斷的可靠性。MSES是一款由美國(guó)Drela教授開(kāi)發(fā)，業(yè)界公認(rèn)較為成熟的軟件，通過(guò)求解Euler方程耦合eN轉(zhuǎn)捩判斷方法進(jìn)行翼型繞流模擬及轉(zhuǎn)捩判斷。計(jì)算狀態(tài)Ma=0.72，Re=2.0×107，升力系數(shù)CL=0.6。圖7是本文方法計(jì)算得到的壓力分布及上下表面轉(zhuǎn)捩位置與MSES軟件結(jié)果的對(duì)比，壓力分布和上下表面轉(zhuǎn)捩位置吻合較好，驗(yàn)證了本文方法在跨聲速流動(dòng)狀態(tài)下對(duì)翼型邊界層轉(zhuǎn)捩預(yù)測(cè)的可靠性。

圖7 NPU-LSC-72613翼型壓力分布和轉(zhuǎn)捩位置比較Fig.7 Comparison of pressure distribution and transition locations on NPU-LSC-72613 airfoil

2.2.3 無(wú)限展長(zhǎng)后掠機(jī)翼NACA642A015繞流自由轉(zhuǎn)捩計(jì)算

以NACA642A015翼型為橫剖面的后掠機(jī)翼在美國(guó)Ames 12英尺的低湍流度風(fēng)洞進(jìn)行了試驗(yàn)[37]。對(duì)于后掠角40°機(jī)翼，采用雙eN方法[38]判斷轉(zhuǎn)捩，對(duì)于風(fēng)洞試驗(yàn)中無(wú)限翼展后掠機(jī)翼邊界層轉(zhuǎn)捩，本文參考文獻(xiàn)[38]中標(biāo)定結(jié)果，TS波和CF波轉(zhuǎn)捩擾動(dòng)放大因子閥值取[NTS,NCF]tr=[10.5,7.5]。圖8給出了Ma=0.27,Re=5.7×106，α=-1.5°狀態(tài)下，計(jì)算壓力分布與試驗(yàn)值對(duì)比。圖9展示了不同雷諾數(shù)下，機(jī)翼上表面計(jì)算轉(zhuǎn)捩位置與試驗(yàn)值比較，可見(jiàn)，計(jì)算值與試驗(yàn)值吻合較好。

本文采用的機(jī)翼轉(zhuǎn)捩判斷方法為雙eN方法，在判斷轉(zhuǎn)捩位置的同時(shí)可以判定轉(zhuǎn)捩是由CF波不穩(wěn)定性導(dǎo)致還是TS波不穩(wěn)定性導(dǎo)致。圖10給出NACA642A015后掠機(jī)翼上表面TS波和CF波擾動(dòng)因子N增長(zhǎng)曲線，從圖中可以看出CF波擾動(dòng)增長(zhǎng)曲線先于TS波達(dá)到對(duì)應(yīng)閥值，邊界層轉(zhuǎn)捩由CF波增長(zhǎng)導(dǎo)致。對(duì)于40°后掠機(jī)翼，由CF波不穩(wěn)定性導(dǎo)致轉(zhuǎn)捩是合理的，轉(zhuǎn)捩位置與試驗(yàn)值也吻合較好，驗(yàn)證了本文RANS求解方法在三維機(jī)翼上轉(zhuǎn)捩預(yù)測(cè)的可靠性。

圖8 NACA642A015后掠機(jī)翼壓力分布Fig.8 Pressure distribution on NACA642A015 swept wing

圖9 NACA642A015后掠機(jī)翼上表面轉(zhuǎn)捩位置Fig.9 Transition locations on NACA642A015 swept wing upper surface

圖10 NACA642A015后掠機(jī)翼上表面TS波和 CF波擾動(dòng)增長(zhǎng)曲線Fig.10 Disturbance growth curves of TS and CF wave on NACA642A015 swept wing upper surface

3 方法驗(yàn)證

3.1 跨聲速自然層流翼型優(yōu)化設(shè)計(jì)

針對(duì)適用于現(xiàn)代中短程民機(jī)的跨聲速層流機(jī)翼設(shè)計(jì)要求，開(kāi)展主翼型氣動(dòng)設(shè)計(jì)研究。機(jī)翼巡航狀態(tài)為Ma=0.75，Re=2.0×107，CL=0.5，則主翼型設(shè)計(jì)狀態(tài)為Ma=0.729，Re=2.0×107，CL=0.556，轉(zhuǎn)捩放大因子閥值取(NTS)tr=9。本節(jié)選取翼型NPU-LSC-72613為基準(zhǔn)外形。

跨聲速層流翼型氣動(dòng)設(shè)計(jì)要求兼顧優(yōu)良的超臨界特性和自然層流特性。在跨聲速狀態(tài)下，前者要求翼型表面激波盡可能弱，以減小激波引起的壓差阻力；后者則要求翼面保持穩(wěn)定的順壓梯度，以抑制邊界層內(nèi)不穩(wěn)定擾動(dòng)波劇烈增長(zhǎng)，保持大范圍的自然層流流動(dòng)。

綜合上述氣動(dòng)設(shè)計(jì)準(zhǔn)則，設(shè)計(jì)了具有穩(wěn)定、較小順壓梯度的翼型上表面局部壓力分布，作為混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)方法中的反設(shè)計(jì)目標(biāo)壓力分布，如圖11 所示。同時(shí)，將翼型總阻力CD作為混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)方法中的直接優(yōu)化設(shè)計(jì)目標(biāo)。

優(yōu)化過(guò)程中，采用8階擾動(dòng)CST參數(shù)化方法，設(shè)計(jì)變量18個(gè)，氣動(dòng)約束為：① 保持升力系數(shù)CL不小于基準(zhǔn)外形；② 力矩系數(shù)Cm不小于基準(zhǔn)外形；幾何約束為翼型相對(duì)厚度t/c不減；優(yōu)化數(shù)學(xué)模型如式(4)所示。優(yōu)化前后翼型氣動(dòng)性能計(jì)算采用C型結(jié)構(gòu)化網(wǎng)格，如圖12 所示。

min obj=w1objinv+w2objdir=

w.r.t.x∈(-0.02,0.02)

s.t.CL≥0.556;Cm≥-0.139;t/c≥12.7%

(4)

圖11 局部目標(biāo)壓力分布與基準(zhǔn)翼型壓力分布Fig.11 Distributions of pressure of local target and baseline airfoil

圖12 C型結(jié)構(gòu)化計(jì)算網(wǎng)格(641×241)Fig.12 C type structured computation grid (641×241)

3.2 混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)

本算例優(yōu)化設(shè)計(jì)進(jìn)行了兩輪，第1輪設(shè)計(jì)中，以基準(zhǔn)翼型為設(shè)計(jì)空間中心，權(quán)重系數(shù)w1、w2分別取0.3和0.7，經(jīng)過(guò)200次CFD計(jì)算后，混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)過(guò)程收斂，收斂歷程曲線如圖13所示?？梢?jiàn)，優(yōu)化設(shè)計(jì)目標(biāo)經(jīng)歷了一次顯著下降，且反設(shè)計(jì)和直接優(yōu)化設(shè)計(jì)目標(biāo)均為下降趨勢(shì)，隨后目標(biāo)函數(shù)值沒(méi)有進(jìn)一步下降。圖14 給出了當(dāng)前最優(yōu)設(shè)計(jì)點(diǎn)的壓力分布曲線，及其與局部目標(biāo)壓力分布和基準(zhǔn)翼型壓力分布的比較。相比于基準(zhǔn)翼型，當(dāng)前最優(yōu)設(shè)計(jì)的吸力面激波減弱，局部壓力分布更接近目標(biāo)分布，但仍然存在明顯差別，需要進(jìn)一步提高反設(shè)計(jì)部分的優(yōu)化效果。經(jīng)校驗(yàn)，當(dāng)前最優(yōu)點(diǎn)處設(shè)計(jì)變量值有77.8%達(dá)到給定設(shè)計(jì)空間邊界，說(shuō)明設(shè)計(jì)空間外可能存在更優(yōu)的氣動(dòng)外形，需調(diào)整當(dāng)前設(shè)計(jì)空間。

圖13 第1輪混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)收斂歷程Fig.13 Convergence history of first round hybrid inverse/direct optimization design

圖14 第1輪設(shè)計(jì)優(yōu)化翼型、基準(zhǔn)翼型和 局部目標(biāo)壓力分布Fig.14 Pressure distributions of optimal airfoil, baseline airfoil, and local target of first round design

通過(guò)將第1輪設(shè)計(jì)最優(yōu)點(diǎn)作為擾動(dòng)中心，改變?cè)O(shè)計(jì)空間，開(kāi)始第2輪混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)過(guò)程。同時(shí)根據(jù)第1輪優(yōu)化結(jié)果分析，權(quán)重系數(shù)w1、w2調(diào)整為0.11和0.89。圖15 給出了第2輪設(shè)計(jì)過(guò)程的優(yōu)化收斂歷程，可見(jiàn)，直接優(yōu)化設(shè)計(jì)目標(biāo)先是小幅上升后又小幅降低，反設(shè)計(jì)目標(biāo)有多次顯著下降。

經(jīng)過(guò)2輪混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)后，圖16給出了優(yōu)化翼型壓力分布，及其與局部目標(biāo)壓力分布和基準(zhǔn)翼型壓力分布的比較?？梢?jiàn)，優(yōu)化翼型的吸力面局部壓力分布與目標(biāo)壓力分布吻合很好，反設(shè)計(jì)效果顯著，而與基準(zhǔn)翼型相比，優(yōu)化翼型吸力面順壓力梯度減弱，激波位置不變，但強(qiáng)度減弱，有利于總阻力降低。圖17給出了優(yōu)化前后翼型幾何形狀比較，優(yōu)化翼型前緣半徑增大，吸力面前緣加速更快，而相對(duì)厚度基本保持不變。圖18 給出了翼型優(yōu)化前后，馬赫數(shù)在流場(chǎng)中的分布云圖對(duì)比，優(yōu)化翼型保持了大范圍的超聲速區(qū)，但激波前馬赫數(shù)小于基準(zhǔn)翼型，激波強(qiáng)度減弱。

圖15 第2輪混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)收斂歷程Fig.15 Convergence history of second round hybrid inverse/direct optimization design

圖16 第2輪設(shè)計(jì)優(yōu)化翼型、基準(zhǔn)翼型及 局部目標(biāo)壓力分布Fig.16 Pressure distributions of optimum airfoil, baseline airfoil and local target of second round design

圖17 優(yōu)化翼型與基準(zhǔn)翼型幾何形狀對(duì)比Fig.17 Comparison of geometry of optimum airfoil and baseline airfoil

圖18 優(yōu)化翼型和基準(zhǔn)翼型馬赫數(shù)云圖Fig.18 Mach number contours of optimum and baseline airfoils

表1給出了優(yōu)化前后翼型氣動(dòng)性能的比較，可見(jiàn)優(yōu)化翼型升力系數(shù)、力矩系數(shù)和相對(duì)厚度嚴(yán)格滿足約束，優(yōu)化翼型總阻力CD下降6 counts(1 count=0.000 1)， 約占基準(zhǔn)翼型阻力的15.5%。 從表中給出優(yōu)化前后翼型阻力分解值及轉(zhuǎn)捩位置可見(jiàn)，優(yōu)化翼型激波阻力CDw降低了5.4 counts， 黏性阻力CDv降低了0.6 count，且上下表面轉(zhuǎn)捩位置(Trans_x)均在55%倍弦長(zhǎng)之后。結(jié)果表明，經(jīng)過(guò)2輪混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)，翼型激波強(qiáng)度顯著降低，具有較好的超臨界特性，并保持了優(yōu)良的自然層流特性，說(shuō)明混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)方法效果顯著。

表1 優(yōu)化翼型與基準(zhǔn)翼型氣動(dòng)性能比較Table 1 Comparison of aerodynamic performances ofoptimum airfoil and baseline airfoil

4 配置基準(zhǔn)翼型和優(yōu)化翼型的跨聲速自然層流機(jī)翼氣動(dòng)性能對(duì)比

為驗(yàn)證所設(shè)計(jì)翼型對(duì)跨聲速自然層流機(jī)翼性能提升的效果，針對(duì)平面幾何形狀如圖19所示的跨聲速層流機(jī)翼進(jìn)行了氣動(dòng)性能對(duì)比計(jì)算。機(jī)翼翼根扭轉(zhuǎn)角為0°，翼尖扭轉(zhuǎn)角為-2°。沿順氣流方向分別配置優(yōu)化翼型和基準(zhǔn)翼型，得到設(shè)計(jì)機(jī)翼和初始機(jī)翼。巡航狀態(tài)Ma=0.75，Re=2.0×107，CL=0.5，對(duì)兩機(jī)翼開(kāi)展考慮自由轉(zhuǎn)捩的繞流數(shù)值模擬計(jì)算。本文設(shè)計(jì)的算例主要針對(duì)于民用客機(jī)，其飛行高度在萬(wàn)米左右，大氣環(huán)境十分安靜，擾動(dòng)較小，擾動(dòng)放大因子臨界值相對(duì)于低空環(huán)境偏高。在歐洲大量飛行試驗(yàn)[39]基礎(chǔ)上，對(duì)試驗(yàn)對(duì)象進(jìn)行穩(wěn)定性分析及擾動(dòng)放大因子臨界值進(jìn)行標(biāo)定，得出了在實(shí)際飛行狀態(tài)(區(qū)別于實(shí)驗(yàn)室環(huán)境)中，層流-湍流轉(zhuǎn)捩的擾動(dòng)放大因子臨界值[NTS,NCF]tr=[10.5,10]，具有參考意義。在優(yōu)化設(shè)計(jì)時(shí)，為了能夠設(shè)計(jì)出對(duì)環(huán)境擾動(dòng)更加魯棒的氣動(dòng)外形，本文選擇擾動(dòng)放大因子臨界值時(shí)乘以相應(yīng)的安全系數(shù)，轉(zhuǎn)捩擾動(dòng)因子閥值取[NTS,NCF]tr=[9,8.5]。取值更低的閥值可以保證設(shè)計(jì)外形更加魯棒，在來(lái)流湍流度擾動(dòng)時(shí)，機(jī)翼表面始終能夠獲得足夠的層流范圍。計(jì)算網(wǎng)格為C-H 型結(jié)構(gòu)化網(wǎng)格，如圖20所示。

圖19 跨聲速層流機(jī)翼平面形狀Fig.19 Planform of transonic laminar flow wing

圖20 機(jī)翼繞流計(jì)算C-H型結(jié)構(gòu)化網(wǎng)格(257×105×73)Fig.20 C-H type structured computation grid for wing flow calculation (257×105×73)

表2給出了設(shè)計(jì)機(jī)翼與初始機(jī)翼氣動(dòng)性能比較，相比于初始機(jī)翼，設(shè)計(jì)機(jī)翼阻力系數(shù)下降8 counts。力矩參考點(diǎn)在平均氣動(dòng)弦的25%弦長(zhǎng)處，力矩系數(shù)絕對(duì)值較小，升阻比由38.88提高到41.46，提升了約6.64%。機(jī)翼優(yōu)化效果沒(méi)有表1中翼型優(yōu)化氣動(dòng)性能的提升明顯，下面進(jìn)行詳細(xì)分析。從阻力組成上看，二維翼型阻力包括黏性壓差阻力、激波阻力和摩擦阻力3種；三維機(jī)翼阻力則包括誘導(dǎo)阻力、黏性壓差阻力、激波阻力和摩擦阻力4種。二維翼型相較于三維機(jī)翼，其激波阻力在總阻力中占比更高。本文優(yōu)化目標(biāo)是保持大范圍層流的基礎(chǔ)上，減弱激波，減少總阻力。從優(yōu)化結(jié)果看，優(yōu)化翼型激波明顯減弱，而層流流動(dòng)保持很好，并沒(méi)有發(fā)生較大變化，總阻力降低的主要貢獻(xiàn)來(lái)自激波阻力的下降。在三維機(jī)翼中，通過(guò)比較優(yōu)化機(jī)翼自然轉(zhuǎn)捩和固定轉(zhuǎn)捩計(jì)算結(jié)果，可知機(jī)翼阻力下降仍然主要源于激波減弱，但由于激波阻力在總阻力中占比不高，減阻百分比稍低。

同時(shí)，該現(xiàn)象也表明為了更好地優(yōu)化機(jī)翼，應(yīng)該考慮下一步將混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)方法直接應(yīng)用于三維機(jī)翼，在設(shè)定局部流動(dòng)特征約束時(shí)，綜合考慮激波阻力、黏性阻力和誘導(dǎo)阻力，達(dá)到總體最優(yōu)效果。

圖21展示了設(shè)計(jì)機(jī)翼和初始機(jī)翼上翼面壓力分布云圖和上下翼面轉(zhuǎn)捩線。相比于初始機(jī)翼，設(shè)計(jì)機(jī)翼上表面順壓梯度減小，激波強(qiáng)度減弱，但仍然保持了大范圍層流流動(dòng)；設(shè)計(jì)機(jī)翼下表面將翼根部分轉(zhuǎn)捩位置從前緣推遲到40%弦長(zhǎng)之后，層流范圍明顯增大?；鶞?zhǔn)機(jī)翼翼根下表面邊界層由于橫流不穩(wěn)定導(dǎo)致前緣轉(zhuǎn)捩，而優(yōu)化機(jī)翼此處橫流不穩(wěn)定得到抑制，從而推遲了轉(zhuǎn)捩的發(fā)生。

表2 設(shè)計(jì)機(jī)翼與初始機(jī)翼氣動(dòng)性能比較Table 2 Comparison of aerodynamic characteristics ofdesigned and baseline wing

順壓梯度不利于抑制橫流CF不穩(wěn)定波的發(fā)展，而機(jī)翼前緣存在劇烈加速區(qū)，是CF波誘導(dǎo)前緣轉(zhuǎn)捩的“高危區(qū)”。比較圖22翼根下表面壓力分布可以看出，優(yōu)化機(jī)翼下表面前緣加速區(qū)更短，有利于穩(wěn)定橫流CF波擾動(dòng)因子在機(jī)翼前緣的劇烈增長(zhǎng)，從而推遲轉(zhuǎn)捩發(fā)生。圖23給出的是初始機(jī)翼和設(shè)計(jì)機(jī)翼下表面翼根處TS波和CF波擾動(dòng)增長(zhǎng)曲線，可以看出轉(zhuǎn)捩由CF波不穩(wěn)定誘導(dǎo)，設(shè)計(jì)機(jī)翼主要是抑制了CF波不穩(wěn)定性推遲了轉(zhuǎn)捩。為了進(jìn)一步分析初始機(jī)翼與設(shè)計(jì)機(jī)翼下表面翼根前緣區(qū)域橫流強(qiáng)度對(duì)比，圖24給出了8.4%展向位置前緣加速區(qū)4個(gè)流向站位上邊界層內(nèi)無(wú)量綱化的橫流速度型的對(duì)比，圖中縱坐標(biāo)Y為距物面無(wú)量綱高度，橫坐標(biāo)W為無(wú)量綱橫流速度。橫流速度型拐點(diǎn)速度越大，拐點(diǎn)與物面距離越小，則橫流強(qiáng)度越大，對(duì)應(yīng)的CF波不穩(wěn)定擾動(dòng)放大率越大。從圖24中可見(jiàn)，在前緣加速區(qū)，隨流向速度的增大，4個(gè)流向站位處的無(wú)量綱橫流強(qiáng)度W沿流向遞減。同時(shí)，設(shè)計(jì)機(jī)翼每個(gè)流向站位處的橫流強(qiáng)度相對(duì)初始機(jī)翼有所減弱。

圖25給出了一定升力系數(shù)范圍內(nèi)，設(shè)計(jì)機(jī)翼和初始機(jī)翼升阻極曲線、力矩系數(shù)和升阻比隨升力系數(shù)的變化曲線，可見(jiàn)，設(shè)計(jì)機(jī)翼氣動(dòng)性能在一定升力系數(shù)范圍內(nèi)均有顯著提高，驗(yàn)證了優(yōu)化翼型優(yōu)良的超臨界特性和自然層流特性，說(shuō)明了在跨聲速自然層流翼型設(shè)計(jì)中混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)方法的有效性。需要說(shuō)明的是，本文考慮的機(jī)翼設(shè)計(jì)狀態(tài)(CL=0.5)的升阻比不在單獨(dú)機(jī)翼的最大升阻比位置，是因?yàn)楫?dāng)考慮機(jī)身阻力后，翼身組合體的最大升阻比在CL=0.5左右。

圖21 設(shè)計(jì)機(jī)翼和初始機(jī)翼上翼面壓力分布云圖及上下翼面轉(zhuǎn)捩線比較Fig.21 Comparisons of pressure contours of upper surface and transition lines on designed and baseline wings

圖22 初始機(jī)翼和設(shè)計(jì)機(jī)翼翼根站位壓力分布比較Fig.22 Comparison of pressure coefficients at wing root station of baseline wing and designed wing

圖23 初始機(jī)翼與設(shè)計(jì)機(jī)翼在翼根站位下表面TS波和 CF波擾動(dòng)增長(zhǎng)曲線Fig.23 Disturbance growth curves of TS and CF wave on lower surfaces at wing root station of baseline wing and designed wing

圖24 初始機(jī)翼和優(yōu)化機(jī)翼翼根站位下 表面邊界層內(nèi)橫流速度型比較Fig.24 Comparisons of cross flow profiles inside boundary-layers on lower surfaces at wing root station of baseline wing and designed wing

圖25 設(shè)計(jì)機(jī)翼和初始機(jī)翼氣動(dòng)性能曲線Fig.25 Aerodynamic performance curves of designed and baseline wings

5 結(jié)論與展望

1) 在基于代理模型的通用優(yōu)化框架下，實(shí)現(xiàn)了混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)方法。將反設(shè)計(jì)目標(biāo)函數(shù)與直接優(yōu)化目標(biāo)函數(shù)通過(guò)權(quán)重系數(shù)組合，構(gòu)造了多目標(biāo)優(yōu)化形式下的混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)問(wèn)題，以實(shí)現(xiàn)局部流場(chǎng)特性反設(shè)計(jì)及總體性能指標(biāo)優(yōu)化。

2) 將混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)方法成功應(yīng)用于跨聲速自然層流翼型設(shè)計(jì)研究。與基準(zhǔn)翼型相比，優(yōu)化翼型實(shí)現(xiàn)目標(biāo)局部壓力分布，總阻力下降15.5%，激波強(qiáng)度減弱，保持了大范圍層流，具有良好的超臨界和層流特性，其優(yōu)良?xì)鈩?dòng)特性在三維機(jī)翼上也得到驗(yàn)證，說(shuō)明了混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)方法處理跨聲速層流翼型設(shè)計(jì)問(wèn)題的有效性。

3) 用多輪優(yōu)化策略應(yīng)對(duì)優(yōu)化設(shè)計(jì)空間選擇問(wèn)題。第1輪優(yōu)化設(shè)計(jì)中，在給定設(shè)計(jì)空間尋找最優(yōu)解，如果優(yōu)化結(jié)果在當(dāng)前設(shè)計(jì)空間邊界，則以當(dāng)前最優(yōu)解為中心，定義新的設(shè)計(jì)空間，開(kāi)展新一輪優(yōu)化設(shè)計(jì)，以避免初始設(shè)計(jì)空間可能不包含最優(yōu)解的問(wèn)題。

下一步擬將混合反設(shè)計(jì)/優(yōu)化設(shè)計(jì)方法直接應(yīng)用于三維機(jī)翼，在設(shè)定局部流動(dòng)特征約束時(shí)，綜合考慮跨聲速層流機(jī)翼的激波阻力、黏性阻力和誘導(dǎo)阻力，達(dá)到總體最優(yōu)效果。