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

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

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

3. 西安交通大學 航天學院，西安 710049

4. 中國空氣動力研究與發(fā)展中心，綿陽 621000

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

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

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

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

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

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

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

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

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

1.1 RANS求解器

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

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

1.2 層流邊界層方程求解

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

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

(1)

(2)

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

(3)

迭代求解，直至解收斂，一般選擇10-15為收斂標準。

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

(4)

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

1.3 基于邊界層特性的轉捩準則

1.3.1 流向失穩(wěn)轉捩預測

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

(5)

式中：A為流向指定位置處的擾動放大率；A0為擾動起始點對應的值；dN/dReδ2和Reδ2cr分別表示為

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

(6)

lgReδ2cr=

(7)

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

對于具有相似解特征的邊界層流動，不同流向站位處x對應的Reδ2是相同的。因而可以給出N因子和流向位置的關系式：

(8)

式中：

從而放大因子可以沿流向從擾動開始開展的位置積分計算，即

(9)

與TS波對應的轉捩閾值進行對比，可獲得對應的轉捩位置。

1.3.2 橫流渦失穩(wěn)轉捩準則

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

圖1 三維機翼擾動示意圖

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

(10)

(11)

通過C1判據(jù)準則，同樣可獲得對應橫流渦誘導的轉捩位置。通過對比，進一步可以獲得最終的轉捩位置，即Tp。

最終，整個轉捩位置計算過程可表示為

Tp=FT(dbl)

(12)

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

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

(13)

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

(14)

式中：xtr為轉捩位置；ltr為轉捩長度，由當?shù)氐睦字Z數(shù)Re(x)和當?shù)氐牧飨蜃鴺藊決定，即

(15)

以某典型機翼為例(圖2)說明所有物面網格點間歇因子的計算方法。圖2中展向中間的實線代表的是不同轉捩計算剖面得到的轉捩點構成的轉捩線。物面上任意網格點x0所處展向站位的轉捩位置及轉捩長度，可通過線性插值表示為

(16)

(17)

式中：下標“1”和“2”分別代表同側翼面、距離點x0最近的兩個轉捩計算剖面位置。將xtr0和ltr0帶入間歇因子函數(shù)，遍歷所有物面網格點，可獲得物面所有網格點對應的間歇因子值。

圖2 轉捩線在物面插值示意圖

求得表面的間歇因子函數(shù)后，在附面層外，間歇因子定義為1。在給定邊界層距離內，法向對應的空間網格點處的間歇因子值g由其最小物面距離對應的表面網格點確定。

1.5 耦合RANS求解器的轉捩計算

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

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

(18)

式中：Q為RANS方程的狀態(tài)變量。轉捩預測問題可以描述為

L(Tr)=Tr-Tp=0

(19)

式中：Tr為不同剖面的強制轉捩位置；Tp為轉捩位置預測值。最終，整個轉捩問題可以描述為以下殘差形式：

(20)

整個轉捩預測系統(tǒng)最終要求滿足該方程的解向量，即Q和Tr。其中，X為幾何設計變量。

圖3 轉捩計算模塊示意圖

1.6 數(shù)值模擬驗證

分別選擇跨聲速風洞試驗(TLFTRM01)[35]和亞聲速飛行試驗(SLFTRM01)[13]進行算例驗證。風洞試驗和飛行試驗模型見圖4和圖5。

圖4 TLFTRM01風洞模型

圖5 SLFTRM01飛行試驗模型

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

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

圖6 TLFTRM01風洞試驗數(shù)值模擬與試驗上翼面Cp對比(Ma=0.7)

圖7 TLFTRM01風洞試驗數(shù)值模擬與試驗數(shù)據(jù)對比

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

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

圖8 SLFTRM01飛行試驗與數(shù)值模擬Cp對比

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

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

2.1 考慮轉捩的耦合伴隨方程構造

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

(21)

式中：為轉捩部分對應的殘差(見式(19))；Tr為轉捩部分對應的狀態(tài)變量;φ為轉捩部分對應的伴隨向量，從而目標函數(shù)關于設計變量X的導數(shù)可寫為

(22)

殘差關于設計變量的導數(shù)進一步表示為

(23)

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

(24)

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

(25)

進而，得到考慮轉捩的耦合伴隨方程:

(26)

相應目標函數(shù)關于設計變量的梯度表示為

(27)

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

(28)

?L/?Q也涉及兩個系統(tǒng)的交互計算，根據(jù)鏈式求導法則，從轉捩模塊的伴隨算子φ開始，即

(29)

(30)

式中：xs為剖面幾何坐標點；xsurf為物面幾何坐標點。同時存在：

(31)

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

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

2.3 梯度信息驗證及效率對比

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

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

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

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

圖10 翼身組合體FFD框

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

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

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

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

圖11 全湍流計算條件下機翼上翼面ρ云圖

圖12 轉捩計算條件下機翼上翼面ρ云圖

表4 CL對100個設計變量的梯度計算花費對比

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

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

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

3.1 優(yōu)化問題數(shù)學描述

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

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

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

翼身組合體構型采用嵌套網格。為了評估網格量對計算結果的影響，基于L0、L1和L2這3套不同疏密程度的網格開展網格收斂性研究，見圖14。L1網格機翼流向353個網格點，法向保證y+<1，總網格數(shù)為814萬。L2網格在L1的基礎上對機翼網格進行粗化。L0網格在L1的基礎上對機翼網格進行細化。

圖14 具有不同網格量的翼身組合體嵌套網格

網格收斂性計算結果如表6所示。定升力系數(shù)條件下，密網格L2阻力預測207.91 counts，L0網格阻力預測218.12 counts，L1網格阻力預測209.97 counts。L1和L2阻力預測結果相差2.06 counts。從工程應用研究角度而言，可以認為L1和L2網格具有相當?shù)挠嬎憔人健＞C合考慮計算效率以及計算精度，選取L1網格進行層流翼減阻優(yōu)化設計研究。

表6 針對初始構型的網格收斂性驗證

3.2 優(yōu)化設計結果

針對翼身組合體構型，分別進行了層流機翼氣動優(yōu)化和全湍流機翼氣動優(yōu)化。圖15給出了優(yōu)化收斂歷程，經過45次主迭代之后，層流機翼優(yōu)化目標函數(shù)達到收斂標準。

表7針對初始構型和優(yōu)化構型給出了不同優(yōu)化設計及數(shù)值模擬條件下各構型的氣動力系數(shù)。初始構型和優(yōu)化構型用“Initial”和“Optimized”代表。表中第1列構型名稱的前綴“LT”和“FT”分別表示層流氣動優(yōu)化和全湍流氣動優(yōu)化。自由轉捩和全湍流數(shù)值模擬條件在表中第2列分別用“LT”和“FT”區(qū)分。對比結果表明自由轉捩數(shù)值模擬條件下，相比初始構型(“Initial-LT”)，層流優(yō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 翼身組合體層流及全湍流機翼優(yōu)化收斂歷程

表7 氣動力系數(shù)優(yōu)化結果對比

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

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

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

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

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

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

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

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

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

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

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

3.3 轉捩抑制機制分析

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

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

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

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

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

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

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

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

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

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

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

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

總之，發(fā)展的層流翼梯度優(yōu)化設計方法有效權衡了TS波和CF渦擾動，顯著推遲了轉捩，反映了邊界層內不同的轉捩機制。

4 結 論

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

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

2) 基于離散伴隨理論，嚴格推導了考慮轉捩的耦合伴隨方程，從理論上保證了相關梯度計算的準確性。對于新增的雅可比矩陣，采用無矩陣存儲技術、鏈式求導及混合自動微分方法進行計算。發(fā)展了基于CK算法的耦合伴隨方程高效求解方法。相比典型LBGS算法，梯度求解效率提高了84.4%。

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

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