基于離散伴隨的層流翼優(yōu)化設(shè)計(jì)方法

楊體浩，王一雯，王雨桐，史亞云，周鑄

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

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

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

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

層流技術(shù)是實(shí)現(xiàn)綠色航空發(fā)展目標(biāo)的核心技術(shù)之一，對(duì)于運(yùn)輸類飛機(jī)，在機(jī)翼、平垂尾、短艙表面維系50%的層流區(qū)，可帶來(lái)全機(jī)總阻力15%的降低[1]。目前層流技術(shù)僅在小型公務(wù)機(jī)、無(wú)人機(jī)以及大型客機(jī)的短艙、翼梢小翼上得到了初步應(yīng)用[2-5]，還未在大型客機(jī)諸如機(jī)翼、尾翼主要升力面部件上得到應(yīng)用。大型客機(jī)的升力面部件后掠角較大、流動(dòng)具有跨聲速、高雷諾數(shù)的特點(diǎn)，導(dǎo)致了具有現(xiàn)象復(fù)雜、多種轉(zhuǎn)捩機(jī)制并存的流動(dòng)物理特征[6]，為層流技術(shù)在大型客機(jī)升力面部件上的應(yīng)用提出了巨大挑戰(zhàn)。構(gòu)建精確、高效、可捕捉多種轉(zhuǎn)捩機(jī)制的轉(zhuǎn)捩預(yù)測(cè)方法，并發(fā)展高效、魯棒的優(yōu)化設(shè)計(jì)方法是層流技術(shù)工程應(yīng)用轉(zhuǎn)換進(jìn)程中關(guān)鍵的一環(huán)。

目前，基于線性穩(wěn)定理論(Linear Stability Theory, LST)的轉(zhuǎn)捩預(yù)測(cè)方法在工程中得到了廣泛應(yīng)用，該方法可以有效預(yù)測(cè)Tollmien-Schlichting(TS)和CrossFlow(CF)不穩(wěn)定性誘導(dǎo)的轉(zhuǎn)捩現(xiàn)象。Krumbein等[7-8]發(fā)展了耦合了基于eN方法及非結(jié)構(gòu)RANS方程求解器的轉(zhuǎn)捩預(yù)測(cè)方法，并對(duì)二維翼型、機(jī)翼以及帶有增升的復(fù)雜三維構(gòu)型進(jìn)行模擬、驗(yàn)證與分析。Liao等[9]基于Fun3d非結(jié)構(gòu)求解器直接提取RANS模擬的邊界層信息，并耦合LST進(jìn)行轉(zhuǎn)捩預(yù)測(cè)研究，并從轉(zhuǎn)捩預(yù)測(cè)精度角度指出對(duì)于小展弦比機(jī)翼考慮邊界層內(nèi)三維效應(yīng)的重要性。Shi等[10]基于結(jié)構(gòu)求解器ADflow耦合基于LST的eN方法，對(duì)低速、跨聲速等構(gòu)型開(kāi)展了大量數(shù)值模擬并與風(fēng)洞試驗(yàn)進(jìn)行對(duì)比，驗(yàn)證了該轉(zhuǎn)捩預(yù)測(cè)模型可滿足工程應(yīng)用研究對(duì)轉(zhuǎn)捩預(yù)測(cè)精度及可靠性的需求?；贚ST的轉(zhuǎn)捩預(yù)測(cè)方法與不同RANS求解器耦合的廣泛應(yīng)用，促進(jìn)了面向工程應(yīng)用的高效、可靠的轉(zhuǎn)捩預(yù)測(cè)數(shù)值方法的發(fā)展。

為了進(jìn)一步提高計(jì)算效率，Arnal[11]基于線性穩(wěn)定理論提出簡(jiǎn)化的eN轉(zhuǎn)捩預(yù)測(cè)方法，基本思想是建立擾動(dòng)放大因子與邊界層信息的關(guān)系，即數(shù)據(jù)庫(kù)方法。Perraud等[12]對(duì)該簡(jiǎn)化的eN方法進(jìn)行了大量數(shù)值驗(yàn)證，結(jié)果表明該方法相比LST可以快速高效的計(jì)算擾動(dòng)放大因子，同時(shí)與LST計(jì)算精度相當(dāng)。Yang等[13]針對(duì)商用客機(jī)真實(shí)運(yùn)營(yíng)環(huán)境中的雷諾數(shù)狀態(tài)，進(jìn)行了飛行試驗(yàn)研究，并對(duì)比驗(yàn)證了基于LST和數(shù)據(jù)庫(kù)的eN方法，結(jié)果表明基于數(shù)據(jù)庫(kù)的eN方法可準(zhǔn)確的預(yù)測(cè)轉(zhuǎn)捩位置。簡(jiǎn)化的eN轉(zhuǎn)捩預(yù)測(cè)方法進(jìn)一步為高效層流翼優(yōu)化方法的建立提供了基礎(chǔ)。

傳統(tǒng)的智能算法在層流翼優(yōu)化中得到了一定應(yīng)用[14-16]，但對(duì)高維設(shè)計(jì)變量問(wèn)題，存在“維度災(zāi)難”問(wèn)題，目前難以滿足具有大規(guī)模設(shè)計(jì)變量的全機(jī)一體化設(shè)計(jì)的發(fā)展趨勢(shì)。近些年，基于伴隨理論的梯度優(yōu)化方法因其高計(jì)算效率，使得該方法在針對(duì)大規(guī)模設(shè)計(jì)變量的數(shù)值優(yōu)化問(wèn)題中得到了快速發(fā)展[17-21]。Kenway等[17]指出，基于離散伴隨的梯度優(yōu)化設(shè)計(jì)是目前針對(duì)大規(guī)模設(shè)計(jì)變量?jī)?yōu)化問(wèn)題最有效的解決途徑?；陔x散伴隨的梯度優(yōu)化在不同知名求解器中得到了應(yīng)用，涌現(xiàn)出大量針對(duì)常規(guī)布局的飛行器，以及包括飛翼、支撐翼等在內(nèi)的復(fù)雜非常規(guī)布局飛行器的成功優(yōu)化實(shí)例[22-24]，極大地促進(jìn)了梯度優(yōu)化設(shè)計(jì)方法在飛行器設(shè)計(jì)中的應(yīng)用。

近年來(lái)，基于伴隨理論的梯度優(yōu)化方法在層流翼優(yōu)化設(shè)計(jì)領(lǐng)域內(nèi)也得到了初步發(fā)展。按照采用的轉(zhuǎn)捩預(yù)測(cè)方法不同，現(xiàn)有基于伴隨理論的層流翼梯度優(yōu)化方法的研究主要分為基于輸運(yùn)模式的梯度優(yōu)化方法和基于eN的梯度優(yōu)化方法。Zhang[25]、Khayatzadeh[26]、Halila[27]等基于Langtry和Menter提出的輸運(yùn)方程[28]構(gòu)建了面向?qū)恿饕硇偷奶荻葍?yōu)化方法，促進(jìn)了層流翼梯度優(yōu)化方法的發(fā)展。但是，相關(guān)研究存在梯度信息求解不精確的不足，如采取凍結(jié)轉(zhuǎn)捩模型的近似處理方法，缺乏伴隨方程精確推導(dǎo)?；跀?shù)據(jù)庫(kù)的eN方法, Lee等[29]發(fā)展了基于離散伴隨的層流翼梯度優(yōu)化方法。但是該方法沒(méi)有精確考慮轉(zhuǎn)捩部分對(duì)梯度信息的影響，導(dǎo)致層流翼優(yōu)化結(jié)果出現(xiàn)轉(zhuǎn)捩位置提前的現(xiàn)象。

目前，可精確考慮轉(zhuǎn)捩部分對(duì)梯度信息影響的層流翼梯度優(yōu)化[30-32]的工作也得到了發(fā)展。Rashad等[30]基于RANS方程和MSES程序的轉(zhuǎn)捩預(yù)測(cè)方法，推導(dǎo)了考慮轉(zhuǎn)捩的離散伴隨方程，借助復(fù)變量微分和有限差分組裝相關(guān)雅可比矩陣，進(jìn)一步求解耦合伴隨方程，獲得精確的梯度信息。利用構(gòu)建的梯度優(yōu)化方法，Rashad等[30]成功進(jìn)行了層流翼型的單點(diǎn)和多點(diǎn)優(yōu)化設(shè)計(jì)研究，顯著降低了翼型摩擦阻力和壓差阻力。Shi等[31]基于簡(jiǎn)化的eN方法和RANS方程，借助鏈?zhǔn)角髮?dǎo)法則及混合自動(dòng)微分，建立了解析形式的考慮轉(zhuǎn)捩的耦合伴隨方程及雅可比矩陣，避免了低效的復(fù)變量微分及精度、效率較低的有限差分的應(yīng)用。同時(shí)，采用了無(wú)矩陣存儲(chǔ)方法，降低了計(jì)算內(nèi)存需求提高了計(jì)算效率，并發(fā)展了高效的CK(Coupled Krylov)方法求解耦合伴隨方程，在層流翼型的氣動(dòng)優(yōu)化上得到了成功應(yīng)用。進(jìn)一步，Shi等[32]將該梯度優(yōu)化設(shè)計(jì)推廣到考慮橫流的準(zhǔn)三維優(yōu)化設(shè)計(jì)中，驗(yàn)證了該梯度優(yōu)化框架具有處理多種轉(zhuǎn)捩機(jī)制并存的復(fù)雜層流翼優(yōu)化問(wèn)題的解決能力。

但是，目前圍繞層流翼梯度優(yōu)化的研究?jī)H僅以簡(jiǎn)單的二維翼型或者無(wú)限展長(zhǎng)后掠翼為研究對(duì)象，尚缺少針對(duì)典型客機(jī)翼身組合體這種復(fù)雜三維構(gòu)型的層流翼梯度優(yōu)化研究。本文基于已有研究基礎(chǔ)，結(jié)合混合自動(dòng)微分、無(wú)矩陣存儲(chǔ)技術(shù)、CK求解算法及伴隨方程理論進(jìn)一步發(fā)展面向三維機(jī)翼的基于離散伴隨的梯度優(yōu)化設(shè)計(jì)方法。以具有20°機(jī)翼前緣后掠角的典型客機(jī)翼身組合體構(gòu)型為研究對(duì)象，開(kāi)展三維層流翼梯度優(yōu)化設(shè)計(jì)研究；驗(yàn)證針對(duì)三維復(fù)雜構(gòu)型發(fā)展的層流翼梯度優(yōu)化設(shè)計(jì)方法的可靠性；探究適用于典型客機(jī)的自然層流超臨界機(jī)翼氣動(dòng)設(shè)計(jì)原理。

1 轉(zhuǎn)捩數(shù)值模擬方法

通過(guò)間歇因子耦合RANS方程和轉(zhuǎn)捩模塊，迭代計(jì)算，直至轉(zhuǎn)捩位置和流場(chǎng)解收斂。轉(zhuǎn)捩模塊主要包含層流邊界層方程和轉(zhuǎn)捩位置計(jì)算準(zhǔn)則兩部分，采用Drela-Giles算法和C1準(zhǔn)則進(jìn)行轉(zhuǎn)捩位置計(jì)算。

1.1 RANS求解器

流動(dòng)控制方程采用結(jié)構(gòu)網(wǎng)格的RANS求解器，該求解器基于MPI并行計(jì)算，采用二階有限體積法進(jìn)行空間離散，支持多種偽時(shí)間推進(jìn)。湍流模型包含Spalart-Allmaras(SA)一方程和k-wSST兩方程模型等，本文使用SA湍流模型。

RANS求解器的關(guān)鍵輸入是固定轉(zhuǎn)捩位置信息。迭代計(jì)算中的轉(zhuǎn)捩位置信息為轉(zhuǎn)捩模塊計(jì)算得到并進(jìn)行松弛處理的轉(zhuǎn)捩位置。RANS求解器流場(chǎng)收斂后(流場(chǎng)殘差小于10-8)，提取機(jī)翼剖面壓力分布作為轉(zhuǎn)捩模塊的輸入。

1.2 層流邊界層方程求解

進(jìn)行穩(wěn)定性分析或者通過(guò)轉(zhuǎn)捩準(zhǔn)則計(jì)算轉(zhuǎn)捩位置都需要獲得邊界層信息，主要包含邊界層速度型、溫度及它們的一階和二階導(dǎo)數(shù)，并進(jìn)一步積分獲得邊界層位移厚度Reδ1、動(dòng)量厚度Reδ2等關(guān)鍵參數(shù)。

兼顧計(jì)算效率和模擬精度，采用基于錐形流假設(shè)的準(zhǔn)三維層流邊界層方程。該非線性層流邊界層方程可以表示為

(1)

(2)

可通過(guò)Lapack數(shù)據(jù)庫(kù)求解該線性方程組，進(jìn)一步得到:

(3)

迭代求解，直至解收斂，一般選擇10-15為收斂標(biāo)準(zhǔn)。

定義邊界層積分參數(shù)向量dbl，其中包含邊界層位移厚度Reδ1，動(dòng)量厚度Reδ2，形狀因子H12，邊界處的密度、馬赫數(shù)、黏性系數(shù)及速度ρe、Mae、μe、Ue。dbl的計(jì)算過(guò)程可表示為

(4)

式中：Cps為剖面壓力系數(shù)分布；xs代表剖面翼型形狀。

1.3 基于邊界層特性的轉(zhuǎn)捩準(zhǔn)則

1.3.1 流向失穩(wěn)轉(zhuǎn)捩預(yù)測(cè)

流向TS擾動(dòng)波失穩(wěn)觸發(fā)轉(zhuǎn)捩的預(yù)測(cè)采用Drela-Giles方法，該模型初始是由Gleyzes等[33]用直線斜率近似N因子和動(dòng)量厚度雷諾數(shù)Reδ2的關(guān)系提出的，即

(5)

式中：A為流向指定位置處的擾動(dòng)放大率；A0為擾動(dòng)起始點(diǎn)對(duì)應(yīng)的值；dN/dReδ2和Reδ2cr分別表示為

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

(6)

lgReδ2cr=

(7)

其中:Hk定義為可壓縮形狀因子，表示為

對(duì)于具有相似解特征的邊界層流動(dòng)，不同流向站位處x對(duì)應(yīng)的Reδ2是相同的。因而可以給出N因子和流向位置的關(guān)系式：

(8)

式中：

從而放大因子可以沿流向從擾動(dòng)開(kāi)始開(kāi)展的位置積分計(jì)算，即

(9)

與TS波對(duì)應(yīng)的轉(zhuǎn)捩閾值進(jìn)行對(duì)比，可獲得對(duì)應(yīng)的轉(zhuǎn)捩位置。

1.3.2 橫流渦失穩(wěn)轉(zhuǎn)捩準(zhǔn)則

圖1左端分別給出了物面流線(Wall streamline)和勢(shì)流流線(External streamline)。后掠角和壓力梯度的影響使得勢(shì)流流線彎曲，且壓力梯度和向心力達(dá)到平衡狀態(tài)。在邊界層內(nèi)，流向速度減小，但是壓力梯度維持不變。因而向心力和壓力梯度的不平衡會(huì)在邊界層內(nèi)部產(chǎn)生第二個(gè)方向的速度型，即橫流速度型，其垂直于勢(shì)流方向。橫流速度在物面和邊界層邊界處都為0，因而存在拐點(diǎn)(圖右端)。對(duì)于橫流擾動(dòng)，分為駐渦和行波。駐渦主要發(fā)生在低湍流度條件下，而行波在高湍流度且表面極其光滑的條件下占主導(dǎo)作用。在航空領(lǐng)域，飛行條件下都是低湍流度環(huán)境，且構(gòu)型表面存在一定的粗糙度，因而一般是駐渦占主導(dǎo)作用。

圖1 三維機(jī)翼擾動(dòng)示意圖

針對(duì)CF渦擾動(dòng)預(yù)測(cè)采用Arnal[12]基于Falkner-Skan-Cooke (FSC)邊界層相似解提出的C1準(zhǔn)則。定義基于橫流位移厚度雷諾數(shù)ReCF：

(10)

(11)

通過(guò)C1判據(jù)準(zhǔn)則，同樣可獲得對(duì)應(yīng)橫流渦誘導(dǎo)的轉(zhuǎn)捩位置。通過(guò)對(duì)比，進(jìn)一步可以獲得最終的轉(zhuǎn)捩位置，即Tp。

最終，整個(gè)轉(zhuǎn)捩位置計(jì)算過(guò)程可表示為

Tp=FT(dbl)

(12)

1.4 匹配RANS求解器的間歇因子函數(shù)計(jì)算

整個(gè)物面的轉(zhuǎn)捩位置，是通過(guò)不同轉(zhuǎn)捩計(jì)算剖面得到的轉(zhuǎn)捩位置返回RANS求解器得到的。轉(zhuǎn)捩位置通過(guò)間歇因子函數(shù)作用在SA湍流模型[34]上，實(shí)現(xiàn)轉(zhuǎn)捩模塊和RANS求解器的耦合，即

(13)

通過(guò)判斷任意物面網(wǎng)格點(diǎn)處的流動(dòng)狀態(tài)，從而確定對(duì)應(yīng)網(wǎng)格點(diǎn)處的間歇因子值。若其處于層流區(qū)，間歇因子值為0；若其處于湍流區(qū)，間歇因子值為1。這種處理方式導(dǎo)致層流/湍流區(qū)的判斷依賴于網(wǎng)格點(diǎn)的疏密，且出現(xiàn)0到1的突變，容易導(dǎo)致轉(zhuǎn)捩迭代不收斂和設(shè)計(jì)空間的不光滑。因而，需要引入間歇因子函數(shù)。采用的間歇因子函數(shù)為

(14)

式中：xtr為轉(zhuǎn)捩位置；ltr為轉(zhuǎn)捩長(zhǎng)度，由當(dāng)?shù)氐睦字Z數(shù)Re(x)和當(dāng)?shù)氐牧飨蜃鴺?biāo)x決定，即

(15)

以某典型機(jī)翼為例(圖2)說(shuō)明所有物面網(wǎng)格點(diǎn)間歇因子的計(jì)算方法。圖2中展向中間的實(shí)線代表的是不同轉(zhuǎn)捩計(jì)算剖面得到的轉(zhuǎn)捩點(diǎn)構(gòu)成的轉(zhuǎn)捩線。物面上任意網(wǎng)格點(diǎn)x0所處展向站位的轉(zhuǎn)捩位置及轉(zhuǎn)捩長(zhǎng)度，可通過(guò)線性插值表示為

(16)

(17)

式中：下標(biāo)“1”和“2”分別代表同側(cè)翼面、距離點(diǎn)x0最近的兩個(gè)轉(zhuǎn)捩計(jì)算剖面位置。將xtr0和ltr0帶入間歇因子函數(shù)，遍歷所有物面網(wǎng)格點(diǎn)，可獲得物面所有網(wǎng)格點(diǎn)對(duì)應(yīng)的間歇因子值。

圖2 轉(zhuǎn)捩線在物面插值示意圖

求得表面的間歇因子函數(shù)后，在附面層外，間歇因子定義為1。在給定邊界層距離內(nèi)，法向?qū)?yīng)的空間網(wǎng)格點(diǎn)處的間歇因子值g由其最小物面距離對(duì)應(yīng)的表面網(wǎng)格點(diǎn)確定。

1.5 耦合RANS求解器的轉(zhuǎn)捩計(jì)算

整個(gè)轉(zhuǎn)捩迭代計(jì)算過(guò)程見(jiàn)圖3。RANS求解器先進(jìn)行固定轉(zhuǎn)捩計(jì)算，流場(chǎng)殘差降低10-8量級(jí)后，將流場(chǎng)解作為轉(zhuǎn)捩模塊的輸入，進(jìn)一步獲得轉(zhuǎn)捩位置。新的轉(zhuǎn)捩位置(松弛因子)返回RANS方程，循環(huán)迭代，直至流場(chǎng)解收斂。

RANS求解器部分求解方程可表示為

(18)

式中：Q為RANS方程的狀態(tài)變量。轉(zhuǎn)捩預(yù)測(cè)問(wèn)題可以描述為

L(Tr)=Tr-Tp=0

(19)

式中：Tr為不同剖面的強(qiáng)制轉(zhuǎn)捩位置；Tp為轉(zhuǎn)捩位置預(yù)測(cè)值。最終，整個(gè)轉(zhuǎn)捩問(wèn)題可以描述為以下殘差形式：

(20)

整個(gè)轉(zhuǎn)捩預(yù)測(cè)系統(tǒng)最終要求滿足該方程的解向量，即Q和Tr。其中，X為幾何設(shè)計(jì)變量。

圖3 轉(zhuǎn)捩計(jì)算模塊示意圖

1.6 數(shù)值模擬驗(yàn)證

分別選擇跨聲速風(fēng)洞試驗(yàn)(TLFTRM01)[35]和亞聲速飛行試驗(yàn)(SLFTRM01)[13]進(jìn)行算例驗(yàn)證。風(fēng)洞試驗(yàn)和飛行試驗(yàn)?zāi)Ｐ鸵?jiàn)圖4和圖5。

圖4 TLFTRM01風(fēng)洞模型

圖5 SLFTRM01飛行試驗(yàn)?zāi)Ｐ?/p>

TLFTRM01風(fēng)洞試驗(yàn)狀態(tài)主要包含自然層流和混合層流試驗(yàn)。自然層流試驗(yàn)工況包含不同攻角，攻角范圍為[-3.69°,3.07°]，馬赫數(shù)為0.70和0.78，雷諾數(shù)為6.5×106。試驗(yàn)?zāi)Ｐ蜋C(jī)翼前緣后掠角為35°。自然層流狀態(tài)，小攻角下為CF渦失穩(wěn)觸發(fā)的層流轉(zhuǎn)捩。隨著攻角的增大，轉(zhuǎn)捩機(jī)制從CF渦轉(zhuǎn)變?yōu)門(mén)S波失穩(wěn)。文獻(xiàn)[35]基于線性穩(wěn)定性理論的轉(zhuǎn)捩預(yù)測(cè)模型結(jié)合試驗(yàn)數(shù)據(jù)標(biāo)定了該風(fēng)洞試驗(yàn)的TS波對(duì)應(yīng)的轉(zhuǎn)捩閾值為NLTS=8.7。本文Drela-Giles方法采用文獻(xiàn)標(biāo)定的TS波轉(zhuǎn)捩閾值。

圖6和圖7分給出了數(shù)值模擬與試驗(yàn)壓力分布和轉(zhuǎn)捩位置的對(duì)比，驗(yàn)證Drela-Giles-C1轉(zhuǎn)捩預(yù)測(cè)方法的預(yù)測(cè)精度。對(duì)比結(jié)果顯示，Ma=0.7不同攻角下上翼面壓力分布吻合度較高，Drela-Giles-C1轉(zhuǎn)捩預(yù)測(cè)結(jié)果與試驗(yàn)值吻合較好，僅有個(gè)別狀態(tài)的轉(zhuǎn)捩預(yù)誤差略大于5%當(dāng)?shù)叵议L(zhǎng)。因此，該方法可準(zhǔn)確預(yù)測(cè)由TS和CF渦擾動(dòng)誘導(dǎo)的轉(zhuǎn)捩現(xiàn)象。

圖6 TLFTRM01風(fēng)洞試驗(yàn)數(shù)值模擬與試驗(yàn)上翼面Cp對(duì)比(Ma=0.7)

圖7 TLFTRM01風(fēng)洞試驗(yàn)數(shù)值模擬與試驗(yàn)數(shù)據(jù)對(duì)比

亞聲速翼套飛行試驗(yàn)涉及自然層流試驗(yàn)和混合層流試驗(yàn)，本文僅關(guān)注自然層流試驗(yàn)。該試驗(yàn)狀態(tài)包含不同雷諾數(shù)和攻角，飛行高度6～7 km。由于模型前緣后掠角5°，因此主要是流向失穩(wěn)主導(dǎo)轉(zhuǎn)捩。Drela-Giles方法采用文獻(xiàn)[13]通過(guò)對(duì)比飛行試驗(yàn)數(shù)據(jù)標(biāo)定的TS波轉(zhuǎn)捩閾值9.0，作為本文TS波失穩(wěn)的轉(zhuǎn)捩判據(jù)進(jìn)行數(shù)值模擬，驗(yàn)證高雷諾數(shù)條件下Drela-Giles-C1的轉(zhuǎn)捩預(yù)測(cè)精度。

圖8給出的典型工況不同截面Cp的對(duì)比結(jié)果顯示，數(shù)值模擬與飛行試驗(yàn)測(cè)量值吻合度較高。典型工況下轉(zhuǎn)捩位置對(duì)比結(jié)果如表1所示，數(shù)值模擬與飛行試驗(yàn)數(shù)據(jù)轉(zhuǎn)捩位置誤差在6%當(dāng)?shù)叵议L(zhǎng)以內(nèi)。因此，高雷諾數(shù)飛行條件下，Drela-Giles-C1可有效預(yù)測(cè)TS波失穩(wěn)觸發(fā)的轉(zhuǎn)捩。

圖8 SLFTRM01飛行試驗(yàn)與數(shù)值模擬Cp對(duì)比

表1 SLFTRM01飛行試驗(yàn)與數(shù)值模擬數(shù)據(jù)對(duì)比

2 層流翼梯度優(yōu)化框架

2.1 考慮轉(zhuǎn)捩的耦合伴隨方程構(gòu)造

針對(duì)全湍流狀態(tài)，其伴隨方程表示為

(21)

式中：為轉(zhuǎn)捩部分對(duì)應(yīng)的殘差(見(jiàn)式(19))；Tr為轉(zhuǎn)捩部分對(duì)應(yīng)的狀態(tài)變量;φ為轉(zhuǎn)捩部分對(duì)應(yīng)的伴隨向量，從而目標(biāo)函數(shù)關(guān)于設(shè)計(jì)變量X的導(dǎo)數(shù)可寫(xiě)為

(22)

殘差關(guān)于設(shè)計(jì)變量的導(dǎo)數(shù)進(jìn)一步表示為

(23)

從而可以得到[dQ/dXdTr/?X]T:

(24)

將式(24)代入式(22)得:

(25)

進(jìn)而，得到考慮轉(zhuǎn)捩的耦合伴隨方程:

(26)

相應(yīng)目標(biāo)函數(shù)關(guān)于設(shè)計(jì)變量的梯度表示為

(27)

2.2 雅克比矩陣組裝及耦合伴隨方程求解

(28)

?L/?Q也涉及兩個(gè)系統(tǒng)的交互計(jì)算，根據(jù)鏈?zhǔn)角髮?dǎo)法則，從轉(zhuǎn)捩模塊的伴隨算子φ開(kāi)始，即

(29)

(30)

式中：xs為剖面幾何坐標(biāo)點(diǎn)；xsurf為物面幾何坐標(biāo)點(diǎn)。同時(shí)存在：

(31)

求解完相關(guān)雅克比矩陣向量后，需要求解耦合伴隨方程。一般有兩種方式求解，一種是傳統(tǒng)的LBGS(Linear Block Gauss-Seidel)算法；另外一種是CK(Coupled Krylov)算法。研究表明CK方法計(jì)算效率遠(yuǎn)高于LBGS方法(詳見(jiàn)3.3節(jié)梯度求解效率對(duì)比)。因而本文采用的CK算法求解耦合伴隨方程，見(jiàn)算法1。

算法1 耦合Krylov算法求解耦合伴隨方程1. function MULT(Z)計(jì)算雅可比向量Z2. (Z,Z)←ZT提取RANS和轉(zhuǎn)捩模塊的雅可比向量3. Y←Q TZ Y←Z并行計(jì)算對(duì)角向量4. Y←Y+Q TZ增加RANS副對(duì)象線向量5. Y←Y+Tr TZ增加轉(zhuǎn)捩模塊副對(duì)角線向量6. Y←(Y,Y)合并RANS和轉(zhuǎn)捩相關(guān)的伴隨向量7. return Y8. end function

2.3 梯度信息驗(yàn)證及效率對(duì)比

選擇翼身組合體構(gòu)型作為梯度信息和計(jì)算效率對(duì)比的研究對(duì)象。機(jī)翼前緣后掠角為20°，后緣后掠角為8.5°。模擬工況為：Re=20×106,Ma=0.78,α=0.2°。翼身組合體構(gòu)型機(jī)翼間歇因子云圖及剖面壓力分布見(jiàn)圖9。

該基準(zhǔn)算例的Free-Form Deformation (FFD)[36]幾何參數(shù)化方法控制框見(jiàn)圖10，圖中機(jī)翼部分幾何設(shè)計(jì)變量有108個(gè)。這里隨機(jī)選擇了10個(gè)幾何設(shè)計(jì)變量進(jìn)行梯度信息求解驗(yàn)證，對(duì)應(yīng)FFD框的標(biāo)號(hào)分別為21、 25、 35、 45、 56、 66、 82、 84、 92和94。將伴隨方程求解的梯度信息與采用最優(yōu)差分步長(zhǎng)的有限差分法(FD)計(jì)算得到的梯度值進(jìn)行對(duì)比。升、阻力系數(shù)(CL、CD)對(duì)設(shè)計(jì)變量的梯度對(duì)比結(jié)果見(jiàn)表2和表3。其中，Var為變量序號(hào)；Adjoint和FD分別為基于耦合伴隨方程與基于有限差分法求得的梯度；Rel err為二者的相對(duì)誤差；hopt表示基于有限差分法求得的梯度對(duì)應(yīng)的最優(yōu)時(shí)間步長(zhǎng)。基于伴隨方程的梯度與有限差分結(jié)果吻合較好，相對(duì)誤差范圍在10-3～10-6。

圖11和圖12分別為全湍流以及自由轉(zhuǎn)捩條件計(jì)算得到的機(jī)翼上翼面伴隨算子ψρ云圖。進(jìn)一步，對(duì)比有限差分、LBGS和CK算法的梯度計(jì)算效率，LBGS算法依賴于松弛因子θ，分別取0.5和1.0。表4給出的對(duì)比結(jié)果表明，針對(duì)具有上百維設(shè)計(jì)變量的求解問(wèn)題，有限差分法遠(yuǎn)無(wú)法滿足梯度優(yōu)化對(duì)梯度信息高效求解的需求。隨著松弛因子取值從0.5增大到1.0，LBGS的梯度求解時(shí)間由5 673.1 s減少到2 641.2 s，計(jì)算效率得到大幅提升。CK算法具有最高的梯度計(jì)算效率，相對(duì)松弛因子取1.0的的LBGS方法，梯度求解時(shí)間減小84.4%。對(duì)比結(jié)果顯示了CK算法求解耦合伴隨方程的高效性。

圖9 翼身組合體構(gòu)型物面間歇因子云圖和剖面壓力分布

圖10 翼身組合體FFD框

2.4 梯度優(yōu)化系統(tǒng)搭建

結(jié)合數(shù)值求解器及伴隨方程，逆距離權(quán)重插值(Inverse Distance Weighting Interpolation, IDW) 網(wǎng)格變形方法[22]、參數(shù)化方法FFD及梯度優(yōu)化器構(gòu)建層流翼梯度優(yōu)化系統(tǒng)見(jiàn)圖13。

表2 伴隨方法和有限差分法計(jì)算得到CL梯度對(duì)比

表3 伴隨方法和有限差分法計(jì)算得到CD梯度對(duì)比

圖11 全湍流計(jì)算條件下機(jī)翼上翼面ρ云圖

圖12 轉(zhuǎn)捩計(jì)算條件下機(jī)翼上翼面ρ云圖

表4 CL對(duì)100個(gè)設(shè)計(jì)變量的梯度計(jì)算花費(fèi)對(duì)比

圖13 層流翼梯度優(yōu)化設(shè)計(jì)框架

3 自然層流超臨界機(jī)翼優(yōu)化設(shè)計(jì)與分析

跨聲速狀態(tài)下，典型客機(jī)翼身組合體通常具有上百維設(shè)計(jì)變量，本節(jié)基于該類構(gòu)型，采用構(gòu)建的基于離散伴隨的梯度優(yōu)化方法，開(kāi)展自然層流超臨界機(jī)翼氣動(dòng)優(yōu)化設(shè)計(jì)研究。

3.1 優(yōu)化問(wèn)題數(shù)學(xué)描述

氣動(dòng)設(shè)計(jì)工況為：Ma=0.78,Re=20×106。機(jī)翼前緣后掠角為20°，后緣后掠角為8.5°。典型跨聲速民機(jī)飛行高度10～12 km，大氣湍流度較低。參考Mack公式[37]以及相似條件下相關(guān)研究[16]，優(yōu)化中TS波轉(zhuǎn)捩閾值取為NLTS=12，CF渦轉(zhuǎn)捩使用C1準(zhǔn)則。

以阻力系數(shù)CD最小為優(yōu)化目標(biāo)的單點(diǎn)優(yōu)化，優(yōu)化問(wèn)題描述見(jiàn)表5。設(shè)計(jì)約束包括油箱容積不減V≥Vinit，翼型厚度不小于0.9倍的初始翼型厚度t≥0.9tinit，俯仰力矩系數(shù)不小于1.28倍的初始翼型俯仰力矩系數(shù)Cm≥1.28Cmy0以及巡航升力系數(shù)約束。優(yōu)化問(wèn)題包含108個(gè)幾何設(shè)計(jì)變量(FFD框控制點(diǎn))。將攻角也作為設(shè)計(jì)變量，總設(shè)計(jì)變量個(gè)數(shù)為109。由于層流-湍流轉(zhuǎn)捩對(duì)頭部幾何尤為敏感，因而前緣FFD控制框的控制點(diǎn)分布更密集，如圖8所示。

表5 翼身組合體優(yōu)化問(wèn)題定義

翼身組合體構(gòu)型采用嵌套網(wǎng)格。為了評(píng)估網(wǎng)格量對(duì)計(jì)算結(jié)果的影響，基于L0、L1和L2這3套不同疏密程度的網(wǎng)格開(kāi)展網(wǎng)格收斂性研究，見(jiàn)圖14。L1網(wǎng)格機(jī)翼流向353個(gè)網(wǎng)格點(diǎn)，法向保證y+<1，總網(wǎng)格數(shù)為814萬(wàn)。L2網(wǎng)格在L1的基礎(chǔ)上對(duì)機(jī)翼網(wǎng)格進(jìn)行粗化。L0網(wǎng)格在L1的基礎(chǔ)上對(duì)機(jī)翼網(wǎng)格進(jìn)行細(xì)化。

圖14 具有不同網(wǎng)格量的翼身組合體嵌套網(wǎng)格

網(wǎng)格收斂性計(jì)算結(jié)果如表6所示。定升力系數(shù)條件下，密網(wǎng)格L2阻力預(yù)測(cè)207.91 counts，L0網(wǎng)格阻力預(yù)測(cè)218.12 counts，L1網(wǎng)格阻力預(yù)測(cè)209.97 counts。L1和L2阻力預(yù)測(cè)結(jié)果相差2.06 counts。從工程應(yīng)用研究角度而言，可以認(rèn)為L(zhǎng)1和L2網(wǎng)格具有相當(dāng)?shù)挠?jì)算精度水平。綜合考慮計(jì)算效率以及計(jì)算精度，選取L1網(wǎng)格進(jìn)行層流翼減阻優(yōu)化設(shè)計(jì)研究。

表6 針對(duì)初始構(gòu)型的網(wǎng)格收斂性驗(yàn)證

3.2 優(yōu)化設(shè)計(jì)結(jié)果

針對(duì)翼身組合體構(gòu)型，分別進(jìn)行了層流機(jī)翼氣動(dòng)優(yōu)化和全湍流機(jī)翼氣動(dòng)優(yōu)化。圖15給出了優(yōu)化收斂歷程，經(jīng)過(guò)45次主迭代之后，層流機(jī)翼優(yōu)化目標(biāo)函數(shù)達(dá)到收斂標(biāo)準(zhǔn)。

表7針對(duì)初始構(gòu)型和優(yōu)化構(gòu)型給出了不同優(yōu)化設(shè)計(jì)及數(shù)值模擬條件下各構(gòu)型的氣動(dòng)力系數(shù)。初始構(gòu)型和優(yōu)化構(gòu)型用“Initial”和“Optimized”代表。表中第1列構(gòu)型名稱的前綴“LT”和“FT”分別表示層流氣動(dòng)優(yōu)化和全湍流氣動(dòng)優(yōu)化。自由轉(zhuǎn)捩和全湍流數(shù)值模擬條件在表中第2列分別用“LT”和“FT”區(qū)分。對(duì)比結(jié)果表明自由轉(zhuǎn)捩數(shù)值模擬條件下，相比初始構(gòu)型(“Initial-LT”)，層流優(yōu)化構(gòu)型(“LT-Optimized-LT”)總阻力降低22.0 counts，減小10.48%。其中壓差阻力CDp降低了13.9 counts，減小10.66%，摩擦阻力CDf降低了8.1 counts，減小10.19%。壓差阻力系數(shù)的降低量占減阻總收益的63%，摩擦阻力系數(shù)的降低量占37%。

圖15 翼身組合體層流及全湍流機(jī)翼優(yōu)化收斂歷程

表7 氣動(dòng)力系數(shù)優(yōu)化結(jié)果對(duì)比

圖16和圖17給出了初始構(gòu)型及層流優(yōu)化構(gòu)型在自由轉(zhuǎn)捩數(shù)值模擬條件下的摩擦阻力系數(shù)分布對(duì)比。圖中摩擦阻力突增的區(qū)域即為轉(zhuǎn)捩線。優(yōu)化后的層流機(jī)翼上、下翼面轉(zhuǎn)捩位置推遲明顯。

圖16 初始構(gòu)型上下翼面摩擦阻力系數(shù)云圖

圖17 層流優(yōu)化構(gòu)型上下翼面摩擦阻力系數(shù)云圖

圖18給出了初始構(gòu)型與層流優(yōu)化構(gòu)型上翼面壓力分布云圖以及選取的6個(gè)剖面站位處壓力系數(shù)Cp分布和剖面翼型對(duì)比。圖中空心圓標(biāo)注了各截面的轉(zhuǎn)捩位置。優(yōu)化后，上翼面轉(zhuǎn)捩位置明顯推遲，尤其在外翼段，直到激波處才觸發(fā)轉(zhuǎn)捩。下翼面轉(zhuǎn)捩位置也得到不同程度的推遲。內(nèi)翼段由于當(dāng)?shù)仄拭胬字Z數(shù)較大，易引起橫流轉(zhuǎn)捩，轉(zhuǎn)捩推遲較為困難。

剖面壓力系數(shù)分布對(duì)比表明，優(yōu)化后層流機(jī)翼內(nèi)翼段激波明顯減弱，剖面A頭部峰值顯著降低，抑制CF渦在頭部過(guò)快失穩(wěn)。除剖面A外，內(nèi)翼段(剖面A-C)上翼面頭部峰值之后的順壓梯度都有不同程度的減弱，抑制CF渦的快速發(fā)展。下翼面前緣負(fù)壓峰值適當(dāng)增大，削弱前緣附近順壓梯度大小，減小橫流擾動(dòng)。外翼段頭部峰值之后的順壓梯度略有降低，在20%c之后，順壓梯度略有增大。剖面E激波位置略有后移，在不顯著增大激波強(qiáng)度的條件下延長(zhǎng)上翼面層流區(qū)。

為了說(shuō)明自然層流超臨界機(jī)翼技術(shù)巨大的減阻能力，將相同設(shè)計(jì)目標(biāo)/約束條件下全湍流優(yōu)化設(shè)計(jì)結(jié)果進(jìn)行對(duì)比，對(duì)比結(jié)果見(jiàn)表7。全湍流數(shù)值模擬條件下，相比初始構(gòu)型(“Initial-FT”)，全湍流優(yōu)化構(gòu)型(“FT-Optimized-FT”)總阻力降低13.2 counts，減小6.17%。剖面壓力分布對(duì)比(圖19)表明，優(yōu)化設(shè)計(jì)構(gòu)型為無(wú)激波構(gòu)型，前緣附近的負(fù)壓峰值略有降低。全湍流優(yōu)化構(gòu)型6.17%的減阻收益，遠(yuǎn)小于層流優(yōu)化設(shè)計(jì)構(gòu)型10.48%的減阻收益。

相同構(gòu)型不同數(shù)值模擬條件(自由轉(zhuǎn)捩模擬/全湍流模擬)對(duì)應(yīng)的阻力系數(shù)不同，見(jiàn)表7。為了得到更準(zhǔn)確的減阻收益，除了全湍流模擬條件外，還需進(jìn)行自由轉(zhuǎn)捩數(shù)值模擬。針對(duì)全湍流優(yōu)化構(gòu)型進(jìn)行自由轉(zhuǎn)捩條件下的數(shù)值模擬(“FT-Optimized-LT”)，摩擦阻力系數(shù)云圖見(jiàn)圖20。由于前緣負(fù)壓峰值略有降低，減小了頭部順壓梯度，一定程度削弱了CF渦發(fā)展，使得上翼面維持的一定的層流區(qū)。自由轉(zhuǎn)捩數(shù)值模擬條件，相比初始構(gòu)型(“Initial-LT”)，全湍流優(yōu)化構(gòu)型(“FT-Optimized-LT”)總阻力系數(shù)降低17.6 counts，其中壓差阻力系數(shù)降低了15.7 counts，摩擦阻力系數(shù)降低了1.9 counts。但是，17.6counts的減阻收益仍然遠(yuǎn)小于層流優(yōu)化構(gòu)型(“LT-Optimized-LT”)22.0 counts的阻力系數(shù)減小量。

圖18 層流機(jī)翼優(yōu)化前后壓力系數(shù)分布對(duì)比

圖19 全湍流機(jī)翼優(yōu)化前后壓力系數(shù)分布對(duì)比

圖20 “FT-Optimized-LT”上下翼面摩擦阻力系數(shù)云圖

總之，優(yōu)化設(shè)計(jì)結(jié)果表明，本文發(fā)展的基于離散伴隨的梯度優(yōu)化方法，能夠有效處理多種轉(zhuǎn)捩機(jī)制并存的復(fù)雜三維層流翼優(yōu)化問(wèn)題。對(duì)于具有中/高雷諾數(shù)的典型民機(jī)構(gòu)型，相比傳統(tǒng)超臨界機(jī)翼，自然層流超臨界機(jī)翼具有更大的減阻潛力。

3.3 轉(zhuǎn)捩抑制機(jī)制分析

對(duì)于20°后掠角、雷諾數(shù)為2 000萬(wàn)的跨聲速構(gòu)型，這種多種轉(zhuǎn)捩機(jī)制并存的復(fù)雜三維層流翼優(yōu)化問(wèn)題，需要借助適宜的壓力分布特征有效平衡TS波和CF渦擾動(dòng)發(fā)展，達(dá)到推遲轉(zhuǎn)捩的目的。

圖18給出的6個(gè)剖面A-F，初始構(gòu)型(“Initial-LT”)上翼面從內(nèi)翼段到外翼段對(duì)應(yīng)的觸發(fā)層流-湍流轉(zhuǎn)捩的類型分別為：CF渦—CF渦—CF渦—CF渦—CF渦—CF渦；下翼面分別為：CF渦—CF渦—CF渦—TS波—TS波—TS波。層流優(yōu)化構(gòu)型(“LT-Optimized-LT”)上翼面對(duì)應(yīng)的從內(nèi)翼段到外翼段的轉(zhuǎn)捩類型分別為：CF渦—CF渦—CF渦—TS波—CF渦—CF渦；而下翼面對(duì)應(yīng)的分別為：CF渦—CF渦—CF渦—CF渦—CF渦—TS波。

選取機(jī)翼典型剖面，從擾動(dòng)波的發(fā)展及邊界層信息角度揭示層流翼優(yōu)化構(gòu)型的轉(zhuǎn)捩抑制機(jī)制。圖21和圖22給出了展向站位z=6.12 m(圖18中截面B)和z=15.67 m(圖18中截面E)對(duì)應(yīng)的TS波和CF渦擾動(dòng)發(fā)展。截面z=6.12 m處的上翼面，初始構(gòu)型和優(yōu)化構(gòu)型TS擾動(dòng)波分別在(x/c)TS=0.061和(x/c)TS=0.5達(dá)到轉(zhuǎn)捩閾值。初始構(gòu)型頭部峰值之后有較強(qiáng)的逆壓梯度，促使TS波擾動(dòng)迅速發(fā)展。相比之下，優(yōu)化構(gòu)型峰值之后為較長(zhǎng)的順壓梯度區(qū)，有效抑制了TS波的發(fā)展。對(duì)于上翼面CF渦，優(yōu)化前對(duì)CF渦達(dá)到轉(zhuǎn)捩閾值的位置分別為(x/c)CF=0.03、(x/c)CF=0.25。初始構(gòu)型負(fù)壓峰值前的大順壓梯度導(dǎo)致CF渦速度失穩(wěn)。優(yōu)化構(gòu)型峰值前順壓梯度減弱，一定程度上抑制了前緣附近CF渦的發(fā)展。綜合上翼面TS波和CF渦的發(fā)展，初始構(gòu)型和優(yōu)化構(gòu)型轉(zhuǎn)捩位置分別為xtr/c=0.03，xtr/c=0.25，均為CF渦失穩(wěn)觸發(fā)的轉(zhuǎn)捩。

圖21 層流翼優(yōu)化前后z=6.12 m剖面轉(zhuǎn)捩位置判斷

圖22 層流翼優(yōu)化前后z=15.67 m剖面轉(zhuǎn)捩位置判斷

截面z=6.12 m處的下翼面，相比初始構(gòu)型，優(yōu)化構(gòu)型以適當(dāng)提高前緣峰值，減小前緣順壓力梯度的方式抑制CF渦在前緣的發(fā)展。順壓梯度的減小減小了對(duì)TS波的抑制能力，導(dǎo)致TS波擾動(dòng)放大率增大。計(jì)算結(jié)果表明，初始構(gòu)型在(x/c)TS=0.6處發(fā)生層流分離，優(yōu)化構(gòu)型TS波在(x/c)TS=0.588處達(dá)到轉(zhuǎn)捩閾值。CF渦分別在初始構(gòu)型和優(yōu)化構(gòu)型的(x/c)CF=0.139和(x/c)CF=0.187達(dá)到轉(zhuǎn)捩閾值。綜合下翼面TS波和CF渦的發(fā)展，初始構(gòu)型和優(yōu)化構(gòu)型轉(zhuǎn)捩位置分別為xtr/c=0.139，xtr/c=0.187，且都是由CF渦失穩(wěn)誘導(dǎo)的轉(zhuǎn)捩。

站位z=15.67 m處下翼面，初始構(gòu)型前緣大逆壓梯度(對(duì)比圖18所示壓力分布)導(dǎo)致TS波在(x/c)TS=0.033變達(dá)到轉(zhuǎn)捩閾值。優(yōu)化構(gòu)型減小了前緣逆壓梯度，顯著削弱了TS波的發(fā)展，最終在(x/c)TS=0.427處達(dá)到轉(zhuǎn)捩閾值。初始構(gòu)型和優(yōu)化構(gòu)型下翼面CF渦達(dá)到轉(zhuǎn)捩閾值的位置分別為(x/c)CF=0.2和(x/c)CF=0.388。

圖22顯示優(yōu)化構(gòu)型的橫流位移厚度雷諾數(shù)在前緣較小，但之后較大，轉(zhuǎn)捩位置推遲的原因是由于壓力梯度的改變，使得對(duì)應(yīng)的CF渦轉(zhuǎn)捩閾值增大，從而CF渦不易失穩(wěn)。最終機(jī)翼下翼面，初始構(gòu)型轉(zhuǎn)捩位置為xtr/c=0.033，優(yōu)化構(gòu)型轉(zhuǎn)捩位置為xtr/c=0.388，均為CF渦失穩(wěn)誘導(dǎo)的轉(zhuǎn)捩。

層流機(jī)翼的優(yōu)化過(guò)程也反應(yīng)了邊界層內(nèi)部流動(dòng)參數(shù)的變化，這里選取展向站位z=6.12 m下翼面邊界層內(nèi)部變化為例進(jìn)行論述。圖23～圖25給出了展向站位z=6.12 m下翼面橫流速度型及其一階和二階導(dǎo)數(shù)。其中，圖23～圖25的橫坐標(biāo)u為橫流速度，縱坐標(biāo)d為距壁面的垂直距離。橫流擾動(dòng)的大小可由邊界層橫流速度型最大值、橫流速度型拐點(diǎn)和橫流速度型一階導(dǎo)數(shù)在拐點(diǎn)處的值表征。圖23顯示，除x/c=0.12站位，優(yōu)化后層流邊界層對(duì)應(yīng)的橫流速度型的最大值整體降低。離開(kāi)機(jī)翼前緣，橫流速度型大小明顯降低。圖24對(duì)比了層流機(jī)翼優(yōu)化前后橫流速度型的二階導(dǎo)數(shù)。優(yōu)化后的拐點(diǎn)基本低于初始構(gòu)型。在圖25對(duì)比了橫流速度型的一階導(dǎo)數(shù)，并用實(shí)心圓點(diǎn)標(biāo)出了一階導(dǎo)數(shù)在拐點(diǎn)處的值。該值表征了切應(yīng)力大小，值越大表明橫流擾動(dòng)越強(qiáng)。對(duì)比結(jié)果顯示，相比初始構(gòu)型，層流優(yōu)化構(gòu)型橫流速度型拐點(diǎn)處的橫流速度梯度略小。結(jié)合邊界層橫流速度型的3個(gè)參數(shù)對(duì)比，可得出優(yōu)化后橫流擾動(dòng)明顯降低的結(jié)論，與橫流轉(zhuǎn)捩判據(jù)參數(shù)和壓力分布變化一致。

圖23 層流翼優(yōu)化前后z=6.12 m剖面下翼面橫流速度型對(duì)比

圖24 層流翼優(yōu)化前后z=6.12 m剖面下翼面橫流速度型二階導(dǎo)數(shù)對(duì)比

圖25 層流翼優(yōu)化前后z=6.12 m剖面下翼面橫流速度型一階導(dǎo)數(shù)對(duì)比

總之，發(fā)展的層流翼梯度優(yōu)化設(shè)計(jì)方法有效權(quán)衡了TS波和CF渦擾動(dòng)，顯著推遲了轉(zhuǎn)捩，反映了邊界層內(nèi)不同的轉(zhuǎn)捩機(jī)制。

4 結(jié) 論

基于Drela-Giles和C1準(zhǔn)則以及離散伴隨理論，發(fā)展了可有效處理多種轉(zhuǎn)捩機(jī)制并存的復(fù)雜三維層流翼優(yōu)化問(wèn)題的高效梯度優(yōu)化方法。針對(duì)具有20°前緣后掠角的典型支線客機(jī)開(kāi)展梯度優(yōu)化設(shè)計(jì)，探究了自然層流超臨界機(jī)翼氣動(dòng)設(shè)計(jì)原理。主要工作總結(jié)如下：

1) 耦合RANS求解器和基于Drela-Giles和C1準(zhǔn)則的轉(zhuǎn)捩模塊，通過(guò)嚴(yán)格線性插值計(jì)算表面間歇因子，實(shí)現(xiàn)氣動(dòng)力的準(zhǔn)確計(jì)算。轉(zhuǎn)捩預(yù)測(cè)結(jié)果與典型風(fēng)洞和飛行試驗(yàn)數(shù)據(jù)吻合度較高，驗(yàn)證了建立的轉(zhuǎn)捩預(yù)測(cè)方法可準(zhǔn)確捕捉TS和CF不穩(wěn)定誘導(dǎo)的轉(zhuǎn)捩現(xiàn)象。

2) 基于離散伴隨理論，嚴(yán)格推導(dǎo)了考慮轉(zhuǎn)捩的耦合伴隨方程，從理論上保證了相關(guān)梯度計(jì)算的準(zhǔn)確性。對(duì)于新增的雅可比矩陣，采用無(wú)矩陣存儲(chǔ)技術(shù)、鏈?zhǔn)角髮?dǎo)及混合自動(dòng)微分方法進(jìn)行計(jì)算。發(fā)展了基于CK算法的耦合伴隨方程高效求解方法。相比典型LBGS算法，梯度求解效率提高了84.4%。

3) 基于建立的層流翼梯度優(yōu)化方法，以具有20°前緣后掠角的典型客機(jī)翼身組合體構(gòu)型，開(kāi)展了層流優(yōu)化設(shè)計(jì)研究。層流翼優(yōu)化結(jié)果表明，通過(guò)梯度優(yōu)化，可實(shí)現(xiàn)減阻10.48%，上下翼面的轉(zhuǎn)捩位置都得到了顯著推遲。相比層流翼，全湍流優(yōu)化結(jié)果僅減阻6.17%，顯示了層流超臨界機(jī)翼巨大的減阻潛力。

4) 優(yōu)化前后擾動(dòng)波的發(fā)展以及邊界層內(nèi)橫流速度型及其一階、二階導(dǎo)數(shù)對(duì)比表明，發(fā)展的層流翼梯度優(yōu)化方法可有效權(quán)衡TS波和CF渦擾動(dòng)，達(dá)到推遲轉(zhuǎn)捩，實(shí)現(xiàn)氣動(dòng)阻力的優(yōu)化目標(biāo)。