高超聲速三維邊界層轉(zhuǎn)捩數(shù)值研究進展及預測軟件

黃章峰，張宇琦

(天津大學 機械工程學院，天津 300072)

0 引言

邊界層轉(zhuǎn)捩是指邊界層流動從層流狀態(tài)轉(zhuǎn)變?yōu)橥牧鳡顟B(tài)的一種自然現(xiàn)象，其物理機理是流動中擾動失穩(wěn)導致的。對于層流，一般存在流動不穩(wěn)定特征，因而在外界因素誘導下很容易產(chǎn)生邊界層擾動，之后擾動會經(jīng)歷線性和非線性增長，在多擾動相互作用產(chǎn)生的三維渦結(jié)構(gòu)破碎后，流動演變?yōu)橥牧鳡顟B(tài)。此時流動中摩擦阻力、熱交換、噪聲和摻混等相比層流狀態(tài)均大幅增加，例如湍流熱流與摩阻大概是層流的3～5 倍[1]。

在航空航天領(lǐng)域，設(shè)計師關(guān)心飛行器的氣動力熱布局，鑒于層流和湍流兩種流動狀態(tài)下熱流與摩阻的巨大差別，超聲速、高超聲速邊界層轉(zhuǎn)捩問題已經(jīng)成為飛行器設(shè)計需要考慮的一個關(guān)鍵因素[2]。李志文等[3]從飛行器總體設(shè)計角度梳理了高超聲速邊界層轉(zhuǎn)捩問題研究的必要性和重要性，總結(jié)了飛行試驗、風洞試驗、理論預示等轉(zhuǎn)捩研究手段的國內(nèi)外研究進展。相比國外，國內(nèi)在轉(zhuǎn)捩研究領(lǐng)域起步晚了近100 年，在高超聲速領(lǐng)域晚了30 多年[4-5]。然而經(jīng)過40 多年的努力，已經(jīng)基本與國際先進水平相當。國內(nèi)轉(zhuǎn)捩研究單位從21 世紀初的以天津大學周恒院士團隊為首的幾家高校，到2010 年以后覆蓋全國高校、研究所，實現(xiàn)了多點開花，特別是2015 年之后，國內(nèi)在高超聲速轉(zhuǎn)捩方面取得了豐碩的研究成果。2016 年，中國空氣動力研究與發(fā)展中心（簡稱氣動中心）聯(lián)合天津大學、中國科學院力學所、清華大學、國防科技大學、北京臨近空間飛行器系統(tǒng)工程研究所等單位共同承擔了科技部“大科學裝置前沿研究”重點研發(fā)計劃中的“高超聲速邊界層轉(zhuǎn)捩機理、預測及控制方法研究”項目。此后，國家自然科學基金資助的轉(zhuǎn)捩研究項目數(shù)也始終保持在兩位數(shù)（見圖1）[6]。

圖1 國家自然科學基金支持轉(zhuǎn)捩研究項目歷年數(shù)量和金額Fig.1 Number and amount of transition research projects supported by NSFC over the years

關(guān)于邊界層轉(zhuǎn)捩方面的研究進展，國內(nèi)多名專家已從不同角度作了回顧與介紹[1,7-19]。從這些綜述可知，轉(zhuǎn)捩研究對象主要集中在如平板、曲板、后掠翼、零攻角圓錐等二維邊界層，在線性增長和弱非線性相互作用兩個階段已獲得較全面的認知。然而實際飛行器表面邊界層以準三維和三維為主，即流向仍然緩慢變化，但展向流場不再均勻，于是不同的展向區(qū)域可能存在截然不同的轉(zhuǎn)捩形態(tài)，且擾動呈現(xiàn)三維和展向非周期特征，這給邊界層流動穩(wěn)定性分析和轉(zhuǎn)捩預測造成了巨大困難。隨著計算能力的增強和分析手段的豐富，三維邊界層的轉(zhuǎn)捩逐漸成為研究的熱點之一。本文將對三維邊界層轉(zhuǎn)捩機理、轉(zhuǎn)捩預測方法和轉(zhuǎn)捩預測軟件等方面的研究進展進行總結(jié)，希望為初學者提供有益借鑒。

1 三維邊界層轉(zhuǎn)捩機理研究進展

目前研究比較多的三維邊界層外形包括后掠翼、三角板、帶攻角圓錐、橢圓錐、BOLT（boundary-layer transition）以及氣動中心的升力體模型HyTRV（hypersonic transition research vehicle）等。研究結(jié)果表明，展向壓力梯度驅(qū)動的橫流、流線匯聚區(qū)的流向渦以及流線發(fā)散區(qū)的接觸線是三維邊界層中常見的流動形態(tài)，與之對應(yīng)的則是橫流不穩(wěn)定性、流向渦不穩(wěn)定性以及接觸線不穩(wěn)定性。

1.1 橫流不穩(wěn)定性及轉(zhuǎn)捩

由于壓力梯度的影響，邊界層內(nèi)存在一個與外部無黏流相垂直的速度分量，即橫流。橫流不穩(wěn)定性屬于無黏拐點失穩(wěn)，擾動通常比黏性不穩(wěn)定性增長更快，因此在三維邊界層中橫流失穩(wěn)機制是主導轉(zhuǎn)捩的機制之一。

牛海波等[20-21]對馬赫數(shù)6 的后掠平板和三角翼進行了風洞試驗，壁面高頻壓力脈動傳感器獲得了10 kHz 的非定常橫流波信號，并通過基于納米粒子示蹤平面激光散射（nano-tracer based planar laser scattering,NPLS）技術(shù)得到了非定常橫流渦結(jié)構(gòu)以及其二次失穩(wěn)結(jié)構(gòu)。趙磊[22]采用改進后的非線性拋物化穩(wěn)定性方程（nonlinear parabolized stability equations,NPSE）準確預測了后掠平板上定常橫流渦的非線性演化過程。Ma 等[23]采用NPSE 和直接數(shù)值模擬（direct numerical simulation,DNS）分析了定常橫流渦和第二模態(tài)的相互作用，并發(fā)現(xiàn)一對同頻率第二模態(tài)擾動相互作用產(chǎn)生定常渦，隨著向下游發(fā)展逐漸比初始定常橫流渦更不穩(wěn)定。Xu 等[24]采用Floquet 理論對后掠翼定常橫流渦進行了二次失穩(wěn)分析，獲得了位于渦頂端的y 模態(tài)和位于渦肩部的z 模態(tài)，發(fā)現(xiàn)了位于渦溝槽部分的新y 模態(tài)。新y 模態(tài)的增長率大于y 模態(tài)，且與z 模態(tài)增長率相當。為了進一步探究新y 模態(tài)在橫流轉(zhuǎn)捩過程中的作用，Chen 等[25]進行了直接數(shù)值模擬（DNS），發(fā)現(xiàn)高超橫流轉(zhuǎn)捩機制與低速情況相似，在二次失穩(wěn)階段由z 模態(tài)主導。Chen 等[26]研究了高焓邊界層橫流后掠翼穩(wěn)定性特征。

后掠板或翼型流場仍存在一個流場均勻的方向，使得采用恰當坐標系可使展向周期條件成立。而有攻角圓錐、橢錐等更復雜的三維外形，流場一般不存在一個展向均勻的方向，穩(wěn)定性分析更加困難。目前學界在首次失穩(wěn)方面主要采用一維穩(wěn)定性分析方法，該方法雖然計算便捷，但需要提前判斷擾動傳播路徑以及展向波數(shù)在該路徑上的變化[27]。全局穩(wěn)定性方法可避免上述困難，Paredes 等[28]和Lakebrink 等[29]應(yīng)用二維特征值穩(wěn)定性分析方法（BiGlobal）計算了橢錐橫流模態(tài)，發(fā)現(xiàn)其模態(tài)形狀函數(shù)分布在整個橫流區(qū)。Chen 等[30]運用約化BiGlobal 和面推進拋物化穩(wěn)定性方法（three-dimensional parabolized stability equations,PSE3D）研究了升力體標模HyTRV 上下表面不同橫流區(qū)的橫流失穩(wěn)特征，發(fā)現(xiàn)每個頻率存在數(shù)量眾多的失穩(wěn)模態(tài)，橫流模態(tài)最大增長頻率在15～30 kHz 之間。Liu 等[31]運用約化方法研究了風洞工況下 6°攻角圓錐橫流失穩(wěn)特征，發(fā)現(xiàn)橫流全局模態(tài)主要集中在背風面，且非定常橫流模態(tài)比定常模態(tài)更不穩(wěn)定。在橫流渦二次失穩(wěn)研究方面，Moyes 等[32]用NPSE 獲得定常橫流模態(tài)沿特定路徑的非線性演化，采用BiGlobal 計算飽和定常橫流模態(tài)的二次失穩(wěn)特征。Choudhari 等[33]通過DNS 計算粗糙元陣列引起的橫流擾動，并用BiGlobal 和PSE3D 研究橫流渦的二次失穩(wěn)擾動，發(fā)現(xiàn)Mack 模態(tài)在圓錐側(cè)面橫流渦影響下可經(jīng)歷最大可能增長。

1.2 流向渦不穩(wěn)定性及轉(zhuǎn)捩

早期針對流向渦問題的全局穩(wěn)定性分析主要集中在HIFiRE-5 模型（2∶1 橢圓錐）的中心線不穩(wěn)定性上。根據(jù)飛行實驗，Choudhari 等[34]最先分析了18 km飛行條件下HIFiRE-5 中心線流向渦不穩(wěn)定性，并發(fā)現(xiàn)了一個位于流向渦頂部的不穩(wěn)定模態(tài)。Paredes 等[28]分析了33 km 處的另一個飛行工況，發(fā)現(xiàn)主導頻率為130 kHz 的對稱模態(tài)。李曉虎等[35]針對風洞工況進行了研究，分析了不穩(wěn)定模態(tài)的無黏不穩(wěn)定性特征。Choudhari 等[36]使用PSE3D 對模態(tài)在流向上的空間演化進行了分析，發(fā)現(xiàn)了PSE3D 和DNS 的結(jié)果符合很好。

對于帶攻角圓錐的分析則是近年來才開展的。陳曦等[37]基于DNS 擾動時均場，對帶攻角小鈍錐的背風流向渦進行了全局穩(wěn)定性分析，根據(jù)基本流速度剖面的剪切層位置和模態(tài)的相速度，將不穩(wěn)定模態(tài)分為兩類，相速度較低的內(nèi)模態(tài)位于近壁剪切層，而相速度較高的外模態(tài)位于流向渦頂部的剪切層。BiGlobal 和PSE3D 分析的結(jié)果與來自DNS 的轉(zhuǎn)捩數(shù)據(jù)之間存在一致性，發(fā)現(xiàn)外模態(tài)不穩(wěn)定性占主導地位，轉(zhuǎn)捩N值約為10。Li 等[38]研究了相同的模型工況，但使用了無入口擾動的層流場，辨識了Mack 模態(tài)，與風洞試驗對比得到的轉(zhuǎn)捩N值約為5。Li 等[39]在背風中心線處布置了一個粗糙元，發(fā)現(xiàn)外模態(tài)增長被顯著延后而內(nèi)模態(tài)則被略微促進。Zhang 等[40]研究了5 mm 鈍度圓錐飛行工況下的背風面流向渦失穩(wěn)機理，發(fā)現(xiàn)內(nèi)模態(tài)具有聲輻射性質(zhì)，而上游的Mack 模態(tài)可在下游轉(zhuǎn)化為剪切模態(tài)并可能引起流向渦破碎。

1.3 附著線不穩(wěn)定性及轉(zhuǎn)捩

附著線附近流場是模型表面邊界層在展向最薄的區(qū)域，此處擾動失穩(wěn)頻率一般遠高于其他展向區(qū)域，對實驗測量和數(shù)值模擬造成巨大困難。Paredes等[28]研究了橢錐附著線失穩(wěn)特征，發(fā)現(xiàn)存在高頻Mack模態(tài)。Xi 等[41]研究了高超聲速來流條件下不同后掠角的后掠圓柱附著線失穩(wěn)特征，發(fā)現(xiàn)大后掠角情形下，在附著線附近呈現(xiàn)典型附著線失穩(wěn)特征，而在遠處則呈現(xiàn)Mack 模態(tài)特征。Chen 等[30]同樣在升力體肩部附著線區(qū)域發(fā)現(xiàn)了高頻Mack 模態(tài)，與前人研究不同的是，此處的附著線流動不再具有展向?qū)ΨQ性，原有的對稱-反對稱模態(tài)交替出現(xiàn)的結(jié)構(gòu)坍縮為一支不穩(wěn)定模態(tài)。他們還利用PSE3D 研究了Mack 模態(tài)沿下游的演化，發(fā)現(xiàn)其擾動局限在附著線附近，從而不太可能與附近橫流區(qū)的擾動發(fā)生相互作用。

圓錐迎風面中心線附近流動也可視為廣義的附著線流動。Paredes 等[42]率先對風洞工況下圓錐迎風面失穩(wěn)開展了全局穩(wěn)定性分析，發(fā)現(xiàn)高頻Mack 模態(tài)。Chen 等[43]在對飛行工況下的圓錐迎風面開展分析后，同樣發(fā)現(xiàn)主導不穩(wěn)定性為Mack 模態(tài)，他們進一步將Mack 模態(tài)分為D 和S 分支，其中D 模態(tài)類似于傳統(tǒng)的附著線模態(tài)，即呈現(xiàn)對稱-反對稱模態(tài)交替出現(xiàn)結(jié)構(gòu)，且對稱模態(tài)最不穩(wěn)定，D 模態(tài)集中在迎風中心線附近，而S 模態(tài)位于側(cè)面。他們還利用PSE3D研究了不同模態(tài)的演化，發(fā)現(xiàn)部分S 模態(tài)會在下游激發(fā)D 模態(tài)。楊鵬等[44]對圓錐迎風中心線轉(zhuǎn)捩開展了DNS 研究，證實了Chen 等[43]的結(jié)果，即轉(zhuǎn)捩最先在迎風中心線附近發(fā)生，且主導失穩(wěn)模態(tài)為Mack 模態(tài)。在非線性演化階段，他們發(fā)現(xiàn)了條帶結(jié)構(gòu)，推測基頻共振是主導非線性機制。

上述轉(zhuǎn)捩機理的研究有助于深刻理解轉(zhuǎn)捩過程，但與準確預測轉(zhuǎn)捩尚有距離。目前主流轉(zhuǎn)捩預測方法是轉(zhuǎn)捩模型和eN方法。

2 三維邊界層轉(zhuǎn)捩模型研究進展

轉(zhuǎn)捩模型是基于雷諾平均Navier-Stokes（Reynoldsaveraged Navier-Stokes,RANS）方程的一種?；椒?，該方法一般建立在現(xiàn)有成熟的湍流模型基礎(chǔ)之上。湍流建模主要是對RANS 方程的雷諾應(yīng)力項進行封閉，可分為雷諾應(yīng)力模型和渦黏模型。轉(zhuǎn)捩模型大都是在渦黏模型基礎(chǔ)上實現(xiàn)對轉(zhuǎn)捩現(xiàn)象相關(guān)統(tǒng)計量的?；?，伴隨著轉(zhuǎn)捩模式理論的發(fā)展不斷在推陳出新。在轉(zhuǎn)捩模型的構(gòu)造與發(fā)展中有兩個重要概念，分別是間歇因子與非湍流脈動。明晰這兩個重要概念的物理意義及其與轉(zhuǎn)捩建模的聯(lián)系，有助于更深入理解各類轉(zhuǎn)捩建模的物理背景與發(fā)展脈絡(luò)。

間歇因子是對實驗觀測到的轉(zhuǎn)捩間歇現(xiàn)象的描述。Emmons[45]經(jīng)過實驗發(fā)現(xiàn)轉(zhuǎn)捩過程是層流狀態(tài)下多個孤立湍斑的產(chǎn)生與發(fā)展的過程，而非普遍認為的層流流動的簡單破碎?；诖耍珽mmons 提出間歇因子的概念來描述轉(zhuǎn)捩過程，即在轉(zhuǎn)捩過程中轉(zhuǎn)捩區(qū)的空間點可以在一段時間內(nèi)是層流狀態(tài)、在另一段時間內(nèi)是湍流狀態(tài)。間歇因子γ表征了該點在統(tǒng)計意義下的流動狀態(tài)：γ=0 對應(yīng)層流，γ=1 對應(yīng)湍流，而轉(zhuǎn)捩介于兩者之間。Dhawan 等[46]給出了間歇因子γ的平板邊界層流向分布實驗擬合式。因此，間歇因子使轉(zhuǎn)捩模型不再局限于“開關(guān)式”的層流/湍流模擬。間歇因子求解主要有代數(shù)求解與輸運方程求解兩種方式。合理的間歇因子求解對于轉(zhuǎn)捩-湍流模型的耦合、轉(zhuǎn)捩區(qū)長度的模擬以及三維轉(zhuǎn)捩流場預測具有重要意義。γ-Reθt和k-ω-γ等轉(zhuǎn)捩建模均離不開間歇因子這一重要物理概念。

非湍流脈動的提出同樣是基于轉(zhuǎn)捩實驗觀測現(xiàn)象。在轉(zhuǎn)捩起始位置之前的層流區(qū)域，存在不同于湍流脈動的擾動形式，稱為層流脈動。Mayle 等[47]提出了非湍流脈動動能這一概念，用以描述上述區(qū)域?qū)恿髅}動的動力學行為。相比前期經(jīng)驗性建模，引入非湍流脈動概念使轉(zhuǎn)捩建模能夠引入更多模態(tài)信息等轉(zhuǎn)捩物理機制，建立更具有轉(zhuǎn)捩機理可解釋性的預測模型，例如k-kL-ε、k-kL-ω、k-ζ、k-ω-γ等轉(zhuǎn)捩模型。

目前見諸文獻的轉(zhuǎn)捩模型已不下20 種，其修正改進模型更不勝枚舉。本文重點選取了早期的低雷諾數(shù)湍流模型以及在學術(shù)界、工業(yè)界應(yīng)用研究相對廣泛的γ-Reθt模型、k-ω-γ模型和層流動能轉(zhuǎn)捩模型4 類典型轉(zhuǎn)捩模型，就轉(zhuǎn)捩模型構(gòu)造、模型適用性研究及模型應(yīng)用情況進行介紹。

2.1 低雷諾數(shù)湍流模型

早期的邊界層轉(zhuǎn)捩模型預測的嘗試，并不是建立相對獨立的轉(zhuǎn)捩輸運方程，也并未在湍流模型中引入類似間歇因子等轉(zhuǎn)捩相關(guān)物理量，而是直接修正低雷諾數(shù)湍流模型來預測轉(zhuǎn)捩流動。主要辦法是通過對低雷諾數(shù)湍流模型中的阻尼函數(shù)進行修正，以臨界雷諾數(shù)為轉(zhuǎn)捩判據(jù)，使渦黏系數(shù)在臨界位置突然增大，造成摩擦系數(shù)及邊界層厚度增長率“陡升”，從而實現(xiàn)轉(zhuǎn)捩現(xiàn)象的流場模擬。Wilcox[48]通過該辦法修正k-ω模型，實現(xiàn)了對槽道、管道湍流以及平板邊界層轉(zhuǎn)捩的有效預測。Schmidt[49]在k-ε模型方面亦有關(guān)于低雷諾數(shù)湍流模型的轉(zhuǎn)捩預測研究。國內(nèi)嚴明等[50]基于低雷諾數(shù)k-ε模型進行了改進，實現(xiàn)了對T3A 平板邊界層的轉(zhuǎn)捩預測。但是相關(guān)研究[51]認為，“低雷諾數(shù)湍流模型對轉(zhuǎn)捩現(xiàn)象的預測只是一種巧合，其阻尼函數(shù)的構(gòu)造準則是模擬黏性次層，而并未考慮轉(zhuǎn)捩的物理機理”。歐洲于20 世紀90 年代啟動的聯(lián)合研究項目“TransPerturb”對大量低雷諾數(shù)湍流模型進行了葉輪機械轉(zhuǎn)捩測試，發(fā)現(xiàn)模型存在的普遍問題是預測的轉(zhuǎn)捩位置過于靠前且轉(zhuǎn)捩區(qū)的長度過短。雖然該類模型未考慮足夠的轉(zhuǎn)捩機理，轉(zhuǎn)捩預測可靠性低，但其構(gòu)造簡單、可實現(xiàn)性強，因此在工程中還是可以基于通用計算流體力學（Computational Fluid Dynamics,CFD）的成熟湍流模型模塊進行修正而實現(xiàn)快速應(yīng)用。

2.2 γ-Reθt 轉(zhuǎn)捩模型

γ-Reθt轉(zhuǎn)捩模型是工程應(yīng)用最為廣泛的轉(zhuǎn)捩模型之一，由Menter 和Langtry 于2004 年[52-53]正式提出，并在2009 年[54]公布了模型的完整構(gòu)造形式和詳細參數(shù)設(shè)置。該模型的基本思想是在SST 湍流模型[55]的基礎(chǔ)上，構(gòu)造關(guān)于間歇因子（γ）和動量厚度雷諾數(shù)（Reθt）的兩個輸運方程。通過求解動量厚度雷諾數(shù)，結(jié)合風洞實驗數(shù)據(jù)建立的經(jīng)驗關(guān)系，來構(gòu)建用以判斷轉(zhuǎn)捩發(fā)生的參數(shù)，即當?shù)鼗瘻u雷諾數(shù)與臨界動量厚度雷諾數(shù)的比值。間歇因子，如前面所述，是對轉(zhuǎn)捩間歇現(xiàn)象的描述，一定程度上可以考慮轉(zhuǎn)捩的物理過程與轉(zhuǎn)捩區(qū)的存在?；陂g歇因子和動量厚度雷諾數(shù)，構(gòu)建適用于現(xiàn)代CFD 的完全當?shù)鼗D(zhuǎn)捩模型。依照該建模思想，基于間歇因子的轉(zhuǎn)捩建模相比于早期的“開關(guān)式”的轉(zhuǎn)捩建模，在層流-轉(zhuǎn)捩-湍流耦合計算以及轉(zhuǎn)捩區(qū)長度預測方面具有重要的進步意義。采用轉(zhuǎn)捩動量厚度雷諾數(shù)輸運方程替代傳統(tǒng)的非當?shù)刈兞縿恿亢穸壤字Z數(shù)來判斷轉(zhuǎn)捩發(fā)生，便于現(xiàn)代CFD 大規(guī)模并行計算，也是γ-Reθt模型的一大特點[56]。

γ-Reθt轉(zhuǎn)捩模型被提出時，還存在很多功能與性能上的不足，需要對模型進行完善與改進以滿足工程需求。Langtry 等[57]通過在動量厚度雷諾數(shù)輸運方程中增加橫流項，以流向渦量強度近似判斷橫流渦強度，完成了模型對橫流功能的拓展，實現(xiàn)了模型對三維邊界層橫流轉(zhuǎn)捩的預測。Menter 等[58]針對伽利略不變性及計算效率低的問題，簡化了模型的構(gòu)造方式，建立了間歇因子單輸運方程轉(zhuǎn)捩模型。此外，模型改進還包括外流改進[59]、分離流動修正[60]、與DES 類方法耦合[61]、可壓縮修正[62]、橫流效應(yīng)[63]、與神經(jīng)網(wǎng)絡(luò)耦合[64]等其他方面。γ-Reθt模型具有完全當?shù)鼗?、可實現(xiàn)性強、易拓展的特點，早期應(yīng)用于對低速平板、二維翼型、渦輪葉片以及典型三維外形的自然轉(zhuǎn)捩、bypass 轉(zhuǎn)捩、分離誘導轉(zhuǎn)捩及橫流轉(zhuǎn)捩的預測。

由于構(gòu)建γ-Reθt模型所依賴的重要經(jīng)驗關(guān)系式是基于低速的風洞試驗與數(shù)值試驗建立的，因此對高超聲速轉(zhuǎn)捩的預測偏差較大，需要對模型進行高超聲速修正。修正方法主要是對模型中決定轉(zhuǎn)捩位置、轉(zhuǎn)捩區(qū)長度的重要經(jīng)驗關(guān)系式進行修改或重新定義，使模型具備高超聲速轉(zhuǎn)捩預測能力。Krause 等[65]將臨界動量厚度雷諾數(shù)和調(diào)整轉(zhuǎn)捩區(qū)長度尺度的函數(shù)直接寫成來流湍流度的函數(shù)，并在雙楔外形上得到了驗證。張曉東等[66]針對高超聲速條件下轉(zhuǎn)捩動量厚度雷諾數(shù)顯著增大的情況，基于高超聲速風洞試驗數(shù)據(jù)，提出了新的轉(zhuǎn)捩動量厚度雷諾數(shù)經(jīng)驗關(guān)系式，同樣在雙楔上得到了驗證。You 等[67]評估了原始γ-Reθt模型對高超聲速邊界層轉(zhuǎn)捩的預測性能，并在模型中分別引入Krause 修正、激波探測器以及壓力梯度，同時評估了它們的影響，完成了雙楔轉(zhuǎn)捩預測的驗證。針對后期改進的間歇因子單輸運方程轉(zhuǎn)捩模型，Aliaga 等[68]將其成功應(yīng)用于平板、旋成體等的高超聲速轉(zhuǎn)捩預測，但并未論述模型及修正細節(jié)。劉清揚等[69]基于γ-Reθt轉(zhuǎn)捩模型，構(gòu)建了考慮可壓縮性的比擬關(guān)系（fRe）輸運方程，發(fā)展了γ-Reθt-fRe轉(zhuǎn)捩模型。陳堅強研究團隊通過對γ-Reθt模型展開系統(tǒng)研究，最大程度發(fā)揮其優(yōu)勢及潛力。他們基于γ-Reθt模型，提出了穩(wěn)定性理論、實驗和DNS 數(shù)據(jù)與轉(zhuǎn)捩準則多層級融合的建模思想，形成了涵蓋馬赫數(shù)等參數(shù)影響的高超轉(zhuǎn)捩數(shù)據(jù)集，構(gòu)造了考慮可壓縮、橫流、粗糙度等因素的高超轉(zhuǎn)捩準則，建立了C-γ-Reθt高超轉(zhuǎn)捩模型[70-74]，在三維復雜高超飛行器（如升力體標模HyTRV[75]）的轉(zhuǎn)捩預測中取得了非常好的效果（見圖2[76]），是目前最有希望應(yīng)用于復雜工程外形的轉(zhuǎn)捩模型。

圖2 HyTRV 飛行器標模（左）和迎風面轉(zhuǎn)捩模型預測的轉(zhuǎn)捩陣面（右）[76]Fig.2 Aircraft standard model of HyTRIV (left) and transition front predicted by transition modeling at windward (right)[76]

2.3 k-ω-γ 轉(zhuǎn)捩模型

Suzen 等[77]建立了一種比較通用的間歇因子輸運方程，方程中包含生成項、耗散項。他們把該間歇因子輸運方程與Menter 的k-ω湍流模型相結(jié)合，在預測低壓渦輪中取得了一定成功。為了便于CFD 計算，需要采用當?shù)刈兞繉δＰ偷氖Х€(wěn)模態(tài)進行?；Ｒ疝D(zhuǎn)捩的第一模態(tài)擾動波和第二模態(tài)擾動波等都有不同的特征時間 τnt（或特征尺度lnt）。Warren 等[78-79]基于特征長度構(gòu)建了層流脈動渦黏系數(shù)模型，并采用當?shù)鼗膮?shù)表示。符松和王亮[80-81]發(fā)展了Warren等方法，把模型推廣至高超聲速轉(zhuǎn)捩情況，命名為kω-γ轉(zhuǎn)捩模型，并引入橫流擾動特征時間，使k-ω-γ轉(zhuǎn)捩模型能在名義上體現(xiàn)第一、第二模態(tài)及橫流失穩(wěn)的影響，初步實現(xiàn)了高超聲速多模態(tài)轉(zhuǎn)捩預測（見圖3）。

圖3 轉(zhuǎn)捩模型預測的高超聲速有攻角圓錐轉(zhuǎn)捩陣面[81]Fig.3 Transition front of hypersonic cone with angle of attack predicted by transition modeling[81]

目前關(guān)于該模型部署于商業(yè)軟件的報道尚未見公開文獻，但該模型仍在不斷發(fā)展以提高其工程適用性。郝子輝等[82]采用k-ω-γ模型在高超聲速圓錐邊界層轉(zhuǎn)捩預測中研究了模型對于頭部鈍度、來流雷諾數(shù)、來流湍流度等的影響規(guī)律。Zhou 等[83]對k-ω-γ模型進行了橫流拓展，構(gòu)造了基于橫流速度和橫流雷諾數(shù)的時間尺度項，拓展了模型的高超聲速橫流轉(zhuǎn)捩預測能力，并應(yīng)用于高超聲速復雜構(gòu)型[84]。Wang 等[85]對k-ω-γ模型的Mack 模態(tài)時間尺度進行了修正，引入了當?shù)鼗拟g度影響，提升了模型在高超聲速圓錐邊界層轉(zhuǎn)捩預測中的性能。Zhao 等[86]開展了k-ωγ模型在平板與直錐算例中的參數(shù)不確定度量化分析。Yang 等[87-88]將粗糙度放大因子（Ar）輸運方程引入模型中，建立了考慮粗糙度影響的k-ω-γ-Ar四方程模型。Zhang 等[89]提出基于流場反演和機器學習的數(shù)據(jù)驅(qū)動框架來提高四方程模型對高超聲速邊界層轉(zhuǎn)捩的預測能力。劉宏康等[90]初步開展了k-ω-γ模型在高超聲速充氣式柔性減速器轉(zhuǎn)捩預測的應(yīng)用研究。

2.4 層流動能轉(zhuǎn)捩模型

層流動能轉(zhuǎn)捩模型一般通過建立層流脈動動能（kL）輸運方程，并與湍流模型輸運方程耦合求解的方式實現(xiàn)轉(zhuǎn)捩建模與轉(zhuǎn)捩預測。具有代表性的層流動能轉(zhuǎn)捩模型有k-kL-ε和k-kL-ω，用來描述非湍流脈動的動力學行為并預測邊界層轉(zhuǎn)捩。

k-kL-ε模型是通過層流脈動動能輸運方程與kε湍流模型進行耦合建立的，是首個考慮轉(zhuǎn)捩前擾動增長以及bypass 轉(zhuǎn)捩機制的完全當?shù)鼗瘑畏匠虒恿鲃幽苻D(zhuǎn)捩模型。該模型由Walters 等[91]提出，并部署于商業(yè)CFD 軟件Fluent 的早期版本中。他們針對邊界層轉(zhuǎn)捩之前的層流脈動沿流向增長的特點，基于層流動能一般性動力學特征的假設(shè)，建立了“湍流動能-層流動能-耗散率”的三個輸運方程，其中層流動能與湍流動能根據(jù)壁面距離與湍流長度的一個截斷尺度來進行區(qū)分，二者之和等于總脈動動能。模型通過對層流動能與湍流動能的建模避免了采用經(jīng)驗公式預測轉(zhuǎn)捩。為提升模型在轉(zhuǎn)捩區(qū)的預測能力，Walters 等[92]延續(xù)上述思想，進一步在k-ω湍流模型的基礎(chǔ)上，完成了k-kL-ω轉(zhuǎn)捩建模，其構(gòu)造形式與kkL-ε模型基本一致。該模型主要部署在較新版本的Fluent 中。

k-kL-ε和k-kL-ω模型在低速平板、渦輪葉片、機翼翼型的邊界層轉(zhuǎn)捩預測中取得了較好的效果，能夠較準確模擬反映轉(zhuǎn)捩位置的壁面摩阻分布、熱流分布等[91-92]。孫潤鵬等[93]將k-kL-ω模型應(yīng)用于NASAMark II 內(nèi)冷葉片轉(zhuǎn)捩預測。陳燦平等[94]在k-kL-ω模型中引入了流線曲率因子與間歇因子以改進其對bypass 轉(zhuǎn)捩預測的性能，在DCA 壓氣機葉珊、S809 翼型、RAE2822 翼型的轉(zhuǎn)捩模型中進行了測試。KO?í?EK 和FüRST 等[95]在開源CFD 平臺OpenFOAM 上實現(xiàn)了k-kL-ω模型，研究了模型在可壓縮流動中的表現(xiàn)性能，將模型應(yīng)用于VKI 和SE-1050 葉柵轉(zhuǎn)捩預測。Salimipour[96]對k-kL-ω模型的分離區(qū)模擬性能進行了改進與提升，將模型應(yīng)用于NACA0012 和Eppler387 翼型的轉(zhuǎn)捩模擬。另外在高超聲速領(lǐng)域中，Dash 和Papp 等[97-98]提出的SSGZkLε-γ的層流動能轉(zhuǎn)捩模型及其改進最具有代表性。Shi 等[99]基于高超音速靜音風洞實驗，提出了一種新的轉(zhuǎn)捩預測工程模型，其預測結(jié)果與HIFiRE-5b 的飛行數(shù)據(jù)非常一致。

3 三維邊界層轉(zhuǎn)捩預測eN方法研究進展

轉(zhuǎn)捩預測eN方法通常指基于線性穩(wěn)定性理論（linear stability theory,LST）[100]的自然轉(zhuǎn)捩預測方法，其理論思想是：邊界層內(nèi)存在各種頻率擾動向下游演化增長；針對層流場通過求解某種簡化形式的線性穩(wěn)定性方程獲得擾動增長率，然后沿擾動波傳播方向積分增長率獲得擾動幅值的指數(shù)增長倍數(shù)N值（又稱放大因子），最后判斷包絡(luò)的N值是否達到臨界的Ntr來確定轉(zhuǎn)捩發(fā)生位置；Ntr的值可通過轉(zhuǎn)捩實驗標定獲取。不過eN方法僅僅考慮了轉(zhuǎn)捩的線性穩(wěn)定性階段，未考慮轉(zhuǎn)捩的感受性和非線性失穩(wěn)階段，所以Ntr值的判據(jù)標準受限于對感受性與非線性失穩(wěn)的認識。因而，eN方法實際上是一個半理論半經(jīng)驗的方法。盡管如此，對于小幅值的環(huán)境擾動或來流擾動誘導的邊界層自然轉(zhuǎn)捩，因為擾動線性失穩(wěn)一般占據(jù)大部分區(qū)域，還因為感受性階段和非線性階段的貢獻可以在N值標定中間接體現(xiàn)，所以eN方法是比較可靠的轉(zhuǎn)捩預測方法。例如高超聲速巡航飛行器，高空大氣環(huán)境擾動通常為小擾動，飛行器表面的轉(zhuǎn)捩現(xiàn)象一般為自然轉(zhuǎn)捩過程，滿足eN方法使用條件。

eN方法早在1956 年[101-102]就已經(jīng)被提出來，對于飛行器大面積處邊界層三維性較弱時，eN方法能夠有效地預測轉(zhuǎn)捩位置。但對于具有復雜外形的三維邊界層，距離工程應(yīng)用還有很大距離[103]。如何面向復雜多塊化、結(jié)構(gòu)與非結(jié)構(gòu)化網(wǎng)格數(shù)據(jù)、有效識別轉(zhuǎn)捩模態(tài)、在三維空間上模態(tài)積分等是eN方法工程應(yīng)用面臨的重要技術(shù)難題。

采用平行流假設(shè)和展向均勻假設(shè)的LST 是常用的理論分析方法。對于橫流區(qū)和附著線區(qū)流動，由于流動在展向變化較緩，LST 能定性甚至定量準確地捕捉擾動波的某些重要信息（如頻率、增長率和相速度等），但對于流向渦區(qū)，流動在展向劇烈變化（展向梯度幅值與法向梯度幅值相當），此時LST 已不能準確預測流場失穩(wěn)行為，而考慮展向和法向流場變化的全局穩(wěn)定性分析方法，則可準確獲得上述三維流場的穩(wěn)定性特征[104]。

根據(jù)流動三維特征的復雜程度，全局穩(wěn)定性分析方法主要可以分為三類：流向基本不變、法向和展向快變的BiGlobal，流向緩變、法向和展向快變的PSE3D，流向、法向和展向均快變的三維特征值穩(wěn)定性分析方法（TriGlobal），下面分別介紹其基本思想。

3.1 二維特征值穩(wěn)定性分析方法（BiGlobal）

在開展全局穩(wěn)定性分析前，需要將不規(guī)則的物理域變換到規(guī)則的計算域下。通常采取的坐標變換形式為 (x,y,z)→(ξ,η,ζ)，其中 (ξ,η,ζ)∈[-1,1]，該坐標系的選取使得邊界層流場沿 ξ方向緩變，于是若采用平行流假設(shè)，則有如下分解：

其中q=(u,v,w,p,T)，u、v和w分別為流向、法向和展向速度，p為壓強，T為溫度。t為時間，代表時均場，為擾動形狀函數(shù)，ω為擾動頻率，α為待求的擾動波數(shù)（實部）和增長率（虛部），c.c.代表共軛。將上式代入Navier-Stokes方程，保留一階小量，可得空間二維特征值（BiGlobal）方程：

其中A0、A1、A2是包含時均場信息的算子，可采用高階有限差分離散，法向和展向根據(jù)計算域和不穩(wěn)定性性質(zhì)的不同分別采用數(shù)目不等的網(wǎng)格數(shù)。法向和展向一般選用齊次邊界條件，在對稱面處采用對稱或反對稱條件。式（2）的特征值采用隱式Arnoldi 迭代求解。

3.2 面推進拋物化穩(wěn)定性方程（PSE3D）

為考慮非平行效應(yīng)的影響，將流場作如下分解：

其中，假設(shè)時均場和擾動形狀函數(shù)均沿流向緩變，則在僅保留其一階流向?qū)?shù)項后，可得面推進拋物化穩(wěn)定性方程（PSE3D）：

其中 L 和 M為包含時均流場信息的微分算子，離散方式與BiGlobal 方程一致，初始條件由BiGlobal 給出，流向推進一般采用一階或二階隱式歐拉格式，為保證形狀函數(shù)沿流向緩變，需在每步迭代α使得

其中ME代表能量權(quán)矩陣，通常采用如下形式：

其中，γ為比熱比，T為溫度，ρ為密度。不同的能量權(quán)矩陣對結(jié)果有一定影響。

上述方法在研究流向渦與接觸線等具有展向局域分布的不穩(wěn)定性時具有較高效率，但在面對展向分布區(qū)域廣而展向波動劇烈的橫流不穩(wěn)定性時則需要很多展向網(wǎng)格點，內(nèi)存消耗大，計算時間長。為此，Chen 等[30]提出約化方法，利用橫流全局模態(tài)具有慢變的幅值和快變的相位的特點，將橫流特征函數(shù)作如下分解：

此時，通過試錯法或者迭代法選取合適的等效波數(shù) β，使得：

即形狀函數(shù)在展向的波動主要吸收進指數(shù)項 exp(iβζ)。將式（7）代入式（2），并消去指數(shù)項，可得約化全局穩(wěn)定性方程：

3.3 三維特征值穩(wěn)定性分析方法（TriGlobal）

當流場在三個維度都非緩變時（如三維凹腔、粗糙元等），流場不穩(wěn)定性需求解三維特征值問題，即TriGlobal。TriGlobal 問題矩陣規(guī)模巨大，顯式寫出系數(shù)矩陣并求解該矩陣特征值這一常規(guī)方法幾乎不可行。為此，Tu 等[105]提出基于隱矩陣投影的無矩陣全局穩(wěn)定性分解方法，為解決強三維穩(wěn)定性分析問題提供可行的思路。該方法的思想是通過反復調(diào)用CFD 求解器構(gòu)造原特征矩陣A（m×m維）的子空間S（n×n維，n?m）實現(xiàn)降維處理，進而找到工程上最關(guān)心的數(shù)個主特征值。具體思路如下：

首先選擇一系列線性無關(guān)的初始擾動：

通過高精度的CFD 工具，讓擾動場演化 ?t時間后得到，

其中矩陣B=exp(A?t)?？梢酝队疤幚順?gòu)造小矩陣

其中Q、R對應(yīng)Um的QR 分解矩陣，滿足Um=QR，且可以證明矩陣S與矩陣A的特征值 ω和特征向量 φ滿足以下關(guān)系，

Tu 等[105]的方法通過CFD 求解和子空間投影，避免了構(gòu)造巨型矩陣A，把全局穩(wěn)定性理論的巨型矩陣特征值問題轉(zhuǎn)化成小矩陣特征值問題，對于高雷諾數(shù)流動，矩陣降維可達4 個量級以上，使三維復雜飛行器的轉(zhuǎn)捩機理分析變得可能。

除上述方法外，基于數(shù)據(jù)驅(qū)動的動態(tài)模態(tài)分解（dynamic mode decomposition,DMD）方法（及其變種）也是常用的全局穩(wěn)定性分析方法[106]。

4 三維邊界層轉(zhuǎn)捩預測軟件研究進展

鑒于eN方法的理論優(yōu)勢與精準效果，21 世紀以來，國內(nèi)外陸續(xù)攻克了eN方法的相關(guān)技術(shù)問題，并匯編成轉(zhuǎn)捩預測軟件，使其在工程實踐中得到了廣泛應(yīng)用。下面介紹幾款基于eN方法開發(fā)的典型轉(zhuǎn)捩預測軟件。

4.1 eMalik 軟件

eMalik 是由NASA 下屬機構(gòu)蘭利研究中心（LaRC）的Malik 主持完成的一款二維穩(wěn)定性分析軟件，并于1991 年底交付給美國空天飛機（NASP）項目組使用，取得了良好的效果。在意識到該軟件對支撐高超聲速飛行器設(shè)計的重要性后，又將該軟件與CFL3D對接，實現(xiàn)層流場計算與轉(zhuǎn)捩分析的功能集成。1992 年底，Malik 還將三維穩(wěn)定性分析代碼的預發(fā)布版本交付給項目組。eMalik 軟件的主要功能是預測典型外形邊界層的轉(zhuǎn)捩模態(tài)，包括第一模態(tài)（TS 波）、第二模態(tài)、橫流模態(tài)、G?rtler 失穩(wěn)模態(tài)等。

eMalik 軟件的產(chǎn)品具體信息并未在互聯(lián)網(wǎng)上公開，文獻提到的兩份關(guān)于eMalik 軟件的報告（No.88-6和No.HTC-8902）國內(nèi)也沒有相關(guān)獲取渠道。不過，追溯Malik 早期發(fā)表的關(guān)于eN方法代碼的文章，可以發(fā)現(xiàn)，名為COSAL[107]的軟件可能與eMalik 軟件存在關(guān)聯(lián)。COSAL 軟件是最早開發(fā)的用于可壓縮邊界層轉(zhuǎn)捩預測的封裝代碼之一，然而并沒有集成到任何CFD 求解軟件上。它需要提供外部層流場，手動輸入所需的波參數(shù)，N值積分方法使用局限，對使用人員專業(yè)水平要求高，因而沒有得到工程上的廣泛應(yīng)用。后來以NASP 計劃為契機，發(fā)展了更為實用的eMalik軟件。還進一步發(fā)展了集成流場計算、穩(wěn)定性分析和轉(zhuǎn)捩預測功能的一系列軟件工具包（graphical transition prediction toolkit,GTPT）[108]。

eMalik 軟件在美國幾家具有軍工背景的單位中得到了使用。從幾篇公開的文獻可以看出，eMalik 軟件應(yīng)用的外形偏簡單，但仍對飛行器發(fā)展、風洞噴管設(shè)計等起到了支撐作用[109]。美國麥克唐納·道格拉斯公司（McDonnell-Douglas Corporation）Elias 等[110]針對NASP 幾何外形測試了eMalik 軟件，成功預測了高超聲速邊界層轉(zhuǎn)捩。普渡大學Schneider[111]利用eMalik 軟件分析了風洞噴管的G?rtler 失穩(wěn)對噴管內(nèi)表面轉(zhuǎn)捩的影響，有效評估了設(shè)計噴管對要求獲得低噪聲層流化噴流的達標情況。美國空軍基地Kimmel 等[112]針對橢圓錐模型利用eMalik 軟件預測了橫流轉(zhuǎn)捩特征，該模型后來成為了飛行試驗HIFiRE-5 的基準外形[113]，也就是美國空軍基地聯(lián)合澳大利亞國防科技組織發(fā)起的高超聲速飛行試驗計劃[114-115]的系列之五，目的是驗證美國新一代高超聲速飛行器發(fā)展的關(guān)鍵技術(shù)。但是，eMalik 軟件作為設(shè)計工具的功能還是很有限，其主要原因是不能分析復雜外形的流動穩(wěn)定性與轉(zhuǎn)捩特征，軟件運行需要大量的人工干預，而且在當時的計算條件下完成分析需要的時間比較長。也許是基于這些原因或局限，后來少有文章對該軟件進行更新報道。

4.2 LASTRAC 軟件

LASTRAC（Langley stability and transition analysis code）軟件是由美國NASA 蘭利研究中心Chang[116]于21 世紀初完成的一套功能強大的轉(zhuǎn)捩預測軟件。雖然它與eMalik 軟件同屬一家單位，開發(fā)時間晚于后者十幾年，但與eMalik 軟件并不是前后傳承關(guān)系。LASTRAC 軟件是從零開始開發(fā)的，基于面向?qū)ο筮M行設(shè)計與編程，以便于后期維護與擴展。LASTRAC軟件在功能上得到了很大拓展，它不僅包含基于LST 的傳統(tǒng)eN方法，還包含考慮非平行效應(yīng)的線性和非線性拋物化穩(wěn)定性分析方法（PSE）。這樣一方面可以采用傳統(tǒng)的基于LST 的eN方法來分析邊界層線性穩(wěn)定性與轉(zhuǎn)捩預測，還可以采用更為高級的基于PSE 的eN方法；另一方面可以利用線性和非線性PSE 模擬邊界層的轉(zhuǎn)捩過程，相較于直接數(shù)值模擬和大渦模擬，大幅降低計算時長。

從適用外形上看，LASTRAC 軟件主要有兩個版本，分別是2D 版[116]和3D 版[117]，于2003 年和2004 年公布。2D 版適用于二維邊界層、旋轉(zhuǎn)體邊界層、無限后掠機翼邊界層、二維射流、軸對稱噴流和渦流等。3D 版進一步引入了三維PSE 方法，將適用范圍拓展到準三維或一般三維邊界層，如有限后掠翼、錐體、有攻角模型等。LASTRAC 軟件基本上為航空航天飛行器的所有典型邊界層的轉(zhuǎn)捩預測與層流控制研究提供了可靠的工程實用工具。2018 年Kline 等[118]進一步考慮了高溫化學反應(yīng)效應(yīng)。2023 年，LASTRAC軟件被集成到了著名的流動仿真軟件OVERFLOW，實現(xiàn)了流動仿真與轉(zhuǎn)捩預測同時進行[119]。此外，LASTRAC 軟件的用戶界面設(shè)計可以減少用戶在識別擾動不穩(wěn)定頻率和波數(shù)方面的人工干預，提高軟件對用戶的友好性。

LASTRAC 軟件一問世就在美國軍工和高校單位推廣起來。NASA 蘭利研究中心的轉(zhuǎn)捩研究團隊（由著名轉(zhuǎn)捩學者Chang、Choudhari、Li Fei、Paredes 等組成）依靠該軟件分析了飛行試驗HIFiRE-1（零攻角和有攻角圓錐外形）和HIFiRE-5（橢圓錐外形）的轉(zhuǎn)捩機理[28,34,120-121]。該團隊還對美國的火星科學實驗（mars science laboratory,MSL）再入飛行器的轉(zhuǎn)捩進行了分析，有效評估了熱載荷分布對精確著陸所需升阻比的影響情況[122]。此外，該團隊還聯(lián)合美國航空航天研究所將LASTRAC軟件與著名非結(jié)構(gòu)有限體積RANS 求解器FUN3D 軟件進行耦合，發(fā)展了用于轉(zhuǎn)捩分析的基本流剖面提取技術(shù)，建立了轉(zhuǎn)捩預測模型，預測了G-3 整機飛行器的N值分布，有效模擬了NLF(1)-0416 翼型和6∶1 標準橢球體轉(zhuǎn)捩前后的熱流摩阻變化（見圖4）[123-124]。美國霍普金斯大學應(yīng)用物理實驗室利用該軟件分析了BOLT 外形的轉(zhuǎn)捩模態(tài)失穩(wěn)機理，支撐了試驗外形的設(shè)計與優(yōu)化以及飛行驗證試驗[125]。該試驗是由美國空軍研究實驗室（AFOSR）牽頭的，主要用于研究具有后掠前緣的凹曲面三維邊界層轉(zhuǎn)捩問題（見圖5）[126-127]。美國波音公司利用該軟件對美國圣母大學馬赫數(shù)6 靜音風洞噴管進行模態(tài)分析，有效支撐了優(yōu)化噴嘴外形和壁溫分布的設(shè)計[128]。

圖4 LASTRAC 和FUN3D 耦合計算得到的G-3 飛行器轉(zhuǎn)捩摩阻等值線[124]Fig.4 Transition friction contours of G-3 aircraft obtained by coupling LASTRAC and FUN3D[124]

圖5 BOLT 飛行器標模（左）和LASTRAC 預測的轉(zhuǎn)捩陣面（右）[126]Fig.5 Aircraft standard model of BOLT (left) and transition front predicted by LASTRAC (right)[126]

4.3 LILO 軟件

LILO 軟件是德國宇航中心（DLR）針對航空飛機轉(zhuǎn)捩問題而設(shè)計的一款工業(yè)分析軟件。該軟件于2004 年9 月首次公開了使用說明書，2006 年7 月又發(fā)布了2.1 版[129]。它是基于Fortran 77 語言開發(fā)的，可分析考慮表面曲率的不可壓/可壓縮準三維邊界層的轉(zhuǎn)捩預測問題。從適用范圍上來說，LILO 軟件相比前面介紹的兩款美國軟件，通用性較弱，在高速航天領(lǐng)域的應(yīng)用方面更是未見到相關(guān)報道。它的優(yōu)勢在于應(yīng)用了數(shù)據(jù)庫方法，包括特征值數(shù)據(jù)庫和轉(zhuǎn)捩N值數(shù)據(jù)庫，與RANS 流場求解器進行耦合，且在計算過程中與流場求解器交互，實現(xiàn)了軟件的高效計算和自動預測功能。特征值數(shù)據(jù)庫是以神經(jīng)網(wǎng)絡(luò)形式表示，可以基于基本流剖面（如Falkner-Skan-Cooke 剖面）通過機器學習訓練獲得[130]。該數(shù)據(jù)庫用于尋找初始模態(tài)值，有助于計算快速收斂到正確解。轉(zhuǎn)捩N值數(shù)據(jù)庫是通過大量實驗標定建立的轉(zhuǎn)捩判據(jù)。當計算的N值曲線達到轉(zhuǎn)捩N值數(shù)據(jù)庫中的臨界值時，即可判斷轉(zhuǎn)捩發(fā)生。該數(shù)據(jù)庫可以將eN方法以全自動的方式應(yīng)用于轉(zhuǎn)捩預測，并且在基本無需提前認識轉(zhuǎn)捩特征、無需用戶干預的情況下，完成轉(zhuǎn)捩位置的迭代分析[131]。

LILO 軟件主要用于指導機翼、機身、短艙、尾翼等航空飛機部件的自然層流（natural laminar flow,NLF）設(shè)計和混合層流控制（hybrid laminar flow control,HLFC）設(shè)計，是一款專用性很強的轉(zhuǎn)捩預測軟件。它一般不單獨使用，而是集成到一些流場模擬軟件上，實現(xiàn)轉(zhuǎn)捩預測分析與轉(zhuǎn)捩場模擬的同步計算。例如，DLR 將LILO 軟件集成到自身開發(fā)的兩款三維RANS求解器FLOWer[130-131]和TAU[132]上，主要服務(wù)于歐洲空客公司。

FLOWer 是一個基于有限體積法開發(fā)的三維可壓縮RANS 軟件，適用于多塊結(jié)構(gòu)化網(wǎng)格，可求解定?；蚍嵌ǔＡ鲌鰡栴}。當轉(zhuǎn)捩預測模塊耦合到FLOWer求解器時，耦合系統(tǒng)成為一個有效預測轉(zhuǎn)捩位置、快速模擬從層流到湍流的轉(zhuǎn)捩流場的CFD 工具，在已知模態(tài)（T-S 波和橫流模態(tài)）轉(zhuǎn)捩判據(jù)的情況下，全程自動化運行，基本無需用戶干預。該套耦合系統(tǒng)的實用性在ONERA 設(shè)計的M6 機翼[133]上得到了檢驗。

TAU 是由DLR 于20 個世紀末開發(fā)完成的一款非結(jié)構(gòu)/混合網(wǎng)格三維RANS 求解器。該求解器可實現(xiàn)高精度層流邊界層流動剖面的計算或提取，更適用于邊界層轉(zhuǎn)捩預測分析。布倫瑞克大學流體力學研究所（ISM）在TAU 軟件中引入耦合了LILO 轉(zhuǎn)捩預測模塊，實現(xiàn)了對三維邊界層轉(zhuǎn)捩的快速自動化預測功能[134]。后來兩家單位不斷合作，將TAU 軟件的轉(zhuǎn)捩預測功能進一步優(yōu)化，對標應(yīng)用于工業(yè)飛機中機翼、尾翼、機身和短艙等復雜配置[135-139]，如荷蘭航空航天實驗室設(shè)計的二段NLR-7301 翼型（圖6 左）[136]、DLR 設(shè)計的DLR-F11 翼身組合體（圖6 右）[136]、NASA的CRM-NLF 機翼（圖7）[139]。

圖6 LILO 預測的NLR-7 301 翼型（左）和DLR-F11 翼身組合體（右）的轉(zhuǎn)捩結(jié)果[136]Fig.6 Transition result of NLR-7 301 Airfoil (left) and DLR-F11 wing-body (right) predicted by LILO[136]

圖7 LILO 預測的CRM-NLF 翼型轉(zhuǎn)捩陣面[139]Fig.7 Transition front of CRM-NLF Airfoil predicted by LILO[139]

4.4 HyTEN 軟件

HyTEN 軟件[140]是一款用于高超聲速邊界層穩(wěn)定性分析與轉(zhuǎn)捩預測的軟件，由氣動中心聯(lián)合天津大學在國家數(shù)值風洞NNW[141]（National Numerical Windtunnel,NNW）工程資助下開發(fā)的。該軟件最早是由天津大學匯編集成LST 代碼和轉(zhuǎn)捩預測代碼，經(jīng)過多個版本迭代更新，成功實現(xiàn)從二維、準三維外形到三維外形的應(yīng)用[142-146]。其中較為典型的是天津大學黃章峰[147]發(fā)展的流動穩(wěn)定性分析和轉(zhuǎn)捩預測（flow instability and transition,FIT）軟件，該軟件實現(xiàn)了飛行器大面積處轉(zhuǎn)捩的自動、高效、準確預測，且在工程單位得到了應(yīng)用。后來以NNW 為契機，氣動中心聯(lián)合天津大學，瞄準工程化、實用化需求，在原來基礎(chǔ)上增加了全局穩(wěn)定性分析模塊、感受性模塊和天地轉(zhuǎn)捩數(shù)據(jù)模塊等，提升了代碼的適用性，并對轉(zhuǎn)捩預測代碼進行優(yōu)化、測試、集成，形成了HyTEN 軟件。此外，HyTEN 軟件也在國內(nèi)多個飛行器設(shè)計部門得到了應(yīng)用，并取得了良好效果。

HyTEN 軟件是基于面向?qū)ο笏枷朐O(shè)計的，軟件核心部分采用Fortran 語言編寫，繼承了Fortran 語言計算效率高的特點，同時兼顧軟件的代碼維護、功能擴展和數(shù)據(jù)對接等。HyTEN 相較于美國eMalik 和LASTRAC，優(yōu)勢包括理論更完備（如全局穩(wěn)定性理論），通用性更強，適用更復雜三維外形，靈活匹配多塊結(jié)構(gòu)和非結(jié)構(gòu)網(wǎng)格流場等。HyTEN 還能自動識別和處理給定的網(wǎng)格和流場，不需要輸入其他網(wǎng)格信息、邊界信息、流場信息等，然后利用插值技術(shù)自動提取邊界層流場剖面信息，從而自動完成流動穩(wěn)定性與轉(zhuǎn)捩預測的前處理準備工作，整個過程無需人工干預。因此，該軟件可以作為獨立軟件分析不同CFD 求解器計算的流場數(shù)據(jù)，適合在行業(yè)內(nèi)推廣應(yīng)用。軟件采用了OpenMP 并行，大大提高計算效率；可調(diào)用感受性數(shù)據(jù)庫，提供N值積分初始幅值，提高轉(zhuǎn)捩預測準度；可調(diào)用轉(zhuǎn)捩的數(shù)據(jù)庫，加快計算收斂速度，實現(xiàn)轉(zhuǎn)捩高效預測。另外，除了一維線性穩(wěn)定性理論，軟件還可以選用全局穩(wěn)定性理論，實現(xiàn)一般三維邊界層的流動穩(wěn)定性分析與轉(zhuǎn)捩預測。

HyTEN 軟件因其通用性和友好性，被廣泛應(yīng)用。目前公開的有：Wan 等[148]分析了有攻角飛行試驗的轉(zhuǎn)捩問題，證實了飛行試驗中測量的5～20 kHz 低頻信號為橫流行波，首次澄清了困擾學界多年的關(guān)于定常-非定常橫流模態(tài)主導性的爭議；涂國華等[149-150]將風洞試驗與飛行試驗數(shù)據(jù)進行關(guān)聯(lián)，建立了風洞試驗轉(zhuǎn)捩數(shù)據(jù)關(guān)聯(lián)方法和轉(zhuǎn)捩天地相關(guān)性，為工程上飛行轉(zhuǎn)捩預測提供了有效方法；陳堅強等[151]分析了HyTRV 流場特征與邊界層穩(wěn)定性特征；對于復雜三維外形，黃章峰等[12]對HIFiRE-5 橢圓錐（圖8 左）、Xiang 等[76]對HyTRV 升力體（圖8 右）的轉(zhuǎn)捩預測結(jié)果（線條）與實驗結(jié)果（背景顏色）對比都有非常好的吻合度。

圖8 HyTEN 預測的HIFiRE-5（左）[12] 和HyTRV（右）[76]的轉(zhuǎn)捩陣面Fig.8 Transition front of HIFiRE-5 (left)[12] and HyTRV (right)[76] predicted by HyTEN

4.5 其他轉(zhuǎn)捩軟件

以上介紹了4 款典型軟件，但基于eN方法的轉(zhuǎn)捩預測軟件不止于此，且在工程應(yīng)用與科學研究中都發(fā)揮了重要作用。例如，美國明尼蘇達大學受美國空軍研究實驗室贊助開發(fā)了一款轉(zhuǎn)捩分析軟件（STABL）[152]，早于LASTRAC 軟件考慮了高溫熱化學平衡或非平衡效應(yīng)，主要用于解決高馬赫數(shù)飛行器的轉(zhuǎn)捩問題，如HTV-2 飛行器、X-51 飛行器、Reentry-F 再入飛行器等。進一步地，他們將STABL 軟件耦合到他們新發(fā)展的非結(jié)構(gòu)網(wǎng)格求解器US3D 軟件[153]，可用于分析飛行器表面燒蝕等物理現(xiàn)象。該軟件已經(jīng)在美國空軍、NASA、桑迪亞實驗室、各大飛機制造商和高校等推廣應(yīng)用[154]。

此外，美國得克薩斯州農(nóng)工大學基于PSE 方法開發(fā)了EPIC 軟件[155]，可解決有攻角圓錐、裙錐、HIFiRE-5 橢圓錐等三維邊界層的線性和非線性穩(wěn)定性問題。ONERA 的Arnal 等[156]根據(jù)自相似流動特征，通過穩(wěn)定性分析形成數(shù)據(jù)集或關(guān)系式，建立了轉(zhuǎn)捩預測的數(shù)據(jù)庫方法。他們前后工作三十多年完成了數(shù)據(jù)庫方法從二維到三維邊界層甚至復雜構(gòu)型的應(yīng)用發(fā)展，可以說是該方面工作的集成者。前面提到的LILO 與FLOWer 耦合系統(tǒng)、LILO 與TAU 耦合系統(tǒng)引入的數(shù)據(jù)庫方法都是基于該成果發(fā)展的。后來他們將數(shù)據(jù)庫方法耦合到ONERA 開發(fā)的elsA 軟件中，并在空客A310 翼型上得到了檢驗與應(yīng)用[157]。瑞典國防研究局將轉(zhuǎn)捩數(shù)據(jù)庫方法耦合到自主開發(fā)的FOI-Edge 求解器，成功應(yīng)用于美國HiLiftPW-1 中的大弦三段梯形翼型的轉(zhuǎn)捩預測[158]。俄羅斯的Boiko 等將eN轉(zhuǎn)捩預測模塊耦合集成到LOGOS 軟件上，形成轉(zhuǎn)捩的自動分析軟件LOTRAN[159-160]，在機翼、長橢球體和發(fā)動機短艙成功應(yīng)用。

在國內(nèi)，西北工業(yè)大學將eN方法耦合到三維RANS 求解器PMNS3D，實現(xiàn)對DLR-F4 翼型[161-162]、空客A320 翼型[163-164]、自然層流機翼[165]、翼型優(yōu)化[166]的轉(zhuǎn)捩自動預測。中國航空工業(yè)空氣動力研究院將eN方法先后耦合到UNSMB[167]和AIR_OVERSET[168]數(shù)值模擬平臺上，并采用法國航空航天研究中心ONERA-M6 翼型、麥道航空公司30P30N 三段翼型和德國宇航中心DLR-F4 翼型等算例進行了驗證。

5 結(jié)論與展望

高超聲速邊界層轉(zhuǎn)捩預測是飛行器熱防護設(shè)計與減阻增升設(shè)計需要考慮的重要問題，本文主要針對三維邊界層的轉(zhuǎn)捩機理、轉(zhuǎn)捩模型、預測方法、預測軟件等四方面研究現(xiàn)狀進行了梳理與討論，主要結(jié)論與認識如下：

1）相比二維邊界層，三維邊界層呈現(xiàn)三類特有的流動結(jié)構(gòu)，即橫流、流向渦和附著線流動，每類結(jié)構(gòu)具有不同的失穩(wěn)特征和轉(zhuǎn)捩機理，需要分別加以深入研究。目前，三維邊界層的線性失穩(wěn)模態(tài)特征已逐漸清晰，但非模態(tài)擾動演化和模態(tài)非線性相互作用機制還不清楚，非線性轉(zhuǎn)捩機制的研究近乎空白，亟待后續(xù)大量持續(xù)的工作。

2）轉(zhuǎn)捩模型的主要實現(xiàn)方式是通過在湍流模型輸運方程基礎(chǔ)上，耦合求解轉(zhuǎn)捩模型輸運方程與湍流模型輸運方程來實現(xiàn)的。伴隨轉(zhuǎn)捩模式理論的發(fā)展，轉(zhuǎn)捩模型處于不斷發(fā)展與改進的過程之中。本文主要介紹了低雷諾數(shù)湍流、γ-Reθt、k-ω-γ、層流動能4 種常用轉(zhuǎn)捩模型的構(gòu)造、適用性及應(yīng)用情況。目前在商業(yè)軟件上部署的轉(zhuǎn)捩模型數(shù)量相對較少，主要是γ-Reθt轉(zhuǎn)捩模型及其改進模型。轉(zhuǎn)捩模型應(yīng)用范圍包括內(nèi)外流條件下各類外形的二維、三維的自然轉(zhuǎn)捩、旁路轉(zhuǎn)捩、分離誘導轉(zhuǎn)捩，速度域涵蓋從低速到高超聲速。但由于轉(zhuǎn)捩機制的復雜性，單一轉(zhuǎn)捩模型在面對來流條件和幾何外形的巨大差異時難以表現(xiàn)出一致的普適性。這對未來發(fā)展具有通用、普適、可靠的轉(zhuǎn)捩模型提出了更高要求。

3）eN方法是基于流動穩(wěn)定性理論的自然轉(zhuǎn)捩預測方法，是目前理論依據(jù)相對充分的一種半理論半經(jīng)驗方法。盡管該方法的理論基礎(chǔ)相當成熟，但如何將科學方法上升到實用的工程技術(shù)是國內(nèi)外一直不斷攻關(guān)的技術(shù)問題。本文梳理了基于eN方法的轉(zhuǎn)捩預測軟件的開發(fā)與工程應(yīng)用情況，重點介紹了eMalik、LASTRAC、LILO、HyTEN 四款典型轉(zhuǎn)捩軟件的優(yōu)缺點。國內(nèi)軟件的發(fā)展相比歐美較晚，但經(jīng)過十幾年的努力實現(xiàn)局部趕超，并且開始推廣到相關(guān)工業(yè)部門。

基于eN方法的轉(zhuǎn)捩預測軟件可實現(xiàn)對航空航天飛行器的橫流、流向渦、附著流等典型三維邊界層的自然轉(zhuǎn)捩預測，特別是在高超聲速領(lǐng)域取得了理論與技術(shù)上的突破，包括實現(xiàn)三維邊界層轉(zhuǎn)捩預測、耦合高溫氣體效應(yīng)、轉(zhuǎn)捩的自動化預測等，深受國防工業(yè)部門的青睞。但是，面對多模態(tài)轉(zhuǎn)捩現(xiàn)象，仍缺乏統(tǒng)一的轉(zhuǎn)捩判據(jù)，需要加強轉(zhuǎn)捩判據(jù)數(shù)據(jù)庫的建設(shè)。