0 引 言

流動由層流轉(zhuǎn)變成湍流的過程稱為轉(zhuǎn)捩，轉(zhuǎn)捩現(xiàn)象普遍存在于流體流動中[1]。例如，大型民機、高空無人機、高空螺旋槳、風力機葉片及微型飛行器表面都存在流動轉(zhuǎn)捩現(xiàn)象。由于影響因素眾多、物理機理復(fù)雜，轉(zhuǎn)捩問題一直是流體力學中尚未完全解決的前沿問題之一[2]。層流和湍流兩種流態(tài)下物體表面的摩擦阻力、熱傳導(dǎo)速率和流動分離位置等特性大為不同，故轉(zhuǎn)捩位置直接影響飛行器氣動力計算的精確性，其中對阻力預(yù)測的影響尤為顯著?？罩锌蛙嚬緦＜抑赋鑫磥硪岣叽笮惋w機的空氣動力學效率，大幅度降低油耗和污染排放，層流減阻技術(shù)是一個非常有前途的發(fā)展方向[3]。因此，發(fā)展高效可靠的轉(zhuǎn)捩預(yù)測方法是飛行器減阻設(shè)計的關(guān)鍵。

隨著計算機技術(shù)的高速發(fā)展，大渦模擬(LES)[4]和直接數(shù)值模擬(DNS)[5]方法已經(jīng)能夠?qū)α鲃愚D(zhuǎn)捩進行精細化研究。但由于其計算量過大，還不適合廣泛的工程應(yīng)用。目前針對層流減阻設(shè)計常用的轉(zhuǎn)捩預(yù)測方法為：基于線性穩(wěn)定性理論的eN方法[6]和基于流場當?shù)刈兞康摩?Reθ轉(zhuǎn)捩模型[7]。最近一種基于動模態(tài)分解的DMD/eN轉(zhuǎn)捩預(yù)測方法[8]也被提出，并展現(xiàn)出良好的發(fā)展前景。下面分別對這三種方法進行詳細介紹。

基于線性穩(wěn)定性理論的eN方法是目前在工程實踐中得到最廣泛認可的轉(zhuǎn)捩預(yù)測方法[6]，eN方法由Smith[9]和Ingen[10]首次提出，目前國外研究機構(gòu)已經(jīng)圍繞翼型、機翼邊界層流動轉(zhuǎn)捩預(yù)測開展了大量的研究工作[11-26]。而針對翼身組合體等復(fù)雜三維外形，國外開展的邊界層轉(zhuǎn)捩預(yù)測方法研究相對較少，且主要使用簡化的eN數(shù)據(jù)庫轉(zhuǎn)捩預(yù)測方法。如法國宇航院(ONERA)在其雷諾平均Navier-Stokes(RANS)方程求解器elsA中加入了Arnal[15]和Casalis[18]發(fā)展的eN數(shù)據(jù)庫方法，分別用來計算TS波和CF波誘導(dǎo)的轉(zhuǎn)捩，并開展了帶增升裝置的翼身組合體外形的轉(zhuǎn)捩預(yù)測研究[27]。德國宇航院(DLR)在其結(jié)構(gòu)化網(wǎng)格RANS求解器FLOWer[28]和非結(jié)構(gòu)/混合網(wǎng)格RANS求解器TAU[29]中加入了eN數(shù)據(jù)庫方法，開展了三維復(fù)雜外形的轉(zhuǎn)捩預(yù)測研究。美國斯坦福大學的Lee和Jameson[30]開展了基于耦合eN數(shù)據(jù)庫方法的RANS求解器的層流機翼優(yōu)化設(shè)計工作。eN數(shù)據(jù)庫方法并不求解線性穩(wěn)定性方程，而只是根據(jù)邊界層的形狀因子、速度型拐點信息等部分特征參數(shù)，在數(shù)據(jù)庫中查詢事先求解的平板邊界層不穩(wěn)定擾動累積放大因子的包絡(luò)線。目前，國內(nèi)已開展的轉(zhuǎn)捩預(yù)測研究工作主要針對翼型、機翼繞流進行，也鮮有使用eN方法對翼身組合體等復(fù)雜三維外形開展轉(zhuǎn)捩預(yù)測的研究工作。在使用穩(wěn)定性理論和eN方法進行后掠機翼邊界層穩(wěn)定性及轉(zhuǎn)捩預(yù)測研究方面，張坤[31]通過線性穩(wěn)定性理論(Linear Stability Theory，LST)，左歲寒[32]通過線性拋物化穩(wěn)定性方程(Linear Parabolized Stability Equations，LPSE)開展了機翼邊界層轉(zhuǎn)捩預(yù)測研究。孫朋朋[33]、靖振榮[34]和黃章峰[35]等通過LST開展了一系列的后掠機翼邊界層穩(wěn)定性及轉(zhuǎn)捩研究，以及使用擴展的O-S方程(Extended Orr-Sommerfeld equation，EOS)和LPSE進行了機翼邊界層的橫流穩(wěn)定性分析和轉(zhuǎn)捩預(yù)測研究。韓忠華等[36]開展了基于線性穩(wěn)定性理論的雙eN方法耦合RANS求解器進行自然層流后掠機翼優(yōu)化設(shè)計的工作。

基于流場當?shù)刈兞康摩?Reθ轉(zhuǎn)捩模型是由Menter和Langtry等[7]提出的經(jīng)驗預(yù)測模型。該模型嚴格基于流場當?shù)刈兞繕?gòu)建，與現(xiàn)代CFD方法兼容，在邊界層流動轉(zhuǎn)捩模擬中獲得了廣泛的應(yīng)用。經(jīng)典的γ-Reθ模型只能預(yù)測流向TS波不穩(wěn)定性導(dǎo)致的轉(zhuǎn)捩，為了能夠預(yù)測后掠機翼上橫流CF波不穩(wěn)定性主導(dǎo)的轉(zhuǎn)捩，需要在現(xiàn)有轉(zhuǎn)捩模型的基礎(chǔ)上添加橫流轉(zhuǎn)捩預(yù)測準則。德國宇航院的Krumbein等[37-38]提出了橫流雷諾數(shù)和當?shù)睾舐咏堑母拍?，實現(xiàn)了橫流雷諾數(shù)的當?shù)厍蠼?，并進行了C1準則和γ-Reθ模型的耦合研究。韓國的Choi和Kwon等[39]對當?shù)貦M流雷諾數(shù)和當?shù)睾舐咏堑母拍钸M行了討論，將一些非當?shù)刈兞坑猛耆數(shù)刈兞看?，對模型做了進一步改進。國內(nèi)西北工業(yè)大學的徐家寬等[40]開展了基于γ-Reθ模型和C1準則的橫流轉(zhuǎn)捩研究。

基于動模態(tài)分解的DMD/eN方法是一種新的轉(zhuǎn)捩預(yù)測方法，由韓忠華等[8,41-42]于2017年最新提出。該方法通過DMD[43-46]來進行穩(wěn)定性分析，并結(jié)合eN方法實現(xiàn)轉(zhuǎn)捩預(yù)測。DMD/eN方法用于翼型轉(zhuǎn)捩預(yù)測與其他轉(zhuǎn)捩預(yù)測方法的主要區(qū)別是，不需要平行流假設(shè)，也不需要求解穩(wěn)定性或經(jīng)驗關(guān)系式方程，而是直接從CFD數(shù)值模擬的結(jié)果中提取出流動穩(wěn)定性信息。因此該方法是一種數(shù)據(jù)驅(qū)動的方法。目前該方法已被應(yīng)用于若干翼型流動轉(zhuǎn)捩預(yù)測算例，將計算結(jié)果與實驗值和其他轉(zhuǎn)捩預(yù)測方法的計算結(jié)果進行了比較，結(jié)果吻合良好，表明該方法具備預(yù)測翼型流動轉(zhuǎn)捩的能力。

本文對這三種實用的轉(zhuǎn)捩預(yù)測方法進行了詳細介紹，展示了所發(fā)展方法的有效性，并在原有研究的基礎(chǔ)上進一步開展了改進研究，主要工作體現(xiàn)在：1) 針對民機后掠機翼的橫流轉(zhuǎn)捩問題，對經(jīng)典γ-Reθ模型添加了當?shù)谻1橫流轉(zhuǎn)捩預(yù)測準則，形成了可以預(yù)測流向和橫流轉(zhuǎn)捩的三維γ-Reθ-C1模型；2) 針對DMD/eN方法中的DMD分析進行了完善，避免了數(shù)值誤差對轉(zhuǎn)捩位置的影響，進一步增強了該方法的魯棒性，并通過NLF0416自然層流翼型轉(zhuǎn)捩預(yù)測算例進行了驗證；3) 基于雙eN轉(zhuǎn)捩預(yù)測方法，開展了設(shè)計升力系數(shù)為0.5的中短程民機的自然層流機翼減阻優(yōu)化設(shè)計工作，充分展示了該方法的工程實用價值。

1 基于線性穩(wěn)定性理論的雙eN轉(zhuǎn)捩預(yù)測方法

本文使用多塊結(jié)構(gòu)化網(wǎng)格和三維雷諾平均Navier-Stokes方程求解器，耦合邊界層方程求解和基于線性穩(wěn)定性理論的雙eN方法，發(fā)展了一套可同時計及Tollmien-Schlichting波和橫流不穩(wěn)定性擾動誘導(dǎo)轉(zhuǎn)捩的翼身組合體流動轉(zhuǎn)捩自動判斷方法。eN方法的思想是對邊界層流動中的小擾動進行線性穩(wěn)定性分析，如果擾動逐漸衰減，則是穩(wěn)定的，如果擾動被放大，則是不穩(wěn)定的。對于不穩(wěn)定的擾動，從其開始放大處起，沿下游方向計算其累積的線性放大倍數(shù)，當累積放大倍數(shù)到達擾動開始放大處振幅的eN倍時，認為轉(zhuǎn)捩發(fā)生?；谄叫辛骷僭O(shè)的三維小擾動形式為:

(1)

圖1 空間放大理論中擾動波數(shù)矢量k與擾動增長矢量AFig.1 Disturbance wave number vector and disturbanceamplification vector in spatial amplification theory

(2)

基于線性穩(wěn)定性理論的eN方法，對機翼邊界層信息求解線性穩(wěn)定性方程，進行擾動累積放大因子的計算，而不對數(shù)據(jù)庫進行查詢，因此不存在近似和簡化。eN方法進行邊界層轉(zhuǎn)捩位置的判斷分兩步進行：首先，計算中性曲線，即擾動的穩(wěn)定區(qū)域和不穩(wěn)定區(qū)域的分界線。通過中性曲線就可以確定不穩(wěn)定擾動的頻率范圍、波數(shù)范圍和不同擾動各自開始被放大的位置；然后，按照適當?shù)姆e分策略，對這些頻率的擾動放大率從中性曲線下半支出發(fā)，沿著擾動傳播路徑進行積分，得到擾動累積放大因子N。這兩步都需要求解三維可壓縮線性穩(wěn)定性方程，其數(shù)值求解方法為：(a) 采用連續(xù)法獲得具有較高精確度的特征變量初值；(b) 將高階的線性穩(wěn)定性控制方程組轉(zhuǎn)化為一階的控制方程組；(c) 采用中心差分方法離散控制方程組；(d) 采用牛頓法線化離散后的控制方程組；(e) 采用塊矩陣消去法求解方程組。

1.1 RANS求解器與轉(zhuǎn)捩預(yù)測方法耦合流程

采用耦合雙eN方法的三維RANS求解器PMNS3D開展面向翼身組合體外形的機翼邊界層轉(zhuǎn)捩自動判斷方法研究。求解器主要由3個模塊構(gòu)成：1) 三維RANS求解器模塊。本模塊采用多塊結(jié)構(gòu)化網(wǎng)格、有限體積法、LU-SGS(Lower-Upper Symmetric Gauss-Seidel)時間推進、多重網(wǎng)格加速收斂技術(shù)。湍流流動模擬采用Spalart-Allmaras(S-A)湍流模型。此模塊對翼身組合體黏性繞流進行數(shù)值模擬，并從計算結(jié)果中提取機翼邊界層外邊界流向速度ue、展向速度we、溫度Te等流動參數(shù)，作為三維層流邊界層方程求解的邊界條件；2) 三維層流邊界層方程求解模塊。因直接從RANS方程解提取的層流邊界層信息難以滿足線性穩(wěn)定性分析的精度要求，此模塊通過非正交貼體坐標系對機翼邊界層方程進行求解，可在對網(wǎng)格量沒有過高要求的條件下，得到滿足三維邊界層穩(wěn)定性分析所需的高精度速度型u、溫度型T以及它們在物面法向的一階導(dǎo)數(shù)u′、T′、二階導(dǎo)數(shù)u″、T″；3) 基于線性穩(wěn)定性理論的雙eN轉(zhuǎn)捩預(yù)測模塊。通過對邊界層內(nèi)部的速度型、溫度型進行線性穩(wěn)定性分析，使用雙eN方法得到邊界層轉(zhuǎn)捩的位置xtr，并反饋給三維RANS求解模塊。重復(fù)迭代至轉(zhuǎn)捩位置收斂，依據(jù)收斂的轉(zhuǎn)捩信息繼續(xù)RANS方程的迭代求解，直至最終達到流場收斂標準。求解器中不同模塊耦合流程如圖2所示。

圖2 RANS求解器耦合轉(zhuǎn)捩預(yù)測模型的流程Fig.2 Coupling structure of RANS solverand transition prediction model

1.2 固定波角方法

對于TS波不穩(wěn)定(TSI)擾動，原有的包絡(luò)線方法[47]始終尋找擾動放大率最大的波角方向進行積分，如果存在較強的CF波，可能會錯誤的把CF波擾動積分到NTS中。本文采用固定波角方法計算TS波擾動累積放大因子，即

(3)

式中：NTS為TS波擾動累積放大因子；f為TS波頻率；x0為中性曲線計算出的小擾動開始放大的x位置。本文采用的固定波角方法，選擇擾動波波角ψ= 0° 方向即邊界層外邊界速度方向，計算不同頻率f擾動的放大率并進行積分，因此能夠保證始終對TS波擾動進行積分，不會混入CF波擾動。

通過后掠角10°和20°的無限翼展機翼邊界層轉(zhuǎn)捩預(yù)測算例，驗證固定波角方法對TS波不穩(wěn)定擾動主導(dǎo)轉(zhuǎn)捩預(yù)測的正確性。計算模型為以NACA642A015翼型為剖面翼型并垂直前緣布置的無限翼展機翼，計算狀態(tài)和轉(zhuǎn)捩位置所參考的實驗數(shù)據(jù)來自美國NASA Ams研究中心在其12英尺低湍流風洞中開展的相關(guān)實驗[48]。計算采用結(jié)構(gòu)化C-H網(wǎng)格，網(wǎng)格量為353×105×9。計算馬赫數(shù)為0.27，雷諾數(shù)范圍從5×106到25×106，迎角分別固定為0.5°和1°，TS波擾動放大因子轉(zhuǎn)捩閥值取為10.5。圖3～圖6分別給出了后掠角分別為10°和20°時機翼上表面轉(zhuǎn)捩位置計算值與實驗值的比較，以及計算壓力分布和TS波擾動累積放大因子曲線?？梢?，采用固定波角方法計算的轉(zhuǎn)捩位置與實驗吻合較好，且從擾動累積放大因子曲線可以看出，當機翼上表面出現(xiàn)逆壓梯度時，TS波劇烈增長，最終其擾動放大因子達到預(yù)設(shè)的轉(zhuǎn)捩放大因子10.5，判定轉(zhuǎn)捩發(fā)生。通過以上算例可以證明本文將固定波角方法用于TS波誘導(dǎo)邊界層流動轉(zhuǎn)捩的預(yù)測是十分可靠的。

圖3 10°后掠無限翼展NACA642A015機翼上表面計算轉(zhuǎn)捩位置與實驗值比較Fig.3 Calculated transition locations on theupper surface of infinite swept wing(NACA642A015, λ=10°) compared with experiments

圖4 20°后掠無限翼展NACA642A015機翼上表面計算轉(zhuǎn)捩位置與實驗值比較Fig.4 Calculated transition locations on theupper surface of infinite swept wing(NACA642A015,λ=20°) compared with experiments

圖5 10°后掠無限翼展NACA642A015機翼上表面計算壓力分布與TS波擾動累積放大曲線Fig.5 Pressure coefficients and N factors ofTSI on the upper surface of infinite swept wing(NACA642A015,λ=10°)

1.3 固定βr/固定f方法

對于CF波不穩(wěn)定(CFI)擾動誘導(dǎo)的轉(zhuǎn)捩，本文采用固定βr/固定f方法[49]計算CF波擾動累積放大因子，即

(4)

圖6 20°后掠無限翼展NACA642A015機翼上表面計算壓力分布與TS波擾動累積放大曲線Fig.6 Pressure coefficients and N factors of TSI on the uppersurface of infinite swept wing(NACA642A015,λ=20°)

式中：NCF為CF波擾動累積放大因子；βr即橫流展向波數(shù)；f為橫流行波頻率，f= 0Hz 時則為橫流駐波；x0為中性曲線計算出的小擾動開始放大的x位置；αi為x方向的擾動放大率。針對本文中短程民機巡航狀態(tài)，選擇對不同展向波數(shù)βr的橫流駐波進行擾動放大率的積分，計算橫流擾動累積放大因子NCF。對于CF波不穩(wěn)定擾動，基于平行流假設(shè)的線性穩(wěn)定性理論沒有考慮擾動幅值本身的變化，存在一定偏差。根據(jù)黃章峰等[35]的研究，線性穩(wěn)定性方程計算得到的擾動放大增長因子與考慮擾動幅值本身變化的拋物化穩(wěn)定性方程相比，形態(tài)趨勢一致但相對偏小。參考此結(jié)論，本文采用基于線性穩(wěn)定性方程預(yù)測CF波不穩(wěn)定誘導(dǎo)的轉(zhuǎn)捩時，采用較小的轉(zhuǎn)捩閥值。

通過后掠角40°和50°的無限翼展機翼邊界層轉(zhuǎn)捩預(yù)測算例，驗證固定βr/固定f方法對CF波不穩(wěn)定擾動主導(dǎo)的轉(zhuǎn)捩預(yù)測的正確性。機翼剖面仍為NACA642A015翼型，計算狀態(tài)和轉(zhuǎn)捩位置實驗值也來自參考文獻[48]。計算采用結(jié)構(gòu)化C-H網(wǎng)格，網(wǎng)格量為353×105×9，計算馬赫數(shù)為0.27，雷諾數(shù)范圍從4.3×106到7.3×106，迎角為-1.5°，CF波擾動放大因子轉(zhuǎn)捩閥值取為7.0。圖7和圖8給出了計算轉(zhuǎn)捩位置和實驗值的比較，吻合較好。本文在CF波擾動計算時僅考慮CF駐波，因此CF波頻率固定為0，分別計算不同波數(shù)的擾動累積放大因子，圖9和圖10分別給出了40°和50°后掠角計算模型在雷諾數(shù)6.1×106時上表面壓力分布和CF波擾動累積放大因子曲線，50°后掠角時計算出的CF波擾動增長更快，相對于40°后掠模型，轉(zhuǎn)捩更加靠前。通過以上算例證明將固定βr/固定f方法用于CF波主導(dǎo)的邊界層流動轉(zhuǎn)捩預(yù)測是可靠的。

圖7 40°后掠無限翼展NACA642A015機翼上表面計算轉(zhuǎn)捩位置與實驗值比較Fig.7 Calculated transition locations on theupper surface of infinite swept wing

圖8 50°后掠無限翼展NACA642A015機翼上表面計算轉(zhuǎn)捩位置與實驗值比較Fig.8 Calculated transition locations on theupper surface of infinite swept wing

圖9 40°后掠無限翼展NACA642A015機翼上表面計算壓力分布與CF波擾動累積放大曲線Fig.9 Pressure coefficients and N factors of CFI on the uppersurface of infinite swept wing(NACA642A015, λ=40°)

圖10 50°后掠無限翼展NACA642A015機翼上表面計算壓力分布與CF波擾動累積放大曲線Fig.10 Pressure coefficients and N factors of CFI on the uppersurface of infinite swept wing(NACA642A015,λ=50°)

1.4 雙eN轉(zhuǎn)捩預(yù)測方法

如前所述，基于線性穩(wěn)定性理論的eN方法是一種半經(jīng)驗的轉(zhuǎn)捩判斷方法，需要與實驗轉(zhuǎn)捩位置對比標定N值。研究表明，機翼后掠角不超過20°時，轉(zhuǎn)捩一般為TS波不穩(wěn)定(TSI)擾動誘導(dǎo)。當機翼后掠角增大時，機翼展向流動增強，產(chǎn)生CF波不穩(wěn)定(CFI)誘導(dǎo)的橫流轉(zhuǎn)捩。一般機翼后掠角大于30°時，邊界層發(fā)生橫流轉(zhuǎn)捩。但是對于實際問題，并不能事先預(yù)知是TSI誘導(dǎo)的轉(zhuǎn)捩，還是CFI誘導(dǎo)的轉(zhuǎn)捩。因此，本文采用既能判斷TS波不穩(wěn)定誘導(dǎo)轉(zhuǎn)捩，也能判斷CF波不穩(wěn)定誘導(dǎo)轉(zhuǎn)捩的雙eN方法。雙eN方法是一種轉(zhuǎn)捩預(yù)測方法的構(gòu)架，可以針對TS波和CF波選擇不同的擾動放大率積分策略構(gòu)建不同的轉(zhuǎn)捩預(yù)測方法。對于TS波的擾動放大增長因子NTS，采用固定波角方法計算；對于CF波的擾動放大增長因子NCF，采用固定βr/固定f方法計算。該方法可由公式(5)表示：

xtr=min(xTS,xCF)，

(5)

根據(jù)事先給定的轉(zhuǎn)捩閥值[(NTS)tr, (NCF)tr]，確定機翼邊界層的轉(zhuǎn)捩位置。如果NTS先到達其對應(yīng)的閥值(NTS)tr，就認為TS波在xTS處誘導(dǎo)轉(zhuǎn)捩發(fā)生；如果NCF先到達其對應(yīng)的閥值(NCF)tr，則認為邊界層流動轉(zhuǎn)捩由CF波在xCF處誘導(dǎo)轉(zhuǎn)捩發(fā)生。

本文對原雙eN方法[47]進行改進。原雙eN方法采用包絡(luò)線方法計算NTS，改進后的雙eN方法采用固定波角方法計算NTS，克服了原雙eN方法在存在較強的CF波時，會錯誤的把CF波擾動積分到NTS中的問題。下面給出一個后掠角為40°的無限翼展機翼算例，說明改進的有效性。計算模型為垂直前緣剖面是NACA642A015翼型的無限翼展機翼，后掠角為40°，馬赫數(shù)為0.27，雷諾數(shù)為7×106，迎角為-1.0°，計算狀態(tài)和轉(zhuǎn)捩位置實驗值也來自參考文獻[48]。圖11展示的是使用改進前后的雙eN方法計算得到的機翼上表面TS波(虛線)和CF波(實線)不穩(wěn)定性擾動放大積分因子發(fā)展規(guī)律。其中，圖11(a)使用原有的雙eN方法，判斷出的機翼上表面轉(zhuǎn)捩位置在流向19%處，而實驗轉(zhuǎn)捩位置為33%，原因是大后掠角時包絡(luò)線方法對NTS的計算是錯誤的。圖11(b)使用改進后的雙eN方法，判斷出的轉(zhuǎn)捩位置在流向37%處，與實驗結(jié)果吻合較好，證明了改進的有效性。

(a) 原有的雙eN方法

(b) 改進后的雙eN方法

采用本文方法，對DLR-F4翼身組合體外形機翼邊界層進行自由轉(zhuǎn)捩數(shù)值模擬，以驗證本文方法對復(fù)雜外形繞流轉(zhuǎn)捩判斷的正確性。對于翼身組合體等三維復(fù)雜外形，如果使用非結(jié)構(gòu)網(wǎng)格，雖沒有嚴格的網(wǎng)格拓撲結(jié)構(gòu)要求，但在計算擾動累積放大因子時，因其相鄰網(wǎng)格單元的無序性，需要采取專門的積分路徑確定方法[50]。本文針對復(fù)雜外形流動，建立了一種便于實現(xiàn)轉(zhuǎn)捩自動判斷的，且通用的翼身組合體結(jié)構(gòu)化網(wǎng)格拓撲結(jié)構(gòu)處理方法，對所有翼面類部件均適用，在計算擾動累積放大因子時，積分路徑可以直接使用機翼表面網(wǎng)格線，相對非結(jié)構(gòu)網(wǎng)格較為簡單便利。具體方法如下：包裹整個機翼(不包含翼梢端面和機翼后緣厚度端面)邊界層生成一塊完整的C型單塊網(wǎng)格，網(wǎng)格塊示意圖如圖12中線框所示。

對于這塊網(wǎng)格，執(zhí)行圖2中右側(cè)虛線框內(nèi)的轉(zhuǎn)捩自動判斷流程。因為機翼表面流向網(wǎng)格線完整連續(xù)，這樣在計算擾動累積放大因子時，積分路徑可以直接使用機翼表面網(wǎng)格線。如果需要對平尾、垂尾等翼面的邊界層進行轉(zhuǎn)捩自動判斷，上述方法同樣適用。而機身部分則不存在網(wǎng)格拓撲結(jié)構(gòu)的特殊要求，可按照常規(guī)方法自由生成網(wǎng)格。

圖12 DLR-F4翼身組合機外形機翼邊界層轉(zhuǎn)捩預(yù)測結(jié)構(gòu)網(wǎng)格拓撲結(jié)構(gòu)Fig.12 Topology of structured grid for boundary layertransition prediction of wing of DLR-F4wing-body configuration

計算狀態(tài)參考DLR在歐洲跨聲速風洞(European Transonic Wind Tunnel，ETW)使用溫敏漆(Temperature Sensitive Paint，TSP)技術(shù)進行的翼身組合體機翼邊界層轉(zhuǎn)捩測量試驗[51]：Ma= 0.785，Re= 6.0×106，α=-0.87°。計算網(wǎng)格為多塊結(jié)構(gòu)化網(wǎng)格，y+約為0.8，整個翼身組合體網(wǎng)格量為420萬，其中包裹機翼的網(wǎng)格塊機翼流向分布185個網(wǎng)格點，展向分布49個網(wǎng)格點，法向分布105個網(wǎng)格點。采用S-A湍流模型模擬湍流流動。放大因子轉(zhuǎn)捩閥值[(NTS)tr, (NCF)tr]根據(jù)經(jīng)驗取[10.5, 7.5]。

使用本文方法預(yù)測的機翼邊界層轉(zhuǎn)捩位置與試驗對比，如圖13所示，其中云圖為TSP試驗結(jié)果[51]，淺色部分為層流區(qū)域，深色部分為湍流區(qū)域，其中，從機翼前緣拖出的三角形湍流區(qū)域是因為機翼表面不光滑造成的湍流楔，并非小擾動引發(fā)的自然轉(zhuǎn)捩；黑色實線為本文方法計算得到的機翼上表面轉(zhuǎn)捩線，和試驗轉(zhuǎn)捩位置基本一致。Fey等[51]分析試驗結(jié)果TSP照片認為機翼內(nèi)段可能為橫流駐波不穩(wěn)定性誘發(fā)的轉(zhuǎn)捩。而本文計算結(jié)果同樣表明，翼根區(qū)域為橫流駐波不穩(wěn)定性擾動導(dǎo)致轉(zhuǎn)捩：圖14為展向25%處站位A和47%處站位B的機翼上表面TS波(虛線)和CF波(實線)不穩(wěn)定性擾動放大積分因子發(fā)展規(guī)律。由圖14(a)可見，展向站位A在流向x/c=0.34位置處，CF波擾動累積放大因子先于TS波達到其對應(yīng)閥值，即CF波不穩(wěn)定性在此站位誘發(fā)轉(zhuǎn)捩。而圖14(b)的展向站位B處為TS波不穩(wěn)定性誘發(fā)轉(zhuǎn)捩，轉(zhuǎn)捩位置在流向x/c=0.45位置處。

此算例表明，本文發(fā)展的雙eN方法不僅能夠?qū)C翼邊界層轉(zhuǎn)捩位置進行較為準確的預(yù)測，同時能夠判斷出機翼內(nèi)段轉(zhuǎn)捩是由橫流不穩(wěn)定性主導(dǎo)的機理，與TSP試驗結(jié)果分析[51]研究結(jié)論一致。驗證了本文方法的正確性，說明本文發(fā)展的方法，針對翼身組合體復(fù)雜外形，能夠?qū)崿F(xiàn)對機翼邊界層的轉(zhuǎn)捩自動判斷。

圖13 本文方法對DLR-F4機翼邊界層轉(zhuǎn)捩預(yù)測結(jié)果與TSP試驗照片對比Fig.13 Transition line on upper surface of DLR-F4wing predicted by proposed method and TSP image

(a) 展向站位A (2Z/B=25%)

(b) 展向站位B (2Z/B=47%)

2 基于流場當?shù)刈兞康摩?Reθ轉(zhuǎn)捩預(yù)測模型

經(jīng)典γ-Reθ轉(zhuǎn)捩模型是Menter和Lantry等[7]于2004年提出的基于流場當?shù)刈兞亢徒?jīng)驗關(guān)系式的轉(zhuǎn)捩預(yù)測模型，包含兩個輸運方程，如下：

(7)

針對經(jīng)典γ-Reθ模型無法預(yù)測橫流轉(zhuǎn)捩的問題，本文做出以下改進：在經(jīng)典γ-Reθ模型的基礎(chǔ)上引入橫流C1準則，發(fā)展了一種可以預(yù)測流向和橫流轉(zhuǎn)捩的γ-Reθ-C1三維轉(zhuǎn)捩模型。

2.1 基于當?shù)谻1準則的橫流轉(zhuǎn)捩預(yù)測模型

橫流C1準則是由法國宇航院的Arnal等[52]基于大量的實驗數(shù)據(jù)，歸納總結(jié)提出的關(guān)于臨界橫流位移厚度雷諾數(shù)Reδ2t和邊界層形狀因子H12的經(jīng)驗準則。C1準則可以表示為：

(8)

當橫流位移厚度雷諾數(shù)Reδ2>Reδ2t時，認為轉(zhuǎn)捩發(fā)生。Reδ2的表達式為：

(9)

其中δ2是邊界層橫流位移厚度，Qe是三維邊界層外邊界的勢流速度，w是橫流速度。

(10)

將當?shù)谻1準則作為橫流轉(zhuǎn)捩判據(jù)加入γ-Reθ模型間歇因子輸送方程的源項，形成γ-Reθ-C1三維轉(zhuǎn)捩模型。間歇因子輸送方程的源項修改為：

Pγ_3D=Flengthca1ρS(γFonset_3D)0.5(1-ce1γ)

(11)

其中Fonset_3D的表達式為：

(12)

圖15說明了γ-Reθ-C1轉(zhuǎn)捩模型和RANS方程求解的耦合過程。

2.2 算例驗證

選取了鐮刀形機翼算例來驗證耦合橫流轉(zhuǎn)捩準則的γ-Reθ-C1轉(zhuǎn)捩模型的正確性。鐮刀形機翼[54]是專門設(shè)計用來驗證橫流轉(zhuǎn)捩預(yù)測方法的一副機翼，分為不同后掠角的幾段，從翼根到翼梢的前緣后掠角分別是0°、30°、45°和55°。由于在每一段的前緣后掠角保持不變，鐮刀形機翼很適合對基于Falkner-Skan-Cooke相似性方程解的當?shù)谻1準則進行驗證。不僅如此，在機翼各分段之間的拐折處的橫流效應(yīng)明顯，會導(dǎo)致由橫流不穩(wěn)定造成的前緣轉(zhuǎn)捩，對橫流轉(zhuǎn)捩預(yù)測方法的適用性要求很高。實驗在布倫瑞克德國宇航院的低湍流度閉環(huán)風洞中進行[54]，風洞湍流度為0.17%，自由來流風速55 m/s，單位弦長雷諾數(shù)為2.75×106。機翼無后掠部分的弦長為0.8m，迎角-2.6°，平均表面粗糙度為9.78μm。

圖15 RANS求解器耦合γ-Reθ-C1轉(zhuǎn)捩模型的流程Fig.15 Coupling structure of RANS solver and γ-Reθ-C1 transition prediction model

采用γ-Reθ-C1模型對鐮刀機翼上表面的橫流轉(zhuǎn)捩進行模擬，計算狀態(tài)為Ma=0.162,Re=2.2×106,α=-2.6°,Tu=0.3%,RT=3。表面網(wǎng)格如圖16所示，計算網(wǎng)格的網(wǎng)格量為400萬，壁面第一層網(wǎng)格高度5×10-6，保證y+小于1。

圖16 鐮刀機翼幾何模型和表面網(wǎng)格Fig.16 Sickle wing geometric model and surface grid

計算得到的表面摩擦系數(shù)云圖如圖17所示，可以看出，相比于不能預(yù)測橫流轉(zhuǎn)捩的經(jīng)典γ-Reθ模型，γ-Reθ-C1模型能夠捕捉到后掠段的橫流轉(zhuǎn)捩，也能夠捕捉到拐折處由于橫流不穩(wěn)定導(dǎo)致的前緣轉(zhuǎn)捩尖角，說明γ-Reθ-C1模型在進行復(fù)雜機翼外形的橫流轉(zhuǎn)捩預(yù)測時適用性較強。但γ-Reθ-C1模型對于后掠角30°分段的轉(zhuǎn)捩位置預(yù)測仍不夠精確，預(yù)測的轉(zhuǎn)捩位置和實驗轉(zhuǎn)捩位置相比稍微靠前，需要進一步的研究和改進?？傮w來說，γ-Reθ-C1模型對機翼外形的橫流流動轉(zhuǎn)捩能夠進行較好的模擬。

圖17 γ-Reθ-C1模型計算得到的鐮刀機翼上表面的表面摩擦系數(shù)云圖Fig.17 Calculated surface friction coefficientcontour of the upper surface over sickle wing

3 基于動模態(tài)分解的DMD/eN轉(zhuǎn)捩預(yù)測新方法

3.1 翼型流場數(shù)據(jù)的動模態(tài)分解方法(DMD)

該部分對DMD方法的實現(xiàn)過程進行簡要介紹。DMD理論的詳細介紹可參見參考文獻[43]。

將從數(shù)值計算中獲得的流場數(shù)據(jù)表示成下面快照序列的形式：

(13)

假設(shè)整個抽樣空間中兩個相鄰的流場快照vi和vi+1可以由線性變化A聯(lián)系起來：

vi+1=Avi

(14)

因此，流動序列表示為如下矩陣形式：

(15)

(16)

將式(16)代入式(15)，并在式(16)兩端左乘UT和右乘以WΣ-1得到：

(17)

對DMD得到的A的特征值進行對數(shù)映射，得到空間模式下的特征值：

λ=lgμ/Δx=λre+λimi

(18)

其中，λre是指數(shù)放大率或衰減率，λim是空間波數(shù)，i表示虛數(shù)(i2=-1)。

3.2 DMD/eN轉(zhuǎn)捩預(yù)測方法

本文提出將DMD方法和eN轉(zhuǎn)捩方法結(jié)合的過程如下：

1) 以駐點為起點，沿弧長s對上、下表面的流場數(shù)據(jù)進行劃分。我們將翼型表面分成L段，如圖18所示。根據(jù)理論分析和經(jīng)驗，我們?nèi)=10。則(0～0.1)s/c為第1段，(0～0.2)s/c為第2段，以此類推(這里的c為無量綱化為1的翼型弦長)；

圖18 翼型表面分段示意圖Fig.18 Schematic of segments to perform DMD analysis

2) 沿翼型流向進行推進求解，對翼型表面每一段的流場數(shù)據(jù)進行DMD計算，得到該段的穩(wěn)定性特征值；

3) 確定有效特征值。本文采用兩種快照數(shù)對某段流場進行DMD分析，選擇不穩(wěn)定特征值(λi,λr)中重疊性最好的特征值分別作為各自流場快照數(shù)在本段的有效特征值，對應(yīng)的模態(tài)即是本段的特征模態(tài)。研究中發(fā)現(xiàn)數(shù)值噪聲對DMD分析的結(jié)果影響較大，在這里選取重疊性最好的特征值的目的是為了避免數(shù)值噪聲干擾。由于不同快照的特征值中重疊最好的兩個也存在些微小差異(并不是完全重合)，故我們以流場快照數(shù)多的結(jié)果為準；

4)N因子計算方法如公式(19)所示。即對DMD分析得到的空間放大率乘以對應(yīng)段的長度，得到翼型表面每段對應(yīng)的擾動放大因子Ni(下標i為段標號)；

Ni=si×λri

(19)

5) 找出Ni超過轉(zhuǎn)捩閥值Ncrit時對應(yīng)的段標號，記為s；連接該段及以前每段末的擾動放大因子N1、N2、……、Ns，得到一條類似包絡(luò)線的N值增長曲線，如圖19所示。在該曲線上線性插值得到轉(zhuǎn)捩閥值對應(yīng)的弧長位置s/c，再通過比較翼型表面網(wǎng)格點弧長與弦長的對應(yīng)關(guān)系，插值得到對應(yīng)的轉(zhuǎn)捩點弦長位置x/c。

圖19 N因子積分示意圖Fig.19 Diagram of N factor integration

以上為DMD/eN方法[8,41-42]的實現(xiàn)步驟，下面將該方法與RANS求解器進行耦合。圖20為繞翼型定常流動的轉(zhuǎn)捩預(yù)測過程。首先，采用RANS求解器對繞翼型的流動進行固定轉(zhuǎn)捩數(shù)值模擬，需選取足夠靠后的初始轉(zhuǎn)捩位置，保證流場充分收斂；其次，提取層流邊界層內(nèi)部的速度信息，形成用于DMD分析的流場快照；而后，采用DMD/eN方法進行流動轉(zhuǎn)捩預(yù)測；最后，將判斷的轉(zhuǎn)捩點位置回帶到RANS求解器進行迭代，直至得到收斂的結(jié)果。

圖20 N-S求解器與DMD/eN方法耦合轉(zhuǎn)捩預(yù)測示意圖Fig.20 Diagram of transition prediction using N-Ssolver coupling with DMD/eN method

3.3 算例驗證

開展不同迎角下NLF0416翼型的轉(zhuǎn)捩預(yù)測計算對該方法進行驗證。計算狀態(tài)是Ma=0.1,Re=4×106。轉(zhuǎn)捩后的湍流區(qū)域采用SA湍流模型。NLF0416翼型計算網(wǎng)格如圖21所示，計算網(wǎng)格量為641×241，第一層網(wǎng)格高度1×10-6c。下面以迎角α=1°為例，較為詳細地展示計算結(jié)果。

圖21 NLF0416翼型計算網(wǎng)格Fig.21 Computational grids of NFL0416 airfoil

采用總快照數(shù)100和200對翼型表面進行DMD分析，分別計算出上、下表面的N值增長。將推進求解中N值達到閥值及之前段的N值表示為隨x方向變化的曲線，并與LST/eN方法的N值增長曲線進行對比，如圖22和圖23所示。根據(jù)以往的經(jīng)驗，在利用eN方法進行NLF0416翼型流動的轉(zhuǎn)捩預(yù)測時，閥值N取為6較為合適。由圖可見，上表面 DMD/eN方法推進求解得到的N值增長曲線與LST/eN方法的包絡(luò)線吻合很好，而下表面DMD/eN方法的N值比LST/eN方法增長要早，但總體趨勢是可比的。

表1為DMD/eN方法與LST/eN方法計算得到的上、下表面轉(zhuǎn)捩點位置與實驗值的對比。實驗測出了翼型表面層流與湍流區(qū)的位置，轉(zhuǎn)捩發(fā)生在兩者之間。由表可知，上表面DMD/eN方法計算的轉(zhuǎn)捩點位置比LST/eN方法更接近實驗值，但兩者都比實驗值稍靠前。下表面兩者均在層流點到湍流點之間，與實驗值吻合良好。

圖22 NLF0416翼型上表面DMD/eN與LST/eN方法的N值增長對比圖(α=1°)Fig.22 Comparison of N growth between DMD/eN andLST/eN method on upper surface of NLF0416 airfoil

由于層流粘性系數(shù)遠小于湍流粘性系數(shù)，故這兩種流態(tài)下翼型表面的摩擦阻力差異很大。圖24為NLF0416翼型全湍流和自由轉(zhuǎn)捩狀態(tài)下的物面摩擦力系數(shù)分布圖。其中，黑色虛線是全湍流模擬的摩阻系數(shù)，藍色點畫線和紅色實線是考慮轉(zhuǎn)捩的分別采用LST/eN和 DMD/eN方法計算得出的結(jié)果?？梢?，LST/eN和DMD/eN兩種方法得到的摩阻系數(shù)分布吻合良好。翼型上、下表面對應(yīng)的摩阻系數(shù)分別在各自轉(zhuǎn)捩點附近有一個突增，對應(yīng)位置在0.33和0.5左右。轉(zhuǎn)捩點之前為層流，全湍流模擬的Cf遠大于自由轉(zhuǎn)捩的結(jié)果，而轉(zhuǎn)捩點之后進入湍流，兩者Cf的值開始逐步接近。

圖23 NLF0416翼型下表面DMD/eN與LST/eN方法的N值增長對比圖(α=1°)Fig.23 Comparison of N growth between DMD/eN andLST/eN method on lower surface of NLF0416 airfoil

MethodXtr_upXtr_lowDMD/eN0.3430.555LST/eN0.3260.558Exp.0.35～0.40.55～0.6

圖24 NLF0416翼型全湍流與自由轉(zhuǎn)捩表面摩擦力系數(shù)對比圖Fig.24 Comparison of skin friction (Cf) for NLF0416airfoil between full turbulence and free transition

以上迎角1°的算例較為詳細的展示了DMD/eN方法對NLF0416翼型的轉(zhuǎn)捩預(yù)測結(jié)果，驗證了本文方法的正確性。下面，將對-6°～6°中的多個迎角狀態(tài)的NLF0416翼型繞流進行轉(zhuǎn)捩點預(yù)測，并與實驗值、文獻[55]及LST/eN方法的計算結(jié)果進行對比。

圖25和圖26給出了NLF0416翼型上、下表面采用以上三種方法計算得到的轉(zhuǎn)捩點隨迎角的變化曲線。其中，黑色空心點和實心點分別為實驗測得的層流和湍流位置；藍色、紅色、綠色線分別代表文獻[55]、DMD/eN方法和LST/eN方法計算得到的轉(zhuǎn)捩點位置所連成的曲線；Xtr_up和Xtr_low分別代表翼型上、下表面的轉(zhuǎn)捩點位置。由圖25可見，隨著迎角的增大，上表面轉(zhuǎn)捩點逐步向翼型前緣移動。DMD/eN的計算結(jié)果基本位于層流和湍流點之間，且和文獻結(jié)果符合較好，而LST/eN方法計算的結(jié)果比實驗值略靠前。由圖26可知，不同迎角下三種方法的轉(zhuǎn)捩點計算結(jié)果均符合較好，且變化趨勢一致，進一步驗證了本文方法的正確性?？梢婋S著迎角的增大，下表面轉(zhuǎn)捩點向翼型后緣移動。這是由于隨著迎角的增加，下表面作為壓力面，順壓梯度增加，更易保持層流，故轉(zhuǎn)捩點后移。上表面的規(guī)律與之相反，但原理是類似的。

圖25 NLF0416翼型上表面轉(zhuǎn)捩點位置隨迎角的變化Fig.25 Variations of transition location along with theangles of attack on upper surface of NLF0416 airfoil

圖26 NLF0416翼型下表面轉(zhuǎn)捩點位置隨迎角的變化Fig.26 Variations of transition location along with theangles of attack on lower surface of NLF0416 airfoil

綜上，NLF0416翼型不同迎角下的轉(zhuǎn)捩預(yù)測結(jié)果與實驗結(jié)果吻合良好，驗證了本文發(fā)展DMD/eN轉(zhuǎn)捩預(yù)測方法的正確性。

4 轉(zhuǎn)捩預(yù)測方法在中短程民機機翼減阻設(shè)計中的應(yīng)用

為了證明本文所述轉(zhuǎn)捩預(yù)測方法的實用性，以雙eN轉(zhuǎn)捩預(yù)測方法為例，將該轉(zhuǎn)捩預(yù)測方法添加到RANS方程求解器中，并在此基礎(chǔ)上利用自主發(fā)展的優(yōu)化平臺SurroOpt[56]開展跨聲速自然層流機翼優(yōu)化設(shè)計研究工作。

目前，在現(xiàn)代民用客機跨聲速高雷諾數(shù)的典型飛行狀態(tài)下，為保持翼面大范圍層流流動，既要求機翼前緣后掠角盡量小于20°，又要求上翼面前部分具有適當順壓梯度。而在相同飛行速度下，以上兩方面設(shè)計均會導(dǎo)致翼面激波的增強，從而部分抵消增大層流范圍帶來的減阻效果。由此可見，應(yīng)用于現(xiàn)代客機的跨聲速自然層流機翼氣動設(shè)計難點在于如何在保持大范圍翼面層流面積的同時保持較小的激波強度。

經(jīng)過仔細考慮，本節(jié)首先配置了適用于現(xiàn)代中短程民機的自然層流機翼，該機翼的剖面翼型為順氣流布置的自然層流超臨界翼型NPU-LSC-72613，該翼型設(shè)計狀態(tài)為Ma=0.72，Re=2.0×107，CL=0.6，該狀態(tài)下翼型壓力分布如圖27。機翼幾何尺寸和設(shè)計狀態(tài)參考了空客A320相關(guān)數(shù)據(jù)，機翼面積123 m2，展長為35 m，展弦比為10.5，梢根比為0.3，前緣后掠角為19°，翼根和翼尖扭轉(zhuǎn)分別為0°和-2°，基于平均氣動弦長的無量綱機翼幾何形狀由圖28給出。初始機翼的設(shè)計巡航狀態(tài)為Ma=0.75，Re=2.0×107，CL=0.5。本節(jié)首先在設(shè)計狀態(tài)下對初始機翼進行考慮轉(zhuǎn)捩影響的繞流數(shù)值模擬計算，計算網(wǎng)格為C-H型結(jié)構(gòu)化網(wǎng)格，大小為257×105×73，其中機翼表面流向布置200個網(wǎng)格，展向布置48個網(wǎng)格，如圖29所示。擾動放大因子轉(zhuǎn)捩閥值取[(NTS)tr, (NCF)tr]=[9.0, 8.5]。

圖30(a)給出了設(shè)計狀態(tài)下初始機翼的壓力分布和上/下表面轉(zhuǎn)捩線?？梢钥闯?，機翼上表面雖然保持了大范圍的層流流動，但是產(chǎn)生了很強的激波；而機翼下表面靠近翼根處層流范圍急劇減小。

圖27 NPU-LSC-72613翼型在設(shè)計狀態(tài)下的壓力分布

圖28 適用于中短程民機的機翼平面形狀Fig.28 Wing geometry for the short- to mid-range aircraft

圖29 針對中短程民機的機翼C-H結(jié)構(gòu)化網(wǎng)格示意圖Fig.29 C-H structured grid for the wing geometryof the short- to mid-range aircraft

為了減小機翼阻力，提升氣動性能，在設(shè)計狀態(tài)Ma=0.75，Re=2.0×107，α=1.76°下開展機翼的優(yōu)化設(shè)計工作。公式(20)給出了具體的優(yōu)化目標和約束表達式，優(yōu)化過程中保持初始迎角不變，優(yōu)化目標是最小化機翼總阻力，約束為保持機翼的升力系數(shù)不小于原升力，力矩絕對值不增大，同時，為保持機翼內(nèi)部容積，翼根(標記為r)和35%展向位置(標記為k)處翼型相對厚度不減，翼尖位置(標記為t)翼型相對厚度不小于12%。

MinCD

W.r.t.x∈[-0.2,0.2]

S.t.CL≥CL0

|CM|≤|CM0|

(t/c)r,k≥(t/c)r0

(t/c)t≥12%

(20)

經(jīng)過400次CFD計算，優(yōu)化收斂。最終機翼阻力下降了11.3框(counts)，且在優(yōu)化機翼處，阻力系數(shù)的代理模型預(yù)測值與真實CFD計算值的相對誤差僅有0.05%，說明代理模型在最優(yōu)點附近空間內(nèi)已經(jīng)十分精確。圖31給出了優(yōu)化收斂歷程，可見優(yōu)化目標阻力系數(shù)為117.3框，而初始機翼在該狀態(tài)下阻力系數(shù)為128.6框，阻力降低了8.8%。表2和表3給出了優(yōu)化前后機翼力系數(shù)和控制剖面相對厚度的比較，可見相比于初始機翼，優(yōu)化機翼的阻力系數(shù)下降了11.3框，其中摩擦力系數(shù)下降1.6框，壓差阻力下降9.7框，而升力保持不變，力矩系數(shù)絕對值減小，機翼相對厚度約束嚴格滿足。

(a) 初始機翼

(b) 優(yōu)化機翼

圖31 適用于中短程民機的機翼優(yōu)化設(shè)計收斂曲線Fig.31 Convergence history of the optimization process

表2 優(yōu)化前后機翼氣動力系數(shù)比較Table 2 Comparisons of aerodynamic coefficientsbetween the baseline and the optimum wing

表3 優(yōu)化前后機翼不同展向位置處相對厚度比較Table 3 Comparisons of relative thicknesses betweenthe baseline and the optimum wing

圖30給出了優(yōu)化前后機翼上翼面壓力分布云圖以及上/下翼面轉(zhuǎn)捩線位置。在上翼面，初始機翼已經(jīng)有大范圍層流流動(占機翼上表面約62.5%)，而優(yōu)化機翼較好地保持了這一優(yōu)良特性，層流流動覆蓋57.8%，且很大程度上減弱了激波強度；在下翼面，初始機翼翼根部分出現(xiàn)了由三維邊界層內(nèi)橫流擾動誘導(dǎo)的前緣轉(zhuǎn)捩，而優(yōu)化機翼有效地抑制了CF波誘導(dǎo)的翼根位置提前轉(zhuǎn)捩，從而顯著推遲層流向湍流的自然轉(zhuǎn)捩，將機翼下翼面翼根位置層流范圍從初始機翼靠近前緣點處擴大到54.8%。圖32給出了優(yōu)化前后機翼典型展向站位截面壓力分布及轉(zhuǎn)捩位置比較?？梢姡跏紮C翼設(shè)計中為穩(wěn)定邊界層內(nèi)TS波，采用了較大的順壓梯度，導(dǎo)致了較高波前馬赫數(shù)，從而產(chǎn)生較強激波，所以盡管實現(xiàn)了翼面大范圍的層流流動，激波阻力的增加抵消了部分層流減阻的效果。優(yōu)化機翼則將上翼面順壓梯度保持在合適的強度，這樣既能穩(wěn)定邊界層內(nèi)TS波的發(fā)展，又不至于產(chǎn)生較強的激波，從而實現(xiàn)機翼總阻力的降低。而在展向15%站位處，可觀察到優(yōu)化機翼下表面在前緣加速區(qū)后保持一段變化緩慢的壓力區(qū)，實現(xiàn)了前緣轉(zhuǎn)捩的推遲。

以上不難看出，將雙eN轉(zhuǎn)捩預(yù)測方法與RANS方程求解器進行耦合，可以有效地預(yù)測TS波和CF波導(dǎo)致的轉(zhuǎn)捩，進而推動層流機翼設(shè)計工作的進行，同時證明了發(fā)展的轉(zhuǎn)捩預(yù)測方法有很好的工程實用價值。

圖32 初始和優(yōu)化機翼典型展向位置截面壓力分布對比Fig.32 Pressure coefficient comparisons of the baseline and the optimum wing at typical spanwise stations

5 結(jié) 論

本文對面向?qū)恿鳒p阻設(shè)計的轉(zhuǎn)捩預(yù)測方法進行了研究，得出如下結(jié)論：

1) 將求解線性穩(wěn)定性方程的雙eN轉(zhuǎn)捩預(yù)測方法推廣到了三維復(fù)雜外形。DLR-F4翼身組合體算例表明，本文發(fā)展的針對三維復(fù)雜外形的轉(zhuǎn)捩自動判斷方法可同時計及TS波和CF波擾動誘導(dǎo)的轉(zhuǎn)捩，較為準確的判斷轉(zhuǎn)捩位置，且能夠正確判斷轉(zhuǎn)捩誘發(fā)機制。

2) 實現(xiàn)了當?shù)谻1橫流轉(zhuǎn)捩準則和經(jīng)典γ-Reθ模型的耦合，形成了γ-Reθ-C1三維轉(zhuǎn)捩模型。利用該模型較為準確的預(yù)測了鐮刀形機翼表面的轉(zhuǎn)捩位置，表明其具有良好的工程應(yīng)用前景。

3) 改進的DMD/eN轉(zhuǎn)捩預(yù)測方法可用于翼型流動轉(zhuǎn)捩判斷，計算的轉(zhuǎn)捩位置與實驗值吻合良好，展現(xiàn)出較好的工程應(yīng)用價值，可作為層流減阻設(shè)計的有力工具。

4) 以雙eN方法的工程應(yīng)用為例，利用該轉(zhuǎn)捩預(yù)測方法，同時結(jié)合代理優(yōu)化方法，開展了針對中短程民機的跨聲速層流機翼優(yōu)化設(shè)計研究。最終結(jié)果表明，相對于基準機翼，優(yōu)化機翼減阻效果明顯，設(shè)計點氣動性能有很大提升，證明了本文發(fā)展的轉(zhuǎn)捩預(yù)測方法魯棒可靠，有很高的工程應(yīng)用價值。

參考文獻：

[1]Fu S, Wang L.Progress in turbulence/transition modelling[J].Advances in Mechacics, 2007, 37(3): 409-416.(in Chinese)符松, 王亮.湍流轉(zhuǎn)捩模式研究進展[J].力學進展, 2007, 37(3): 409-416.

[2]Zhou H.Studies on transition and turbulence[C]//2003 Advanced Research Papers on Aerodynamics.Beijing: Aerodynamic society of China, 2003: 87-93.(in Chinese)周恒.關(guān)于轉(zhuǎn)捩和湍流的研究[C]//2003空氣動力學前沿研究論文集.北京: 中國空氣動力學會, 2003: 87-93.

[3]Seitz A, Kruse M, Wunderlich T, et al.The DLR Project LamAiR: design of a NLF forward sweep wing for short and medium range transport application[C]//29th AIAA Applied Aerodynamics Conferences, AIAA-2011-3526, 2011.

[4]Urbin G, Knight D.Large-eddy simulation of a supersonic boundary layer using an unstructured grid[J].AIAA Journal , 2001, 39(7): 1288-1295.

[5]Schlatter S.Assessment of direct numerical simulation data of turbulent boundary layers[J].Journal of Fluid Mechanics, 2010, 659: 116-126.

[6]Krumbein A, Krimmelbein N, Grabe C.Streamline-based transition prediction techniques in an unstructured computational fluid dynamics code[J].AIAA Journal, 2017, 55(5): 1548-1564.

[7]Menter F R, Langtry R B, Likki S R, et al.A correlation-based transition model using local variables: part I-model formulation[J].Journal of Turbomachinery, 2004, 128(3): 413-422.

[8]Han z h, Wang S N, Han L, et al.Transition prediction based ondynamic mode decomposition for flows over airfoils[J].Acta Aeronautica et Astronautica Sinica, 2017, 38(1): 120034-1- 120034-17.(in Chinese)韓忠華, 王紹楠, 韓莉, 等.一種基于動模態(tài)分解的翼型轉(zhuǎn)捩預(yù)測新方法[J].航空學報, 2017, 38(1): 120034-1-120034-17.

[9]Smith A M O, Gamberoni N.Transition, pressure gradient and stability theory[M].Long Beach: Douglas Aircraft Company, 1956.

[10]Ingen J L V.A suggested semi-empirical method for the calculation of the boundary layer transition region.VTH-74[R].Delft: Delft University of Technology, 1956.

[11]Drela M.Newton solution of coupled viscous/inviscid multielement airfoil flows[C]//21st Fluid Dynamics, Plasma Dynamics and Lasers Conference.1990: 1470.

[12]Cebeci T, Stewartson K.On stability and transition in three-dimensional flows[J].AIAA Journal, 1980, 18(4): 398-405.

[13]Malik M R.COSAL: A black-box compressible stability analysis code for transition prediction in three-dimensional boundary layers.NASA-CR-165925[R].Washington D C: NASA, 1982[14]Mack L M.Stability of three-dimensional boundary layers on swept wings at transonic speeds[C]//IUTAM.Symposium Transsonicum III.Berlin Heidelberg: Springer, 1989: 209-223.

[15]Arnal D.Transition prediction in transonic flow[C]//IUTAM.Symposium Transsonicum III.Berlin Heidelberg: Springer, 1989: 253-262.

[16]Arnal D, Casalis G, Juillen J C.Experimental and theoretical analysis of natural transition on “infinite” swept wing[C]//IUTAM.Laminar-Turbulent Transition.Berlin Heidelberg: Springer, 1990: 311-325.

[17]Radespiel R, Graage K, Brodersen O.Transition predictions using Reynolds-averaged Navier-Stokes and linear stability analysis methods.AIAA-1991-1641[R].Reston: AIAA, 1991.

[18]Casalis G, Arnal D.ELFIN II subtask 2.3: Database method-Development and validation of the simplified method for pure cross-flow instability at low speed: ELFIN II-145[R].Toulouse: ONERA-CERT, 1996.

[19]Stock H W, Haase W.Feasibility study ofeNtransition prediction in Navier-Stokes methods for airfoils[J].AIAA Journal, 1999, 37(10): 1187-1196.

[20]Nebel C, Radespiel R, Wolf T.Transition prediction for 3D flows using a Reynolds-averaged Navier-Stokes code and N-factor methods.AIAA-2003-3593[R].Reston: AIAA, 2003.

[21]Krumbein A.Automatic transition prediction and application to high-lift multi-element configurations.AIAA-2004-2543[R].Reston: AIAA, 2004.

[22]Stock H W.Infinite swept wingNavier-Stokes computations with eNtransition prediction[J].AIAA Journal, 2005, 43(6): 1221-1229.

[23]Krumbein A.Automatic transition prediction and application to 3D wing configurations.AIAA-2006-0914[R].Reston: AIAA, 2006.

[24]Cliquetj, Houdeville R, Arnal D.Application of laminar-turbulent transition criteria in Navier-Stokes computations.AIAA-2007-515[R].Reston: AIAA, 2007.

[25]Krumbein A.eNtransition prediction for 3D wing configurations using database methods and a local, linear stability code[J].Aerospace Science and Technology, 2008, 12(8): 592-598.

[26]Perraud J, Arnal D, Casalis G, Archambaud J P, Donelli R.Automatic transition predictions using simplified methods[J].AIAA Journal, 2009, 47(11): 2676-2684.

[27]Perraudj, Cliquet J, Houdeville R, et al.Transport aircraft three-dimensional high-lift wing numerical transition prediction[J].Journal of Aircraft, 2008, 45(5): 1554-1563.

[28]Krumbein A.Automatic transition prediction and application to 3D high-lift configurations.AIAA-2006-3164[R].Reston: AIAA, 2006.

[29]Krimmelbein N, Radespiel R, Nebel C.Numerical aspects of transition prediction for three-dimensional configurations.AIAA-2005-4764[R].Reston: AIAA, 2005.

[30]Lee J D, Jameson A.Natural-laminar-flow airfoil and wing design byadjoint method and automatic transition prediction.AIAA-2009-897[R].Reston: AIAA, 2009.

[31]Zhang K, Song W P.Application of the full eNtransition method to the infinite swept-wing's transition prediction[J].Journal of Northwestern Polytechnical University, 2011, 29(1): 142-147.(in Chinese)張坤, 宋文萍.eN方法在無限展長后掠翼邊界層轉(zhuǎn)捩預(yù)測中的初步應(yīng)用[J].西北工業(yè)大學學報, 2011, 29(1): 142-147.

[32]Zuo S H, Yang Y, Li D.Investigation on cross-flow instabilities in swept-wing boundary layers with linear parabolized stability equations[J].Chinese Journal of Computational Physics, 2010, 27(5): 665-670.(in Chinese)左歲寒, 楊永, 李棟.基于線性拋物化穩(wěn)定性方程的后掠翼邊界層內(nèi)橫流穩(wěn)定性研究[J].計算物理, 2010, 27(5): 665-670.

[33]Sun P P, Huang Z F.Effect of sweep angle on stability and transition in a swept-wing boundary layer[J].Journal of Beijing University of Aeronautics and Astronautics, 2015, 41(7): 1313-1321.(in Chinese)孫朋朋, 黃章峰.后掠角對后掠機翼邊界層穩(wěn)定性及轉(zhuǎn)捩的影響[J].北京航空航天大學學報, 2015, 41(7): 1313-1321.

[34]Jing Z R, Sun P P, Huang Z F.Effect of attack angle on stability and transition in a swept-wing boundary layer[J].Journal of Beijing University of Aeronautics and Astronautics, 2015, 41(11): 2177-2183.(in Chinese)靖振榮, 孫朋朋, 黃章峰.小迎角對后掠機翼邊界層穩(wěn)定性及轉(zhuǎn)捩的影響[J].北京航空航天大學學報, 2015, 41(11): 2177-2183.

[35]Huang Z F, Lu X Z, Yu G T.Cross-flow instability analysis and transition prediction of airfoil boundary layer[J].Acta Aerodynamica Sinica, 2014, 32(1): 14-20.(in Chinese)黃章峰, 逯學志, 于高通.機翼邊界層的橫流穩(wěn)定性分析和轉(zhuǎn)捩預(yù)測[J].空氣動力學學報, 2014, 32(1): 14-20.

[36]Han Z H, Chen J, Zhu Z, et al.Aerodynamic design of transonic natural-laminar-flow (NLF) wing via surrogate-based global optimization.AIAA 2016-2041[R].Reston: AIAA, 2016.

[37]Grabe C, Krumbein A.Correlation-based transition transport modeling for three-dimensional aerodynamic configurations[J].Journal of Aircraft, 2013.

[38]Grabe C, Krumbein A.Extension of theγ-Reθtmodel for prediction of crossflow transition.AIAA 2014-1269[R].Reston: AIAA, 2014.

[39]Choi J H, Kwon O J.Enhancement of a correlation-based transition turbulence model for simulating crossflow instability[J].AIAA Journal, 2015.

[40]Xu J K, Bai J Q, Qiao L, et al.Transition model for predicting crossflow instabilities[J].Acta Aeronautica et Astronautica Sinica, 2015, 36(6): 1814-1822.(in Chinese)徐家寬, 白俊強, 喬磊, 等.橫流不穩(wěn)定性轉(zhuǎn)捩預(yù)測模型[J].航空學報, 2015, 36(6): 1814-1822.

[41]Wang S N, Han Z H, Song W P.A data-driven transition prediction method for flows over airfoils[J].Physics of Gases, 2017, 2(3): 5-16.(in Chinese)王紹楠, 韓忠華, 宋文萍.一種數(shù)據(jù)驅(qū)動的翼型流動轉(zhuǎn)捩預(yù)測方法[J].氣體物理, 2017, 2(3): 5-16.

[42]WU M M, Han Z H, Wang S N, et al.A DMD-based automatic transition prediction method for flows over airfoils.AIAA 2017-4303[R].Reston: AIAA, 2017.

[43]Schmid P J.Dynamic mode decomposition of numerical and experimental data[J].Journal of Fluid Mechanics, 2010, (656): 5-28.

[44]Rowley C W, Mezic I, Bagheri S, et al.Spectral analysis of nonlinear flows[J].Journal of Fluid Mechanics, 2009, (641): 115-127.

[45]Jovanovic M R, Schmid P J, Nichols J W.Sparsity-promoting dynamic mode decomposition.physics of fluids[J].2014, 26(2): 024103-1-024103-22.

[46]Seena A, Sung H J.Dynamic mode decomposition of turbulent cavity flows for self-sustained oscillations[J].International Journal of Heat and Fluid Flow, 2011, 32: 1098-1110.

[47]Zhang K, Song W P.Application of the full en transition method to the infinite swept-wing′s transition prediction[J].Journal of Northwestern Polytechnical University, 2011, 29(1): 142-147.(in Chinese)張坤, 宋文萍.eN法在無限展長后掠翼邊界層轉(zhuǎn)捩判斷中的初步應(yīng)用[J].西北工業(yè)大學學報, 2011, 29(1): 142-147.

[48]Boltz F W, Kenyon G, Allen C Q.Effects of sweep angle on the boundary-layer stability characteristics of an untapered wing at low speeds[R].NASA 19980227185 (1960).

[49]Zhang K, Song W P.VariedeNmethods used in the three dimensional boundary layer flows[C]//11th Russian-Chinese Conference on Fundamental problems of Aircraft Aerodynamics, Flight Dynamics, Strength and Flight Safety 2011.Zhukovsky: TsAGI, 2011.

[50]Krumbein A, Krimmelbein N, Grabe C.Streamline-based transition prediction techniques in an unstructured computational fluid dynamics code[J].AIAA Journal, 2017, 55(5): 1548-1564.

[51]Fey U, Egami Y, Engler R.High Reynolds number transition detection by means of temperature sensitive paint.AIAA 2006-0514[R].Reston: AIAA, 2006.

[52]Arnal D, Coustols E, Juillen J C.Experimental and theoretical study of transition phenomena on an infinite swept wing[J].La Recherche Aerospatiale(English Edition), 1984 (4): 39-54.

[53]Cooke J C, Howarth L.The boundary layer of a class of infinite yawed cylinders[J].Mathematical Proceedings of the Cambridge Philosophical Society, 2008, 46(4): 645-648..

[54]Petzold R, Radespiel R.Transition on a wing with spanwise varying crossflow evaluated with linear stability theory[C]//43rd AIAA Fluid Dynamics Conference.2013: 2466.

[55]Basha W.Accurate drag prediction for transitional external flow over airfoils[D].Concordia University, Canada, 2006: 56-85.

[56]Han Z H.SURROOPT: A generic surrogate-based optimization code for aerodynamic and multidisciplinary design.In: 30th congress of the international council of the aeronautical sciences.Korea, 2016.