王一雯，蘭夏毓，史亞云，華俊，白俊強(qiáng)，周鑄

1.西北工業(yè)大學(xué) 無(wú)人系統(tǒng)技術(shù)研究院，西安 710072

2.西北工業(yè)大學(xué) 航空學(xué)院，西安 710072

3. 西安交通大學(xué) 航天航空學(xué)院，西安 710049

4. 中國(guó)航空研究院，北京 100012

5.中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000

為了提高經(jīng)濟(jì)效益、緩解環(huán)境污染，實(shí)現(xiàn)“綠色航空”長(zhǎng)遠(yuǎn)發(fā)展目標(biāo)，航空領(lǐng)域提出了各項(xiàng)技術(shù)以改善飛機(jī)性能[1-7]。研究表明，如果能夠在飛機(jī)機(jī)翼、短艙以及尾翼等區(qū)域保持40%的層流區(qū)域，相比于全湍流流動(dòng)，可減小16%的全機(jī)阻力，直接節(jié)省10%的燃油消耗。因而，美國(guó)國(guó)家航空航天局、德國(guó)宇航中心、中國(guó)商用飛機(jī)有限責(zé)任公司等國(guó)際重要航空研究機(jī)構(gòu)將層流技術(shù)視為未來(lái)最有前景的減阻技術(shù)之一。

層流減阻控制技術(shù)分為被動(dòng)和主動(dòng)控制。自然層流技術(shù)(Natural Laminar Flow，NLF)是一種重要的被動(dòng)層流控制技術(shù)。主動(dòng)層流控制技術(shù)包含了全層流控制(Laminar Flow Control，LFC)與混合層流控制(Hybrid Laminar Flow Control，HLFC)等[8]。相比于LFC，HLFC僅在機(jī)翼前緣采用吸氣控制，在弦向中部通過(guò)形面設(shè)計(jì)維持有利的壓力分布特征。HLFC技術(shù)不會(huì)破壞現(xiàn)有機(jī)翼的翼盒結(jié)構(gòu)，在有效推遲層流-湍流轉(zhuǎn)捩的同時(shí)可極大減小需用吸氣體積流量。因此，HLFC技術(shù)是目前最有希望應(yīng)用于工程中的主動(dòng)層流控制技術(shù)[9]。波音和空客分別基于B757和A320開(kāi)展了大量HLFC飛行試驗(yàn)研究[10-12]。試驗(yàn)結(jié)果表明，HLFC可在大型客機(jī)的機(jī)翼及尾翼表面維持可觀的層流區(qū)。

隨著計(jì)算流體力學(xué)(Computational Fluid Dynamics，CFD)方法和計(jì)算機(jī)的發(fā)展，基于CFD的氣動(dòng)優(yōu)化設(shè)計(jì)在航空領(lǐng)域得到了廣泛應(yīng)用。針對(duì)層流翼的氣動(dòng)優(yōu)化方法主要分為非梯度優(yōu)化設(shè)計(jì)方法和梯度優(yōu)化設(shè)計(jì)方法。非梯度類(lèi)算法[13-15]通常計(jì)算量過(guò)大、成本過(guò)高。雖然借助先進(jìn)的代理模型能夠顯著縮減時(shí)間成本，但是該方法僅能高效處理具有百維以?xún)?nèi)設(shè)計(jì)變量的優(yōu)化問(wèn)題[16]。一旦設(shè)計(jì)變量過(guò)多，基于非梯度算法的優(yōu)化方法將面臨“維度災(zāi)難”。同時(shí)，代理模型普遍存在泛化能力不足的情況。與非梯度優(yōu)化方法相比，梯度優(yōu)化方法具有更高的收斂速度和計(jì)算效率?；陔x散伴隨理論的梯度求解方法的計(jì)算成本與設(shè)計(jì)變量規(guī)模幾乎無(wú)關(guān)。因此，基于離散伴隨的梯度優(yōu)化方法是解決具有大規(guī)模設(shè)計(jì)變量?jī)?yōu)化問(wèn)題最有效的方法之一，已被廣泛應(yīng)用于全湍流工程優(yōu)化[17-18]，并在層流翼優(yōu)化方面也得到了初步發(fā)展。

針對(duì)層流翼優(yōu)化設(shè)計(jì)問(wèn)題，Driver和Zingg[19]將二維雷諾平均Navier-Stokes(Reynolds Average Navier-Stockes,RANS)求解器與MSES(Method of Streamline-based Euler Simulation)結(jié)合，并對(duì)二維翼型的轉(zhuǎn)捩位置進(jìn)行預(yù)測(cè)。Rashad和Zingg[20]在Driver的研究基礎(chǔ)上，基于全湍流流動(dòng)的離散伴隨方程，推導(dǎo)了考慮層流轉(zhuǎn)捩的耦合伴隨方程，并將建立的梯度優(yōu)化應(yīng)用于層流翼型的優(yōu)化設(shè)計(jì)中。史亞云等[21-23]將RANS求解器與簡(jiǎn)化的轉(zhuǎn)捩預(yù)測(cè)模塊相結(jié)合，開(kāi)發(fā)了一種基于離散伴隨的層流翼梯度優(yōu)化框架。Kaya和Tuncer[24]采用基于離散伴隨的自然層流氣動(dòng)優(yōu)化框架，對(duì)NLF(1)-0416翼型進(jìn)行優(yōu)化設(shè)計(jì)。相關(guān)研究證明，基于伴隨優(yōu)化的梯度優(yōu)化設(shè)計(jì)方法能夠高效、可靠地解決層流翼氣動(dòng)外形優(yōu)化問(wèn)題。

近年來(lái)，基于離散伴隨理論的層流翼梯度優(yōu)化方法得到了不斷發(fā)展。但是，相關(guān)研究主要集中在自然層流翼優(yōu)化領(lǐng)域，對(duì)于考慮吸氣控制的混合層流翼的梯度優(yōu)化的研究很少。本文針對(duì)混合層流翼氣動(dòng)優(yōu)化設(shè)計(jì)問(wèn)題，基于離散伴隨理論，構(gòu)建一種可考慮吸氣影響的層流翼梯度優(yōu)化設(shè)計(jì)方法。以二維翼型為研究對(duì)象，開(kāi)展自然層流翼型多點(diǎn)優(yōu)化設(shè)計(jì)研究，并從優(yōu)化結(jié)果角度出發(fā)對(duì)比層流翼梯度優(yōu)化和非梯度優(yōu)化方法。進(jìn)一步，采用構(gòu)建的層流翼梯度優(yōu)化方法，開(kāi)展混合層流翼型多點(diǎn)優(yōu)化設(shè)計(jì)，并與自然層流翼型優(yōu)化設(shè)計(jì)結(jié)果進(jìn)行對(duì)比，驗(yàn)證發(fā)展的可考慮吸氣影響的層流翼梯度優(yōu)化設(shè)計(jì)方法的可靠性，探究混合層流翼型減阻原理。

1 考慮邊界層吸氣控制的轉(zhuǎn)捩預(yù)測(cè)方法

采用基于線性穩(wěn)定性理論的簡(jiǎn)化方法計(jì)算擾動(dòng)放大因子N，并結(jié)合eN方法進(jìn)行轉(zhuǎn)捩預(yù)測(cè)，在此基礎(chǔ)上，將吸氣參數(shù)引入邊界層方程，實(shí)現(xiàn)考慮吸氣控制的轉(zhuǎn)捩預(yù)測(cè)。進(jìn)一步，將轉(zhuǎn)捩模塊與RANS流場(chǎng)求解進(jìn)行耦合，建立基于RANS方程的轉(zhuǎn)捩預(yù)測(cè)方法。

1.1 轉(zhuǎn)捩預(yù)測(cè)方法

構(gòu)建的轉(zhuǎn)捩預(yù)測(cè)方法主要包含：層流邊界層方程、穩(wěn)定性分析和eN方法。采用基于錐形流假設(shè)[25]的準(zhǔn)三維層流邊界層方程，具體方程推導(dǎo)細(xì)節(jié)可參考文獻(xiàn)[26]。將非定常守恒變量表示為時(shí)均量和擾動(dòng)量，即

(1)

(2)

根據(jù)式(2)，流向擾動(dòng)的增長(zhǎng)率為

(3)

對(duì)式(3)取對(duì)數(shù)，并進(jìn)行積分，可得

(4)

式中：A為波幅；A0=e-αi(x0)為初始波幅；N為擾動(dòng)放大因子，由擾動(dòng)放大率αi積分得到；x0為擾動(dòng)開(kāi)始的位置。eN方法通過(guò)對(duì)擾動(dòng)放大率αi沿流向積分，得到不同流向位置的擾動(dòng)放大因子，當(dāng)某一位置的N達(dá)到轉(zhuǎn)捩閾值時(shí)，即認(rèn)為發(fā)生轉(zhuǎn)捩。若轉(zhuǎn)捩發(fā)生時(shí)，擾動(dòng)放大因子N未達(dá)到轉(zhuǎn)捩閾值，則認(rèn)為是由邊界層分離流誘導(dǎo)發(fā)生轉(zhuǎn)捩[22]。

為了提高正向計(jì)算和梯度求解效率，采用基于線性穩(wěn)定理論提出的擾動(dòng)放大因子模型(Amplification Factor Model，AFM)[27-28]，利用邊界層參數(shù)直接計(jì)算T-S波的放大因子。AFM方法最初由Gelyzes等[27]根據(jù)直線斜率近似N因子和動(dòng)量厚度雷諾數(shù)Reδ2的關(guān)系提出，公式為

(5)

式中：H12為形狀因子；Reδ2cr為臨界動(dòng)量厚度雷諾數(shù)。在此基礎(chǔ)上，Drela和Giles[28]提出

2.5tanh(1.5(Hk-3.1)]}2+0.25}0.5

(6)

式中：Hk為可壓縮形狀因子。對(duì)于具有相似解特征的邊界層流動(dòng)，不同流向站位對(duì)應(yīng)的動(dòng)量厚度雷諾數(shù)是相同的。因此放大因子N與流向站位x的關(guān)系可以表示為

(7)

式中：Ue、μe分別為邊界層邊界處的合速度和運(yùn)動(dòng)學(xué)黏性系數(shù)；δ2為邊界層動(dòng)量厚度。

(8)

(9)

將經(jīng)驗(yàn)關(guān)系式(8)、式(9)代入式(7)，可得到式(10)。

(10)

式中：ρ為密度；m(Hk)、l(Hk)均指與Hk相關(guān)的經(jīng)驗(yàn)公式。對(duì)式(10)沿流向積分，即

(11)

式中：x0為不穩(wěn)定性擾動(dòng)開(kāi)始發(fā)展的臨界點(diǎn)，由Reδ2cr決定。根據(jù)式(11)和轉(zhuǎn)捩閾值，可求得轉(zhuǎn)捩位置。

1.2 邊界層吸氣控制

求解邊界層方程可得到邊界層信息，進(jìn)一步可進(jìn)行穩(wěn)定性分析。在層流邊界層方程求解中，壁面邊界條件為

(12)

遠(yuǎn)場(chǎng)邊界條件為

(13)

吸氣控制能夠通過(guò)改變基流條件，進(jìn)而有效抑制邊界層不穩(wěn)定性。當(dāng)引入吸氣控制時(shí)，壁面法向邊界條件為

v=vs

(14)

式中：vs為垂直于物面的連續(xù)吸氣速度，可表示為

(15)

在實(shí)際應(yīng)用中，通常在機(jī)翼表面鋪設(shè)有孔蒙皮進(jìn)行吸氣控制，蒙皮表面的孔徑約為50 μm。由于垂直于機(jī)翼表面的吸力速度的量值很小，屬于微吸氣，因而對(duì)壓力分布的影響可以被忽略。在工程應(yīng)用中，一般采用面吸氣速度替代孔吸氣速度，即考慮整個(gè)吸氣區(qū)域?yàn)榫鶆蛭鼩猓⒃诿總€(gè)網(wǎng)格點(diǎn)處都設(shè)置對(duì)應(yīng)的面吸氣速度，最終保證其積分得到的吸氣質(zhì)量流量與風(fēng)洞試驗(yàn)中的吸氣質(zhì)量流量一致即可[29-31]。

為了便于求解面吸氣速度vs，通常定義其與孔吸氣速度vh之間的關(guān)系為

vs=vhσ

(16)

(17)

式中：σ為孔隙率;dh為吸氣孔直徑;S1為相鄰的4個(gè)孔形成的表面積。通過(guò)式(16)，可將單個(gè)孔吸氣速度vh轉(zhuǎn)換為整個(gè)面積為S1的吸力面的吸氣速度vs。

各吸氣孔對(duì)應(yīng)的孔吸氣速度與有孔蒙皮內(nèi)外的壓差相關(guān)，即

(18)

式中：Pext為有孔蒙皮外部壓強(qiáng)；Pplenum為有孔蒙皮內(nèi)部壓強(qiáng)；ΔP為壓差；δ為蒙皮厚度；A′、B′為實(shí)驗(yàn)標(biāo)定參數(shù)。結(jié)合式(16)～式(18)，即可求得機(jī)翼表面吸氣速度vs。

在此基礎(chǔ)上，定義無(wú)量綱吸氣參數(shù)Cq：

(19)

式中：U∞為自由來(lái)流速度。將式(19)代入式(14)，進(jìn)一步結(jié)合轉(zhuǎn)捩預(yù)測(cè)可獲得考慮吸氣控制的轉(zhuǎn)捩位置。

1.3 耦合轉(zhuǎn)捩預(yù)測(cè)和RANS方程

轉(zhuǎn)捩預(yù)測(cè)模塊與RANS方程的耦合方式如圖1所示。首先，給定固定轉(zhuǎn)捩位置，進(jìn)行RANS計(jì)算，直至流場(chǎng)殘差降至10-8量級(jí)，得到流場(chǎng)解。其次，將計(jì)算得到的流場(chǎng)解(流場(chǎng)壓力分布Cp)與對(duì)應(yīng)的不同機(jī)翼展向站位翼型的幾何信息作為轉(zhuǎn)捩模塊的輸入，通過(guò)引入吸氣控制的層流邊界層方程和AFM方法可以得到新的轉(zhuǎn)捩位置。最后，將新的轉(zhuǎn)捩位置與固定轉(zhuǎn)捩位置進(jìn)行對(duì)比，當(dāng)轉(zhuǎn)捩位置殘差降低到10-6～10-8量級(jí)視為流動(dòng)收斂，否則將求出的轉(zhuǎn)捩位置作為新的輸入，繼續(xù)求解RANS方程。

圖1 轉(zhuǎn)捩預(yù)測(cè)與RANS方程耦合流程圖

1.4 轉(zhuǎn)捩預(yù)測(cè)方法驗(yàn)證

選取基于某型公務(wù)機(jī)的自然/混合層流翼套飛行試驗(yàn)[31]對(duì)構(gòu)建的轉(zhuǎn)捩預(yù)測(cè)方法的預(yù)測(cè)精度進(jìn)行驗(yàn)證,翼套試驗(yàn)段模型如圖2所示。翼套前緣后掠5°，因此主要是流向失穩(wěn)主導(dǎo)轉(zhuǎn)捩。翼套前緣橙色框?yàn)榭筛鼡Q部件。依據(jù)具體試驗(yàn)要求可安裝自然層流和混合層流試驗(yàn)段。對(duì)于混合層流試驗(yàn)，上翼面前緣16%范圍為由4個(gè)獨(dú)立的弦向等弧長(zhǎng)吸氣控制腔體構(gòu)成的吸氣控制區(qū)域。采用文獻(xiàn)[32]中飛行試驗(yàn)數(shù)據(jù)標(biāo)定的T-S波轉(zhuǎn)捩閾值9.0作為本文T-S波失穩(wěn)的轉(zhuǎn)捩判據(jù)進(jìn)行數(shù)值模擬。

圖2 翼套模型

1.4.1 自然層流轉(zhuǎn)捩預(yù)測(cè)方法驗(yàn)證

對(duì)于自然層流飛行試驗(yàn)，選取的典型驗(yàn)證狀態(tài)為：飛行高度H=7 km，飛行馬赫數(shù)Ma=0.458，基于平均氣動(dòng)弦長(zhǎng)的雷諾數(shù)Rec=12.22×106，飛行攻角AoA=2°。圖3為CFD計(jì)算物面網(wǎng)格，其中機(jī)翼展向布置113個(gè)點(diǎn)，弦向布置201個(gè)點(diǎn),并對(duì)翼套進(jìn)行加密,全機(jī)網(wǎng)格量為1 100萬(wàn)。

圖3 CFD計(jì)算網(wǎng)格

圖4對(duì)比了數(shù)值模擬結(jié)果與飛行試驗(yàn)結(jié)果，由圖4(a)可以看出，二者的壓力分布整體吻合較好，僅在頭部附近有所區(qū)別，圖4(b)、圖4(c)分別給出了數(shù)值模擬得到的NTS增長(zhǎng)云圖和對(duì)應(yīng)試驗(yàn)段上表面的紅外線(Infrared Ray,IR)圖像,當(dāng)NTS達(dá)到試驗(yàn)標(biāo)定的轉(zhuǎn)捩閾值時(shí)(即NTS=9)，采用黃色實(shí)線標(biāo)出其數(shù)值模擬轉(zhuǎn)捩線。由圖4(c)可以看出，數(shù)值模擬得到的轉(zhuǎn)捩線與紅外測(cè)量得到的轉(zhuǎn)捩線基本重合。該結(jié)果驗(yàn)證了本文建立的轉(zhuǎn)捩預(yù)測(cè)方法能夠準(zhǔn)確預(yù)測(cè)NTS增長(zhǎng)趨勢(shì)，從而準(zhǔn)確預(yù)測(cè)轉(zhuǎn)捩位置。

圖4 自然層流數(shù)值計(jì)算與飛行試驗(yàn)結(jié)果對(duì)比

1.4.2 混合層流轉(zhuǎn)捩預(yù)測(cè)方法驗(yàn)證

對(duì)于混合層流飛行試驗(yàn)，選取的典型驗(yàn)證狀態(tài)為：飛行高度H=7 km，飛行馬赫數(shù)Ma=0.458，基于平均氣動(dòng)弦長(zhǎng)的雷諾數(shù)Rec=12.22×106，飛行攻角AoA=2.2°。4個(gè)獨(dú)立吸氣控制腔體的吸氣控制體積流量分別為0.017 808、0.003 575、0.003 542、0.005 069 m3/s。此外，用于獲得面吸氣速度vs的式(18)中的參數(shù)A′、B′采用文獻(xiàn)[31]中的標(biāo)定結(jié)果，如表1所示。

表1 不同吸氣腔的標(biāo)定參數(shù)[31]

圖5(a)給出了數(shù)值模擬和試驗(yàn)壓力分布對(duì)比，可以看出二者吻合較好?；谠囼?yàn)參數(shù)，得到如圖5(b)所示的吸氣參數(shù)云圖。圖5(c)、圖5(d)分別給出了數(shù)值模擬得到的NTS增長(zhǎng)云圖及飛行試驗(yàn)得到的試驗(yàn)段IR圖像。由圖5(d)可以看出，數(shù)值模擬的轉(zhuǎn)捩線與IR圖像的轉(zhuǎn)捩線基本重合。因此采用的轉(zhuǎn)捩預(yù)測(cè)方法能夠準(zhǔn)確預(yù)測(cè)考慮吸氣影響的翼型轉(zhuǎn)捩位置。

圖5 混合層流數(shù)值計(jì)算與飛行試驗(yàn)結(jié)果對(duì)比

2 基于離散伴隨的層流翼梯度優(yōu)化方法

2.1 考慮轉(zhuǎn)捩的耦合伴隨方程推導(dǎo)與求解

2.1.1 耦合伴隨方法推導(dǎo)

在全湍流伴隨方程的基礎(chǔ)上，開(kāi)展考慮轉(zhuǎn)捩的耦合伴隨方程研究。對(duì)于層流-湍流轉(zhuǎn)捩預(yù)測(cè)，相比全湍流(狀態(tài)變量為Q)，會(huì)新增狀態(tài)變量，即轉(zhuǎn)捩位置/轉(zhuǎn)捩長(zhǎng)度Tr。進(jìn)一步，將增廣的殘差變量、狀態(tài)變量Y及對(duì)應(yīng)的伴隨向量Ψ表示為

(20)

(21)

式(21)可改寫(xiě)為

(22)

目標(biāo)函數(shù)I關(guān)于設(shè)計(jì)變量的全導(dǎo)數(shù)可表示為

(23)

將式(22)代入式(23)，可得到

(24)

可將考慮轉(zhuǎn)捩的的耦合伴隨方程寫(xiě)為

(25)

式中：?I/?Q、?I/?Tr分別為目標(biāo)函數(shù)對(duì)RANS方程狀態(tài)變量和轉(zhuǎn)捩狀態(tài)變量的偏導(dǎo)數(shù)。進(jìn)一步，可將目標(biāo)函數(shù)的全導(dǎo)數(shù)表示為

(26)

2.1.2 耦合伴隨方程求解

采用Copuled Krylov方法(CK)方法求解耦合伴隨方程式(25)。其中，考慮轉(zhuǎn)捩的耦合伴隨方程中涉及的偏導(dǎo)數(shù)雅可比矩陣采用無(wú)矩陣形式進(jìn)行存儲(chǔ)，結(jié)合鏈?zhǔn)角髮?dǎo)法則與反向自動(dòng)微分進(jìn)行計(jì)算?？紤]邊界層吸氣控制效應(yīng)的邊界條件的相關(guān)偏導(dǎo)數(shù)通過(guò)反向自動(dòng)微分獲得。

采用CK方法求解式(25)的主要步驟如下：

1) 計(jì)算相關(guān)雅可比矩陣與向量點(diǎn)積后的向量Z。

2.2 基于耦合伴隨方程的梯度驗(yàn)證

采用如圖6所示的翼型進(jìn)行梯度驗(yàn)證，翼型上表面吸氣位置設(shè)置在20%c處，吸氣強(qiáng)度為0.02%，上下翼面各12個(gè)FFD控制點(diǎn)，其中翼型前緣上下控制點(diǎn)實(shí)際為1個(gè)設(shè)計(jì)變量，后緣上下控制點(diǎn)保持不變，因此實(shí)際控制變量為21個(gè)。

圖6 翼型幾何、FFD控制點(diǎn)及轉(zhuǎn)捩位置示意圖

由于采用的翼型優(yōu)化設(shè)計(jì)均為定升力優(yōu)化，因此僅對(duì)比伴隨算法和有限差分算法求得的Cd關(guān)于設(shè)計(jì)變量的梯度對(duì)比，如表2所示。由表2可知，基于耦合伴隨方程與基于有限差分法求得的梯度符合較好，其相對(duì)誤差范圍為10-3～10-7，符合有限差分方法的截?cái)嗾`差精度。該梯度驗(yàn)證結(jié)果證明在考慮邊界層吸氣控制的情況下，本文采用的耦合伴隨方法可精確計(jì)算梯度信息。

表2 伴隨算法和有限差分算法求得目標(biāo)函數(shù)CD關(guān)于設(shè)計(jì)變量的梯度對(duì)比

2.3 梯度優(yōu)化設(shè)計(jì)框架

SNOPT(Sparse Nonlinear Optimizer)梯度優(yōu)化工具采用序列二次規(guī)劃(Sequential Quadratic Programming，SQP)算法，可以處理設(shè)計(jì)變量和約束數(shù)量上千的大規(guī)模非線性?xún)?yōu)化問(wèn)題。在優(yōu)化過(guò)程中，需要將目標(biāo)函數(shù)及梯度作為SNOPT的輸入，SNOPT最終返回設(shè)計(jì)變量。結(jié)合自由變形(Free-Form Deformation，F(xiàn)FD)參數(shù)化方法、逆距離權(quán)重(Inverse Distance Weighted，IDW)網(wǎng)格變形技術(shù)、CFD流場(chǎng)及伴隨方程求解方法、梯度優(yōu)化算法，建立如圖7所示的氣動(dòng)外形優(yōu)化設(shè)計(jì)框架。

圖7 梯度優(yōu)化設(shè)計(jì)框架

3 層流翼型優(yōu)化設(shè)計(jì)與分析

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

參考Honda Jet的設(shè)計(jì)狀態(tài)[33]，進(jìn)行多點(diǎn)自然層流翼型梯度優(yōu)化，在此基礎(chǔ)上，將梯度優(yōu)化結(jié)果與非梯度優(yōu)化(采用NSGAII(Non-dominated Sorting Genetic Algorithms II)算法)結(jié)果進(jìn)行對(duì)比。

初始翼型為一個(gè)典型層流翼型，如圖8所示，上、下翼面各12個(gè)FFD控制點(diǎn)，翼型后緣最后一排FFD控制點(diǎn)(如圖8后緣圓內(nèi)設(shè)計(jì)點(diǎn))不作為設(shè)計(jì)變量；翼型前緣第1排FFD控制點(diǎn)(如圖8前緣圓內(nèi)設(shè)計(jì)點(diǎn))采用約束移動(dòng)，即上、下控制點(diǎn)在y向以相同移動(dòng)量反向移動(dòng)，用于控制翼型前緣半徑；其他FFD控制點(diǎn)均選為設(shè)計(jì)變量。最終整個(gè)優(yōu)化問(wèn)題具有21個(gè)設(shè)計(jì)變量。計(jì)算網(wǎng)格見(jiàn)圖9，轉(zhuǎn)捩閾值選為9.0。

圖8 初始翼型及FFD設(shè)計(jì)變量

圖9 初始翼型空間網(wǎng)格

優(yōu)化目標(biāo)及約束的數(shù)學(xué)描述為

(27)

式中：CDi為每個(gè)設(shè)計(jì)狀態(tài)的阻力系數(shù)；wi為對(duì)應(yīng)的權(quán)重系數(shù)。幾何優(yōu)化約束為翼型面積Sf不小于初始翼型面積Sinit，翼型厚度ty不小于初始翼型厚度(tyinit)的95%，氣動(dòng)優(yōu)化約束為第4個(gè)設(shè)計(jì)工況下的力矩系數(shù)Cm4不小于-0.04。

如表3所示，設(shè)計(jì)狀態(tài)1為巡航狀態(tài)，設(shè)計(jì)狀態(tài)3為爬升狀態(tài)，此外，為了降低優(yōu)化結(jié)果對(duì)飛行條件變化的敏感性，增加針對(duì)設(shè)計(jì)狀態(tài)2、設(shè)計(jì)狀態(tài)4的優(yōu)化設(shè)計(jì)。其中，考慮設(shè)計(jì)狀態(tài)2(巡航)以降低翼型型阻，并且為了在低雷諾數(shù)、高馬赫數(shù)狀態(tài)下最小化配平阻力，增加設(shè)計(jì)狀態(tài)4(巡航)。

表3 Honda Jet多點(diǎn)優(yōu)化設(shè)計(jì)狀態(tài)

圖10給出的初始翼型和優(yōu)化翼型對(duì)比顯示，不同優(yōu)化算法得到的設(shè)計(jì)結(jié)果呈現(xiàn)很強(qiáng)的一致性。相比于初始翼型，優(yōu)化翼型40%c～60%c區(qū)域上翼面彎度增大，最大厚度點(diǎn)位置適當(dāng)后移，上翼面60%c～90%c區(qū)域彎度減小。下翼面最大厚度點(diǎn)位置顯著后移。在翼型前緣附近，優(yōu)化翼型均呈現(xiàn)出上翼前緣附近彎度適當(dāng)增大，下翼面前緣附近彎度適當(dāng)縮小的變化趨勢(shì)。

圖10 初始翼型與優(yōu)化翼型幾何對(duì)比

表4分別給出4種設(shè)計(jì)狀態(tài)下，初始翼型、非梯度優(yōu)化翼型(NSGAII)、梯度優(yōu)化翼型(SNOPT)的阻力系數(shù)(CD)、上下表面轉(zhuǎn)捩位置(xU,tr/c、xL,tr/c)。初始翼型和2類(lèi)優(yōu)化翼型的壓力分布如圖11所示，其中，圖11(a)～圖11(d)分別對(duì)應(yīng)設(shè)計(jì)狀態(tài)1～4的壓力分布對(duì)比結(jié)果。圖12給出了不同設(shè)計(jì)狀態(tài)下的T-S擾動(dòng)波對(duì)比，其中設(shè)計(jì)狀態(tài)1～ 4分別對(duì)應(yīng)Point 1～Point 4。

表4顯示，在設(shè)計(jì)狀態(tài)1，梯度優(yōu)化翼型和非梯度優(yōu)化翼型分別減阻15.9 counts、11.6 counts；在設(shè)計(jì)狀態(tài)2，SNOPT、NSGAII優(yōu)化翼型分別減阻34.1、29.6 counts(1 count=1×10-4)。2類(lèi)優(yōu)化翼型上、下翼面轉(zhuǎn)捩位置相近，且比初始翼型均有明顯推遲。設(shè)計(jì)狀態(tài)1、設(shè)計(jì)狀態(tài)2上翼面轉(zhuǎn)捩點(diǎn)推遲到55%c以后，下翼面轉(zhuǎn)捩點(diǎn)推遲到68%c以后。圖11(a)、圖11(b)的壓力分布對(duì)比顯示，相比于初始翼型，2類(lèi)優(yōu)化翼型均通過(guò)適當(dāng)減小上翼面順壓力梯度大小的方式，增大了順壓梯度范圍，實(shí)現(xiàn)了轉(zhuǎn)捩位置的有效推遲。在下翼面，優(yōu)化翼型維持了更長(zhǎng)的順壓梯度區(qū)。

對(duì)于設(shè)計(jì)狀態(tài)3，梯度優(yōu)化翼型和非梯度優(yōu)化翼型分別減阻2.9、9.7 counts，如表4所示。相比于梯度優(yōu)化翼型，非梯度優(yōu)化翼型具有更大的減阻收益。圖11(c)顯示，梯度優(yōu)化翼型上翼面負(fù)壓峰值較高，峰值之后維持了一定弱逆壓力梯度，導(dǎo)致優(yōu)化翼型上表面NTS發(fā)展迅速，在37.5%c處達(dá)到轉(zhuǎn)捩閾值(如圖12的Point 3所示)，從而使翼型上表面轉(zhuǎn)捩提前。而非梯度優(yōu)化翼型上翼面負(fù)壓峰值較低，峰值之后維持較長(zhǎng)的壓力平臺(tái)區(qū)，使得上翼面轉(zhuǎn)捩位置比初始翼型推遲了14.2%，相比于梯度優(yōu)化翼型，上表面層流區(qū)延長(zhǎng)近20%)。此外，2類(lèi)優(yōu)化翼型的下翼面均保持較長(zhǎng)的順壓梯度，從而顯著抑制了T-S波的發(fā)展。

表4 初始翼型與優(yōu)化翼型氣動(dòng)數(shù)據(jù)對(duì)比

圖11(d)顯示，初始翼型在設(shè)計(jì)狀態(tài)4具有較強(qiáng)的激波，梯度優(yōu)化翼型通過(guò)提高頭部負(fù)壓峰值，消除初始翼型的大順壓梯度，顯著削弱激波強(qiáng)度。而非梯度優(yōu)化翼型通過(guò)減小40%c之后區(qū)域的順壓梯度，適當(dāng)削弱了激波強(qiáng)度。在層流保持方面，壓力分布(圖11(d))和圖12給出的NTS增長(zhǎng)曲線顯示，初始翼型上翼面為激波引起的轉(zhuǎn)捩。優(yōu)化翼型通過(guò)將激波位置后移，推遲了轉(zhuǎn)捩的發(fā)生。由于梯度優(yōu)化翼型維持一定弱逆壓力梯度來(lái)顯著減小激波阻力，因此其上翼面層流區(qū)范圍略小于非梯度優(yōu)化翼型。但是在下翼面，梯度優(yōu)化翼型維持更長(zhǎng)的層流區(qū)范圍。最終，相比于初始翼型，梯度優(yōu)化翼型、非梯度優(yōu)化翼型分別減阻21.7、19.1 counts。

圖11 初始翼型與優(yōu)化翼型壓力分布對(duì)比

圖12 多點(diǎn)優(yōu)化不同設(shè)計(jì)狀態(tài)放大因子曲線對(duì)比

自然層流翼型多點(diǎn)優(yōu)化設(shè)計(jì)結(jié)果表明，梯度優(yōu)化翼型與非梯度優(yōu)化翼型外形及壓力分布形態(tài)特征相似，但仍存在一定差異。梯度優(yōu)化結(jié)果更加追求巡航狀態(tài)(設(shè)計(jì)狀態(tài)1、2、4)的氣動(dòng)性能，而非梯度優(yōu)化結(jié)果更加追求爬升狀態(tài)(設(shè)計(jì)狀態(tài)3)的氣動(dòng)性能。導(dǎo)致這一現(xiàn)象的原因可能與層流翼多點(diǎn)優(yōu)化問(wèn)題存在多極值特性相關(guān)。總之對(duì)比結(jié)果表明，構(gòu)建的基于離散伴隨的層流翼梯度優(yōu)化方法可有效實(shí)現(xiàn)層流翼氣動(dòng)性能的大幅提升。

3.2 考慮吸氣控制的混合層流翼型優(yōu)化設(shè)計(jì)

參考Honda Jet的設(shè)計(jì)狀態(tài)[33]，在自然層流翼型優(yōu)化的基礎(chǔ)上，引入邊界層流吸氣控制，進(jìn)行混合層流翼型的多點(diǎn)梯度優(yōu)化設(shè)計(jì)。其中，翼型上表面前20%c區(qū)域?yàn)槲鼩饪刂茀^(qū)域，如圖13所示，整個(gè)吸氣控制區(qū)域采用均勻吸氣，其吸氣強(qiáng)度Cq=-0.000 2。初始翼型及FFD設(shè)計(jì)變量、空間網(wǎng)格、設(shè)計(jì)狀態(tài)、優(yōu)化目標(biāo)/約束及不同目標(biāo)目標(biāo)間的權(quán)重系數(shù)的選取均與3.1節(jié)相同。

圖13 初始翼型吸氣控制區(qū)域示意圖

表5給出了初始翼型和自然層流、混合層流優(yōu)化翼型的氣動(dòng)數(shù)據(jù)(CD)及上下表面轉(zhuǎn)捩位置(xU,tr/c、xL,tr/c)對(duì)比。對(duì)于初始翼型，吸氣控制的引入僅對(duì)設(shè)計(jì)狀態(tài)3的氣動(dòng)性能有較為明顯的提升，其上翼面轉(zhuǎn)捩位置推遲了5.2%c，氣動(dòng)減阻2.9 counts。但是，其余3個(gè)設(shè)計(jì)狀態(tài)的氣動(dòng)性能無(wú)明顯改善。這是由于吸氣控制難以推遲強(qiáng)逆壓力梯度引起的轉(zhuǎn)捩，而初始構(gòu)型的轉(zhuǎn)捩主要由強(qiáng)逆壓力梯度甚至激波導(dǎo)致。

表5 初始翼型與NLF、HLFC優(yōu)化翼型氣動(dòng)數(shù)據(jù)對(duì)比

相比于自然層流優(yōu)化翼型，混合層流優(yōu)化翼型上表面轉(zhuǎn)捩位置分別推遲6.1%c、7%c、26.4%c、8.6%c，其氣動(dòng)阻力分別減小2.0、1.9、14.3、3.8 counts,其相對(duì)減阻收益分別為6.1%、5.9%、33.3%、9.5%。在4種設(shè)計(jì)狀態(tài)下，引入吸氣控制均能顯著推遲優(yōu)化翼型上表面轉(zhuǎn)捩位置，從而實(shí)現(xiàn)其氣動(dòng)阻力的減小。

圖14、圖15分別給出了初始翼型與自然層流、混合層流優(yōu)化翼型及剖面壓力分布對(duì)比，圖16 給出了不同設(shè)計(jì)狀態(tài)下的T-S擾動(dòng)波對(duì)比。

優(yōu)化前后翼型對(duì)比(圖14)顯示，相比自然層流優(yōu)化翼型，混合層流優(yōu)化翼型40%c～70%c區(qū)域上翼面彎度進(jìn)一步增大，最大厚度位置后移明顯。上翼面前緣15%c區(qū)域彎度也有一定程度增大。此外，翼型后緣區(qū)域整體彎度減小。

圖14 初始翼型與優(yōu)化翼型對(duì)比

圖15給出的剖面壓力分布對(duì)比圖顯示，相比于自然層流翼型，各個(gè)設(shè)計(jì)狀態(tài)下，混合層流優(yōu)化翼型上翼面強(qiáng)逆壓力梯度起始位置后移8%c～10%c。導(dǎo)致這一現(xiàn)象的主要原因在于混合層流優(yōu)化翼型最大厚度位置后移，上翼面局部彎度增大。

對(duì)于設(shè)計(jì)狀態(tài)1、2，在上翼面前20%c區(qū)域(吸氣控制區(qū)域)，由于吸氣控制完全抑制了該區(qū)域T-S擾動(dòng)波的發(fā)展(見(jiàn)圖16(a)、圖16(b))，混合層流優(yōu)化翼型維持了比自然層流優(yōu)化翼型更小的順壓梯度。在20%c～40%c區(qū)域內(nèi)，混合層流優(yōu)化翼型順壓梯度迅速增大，以抑制T-S波的快速發(fā)展。之后，40%c～60%c區(qū)域內(nèi)維持弱逆壓力梯度，T-S波迅速發(fā)展，在接近基于強(qiáng)逆壓力梯度的壓力恢復(fù)區(qū)時(shí)，NTS值快速接近轉(zhuǎn)捩閾值，如圖16(a)、圖16(b)所示。

對(duì)于設(shè)計(jì)狀態(tài)3，除吸氣控制區(qū)域維持弱逆壓梯度外，上翼面前60%c區(qū)域維持壓力平頂，T-S 擾動(dòng)波得到完全抑制，如圖15(c)、圖16(c)所示。對(duì)于設(shè)計(jì)狀態(tài)4，相比自然層流優(yōu)化翼型，混合層流優(yōu)化翼型在40%c以后維持了更長(zhǎng)的逆壓力梯度區(qū)，以更弱的激波進(jìn)行壓力恢復(fù)。過(guò)大的逆壓力梯度范圍會(huì)加速T-S擾動(dòng)波的發(fā)展。但是吸氣控制的引入以及在20%c～40%c范圍內(nèi)維持更大的順壓力梯度，使得前40%c區(qū)域內(nèi)NTS值的增長(zhǎng)非常緩慢，如圖15(d)、圖16(d)所示。最終，相比于自然層流優(yōu)化翼型，混合層流優(yōu)化翼型在上翼面維持了更長(zhǎng)的層流區(qū)范圍。

圖15 初始翼型與HLFC優(yōu)化翼型壓力分布對(duì)比

圖16 多點(diǎn)優(yōu)化不同設(shè)計(jì)狀態(tài)放大因子曲線對(duì)比

為了進(jìn)一步說(shuō)明吸氣控制對(duì)翼型邊界層的影響，對(duì)比了初始翼型在吸氣控制前后的近壁區(qū)速度型，如圖17所示。其中，橫坐標(biāo)u為流向速度，d為距壁面的垂直距離。相比于無(wú)吸氣控制的初始翼型上翼面近壁區(qū)速度型，吸氣控制使得上翼面的速度型更加“飽滿(mǎn)”，進(jìn)而減小其邊界層厚度，以達(dá)到抑制TS波發(fā)展的目的。

圖17 考慮吸氣控制影響的初始翼型上翼面沿流向不同站位處的速度型對(duì)比

4 結(jié) 論

基于離散伴隨理論，發(fā)展了一種可考慮吸氣影響的層流翼型梯度優(yōu)化設(shè)計(jì)方法。在此基礎(chǔ)上，分別進(jìn)行了自然層流和混合層流翼型優(yōu)化設(shè)計(jì)研究，得到如下主要結(jié)論：

1) 利用自然/混合層流翼套飛行試驗(yàn)對(duì)轉(zhuǎn)捩預(yù)測(cè)方法的可靠性進(jìn)行驗(yàn)證。數(shù)值模擬和試驗(yàn)對(duì)比結(jié)果顯示，基于AFM方法的轉(zhuǎn)捩預(yù)測(cè)方法能夠有效捕捉T-S擾動(dòng)波的發(fā)展，可滿(mǎn)足自然/混合層流翼優(yōu)化設(shè)計(jì)研究對(duì)轉(zhuǎn)捩預(yù)測(cè)精度的需求。

2) 基于離散伴隨理論，推導(dǎo)并求解了考慮邊界層吸氣控制效應(yīng)的轉(zhuǎn)捩耦合伴隨方程?？紤]吸氣控制的影響下，構(gòu)建的基于耦合伴隨方程的梯度求解與基于有限差分的梯度求解的絕對(duì)誤差在10-9～10-5量級(jí)范圍。

3) 利用建立的梯度優(yōu)化方法及基于NSGAII的非梯度優(yōu)化方法，開(kāi)展了自然層流翼型多點(diǎn)優(yōu)化設(shè)計(jì)研究。優(yōu)化結(jié)果顯示，2種不同優(yōu)化方法得到的設(shè)計(jì)結(jié)果呈現(xiàn)出較強(qiáng)的一致性。對(duì)比結(jié)果表明構(gòu)建的基于離散伴隨的層流翼梯度優(yōu)化方法可有效實(shí)現(xiàn)層流翼氣動(dòng)性能的大幅提升。

4) 利用建立的梯度優(yōu)化方法，開(kāi)展了混合層流翼型多點(diǎn)優(yōu)化設(shè)計(jì)研究。相比于自然層流優(yōu)化翼型，混合層流優(yōu)化翼型氣動(dòng)阻力分別減小了2.0、1.9、14.3、3.8 counts,其相對(duì)減阻收益為6.1%、5.9%、33.3%、9.5%。優(yōu)化結(jié)果表明，混合層流翼型具有更強(qiáng)的層流保持能力，可更有效地實(shí)現(xiàn)激波阻力與摩擦阻力間的權(quán)衡。