航空航天CFD物理模型和計算方法的述評與挑戰

閻 超, 屈 峰, 趙雅甜, 于 劍,*, 武從海, 張樹海

(1. 北京航空航天大學 航空科學與工程學院, 北京 100191；2. 西北工業大學 航空學院, 西安 710072；3. 中國空氣動力研究與發展中心 空氣動力學國家重點實驗室, 綿陽 621000)

0 引 言

CFD在國內外航空航天領域中的地位日益重要，成為飛行器研制中不可或缺的重要手段。這不僅表現在對CFD越來越多的依賴，一些時候，如在極端氣動力/熱環境、流動結構或流動機理剖析、緊急任務或復雜事故分析等情況下，CFD是唯一可以依賴的手段。

CFD歷經百年，但真正實用化只有約40年。40年來，CFD在應用上的普及和發展速度如此之快，幾乎超出了所有人的意料。但CFD在計算方法和物理模型等核心理論上的進步卻是步履蹣跚、甚至止步不前，這也幾乎超出了所有人的意料。例如，1979年, 美國工程院院士、NASA Ames研究中心航天部主任Chapman預測[1]LES將于20世紀90年代應用于飛行器全機的計算，可事實是：距他預測的時間過去近30年了，這個預測至今還沒有實現，而且在可預見的未來也不太可能實現。雖然每年都有數以千計的CFD論文等成果發表，但真正被大家接受、甚至走入應用的卻很少。

為了進一步論證上述可能引起爭議的觀點，這里簡述一下國內外的CFD現狀。

在國內，眾所周知，CFD的應用無論是深度還是廣度，都有了飛速的發展，一些十年前我們甚至想都不敢想的復雜CFD計算，也許今天新入門的研究生都可以輕松完成。這里給出一個實例：圖1是筆者碩士研究生在2017年完成的復雜飛行器全機非定常湍流DES模擬，全模7000萬網格，湍流能譜平均誤差2%。十年前，這樣的CFD計算是不可想象的。

圖1 復雜飛行器全機非定常湍流瞬時渦結構的精細模擬(Q=300，Ma染色，Ma=2.5，H=20 km，α=10°)

但是，CFD真的大踏步進步了嗎？不妨簡單考察一下CFD的核心理論也就是它的計算方法和物理模型的現狀：

(1) 計算方法：目前航空航天領域應用最廣的Roe格式是1981年發表的，39年了，但它早被人知的原始缺陷依然存在；

(2) 物理模型：目前應用最廣的湍流模型之一k-ε模型(它也是另一個應用最廣的湍流模型之一——SST模型——的組成部分)出現在1972年，48年了，但它早被人知的原始缺陷也依然存在。

也就是說，在CFD核心方法上，我們已經陷入了長達40年沒有本質進步的困境。更讓我們憂慮的是，目前，仍然很難看到走出這些困境的曙光。

2017年至2019年，中國空氣動力研究與發展中心等多家單位聯合開展了“第一屆航空CFD可信度研討會”[2]，公開發布了計算模型、基礎網格和試驗數據，對多家單位提供的計算結果進行了統計分析，并與相關風洞試驗結果進行了對比。會議認為：國內自主研發的CFD軟件與美國先進CFD軟件計算結果精度相當，CFD計算結果與風洞試驗結果定性吻合較好,定量方面距離精細化設計尚有不小的差距。

為了全面了解和評估航空航天中CFD的現狀、梳理存在的問題，美國、歐洲等的政府部門、科研機構、企業、高校等開展了很多的國際合作、項目研究、專題討論、系列研討會等，比較知名的如：The European Research Community on Flow, Turbulence, and Combustion (ERCOFTAC)、Workshops on CFD Uncertainty Analysis、CFDVAL2004 Workshop、The Drag Prediction Workshop (DPW)、High Lift Prediction Workshop (HiLiftPW) series、A European Project on the Development of Adaptive Higher-Order Variational Methods for Aerospace Applications(ADIGMA)、Advanced Turbulence Simulation for Aerodynamic Application Challenges(ATAAC)。這其中DPW影響最大，DPW參與者幾乎包含了國際上所有知名的航空航天單位，其大量CFD結果的數據統計、分類歸納和物理分析等成果非常珍貴、意義重大，因為本文的部分觀點也基于這些成果，因此這里先簡單介紹DPW。

DPW起源始于1998-1999年AIAA會議的走廊討論，這些走廊討論的幾位成員在AIAA應用空氣動力學技術委員會的贊助下，于2000年1月正式成立了第一屆DPW組織委員會。此后，經過二十年的不懈努力和六次國際研討會，DPW取得了出乎意料的成功。與會者和組織者都認為DPW是絕對的成功，因為它是CFD研究者、軟件開發者、航空航天CFD使用者等“一線基層”集體努力的計算分析、公開誠實的交流討論。DPW通過廣泛、長期、充分、公開、多角度、多層次的國際合作和研討，評估CFD的實用性能(尤其關注阻力)、公正評價CFD軟件的有效性、發現需要進一步研究的問題并探索其解決方法，DPW因此成為航空航天CFD領域影響最大的國際活動。會議發布了大量的幾何構型、各種網格、數據統計、結果對比分析和建議等，這些均在公共可訪問的網站上發布，包括數據庫、研討會記錄、大量公開發表的論文、相關出版物和演示文稿等。

DPW的成果太豐富，本文只能根據筆者自己的閱讀理解，簡單地濃縮其要點(為簡潔，以下用DPW后面加數字N表示第N次DPW會議)[3]：

1) 普遍對CFD結果失望，這主要表現在不同軟件的CFD結果散布較大，超出了預期。雖然通過不斷提高網格質量和數量、強化收斂性等，降低了散布程度，但仍然大于飛機設計者對CFD的期望，而且俯仰力矩總體上并沒有顯示出隨著網格分辨率的提高而趨于收斂的明顯趨勢。整體說來，CFD置信度仍未降到可與實驗相競爭的水平。

2) 網格問題(網格質量、網格分辨率、網格收斂性)一直是DPW研究的焦點，數據統計證實了網格質量和網格分辨率的重要性。網格數量也從最初的250萬平均基準網格、到后來最大24億個網格。但最終并沒有給出明確的或指導性結論，計算結果散布問題仍然存在。比如：對于類似的網格分辨率，DPW6結果的變化通常比DPW5更高，DPW6俯仰力矩的變化也高于DPW3的結果。

3) 會議通過設計各種獨特的網格集合，以盡量消除網格變化對計算結果的影響。但結果是，數據仍然存在較大的散布，所以會議認為湍流模型、軟件、邊界條件和迎風格式等其他因素也很重要。因為工作量大、協調困難等原因，對這些問題沒有開展深入的研究。

4) 統計分析DPW2會議的結果，在網格、程序、湍流模型、轉捩模型、黏性模型等影響CFD計算結果的因素中，發現湍流模型、轉捩模型以及網格對CFD計算結果的影響位居前三強，如圖2所示，分別占到約15%、11%和11%，這也是世界上第一次量化認識了CFD結果的影響因素。

圖2 影響CFD計算精度的因素統計結果[4]Fig.2 Comparison of CFD components on accuracy[4]

5) 從DPW2開始，就有一個共識：關鍵區域的流動分離使得很難對網格收斂性和阻力預測得到有意義的結論。很遺憾，后來的會議雖然針對這一問題設計了多種方法，這一共識最后也沒有得到確認或否定，因為影響因素多且復雜。例如，網格細化后，分離區變大還是變少？對于這一個看似簡單的問題，幾個世界知名的CFD軟件，給出了令人困惑的結論：變小(CFL3D,OVERFLOW)、變大(UPACS,TAS,FUN3D,ONERA-elsA，BCFD)、不變(Edge，DLR-TAU)。

實際上，HiLiftPW和國內“第一屆航空CFD可信度研討會”[2]也有類似的結論：網格規模的增加并沒有明顯提高計算結果之間的數據散布度；應該加強CFD方法、湍流模型和可信度分析方法研究。

CFD在應用上的飛速前進和理論上的步履蹣跚，讓我們有必要認真地審視他的現狀、籌劃他的未來。確實，近些年來，美國、歐洲、日本等陸續發表了一些這方面的研究報告。其中，特別要提到的是，2014年NASA組織了185名有關領域的資深專家，進行了150多次調查和討論，發表了研究報告“CFD Vision 2030 Study: A Path to Revolutionary Computational Aerosciences” (以下簡稱CFD2030)[5]。該報告評估了航空航天領域CFD的成就與不足、需求與差距，對CFD在2030年應具備的能力進行了規劃、給出了實現的具體建議。這是一份在CFD發展史上具有重要意義的指導性文件。顯然，我國也有必要思考類似的問題，很高興國內的一些單位已經開始這樣有意義的工作。

本文試圖對航空航天領域CFD應用的關鍵問題，進行一點梳理和思考。鑒于CFD涉及的面太廣，很難在一篇期刊論文中面面俱到地論述清楚，如NASA的CFD2030共58頁，雖然作者十年前也曾做過這樣的努力[6]，但是十年后的今天，為了論述得更深入，本文欲聚焦一下重點：面向航空航天領域的CFD應用，針對其核心問題即CFD計算方法和物理模型，敘述、評論其發展的現狀、存在的問題以及面臨的挑戰，以期認識現狀、展望未來、推動應用。

1 物理模型

無論是在實際工程應用、還是在科學研究領域，目前廣泛使用的CFD方法是基于物理模型的雷諾平均N-S方程(Reynolds Averaged Navier-Stokes Equations, 簡稱RANS方程)方法。其基本分類如圖3所示。即：雷諾應力模型(Reynolds-Stress Model, RSM)和渦黏性模型(Eddy Viscosity Model, EVM) 兩大陣營組成了RANS CFD的物理模型。目前廣泛使用的物理模型幾乎都是線性EVM，盡管它有一些“致命的弱點”，但因為其簡單、經濟、魯棒性好、容易實現、歷史悠久、性能特點了解較多、使用者眾多等原因，因此本文重點討論線性EVM，主要包括湍流模型和在湍流模型基礎上通過補充經驗關系式或輸運方程發展而來的轉捩模型。

圖3 RANS方法分類Fig.3 Classification of RANS models

1.1 湍流模型

在CFD長期的發展、研究和應用中，越來越認識到湍流模型的重要性。Spalart[7]在湍流模型綜述論文指出“湍流模型是CFD的決定因素”；2018年， Durbin[8]在Annu.Rev.FluidMech.中明確“湍流模型是大量應用CFD的核心”。

但是，湍流模型一般被認為是CFD中最大的誤差源、最薄弱的環節、“最不靠譜”，多次被稱為是CFD的Achilles heel(原指荷馬史詩中英雄阿喀琉斯的腳跟，現在一般指致命的弱點、要害)[9-10]。因此當CFD結果與實驗結果不一致時，湍流模型很容易成為替罪羊[10]。事實上，本文開頭介紹的DPW2的統計分析結果也表明：湍流模型對CFD計算精度的影響最大。因此，湍流模型成為制約CFD發展和應用的主要瓶頸。

即便在目前熱門的RANS/LES混合類方法中，湍流模型仍然是成功的基礎：因為一般湍流主要在近壁產生、輸運、耗散等，因此混合方法強烈依賴近壁區湍流特性模擬的準確性，這些都是由湍流模型決定的。

圖4表示的是在過去五十年中，商業飛行器設計一個循環中所需的試驗模型數隨時間變化。圖中藍色軸線表示在過去五十年中，在商業飛行器設計一個循環中所需的試驗模型數隨時間的變化，紅色的是常用湍流模型發展的時間軸。由圖可見，湍流模型在近二十多年中，沒有取得顯著的進展。結合下面紅色的常用湍流模型發展時間軸來看，模型發展之初，在應用中成就斐然，試驗模型數顯著下降。但是近二十多來來，該數據并沒有發生大的變化，導致這一瓶頸的主要原因是湍流模型的預測精度不足。

圖4 商業飛行器設計一個循環中所需的試驗模型數隨時間變化[11]Fig.4 History of number of wind tunnel test during one design period of commercial airplanes[11]

2010年美國航空航天局(AIAA)下屬的湍流模型基準測試工作小組針對湍流模型的應用和未來發展進行了深刻的討論[12]，結論之一即是：盡管對于許多工程流動的模擬不盡完美，RANS將在未來20～50年內依舊得到廣泛的使用。目前，國內外CFD界普遍認為：在未來的20～50年，基于湍流模型的RANS方法仍然是CFD的主流。

近20多年來，隨著計算機能力的普遍提升，可以通過使用數量更多、質量更好的網格來提高計算格式方面的精度。但是湍流模型的問題來自本身固有的缺陷，改善網格對湍流模型并沒有本質的幫助，因此，湍流模型的問題愈加突出。

首先，這里先簡單討論一下Boussinesq線性渦黏性湍流模型的基本原理和特點：Boussinesq基于渦黏性假設，認為湍流脈動引起的動量交換與分子熱運動引起的動量交換存在可以類比的關系，引入渦黏性系數將雷諾應力與平均速度應變率聯系起來，建立了線性渦黏性模型的本構關系：

(1)

線性渦黏性模型建立了雷諾應力——平均速度應變率之間的線性關系(應力和應變之間用一個標量的渦黏性系數，而不是張量聯系)，這導致了諸多缺陷：

·各向同性：渦黏性系數的各向同性或標量特性；

·坐標不變性：雷諾應力依賴于應變率張量；

·雷諾應力張量與平均速度應變率之間瞬時平衡；

·雷諾應力與平均場應變率方向相同；

……

各向同性是EVM的致命缺點之一，如航空航天領域最關心的壁湍流，由于壁面的抑制作用而具有非常明確的各向異性。

這些缺陷導致了很多后果，如無法正確預測分離流動、角區二次流、旋轉效應、曲率效應等等。大量的理論分析和數值實驗表明：一般說來，渦黏性模型只有對二維剪切主導、正應力不重要的簡單流動才能給出可靠的結果；對于高速、轉捩、流線彎曲、旋轉、逆壓梯度，尤其是分離流動等復雜湍流，渦黏性模型存在很大缺陷。

最常見的湍流是邊界層等壁湍流，在航空航天領域尤其如此。湍流模型的近壁區處理是一個重要、而又常被忽視的問題。由于近壁區流動梯度非常大，而且該區也是湍流產生、輸運、耗散、擴散等效應最強烈的區域，模型中的近壁區處理方法對最終的湍流模擬精度如摩擦系數、熱流等影響很大，因此是一個非常重要、不幸又是一個非常困難的問題，更嚴重的是很多湍流模型的使用者對這個關鍵、復雜而又困難的問題認識不足，因此，它又常常成為湍流模型中的“替罪羊”。

按照近壁區處理方式的不同，一般把湍流模型分為低雷諾數湍流模型和高雷諾數湍流模型，前者可以直接積分到壁面，后者一般使用壁函數等橋接方法。

低雷諾數湍流模型，依賴引入一些以湍流雷諾數和壁面距離等表達的阻尼函數，借此考慮對模化輸運方程中各項的近壁影響，這方面的研究是湍流模型經久不息的研究熱點之一。實際上，目前廣泛使用的k-ε模型自1972年提出后，近50年來的發展主要就是關于黏性阻尼函數和耗散率方程的修正。

目前，低雷諾數模型的使用者較少，這是因為存在以下問題：(1) 一些方程系數和函數需要調整且有一定的任意性，對大壓力梯度流動等缺乏普適性；(2) 近壁大梯度等帶來的數值剛性困難；(3) 網格量多、計算步長約束大、計算量大。

高雷諾數湍流模型不使用壁阻尼或壁修正等方法，而多是使用壁函數。壁函數就是放棄緊鄰壁區域的計算，而在該區域的頂部(壁上第一個網格點處)施加邊界條件。在緊鄰壁區域內，假定湍流和平均速度遵循指定的關系，如對數律。這就避免了對網格要求很高的緊鄰壁區域的求解，大大減少了近壁區的網格量，增大了時間步長，因此有效提高了計算效率。但其缺點也很明顯：所指定的關系缺乏普適性，尤其是復雜流動，因此計算精度和可靠性較差。而且，在應用中，實際生成的網格常常不能滿足壁函數關于y+等較苛刻的要求，這會帶來明顯的誤差、而且不容易被重視和發現。

回顧湍流模型的發展歷史，不難發現，人們在更經濟的壁函數和更精細可靠的直接積分到壁面方法之間徘徊、曲折前進。但目前，大多數計算者都更喜歡使用壁函數，因為其簡單、經濟、魯棒，當然這是以犧牲精度和普適性等為代價的。也因此可以推論：對壁函數進行進一步改進和通用化是當前CFD最急迫的任務之一。

湍流模型給使用者帶來的困難不僅是理論上的缺陷，還包括認識上的困惑。這種困惑一方面來自湍流模型本身的種類很多、假設條件和適用范圍復雜、各種經驗性系數繁多等等，另一方面來自湍流模型本身性能表現的模糊不清甚至是自相矛盾[10]，令人困惑的是，同一種模型在不同的軟件上得出的結論也常不一致[10]、甚至一些知名軟件也是這樣。另外，研究和應用實踐也發現：無論是k-ε、k-ω還是SST湍流模型，其中封閉參數的微小改變可能會導致模式性能的顯著改善或惡化。近年來的不確定度量化分析也證實了這一結論。

準確、系統地認識湍流模型的性能是研究者和使用者都非常關心的問題。但如前所述，雖然長期以來這方面開展了大量的研究、分析、評估等等，但直到今天，即便對于最簡單或最常用的湍流模型，我們也無法給出嚴格、明確的性能結論。盡管如此，一般說來，經過幾十年的不斷實踐、總結和摸索，對各種湍流模型的性能還是有了一些粗略的認識，這里筆者根據自己的理論認識和實踐經驗給出航空航天領域常用湍流模型的一些一般性的認識：

1) SA模型：該模型從經驗和量綱分析出發得到了關于渦黏性的輸運方程，模型中的常數通過大量、廣泛的實驗數據進行了精細標定，是一個強調應用的實用化模型。突出的優點是魯棒性非常好、計算效率高、簡單實用，是航空航天領域應用最廣的模型之一，尤其在航空領域受到普遍歡迎，對附著邊界層等較簡單的流動效果很好。例如在上述的DPW、HiLiftPW等活動中，SA的使用者數量遠超其他模型。但是，一般說來，SA模型的普適性、精度(如渦黏性有時偏大)、對分離等復雜流動的模擬能力等有待進一步提高。

2)k-ε模型：這是RANS 歷史上第一個實用化的、完備的兩方程模型，歷史悠久、發展和演化的類型多、使用者眾、應用廣泛，發展了很多種針對各種情況的修正方法等，但該模型的邊界條件處理困難、數值剛性大、穩定性較差，對逆壓梯度、分離等模擬效果欠佳，需要普適性差、經驗性強的阻尼函數或橋函數等。尤其是在一些較復雜的邊界層中，因為ε的擴散模型性能不佳，會高估湍流長度尺度等，這將導致分離抑制或延遲、回流區短、表面摩擦系數和熱流過大等難以克服的缺陷。相比SA、SST等模型，現在的趨勢是問題較多、應用偏少，在最近的一些大型湍流模型的性能評估中，甚至沒有進入評估名單，在DPW、HiLiftPW等系列活動中，使用者也越來越少。

3)k-ω模型：該模型的發展歷史悠久，種類很多，這其中應用最成功、最受歡迎的是不同時期、不同版本的Wilcoxk-ω兩方程模型。該模型數值穩定性好，壁面處可以使用Dirichlet型邊界條件，處理難度低、精度高、數值剛性小；因為不需要計算到壁面的最小距離，可以直接積分到壁面，不需要顯式壁面衰減函數，因此便于應用在復雜外形流動中；對逆壓梯度、分離等問題的模擬性能良好。因此，總體上說來是一個性能優秀、特點鮮明的模型。但其缺點也較明顯：對ω的自由來流邊界條件過分敏感、湍流近壁漸近特性欠佳等。需要說明的是，近些年該模型的不斷完善也在努力改進這些缺點。

4)SSTk-ω模型：剪切應力輸運(Shear-Stress-Transport，SST) 模型，是目前廣泛使用的EVM模型中評價性能最好一個，無論是模型研究者還是有經驗的實戰使用者對它都比較推崇。SST的基本原理是：在近壁處采用k-ω模型、在邊界層邊緣和自由剪切層采用k-ε模型，并通過Bradshaw假設(剪應力正比于湍動能τ=ρa1k)引入了雷諾剪切應力輸運的影響。因此它巧妙結合了k-ω、k-ε和JK模型的優點。除了具有上述k-ω的穩定性好、精度高等優點外，還可以較好地處理湍流剪切應力在逆壓梯度和分離邊界層內的輸運，故SST模型能更好地預測逆壓梯度和邊界層分離等較復雜的流動情況。因此，SST模型以其精度高、魯棒性好、適用性較廣而成為航空航天領域應用最廣的湍流模型之一。

5) 雷諾應力模型RSM：目前CFD界對RSM認識和使用者都很少，但不少人看好RSM的應用前景，甚至認為是未來的主力模型，因此這里較詳細地討論一下。

RSM通過求解由N-S方程嚴格推導來的、完整的傳輸方程，從而提供額外的湍流動量通量，在RANS框架下是對湍流最全面的描述。因此被認為是最自然、最符合邏輯的模型，它可以克服渦黏性模型的一些嚴重缺陷，如各向同性等。但是長期以來的評估分析認為：雖然一般情況下RSM模型比EVM模型具有更高的精度、更廣的普適性，尤其是對曲率、旋轉、分離等復雜流動情況，但這種優勢并不是全面的、明確的。

RSM的主要問題在于：1) 雷諾應力方程的建模問題，尤其是壓力應變項的建模復雜、困難；2) 數值剛性問題較嚴重；3) 計算量較大。但近十多年來，RSM模型在建模和計算方法等方面都取得了較大的進展：模化精度進步較大，數值方法進步明顯：計算量已不是大問題，魯棒性得到改善。例如：使用RSM模型，成功地計算出幾個復雜的渦輪機械流動，所用的計算時間與k-ε模型相比僅多了30%[13]。

遺憾的是，至今愿意使用RSM模型的人還是不多，分析原因如下：

1) 一般認為，相比EVM，雖然RSM的性能有改善，尤其對復雜流動，但并沒有顯示出普遍性、壓倒性的優勢；

2) CFD界普遍對RSM模型的進展認識不足，習慣用過去的眼光看待今天的RSM；

3) 使用者的習慣或者慣性，也就是說不太愿意接觸新的、更加復雜難用的模型；

4) 軟件開發和推廣困難，因為RSM雖然性能上有優勢，但仍然不如EVM使用方便，也還可能包括商業上利益驅動不足。

長期以來，美國、歐洲等開展了多次針對湍流模型的專題研究、性能評估、分析研討等，近些年來，尤其關注EVM和RSM的對比分析和評估。這是因為一方面EVM在目前的CFD占主導地位、而RSM有希望成為未來的主流；另一方面，RSM在建模方法、數值剛性、計算效率等方面取得了較大的進步，而EVM的表現不盡如人意。

為了對EVM和RSM的性能有一個直觀的認識，也為了介紹上述的性能評估，作為示例，這里給出一個比較嚴謹的湍流模型性能研究。

文獻[14]使用多個經典標模，系統深入地對比研究了4個經典湍流模型：EVM(SA、SST)和RSM(SSG /LRR-ω、JHh-v2)。這4個湍流模型是目前學術界公認的、具有代表性的優秀湍流模型。計算結果分析表明：對于簡單標模如小壓力梯度平板邊界層，4個模型計算結果非常一致；隨著標模復雜性的增加，一般情況下，RSM表現出了有更大的適用性和更高的精度，但這種結果的改善常常不是很明顯，也不是很普遍。

這里較詳細地討論一個計算結果改善明顯的標模算例：ONERA M6機翼。這是一個經典的湍流模型驗證標模，常用來檢驗湍流模型的跨聲速性能，尤其激波誘導分離的模擬水平。使用DLR TAU軟件計算，流動參數：Ma=0.84，Re=11.72×106，0.03°≤α≤ 6.06°。

α=3.06°之前的計算結果表明，不同湍流模型的預測之間幾乎沒有任何差異。但α= 4.08°時，在機翼外側，RSM模型的結果部分同EVM模型結果大相徑庭。圖5顯示了兩個機翼外側展向截面站位分別為80%和95%的壓力分布。可以看到，SA和SST模型的預測與實驗結果明顯不同，而SSG/LRR-ω模型所預測的壓力分布與測量結果都非常吻合，JHh-v2 RSM預測的激波位置在80%展向截面處同實驗更加吻合、但在95%展向截面時偏差明顯。圖6顯示了機翼上表面的摩擦線圖，可以看到，SST模型預測的機翼外半部有很大的分離，而SSG/LRR-ω模型預測的分離僅限于激波根部后面的小區域，根據壓力分布，這個結果與實驗吻合得更好，JHh-v2 模型預測的分離流展向延伸得更小。

圖5 ONERA M6不同展向站位的壓力分布[14]Fig.5 Pressure distribution at different spanwise locations of ONERA M6[14]

圖6 ONERA M6機翼上表面摩擦線和壓力分布(α= 4.08°)SST模型(左)，SSG /LRR-ω模型(中)和JHh-v2模型(右)[14]

進一步的研究還發現了另一個不希望的結果：除了SSG/LRR-ω外，其他3個模型顯示出了對網格質量的高度敏感性，只有SSG/LRR-ω模型在兩套網格上都與實驗結果一致吻合。在上述計算網格的基礎上，保持網格總量基本不變、重點在激波附近加密網格，結果4個湍流模型都給出了非常好的結果，如圖7，對比圖5，網格效應明顯，顯然這個結果不是使用者希望的。

圖7 ONERA M6 η=0.8展向站位的壓力分布[14]Fig.7 Pressure distribution at the spanwise location of η=0.8 for ONERA M6[14]

其他較復雜的經典標模的計算結果表明：在翼尖渦旋、彎管流動等算例中，RSM的優勢更明顯，這是因為線性渦黏性模型無法考慮流線曲率等效應。

1.2 轉捩模型

轉捩是流動從分層穩定的層流狀態演化為混沌紊亂的湍流狀態的連續物理過程。由于內在演化機制高度復雜，轉捩流動被認為是經典物理學遺留的最具有挑戰性的問題之一，具有重要的理論研究意義。同時，轉捩現象與社會經濟發展和科技進步的重大科學問題密切相關。國務院2018年印發的《“十三五”國家科技創新規劃》部署了面向2030年的一批重大科技項目，其中的航空發動機及燃氣輪機項目重點攻關的多類復雜流動與轉捩密不可分(如圖8所示)。以渦輪葉片為例，受自由來流湍流度、來流雷諾數、流動分離和非定常尾流等多類因素影響，渦輪葉片的流動大多涉及邊界層轉捩現象，包括但不限于旁路轉捩、分離誘導轉捩、逆轉捩(再層流化)以及上游葉片尾流導致的周期性轉捩等。早期研究認為尾流誘導的轉捩會增加渦輪葉片流動損失，但隨著認知深入，人們發現上游尾流效應還可以被用來控制葉片吸力側的邊界層轉捩過程，從而抑制流動分離，在不犧牲效率的情況下提升低壓渦輪葉片的升力、節約成本。由此可見，深入開展轉捩研究、準確預測轉捩對現代渦輪設計和性能評估有重要指導意義。此外，邊界層轉捩還將顯著影響高超聲速飛行器的壁面摩擦阻力、氣動熱環境、流動分離位置和進氣道流場品質等，是發展高超聲速飛行器技術必須攻克的瓶頸問題。

圖8 高涵道比渦扇發動機中的湍流/轉捩現象[15]Fig.8 Modern bypass engine section and examples of CFD simulation zones considered[15]

目前針對邊界層轉捩/湍流的數值模擬研究方法主要包括DNS、LES、轉捩準則、穩定性理論和基于RANS的湍流/轉捩模型方法。這其中，基于RANS的轉捩模型方法由于優越的經濟性和魯棒性，是現有條件下最有可能應用于工程復雜外形邊界層轉捩預測的數值模擬方法。

經過二十余年的快速發展，轉捩模型方法的研究取得了長足的進步，截止到目前，文獻公開發表的基于RANS方程的轉捩模型就已多達幾十種，其宗旨都是合理預測轉捩起始位置(即轉捩判據的建立)和轉捩區長度。同時，轉捩模型的發展伴隨著人們對轉捩物理機理認知的日益豐富也在不斷進步。這方面里程碑式的工作是轉捩流動間歇性和層流脈動(或稱非湍流脈動)的發現。這兩個轉捩流動關鍵特征的提出，為轉捩模型的發展和改進提供了重要的理論依據和指導價值，由此發展出一系列面向工程應用背景的轉捩數值預測手段。根據所構造的轉捩模型對這兩個重要概念的模化程度，本文將現有渦黏性RANS轉捩模型分為低雷諾數湍流模型、基于間歇因子的轉捩模型和基于層流動能的轉捩模型三類。

1.2.1 低雷諾數湍流模型

低雷諾數湍流模型既沒有考慮轉捩的間歇現象，也沒有描述層流動能，其核心思想是通過阻尼函數限制黏性子層湍動能的生成，從而在一定程度上達到模擬邊界層轉捩的效果。比較著名的工作是Wilcox的研究，他在k-ω湍流模型的封閉參數和渦黏性中引入了一個關于湍流雷諾數的人工阻尼函數[16]。湍流雷諾數表達式為ReT=ρk/ωμ，通過變換可表示為渦黏性系數μt和分子黏性系數μ的黏性比ReT=μt/μ。流動為層流時，湍流雷諾數遠小于1(ReT?1)。當流動接近轉捩起始位置時ReT近似為1，從而導致湍動能k迅速增長。隨著k增長，渦黏性和壁面摩阻開始增大，ReT隨之變大，一旦ReT?1湍動能比耗散率ω開展增長，直到k方程的生成項和耗散項達到平衡，此時層流到湍流的轉捩過程結束，流動轉捩為湍流。基于ReT構造的低雷諾數湍流模型還有很多，例如Jones和Launders的k-ε模型[17]，Hadzic和Hanjalic的SMC模型[18]等。Langtry和Sjolander的模型也可歸為此類，不同的是，其用渦雷諾數ReV代替ReT作為轉捩起始位置的判定因子[19]。通過上述描述可發現，低雷諾數湍流模型對轉捩現象的模擬純粹由模型的數值特性引起，并無物理背景。Zheng等[20]的研究發現，Wilcox的k-ω低雷諾數模型預測的轉捩起始位置過于提前，且不能正確預測分離誘導轉捩。Rumsey[21]通過深入分析認為低雷諾數模型對轉捩的預測只是一種巧合，并將其稱之為“偽轉捩”。越來越多的研究指出“只應用湍流模型而不考慮轉捩現象的間歇性，對轉捩過程的模擬是脆弱且不可靠的”[22]。因此，盡管低雷諾數湍流模型作為早期的轉捩研究工具，對轉捩模型的發展意義重大，但建立基于轉捩物理機理的轉捩模型刻不容緩。

1.2.2 間歇因子轉捩模型

這類模型考慮了轉捩物理過程中的間歇現象，從而將層流、轉捩和湍流等不同發展階段的求解解耦。轉捩區域流動的間歇現象最早由Emmon[23]發現，自此轉捩問題的理論研究有了實質性突破。所謂間歇現象是指在轉捩過程中，同一空間位置層、湍流交替出現的現象，通常用間歇因子γ來描述。對于γ的求解，目前主要有兩種方法：第一種方法也是最簡單直接的方法，是根據實驗數據擬合出代數經驗公式；第二種方法是根據擬合公式構造γ輸運方程。

首先介紹第一種方法，Dhawan和Narasimha[24]最早根據平板邊界層的實驗數據給出了γ沿流向分布的經驗擬合式：

(2)

其中，n代表展向單位距離的湍斑生成速率，σ表示湍斑擴散參數，xt為轉捩起始點。上式可無量綱化為以下兩種形式：

(3)

γ(y)=0.5[1-erf(ζ)]，

ζ=5(y/δ*-0.78)/8

(4)

其中，δ*為邊界層位移厚度。Addison和Hodson[28-29]提出了Prescribed Unsteady Intermittency Model (PUIM)方法以使γ在流向和法向的分布均吻合實驗數據。PUIM中γ是一個關于壁面Emmon斑(湍斑)發展特征的積分表達式，其只對自然轉捩和旁路轉捩有效。此外，基于代數經驗公式的方法還有Steelant和Dick[30]提出的γ模型，Fürst等[31]建立的模型等，在此不一一贅述。這種通過代數經驗擬合式求解γ的方法應用范圍有限，其嚴重依賴于來流和外形條件，無法反映轉捩區間歇現象的非當地效應和歷史效應。

為克服這一缺陷，人們嘗試引入γ輸運方程來求解三維流動中的γ分布規律。Libby[32-33]在1975年構造了間歇因子輸運方程，Dopazo[34]、Chevray和Tutu[35]在其基礎上對輸運方程進行了改進和發展。Cho和Chung[36]將γ輸運方程與k-ε湍流模型相耦合，開創性地提出了k-ε-γ轉捩模型，并將其成功應用于平面射流、圓射流和尾流等剪切流的轉捩預測。以上輸運方程均針對自由剪切流發展而來，不適用于邊界層流動中的轉捩模擬。Savill[37]在Cho和Chung[36]模型的基礎上發展了適用于邊界層轉捩的γ輸運方程。Steelant和Dick則[38]借鑒Dhawan和Narasimha經驗擬合式[24]構造了新的γ輸運方程。比較代表性的工作還包括Suzen和Huang[39]，Suzen等[40]和Pecnik等[41]發展的邊界層流動γ輸運方程，這些模型在特定三維流場計算中都獲得了滿意的結果。

然而，間歇因子只描述轉捩過程，為實現流動的轉捩預測，往往還需要一個與之分離的轉捩起始位置判據。轉捩起始位置受眾多因素影響，加劇了預測難度[42]。目前轉捩起始位置判據通常基于經驗公式發展而來，以動量厚度雷諾數Reθ最為常見。首先積分求解層流流場中每個流向位置Reθ值，一旦判定當地Reθ超過經驗確定的臨界值，即可認為轉捩發生。然而需要積分求解的非當地變量的引入，限制了這類轉捩模型在非結構并行求解器中的應用。為此，閻超等[43-44]提出了適用于大規模并行計算的網格重排序循環盒子法來識別邊界層參數。Menter和Langtry等[45-46]通過使用應變率為底的雷諾數Rev代替動量厚度雷諾數Reθ，創新性地提供了一個將轉捩判據進行完全當地化處理的框架，進而提出了γ-Reθ轉捩模型。該工作具有十分重要的借鑒意義，在學術研究和工程應用方面受到了廣泛關注。大量算例驗證表明，γ-Reθ模型在渦輪機、風力發電機以及其他工業應用中都可得到滿意的預測結果[46]。圖9給出了γ-Reθ模型計算的T106渦輪葉片在不同尾跡周期的湍動能分布與實驗值的對比，可見計算與實驗結果吻合良好。但需要指出的是，γ-Reθ模型的構造完全依賴于低速流動實驗數據的經驗擬合，并不反映邊界層轉捩過程中的物理機制，導致適用范圍有限，尤其在應用于擬合數據適用范圍之外的其他轉捩問題時，如高速邊界層轉捩，其模型構造往往需要進行重新擬合與標定。郝子輝和閻超等[47]構造了適用于高超聲速流動的渦雷諾數Rev與邊界層動量厚度雷諾數Reθ的擬合關系，實現了γ-Reθt轉捩模型的高超聲速修正。

圖9 T106渦輪葉片不同尾跡周期(t/T)湍動能分布云圖(左:實驗，右:γ-Reθ模型)[46]Fig.9 Contours of turbulent kinetic energy at various wake passing periods (t/T) of T106 turbine blade cascade (left: Experimental data, right: γ-Reθ model)[46]

總體而言，基于間歇因子的轉捩模型可較為合理地描述轉捩過程，但通常需要額外構造一個轉捩起始位置判據，且目前該判據的建立高度依賴于實驗數據的經驗擬合，而對轉捩前區域的物理機制考慮不足，極大地限制了該類轉捩模型的適用流動范圍。

1.2.3 層流動能轉捩模型

隨著對邊界層轉捩機理認知的不斷深入，人們發現通過對轉捩前層流動能的刻畫即可實現轉捩起始位置的自動判斷，為此發展出一系列層流動能轉捩模型。顧名思義，這類模型考慮了轉捩物理過程中層流動能的發展，當邊界層轉捩流動前的層流動能發展到一定程度時判定轉捩發生，因此又被稱為基于物理機理的轉捩模型(Physics-Based Transition Models)。

實驗研究發現，轉捩起始點前的層流區域中的流動并不是穩定的，其存在類似湍流脈動的結構，Mayle和Schulz[51]將其稱為“層流脈動”或“非湍流脈動”。層流脈動與湍流脈動有著顯著的區別，具體表現在以下三個方面：首先，從能量分布的角度來看，層流脈動的能量幾乎全部包含在流向分量上；其次，從動力學特性的角度來看，層流脈動不存在經典湍流脈動中能量由大尺度向小尺度輸運的能量級串過程；最后，層流脈動頻率相對較低，除了壁面附近，其耗散也較低。層流脈動結構的發現對轉捩預測極具吸引力，并由此開辟了一個新的轉捩模型研究思路。就目前而言，邊界層層流脈動發展的物理機制還未被完全揭示，不過主流思想認為：邊界層對自由來流中特定尺度渦的選擇性以及邊界層內低頻擾動被流動剪切作用放大的過程對層流脈動的發展影響顯著。Bradshaw[52]提出了Splat機制來揭示層流脈動的演化發展過程。Splat機制認為，壁面的約束作用強迫邊界層內法向脈動方向發生改變，導致當地局部壓力梯度的出現，進而誘導層流脈動的產生和發展。盡管如此，層流脈動的發展規律和頻率選擇性仍尚未被完全揭示。

目前層流脈動的主流計算方法包括兩種，一是建立層流脈動動能kL的輸運方程，二是比擬湍流渦黏性系數μt對層流脈動黏性系數μnt進行模化。

層流脈動動能kL輸運方程的建模關鍵是層流脈動的產生和能量由層流脈動向湍流脈動的輸運。Mayle和Schulz[51]首次通過類比湍動能kT的輸運方程構造了kL的輸運方程，不同的是將kT方程產生項的應力/應變替換成了壓力脈動。這也反映了層流脈動不同于湍流脈動的產生機制。Walters等基于Splat機制，將邊界層內的脈動能量分為有壁面約束的層流動能kL和無壁面約束的湍動能動能kT，由此構建了kL輸運方程，并分別基于k-ε和k-ω湍流模型發展了kT-kL-ε[53]和kT-kL-ω[54]轉捩模型，成功實現了對轉捩起始點的自動判斷。秦玉培和閻超等將kT-kL-ω轉捩模型拓展至高超聲速邊界層的自然轉捩[55]預測，并針對原始模型對轉捩峰值位置“overshoot”現象預測不足的問題，通過分析湍動能擴散時間尺度與湍動能輸運時間尺度之間的關系構造了新的代數型間歇因子，成功實現了對 “overshoot”現象[56]的準確模擬(如圖10所示)。劉再接和閻超等[57]進一步拓展了kT-kL-ω模型對粗糙單元誘導轉捩現象的預測功能。陳燦平等[58]通過引入分離敏感參數修正了小尺度渦黏性系數，改進了層流動能轉捩模型對于流動分離區域數值預測的準確性。

圖10 超聲速尖錐壁面熱流分布[56]Fig.10 Surface heat flux distribution for hypersonic straight cone[56]

第二種描述層流脈動的方法是直接模化層流動能黏性系數μnt，該思路最初由Warren等[59]提出。他們以k-ζ湍流模型為基礎，將其中的渦黏性系數μt替換為有效渦黏性系數μeff，從而在黏性系數中顯示層流脈動的影響。μeff的表達式為μeff=(1-γ)μnt+γμt，其中μnt的構造利用了穩定性分析的模態理論結果。國內王亮和符松[60-61]借鑒Warren等[59]工作，發展了完全基于當地變量的k-ω-γ轉捩模型。該模型的特色之處是對超/高超聲速邊界層轉捩現象的預測。趙雅甜和閻超等[62-63]通過對比評估k-ω-γ和γ-Reθ、kT-kL-ω轉捩模式在高超聲速充氣式柔性減速器外形轉捩現象預測中的性能發現：γ-Reθ和kT-kL-ω模式預測的轉捩位置嚴重偏離實驗值，而k-ω-γ模型計算得到的轉捩起始位置和轉捩陣面均與實驗值吻合良好(如圖11所示)。周玲和閻超等進一步對k-ω-γ轉捩模型進行了雷諾數效應修正[64]、頭部鈍度修正[64]和橫流拓展[65]，并拓展應用于復雜工程外形的轉捩預測，取得了較為滿意的結果。

圖11 k-ω-γ模型(CFD)計算的高超聲速充氣式柔性減速器壁面熱流分布云圖與實驗值(Exp.)對比[62]

總體而言，層流動能轉捩模型較為充分地考慮了轉捩過程的物理機理，顯示出了較大的應用潛質和優勢。然而，由于目前層流動能的機制尚未被完全揭示，再加上多種不穩定擾動模態相互干擾的復雜性，現有基于層流動能的轉捩模型依然難以解析地給出三維邊界層流動的轉捩起始位置和轉捩過程。同時，無論是kL輸運方程的建立，還是層流動能黏性μnt的構造，都不可避免地增加了轉捩模型中經驗參數的數量，從而導致預測結果的更大不確定性。針對這一現象，趙雅甜和閻超等[66-67]、張金程和符松[68]開展了相應的不確定度量化研究工作。

本文根據轉捩模型構建思路和構造方式的差異將轉捩模型分為三類，但需要注意的是各類方法間并沒有明確清晰的界限。例如，王亮和符松提出的k-ω-γ模型[60-61]既考慮了轉捩前的層流脈動動能，又建立了間歇因子輸運方程，同時該模型以SST湍流模型為基礎框架，因此還會受到低雷諾數模型中阻尼函數的影響。此外近年來還有一些新的方法不斷涌現，他們并不屬于這三類方法的范疇。例如，依據層湍流發展過程中的動力學近似發展的三維唯象V-model模型[69]以及基于機器學習的數據驅動湍流/轉捩建模方法[70-71]等。這些新方法的提出豐富了轉捩的模擬手段和求解思路，具有重要的研究價值。

2 計算方法

在航空航天領域，CFD是通過數值求解流體動力學方程，來獲取各種流動條件下的流場數據以及作用在繞流物體上的氣動載荷的。因此，直接針對控制方程進行離散化求解的計算格式，如通量格式和高階格式，就一直是學界研究的熱點。它們通過模擬離散化控制方程所表征的流動演化，起著反映流動物理性質的作用。而CFD領域中的其他關鍵技術，如湍流模擬、轉捩預測等方面的研究，都需要建立在方程物理性質能夠被準確刻畫的基礎之上。所以，計算方法也一直被視作為CFD的基石之一。而NASA發表的CFD2030報告中也明確指出：計算格式是未來需要大力發展的關鍵技術之一。為此，本節將針對航空航天中通量格式和高階格式的研究進展分別進行評述。

2.1 通量求解方法研究進展

通量格式通常被大致分為兩類：中心類通量格式和迎風類通量格式[72]。而較之于中心類通量格式，迎風通量格式是在精確Riemann解的基礎上，通過近似求解Riemann問題，考慮了實際流動過程中波的傳播方向、物理意義更加明確。此外，這類方法還具有計算效率和計算精度的綜合優勢，因此成為了當前CFD中應用最為廣泛的計算方法。例如，工業界廣泛使用的商業或in-house軟件，如FLUENT、CFD++、CFL3D、FUN3D、OVERFLOW、BCFD等，均采用迎風類通量格式作為自己的核心算法。所以，本文將著重針對迎風類通量格式展開討論。

目前，根據各自構造思想的不同，迎風類通量格式被分為三類。

第一類是基于近似求解Riemann問題的通量差分分裂類格式(Flux Difference Splitting, FDS)。其中，Roe格式因為其對接觸間斷與激波的高分辨率而應用最廣。但是在模擬膨脹波時，Roe格式會因為違反熵條件而產生非物理解，需要引入熵修正。熵修正需要經驗性人工參數的設置，不同的參數設置對流場計算精度的影響較大[73]。

第二類是矢通量分裂格式(Flux Vector Spliting, FVS)。它們基于特征值符號的信息對通量項進行分裂和定向離散，比較有代表性的是Steger and Warming格式和van Leer格式。這類方法有激波捕捉能力強、可靠性高，以及計算效率高等優點。但是該類格式無法精確捕捉接觸間斷，黏性區分辨率較差。

第三類是混合的AUSM(Advection Upstream Splitting Method)類格式。這類格式通過對壓強通量和對流通量分別進行相應的通量分裂，結合了FVS類格式的高效率以及FDS類格式的高間斷分辨率，并且克服了二者的缺點。其中，基于對流通量或者壓強通量構造方式的不同，這類方法包括AUSM+、AUSMDV、AUSMPW+等一系列格式[74]。此外，還有一些方法采用了其它不同的通量分離模式，例如Zha-Bilgen通量分裂方法以及Toro通量分裂方法等[75-76]。

從上述三類迎風類通量格式的簡要介紹可以看出，作為CFD領域所圍繞的發展主線之一，通量格式無論在理論研究方面還是在工程應用方面，近幾十年來都受到了廣泛的關注并取得了一定的進展。然而正如本文引言部分所述，作為航空航天領域所廣泛使用的通量計算方法，它們卻一直沒有取得太多本質上的進步，如下一些不足致使它們距離“完美”還非常遙遠：

1) 在模擬強激波時，容易出現激波異常現象；

2) 在計算膨脹流動時容易出現非物理的“overheating”解；

3) 在模擬橢圓型的低速不可壓流動時，會出現精度下降、收斂困難等問題，無法適用于全速域流動的高精度數值模擬；

4) 在針對流動分離、激波/邊界層干擾這類多維復雜流動進行數值模擬時，會出現精度降低和計算穩定性下降等問題。

下面我們將針對這四個問題，就近年來國內外學者所取得的研究成果以及遇到的挑戰進行簡要概述。

2.1.1 激波異常問題

激波異常現象是指CFD在模擬高馬赫數流動中的強激波問題時，容易因為出現非物理的激波不穩定解而得到錯誤的氣動力/熱載荷。1988年Perry等在使用Roe格式計算高超聲速圓柱繞流問題時首先發現了該問題，并將其稱之為“紅玉”現象(圖12)[77]。在此之后，大量學者針對該問題開展了深入的研究，其中比較有代表性的工作可以參見文獻[78,79]。但是，2010年Roe等通過計算一維/二維靜駐激波問題指出，即使是公認的激波捕捉最為穩定的van Leer格式，在平行于激波方向的網格尺度過小時，也會出現非物理的激波異常現象(圖13)[80]。也就是說，現有的所有通量格式都無法完全避免激波異常現象的出現，而激波處的網格長細比顯著影響著它們的激波穩定性。

圖12 激波異常現象Fig.12 The shock anomaly phenomenon

圖13 1.5D靜駐激波問題中van Leer’s FVS格式在不同網格下計算得到的密度分布圖[80]Fig.13 Density contours of the 1.5D steady shock problem obtained by van Leer’s FVS[80]

2.1.2 “overheating”問題

2011年，Liou在 “Open Problems in Numerical Fluxes: Proposed Resolutions”的文章中首次明確指出，現有通量格式在計算膨脹流動或對撞流動時都會出現非物理的溫度跳躍，即 “overheating”現象(圖14)。而加密網格、增加時間迭代步數或者采用高階重構方法均無法有效解決該問題[81]。針對這一不足，Liou開展了詳細的理論分析并發現，造成“overheating”現象出現的直接原因是當前通量格式在計算無黏流動時無法保證熵守恒。基于該結論，他提出了一種熵守恒的通量計算方法。但遺憾的是，這種方法需要與時間推進強耦合，不適用于現有的CFD軟件體系，因此在航空航天領域沒有得到推廣應用。

圖14 膨脹流動模擬中的“overheating”現象[81]Fig.14 The ‘overheating’ phenomenon in simulating the receding flow problem [81]

2.1.3 全速域流動模擬問題

特征分析表明，描述無黏流動的Euler方程在不同速域中的數學特性差別較大。例如，描述低速近不可壓流動時，Euler方程為橢圓型；描述超聲速流動時，Euler方程為雙曲型。理論上針對不同特性的微分方程組應采用不同的空間離散方法[82]。因此，上述基于可壓縮Euler方程求解進行構造的通量格式難以實現對低速流動的高精度模擬。甚至在很多情況下，即使加密網格或者采用高階重構方法，它們也依舊無法給出令人滿意的結果(圖15)[83]。

圖15 低速無黏NACA0012翼型等密度分布圖(Ma=0.05)[83]Fig.15 Density contours of inviscid flow over the NACA0012 airfoil[83]

針對該問題，Guillard等通過對可壓縮Euler方程開展低速漸近性分析發現，在低速近不可壓流動中，流場中壓強擾動的幅值與馬赫數滿足如下的關系式：

(5)

而傳統的通量格式在模擬低速流動時得到的流場解無法滿足式(5)。因此，它們難以高精度模擬低速流動[84]。為了改善這一缺陷，Weiss與Turkel等在特征分析可壓縮Euler方程的基礎上，通過引入預處理矩陣，有效改善了離散化可壓縮Euler方程在低速近不可壓流動中的剛性問題，提高了通量格式在低速近不可壓流動區域的模擬精度[85]。然而，這類預處理方法存在如下不足：自由來流條件為超聲速時，預處理方法完全失效，如果此時流場存在低速流動區域，該區域的計算精度難以保證；需要設置經驗性人工參數，這種參數應該根據具體的流場變化，難以找到普適的量化準則。針對該問題，李雪松對Roe格式開展了低速漸近性分析。通過量綱分析他發現，Roe格式中壓強耗散項系數過大導致其難以高精度求解低速流動。基于該結論，李雪松將無量綱化離散方程中壓強耗散項的量級修正至O(M)，改進了Roe格式的低速模擬精度[86]。Thornber、Rieper、FillionShima和Kitamura等也開展了類似的工作[87]。

但是，這些研究并沒有關注格式在捕捉強激波時所遇到的穩定性問題。為此，筆者借鑒李雪松等的工作，針對Roe格式同時開展了高速的激波線性小擾動分析及其低速的漸近性分析，揭示了通量格式中各耗散項系數在低速流動模擬以及激波捕捉中應滿足的構造準則是相互補充而非互相矛盾的。基于上述結論，筆者提出了一種適用于全速域流動模擬的RoeMAS格式[88]。此后，Tsoutsanis、Chen、Xie等也開展了類似的工作[89-91]。

下面展示幾個全速域通量格式在航空航天中的應用實例。

圖16和圖17分別對比了低速NACA0012和NACA0042翼型有黏繞流算例中不同格式計算得到的壁面壓強系數分布[92]。其中，AUSMPWM格式是基于AUSMPW+格式，通過抑制低馬赫數情況下的壓強耗散而構造的一種全速域通量格式[93]。可以看到，在低速流動模擬中，全速域通量格式的模擬精度明顯優于傳統通量格式。

圖16 低速有黏NACA0012翼型壁面壓強系數分布對比(Ma=0.15, Re∞=6×106 )[92]

圖17 低速有黏NACA4412翼型壁面壓強系數分布對比(Ma=0.09, Re∞=1.52×106 )[92]

圖18展示了高超聲速鈍錐繞流問題中，不同格式計算得到的駐點熱流值同網格雷諾數的對應關系[94]。從中可以看出，采用全速域通量格式可顯著改善高超聲速氣動熱預測精度對計算網格的依賴性，有利于降低使用者的技術門檻。但是遺憾的是，針對全速域通量格式在知名商業軟件或者in-house工程CFD代碼中的應用推廣，尚沒有相關公開報道。

圖18 高超聲速鈍錐擾流駐點熱流值隨網格雷諾數變化曲線對比[94]Fig.18 Schemes’ stagnation heating values vs. the cell Reynolds number[94]

2.1.4 多維流動模擬問題

2017年Abgrall和Roe在HandbookofNumericalAnalysis一書中先后指出，基于維度分裂思想進行構造的傳統通量格式，在針對流動分離、激波/邊界層干擾這類多維復雜流動進行數值模擬時，僅考慮了沿著網格界面法向傳播的信息而忽略了沿著網格界面橫向傳播的信息，因此無法準確描述空間中不同維度之間的相互影響，容易出現精度降低和計算穩定性下降等問題。為了克服該缺陷，諸多學者開展了關于多維通量格式的相關研究。其中比較有代表性的工作包括文獻[95-97]。但是，這些多維通量格式構造復雜，或針對定常流動求解只能顯式推進、無法結合隱式時間推進方法實現求解加速，所以尚未在復雜流動模擬中得到廣泛應用。為了解決這個問題，Balsara通過在網格角點處構造二維干擾波系，提出了一種形式簡單的多維HLL通量求解方法[98]。圖19所展示的徑向沖擊波算例的計算結果表明，該方法在求解多維波傳播時，可以基本實現流場解與網格分布的相互解耦[99]。基于該工作，Mandal 等在Zha-Bilgen通量分裂思想的基礎上，將流動控制方程中的無黏通量項分解成壓強通量項及對流通量項，并通過直接采用Balsara所構造的多維干擾波模型，實現了可精確捕捉線性接觸間斷的多維GM-HLL-CPS-Z通量求解方法[100]。袁禮等基于TV格式的分裂思想，也提出了類似的多維HLL通量求解方法[101]。筆者基于坐標變化，將上述工作推廣至一般曲線坐標系，提出了一種適用于復雜飛行器流動模擬的多維通量求解方法MULTV格式[102]。高超聲速二維雙橢球繞流算例的計算結果驗證了該方法在多維復雜流場模擬中的精度優勢(圖20、圖21)。然而，當前關于真正多維Riemann通量求解方法的研究仍比較初步，針對三維、寬速域流動問題的算法研究較為匱乏。

圖19 徑向沖擊波問題等密度分布圖[99]Fig.19 Density contours of the spherical blast wave case[99]

圖20 高超聲速二維雙橢球等馬赫數云圖[101]Fig.20 Mach contours of the hypersonic viscous flow over the two-dimensional double-ellipsoid[101]

圖21 不同格式得到的壁面斯坦頓數分布(高超聲速二維雙橢球算例)[102]Fig.21 Distributions of the schemes’ Stanton numbers along the meridians of the windward side[102]

2.2 高階格式

假設網格尺度為h、數值解的誤差為e，如果e∝hk成立，則稱數值格式為k階精度。一般將離散精度高于二階(即k≥3)的計算格式統稱為高階格式[103]。相比于低階格式，高階格式被認為在花費較低計算代價的前提下具備獲得更高精度的巨大潛力，因而引起了CFD研究者和用戶的廣泛關注。國外方面，專注于高階格式的研討會“The International Workshop on High-Order CFD Methods”[104]已分別于2012年、2013年、2015年、2016年和2018年舉辦了五屆。國內方面，已連續舉辦了兩屆(2018年和2019年)“計算流體力學中高精度數值方法及應用”(Trends of High Order Numerical Methods for Computational Fluid Dynamics and Their Applications)研討會。

CFD發展到今天，出現了多種高階格式，主要包括：

1) 有限差分類高階格式： WENO(weighted essentially non-oscillatory)格式[105-106]、TENO(targeted ENO)格式[107]、緊致格式[108]、保單調格式[109]、WCNS格式[110]、ENN格式[111]等。這類方法最突出的優點是高精度、高效率，編程實現相對直觀，較多應用于外形相對簡單的湍流大渦模擬和直接數值模擬。然而，有限差分類格式只適用于結構網格，并且其數值解的質量與度量系數的計算精度相關，因此要求網格要足夠光滑。這一要求很難在復雜構型的網格生成中得到滿足，也限制了此類方法向實際工程問題的推廣。

2)有限體積類高階格式：有限體積WENO格式[112]、高階緊致重構格式[113]等。有限體積格式適用于非結構網格，因此，對于工業級CFD軟件，二階精度的有限體積格式是最常用的數值離散框架。然而有限體積法需要擴展重構模板來實現高階精度和增強穩定性，而由于非結構網格的無序性，使得高階有限體積法無論是算法設計還是程序編寫都較為繁瑣。同時，擴張的模板也給計算效率帶來較為嚴重的負面影響。為應對上述問題，Wang等[113]通過引入相鄰單元導數的信息來有效控制模板的范圍。

3)有限元類高階格式：SD (Spectral Difference)[114]、SV (Spectral Volume)[115]、DG(Discontinuous Galerkin)[116]、FR(Flux Reconstruction)[117]、SUPG (Streamline upwind/Petrov-Galerkin)[118]等。這類格式與有限差分和有限體積類格式的最大差別是：有限差分和有限體積類格式通過擴張插值或重構模板，引入周圍單元更多的信息來提高精度，而有限元類高階格式則通過在一個網格單元內引入更多自由度來提高精度，而無需擴展模板。這一特征決定了有限元類格式可以較容易地在非結構網格上實現任意高階精度，也是近年來獲得較多關注的高階格式。有限元類高階格式在帶來高精度的同時，也在網格生成、時間推進以及激波捕捉等方面迎來了挑戰。本節后面將詳細論述。

4) 其他衍生類型高階格式：PnPm[119]、DG/FV混合格式[120]、BGK(Bhatnagar-Gross-Krook)[121]格式等。

本節以下將從應用現狀、計算精度、時間推進、激波捕捉、計算成本、開源軟件等方面進行評述。因篇幅限制，以下內容將重點關注發展歷史較長、關注較多的兩種高階格式：WENO和DG。

2.2.1 高階格式的應用現狀

高階格式發展至今，已在多個方向得到了不同程度的應用。

2.2.1.1 氣動聲學數值模擬中的應用

1)激波與旋渦的相互作用[122]。激波與旋渦相互作用是激波噪聲的一個簡化模型，通過研究激波與旋渦相互作用，可以深入了解激波噪聲產生機制。張樹海等采用五階精度WENO格式，系統模擬了激波與強旋渦相互作用系統的模擬，發現激波與強旋渦相互作用存在多級干擾。第一級干擾是入射激波和初始旋渦的相互作用，第二級干擾是反射激波和變形旋渦的相互作用，第三級干擾是渦心附近產生的小激波串和變形旋渦的干擾，如圖22所示。與此同時，旋渦的變形也存在多級特征。在第一級干擾中，初始圓形旋渦被壓成橢圓形；在第二級干擾中，橢圓形旋渦被壓回成圓形旋渦；在第三級干擾中，圓形旋渦又被壓縮成橢圓形。此外，這種多級干擾會引起旋渦在水平線附近做振蕩運動。根據激波與強旋渦相互作用過程中的多級干擾過程，當時張樹海等預測，這種多級干擾過程會產生更為復雜的聲波。

圖22 激波與強旋渦的干擾過程 (Ms=1.2, Mv=1.0)Fig.22 Interaction between a shock wave and strong vortex (Ms=1.2, Mv=1.0)

2) 激波與剪切層相互作用[123]。激波與剪切層相互作用是一個典型的問題，也可以看做是研究超聲速噴流激波噪聲的一個簡化模型。為了了解超聲速噴流激波噪聲的產生機制，劉旭亮和張樹海設計了一個激波與剪切層相互作用的模型。一道斜激波打到一個剪切層上，下邊界是一個反射固壁，透射激波打到固壁上后，形成反射激波，反射激波再次與剪切層相互作用發生反射，形成類似于超聲速噴流的激波串結構。通過較為系統的數值模擬，揭示了兩種噪聲產生機制，一是激波與旋渦相互作用，發生在圖23所示的區域1，另一是激波泄露機制，發生在圖23所示的區域2。

圖23 激波穿過剪切層的兩個區域Fig.23 Two zone of shock wave passing through supersonic shear layer

3) 軸對稱超聲速噴流噪聲[124]。超聲速噴流噪聲是一種典型的激波噪聲，與很多工程密切相關。采用五階WENO 格式直接求解軸對稱Navier-Stokes方程，數值模擬了射流馬赫數(完全膨脹)為Mj=1.19(嘯聲呈現軸對稱模態，三維效應不明顯)的軸對稱欠膨脹射流，其對應的射流聲速馬赫數為Ma=1.05，雷諾數為Re=6.216×105，噴管直徑為D=25.4 mm，噴管厚度為0.4D，環境溫度為288.15 K。圖24給出了噴管唇口壁面[0.0,0.642D]和噴管出口平面[0.0,2.0D]上壓力監測點處噪聲信號的頻譜信息,包括頻率和聲壓級。結果顯示，在大于5000 Hz的頻率區間內，有4個尖銳的聲壓級突起：前兩個聲壓級突起的頻率為65 657 Hz (128 dB)和8637 Hz (123 dB)，分別是嘯聲A1模態和嘯聲A2模態；另外兩個聲壓級突起的頻率為12 087 Hz(125 dB) 和14 310 Hz (123 dB)，分別是嘯聲B模態的諧頻和嘯聲A0模態。綜合流場結構和聲場信息可以得到的嘯聲模態的圖像及其產生位置，其中流場結構用數值紋影(密度梯度的模)表示，聲場由脹量場表示。根據嘯聲頻率可以得到嘯聲的波長，利用波長可以在脹量場中清晰地分辨出嘯聲A1模態和A2模態,相應的波長分別為λ=2.04D和λ=1.51D。從圖25中可以清晰地分辨出向上游傳播的嘯聲，并且嘯聲產生位置都是在第三個激波柵格和第五個激波柵格之間的區域。

(a) 噴管唇口壁面

(b) 噴管出口平面

圖25 噴管厚度為0.4D射流馬赫數為1.19的射流嘯聲圖像Fig.25 Screech of supersonic jet with the diameter being 0.4D at Ma=1.19

2.2.1.2 湍流數值模擬中的應用

1) NLR7301 翼型極限環振蕩[125]。采用流固耦合方法結合五階WENO差分格式計算NLR7301翼型極限環振蕩問題。圖26是彈性支撐NLR7301翼型的模型。流動計算中黏性項采用四階中心差分格式，湍流模型為Spalart-Allmaras模型。來流馬赫數為0.753。時均迎角為-0.45°，初始流場分別采用了定常流計算的解和均勻來流。圖27給出了計算結果的極限環振蕩的振幅，可以看到，盡管中間過程有所不同，兩種初始條件達到了幾乎同樣的極限環。

圖26 彈性支撐NLR7301翼型的模型Fig.26 Model of the elastically mounted NLR7301 airfoil

圖27 彈性支撐NLR7301翼型的俯仰運動情況，Ma=0.753Fig.27 Predicted plunge motion for the model of the elastically mounted NLR7301 airfoil

2) 渦流發生器高超聲速繞流直接數值模擬。圖28為基于高階有限差分WENO格式的高超聲速可壓縮流動直接數值模擬，計算構型為微型渦流發生器。這類問題的特點是外形相對簡單，利于保證網格質量，同時對精度和效率要求較高，所以一般采用有限差分類高階格式。

(a) 計算構型

(b) 微型渦流發生器附近網格

(c) 密度瞬時云圖

3) 整車大渦模擬。對于這類復雜幾何構型的精細模擬，則需采用有限元類高階格式，以便在非結構網格上實現高階精度。圖29為采用五階DG格式對整車開展大渦模擬的結果及網格示意圖[128]。這類模擬的目的是對車身部件流場特征進行分析，優化汽車高速行駛的空氣動力學性能，關鍵是要準確刻畫漩渦結構的發展和相互作用，因而必須采用高階格式。

(b) 瞬時流場

2.2.1.3 面向工程問題的定常氣動力計算

目前這類算例多見于運輸機構型(如圖30所示)的氣動力計算，這類構型能夠獲得高質量的結構網格，王運濤等[129]已成功地使用五階WCNS格式對這類構型開展了計算。該問題對于有限元類格式帶來的主要挑戰是如何克服高階離散的剛性使計算收斂到定常解，目前公開發表的文獻中已經出現三階DG格式的計算結果[130]。

圖30 運輸機構型氣動力特性計算模型示意圖Fig.30 Computational model for a transport aircraft

2.2.2 高階格式的計算精度

目前人們普遍的共識是與采用同樣網格規模的二階格式相比，無論是對于定常RANS計算還是多尺度求解問題，高階格式的計算精度都將獲得顯著的提升。圖31對比了典型高超聲速飛行器繞流在有限差分框架下二階格式和五階格式的計算結果，可以看出二階格式得到的“二次渦”結構清晰程度明顯弱于五階格式，二階格式過大的數值耗散抑制了邊界層的不穩定發展，抹平了流場中的細節結構。此外，對于定常氣動力/熱計算，盡管湍流模型的誤差常常起主導作用，高階有限差分格式也體現出了顯著的優勢[131-132]。

(a) 流動模擬結果示意圖

(b) 背風面x方向不同站位截面的密度云圖，二階差分格式

(c) 背風面x方向不同站位截面的密度云圖，五階差分格式

對于給定的網格，DG等有限元類格式提高近似階數帶來的精度提升要遠高于有限差分類格式。這是因為有限元類格式提高階數會增加自由度數，而有限差分離散通過擴展插值模板提升精度，其自由度數是固定的。因此，有限元類格式能夠在較低自由度個數下獲得更高的計算精度。這一結論在更復雜的問題中也得到了驗證。圖32給出了CRM(Common Research Model)模型計算結果對比。從中可以看到，在采用相當甚至更少自由度的前提下，高階SUPG和DG格式所得阻力系數落在了收斂區域內，而此時SUPG和DG所采用的網格要遠稀疏于二階有限體積法。對于P2(三階格式)情況，計算的分散性大大改善，即在最粗的網格上也得到了比較精確的阻力系數值。

圖32 CRM模型阻力系數計算結果[133](HOMA和COFFE為SUPG方法，DG3D和SANS代表DG方法)Fig.32 Drag coefficients of the CRM model[133] (HOMA and COFFE use the SUPG method, DG3D and SANS use the DG method)

如前所述，對于高階有限差分格式，需要解決的問題主要是如何在復雜構型上生成高質量的結構化網格。而對于高階有限元類格式，要實現上述高階精度仍然面臨諸多挑戰。首先，如果近似解的多項式的階數是P次，曲面邊界需要采用P+1次近似才能使得DG類格式達到設計精度。而復雜外形的高階網格生成仍極具挑戰性，相關功能在網格生成軟件中的集成也還處于起步階段。其次，對于多尺度求解問題，高階有限元格式存在因不能精確積分非線性通量而帶來的混淆誤差，導致數值計算不穩定，這一現象在“高階精度+稀疏網格”的組合下尤其嚴重[134]。此外，高階有限元類格式一個顯著的特點是在高波數區域具備很強的耗散，這一般被看作是此類方法適用于隱式大渦模擬的主要依據。但高階有限元類格式的多模態屬性使得控制和調節格式的傅立葉特性變得困難(如DG格式的auto-correction機制[135])，這給格式性能的優化以及亞格子模型的構造都帶來了挑戰。

2.2.3 高階格式的時間求解方法

現代CFD方法一般都采用半離散形式，時間格式的構造不依賴具體的空間離散，有限差分與有限體積法所采用的時間格式一般不需隨格式精度變化。因此，本部分主要關注有限元類高階格式的時間求解。定常CFD計算常用的時間格式主要包括兩大類：近似算子分裂方法和Newton-Krylov迭代法。對于有限元類格式，Block Lower-Upper Symmetric Gauss-Seidel (BLU-SGS) 格式是典型的近似算子分裂方法，通過求解非線化方程，避免了方程線化帶來的誤差，能夠取得更快的收斂，因而被推廣到了諸如SD[136]、FR[137]、CPR[138]等方法中。另外，Nastase 等[139]系統比較了線性 Element Jacobi(EJ)、非線性 Element Jacobi 以及 Linear Gauss-Seidel(LGS)格式。結果表明，LGS格式雖然需要額外存儲非對角矩陣，但卻能夠取得最優收斂效果。

在基于 Newton-Krylov 的迭代法中，最常用的為廣義最小殘差法 (GMRES)。對Newton-Krylov方法的研究首先集中在預處理技術上，Persson 等[140]對比了多種預處理技術，認為線性多重網格以及不完全因子分解(incomplete lower-upper factorization，ILU)技術配合使用效果較佳。Diosady 等[141]提出了基于張量積的預處理技術，并對比了與算子分裂預處理子的性能差異。迭代法涉及大型線性方程組求解，會導致計算量和存儲量急劇增加。例如，采用DG格式求解三維Navier-Stokes方程，假設使用 100 000 個四面體單元，應用 3 次多項式，雙精度計算，存儲雅克比(Jacobian)矩陣非零塊的內存消耗約為102.1 GB，除需要大量存儲空間外，計算Jacobian矩陣本身也將花費大量的計算資源。Luo等[142]探討了避免矩陣存儲的技術，能夠有效降低內存消耗。Xia 等[143]對比了 GMRES 方法多種雅可比矩陣(Jacobian Matrix)計算方法的優劣，并評估了一系列隱式格式的計算效率。

高階有限元類格式空間離散后得到的方程系統剛性往往較強(尤其是對于大長細比網格[144]、跨聲速、復雜外形等情況)，求解定常解的難度和代價都會急劇升高。這一點也已在評估近似算子分裂方法和 GMRES 迭代法兩大類隱式方法的研究工作[145]中得到體現。計算結果表明，即使對于平板、翼型這類簡單問題，BLUSGS也常常計算發散， GMRES類格式收斂性能更好，但也需要謹慎選擇預處理方法。高階格式定常計算的求解策略往往需要先用低階格式計算獲得初始解，隨后再用高階格式續算直到計算收斂。此外，此類方法的收斂速率嚴重依賴于網格，即隨著網格的加密收斂性能迅速惡化(如圖33所示)。上述困難是目前在復雜構型問題中少有四階或更高階有限元類格式計算結果出現的主要原因。

(a) P=1 (二階格式)

(b) P=2 (三階格式)

對于非定常計算，高階時間離散格式如隱式 Runge-Kutta(Implicit Runge-Kutta, IRK)方法，已被推廣應用于 DG 非定常算法的構造中[146]。Wang 等[147]研究了從低階到高階一系列時間離散格式，發現在同樣的誤差要求下，高階時間離散方法比低階時間離散方法的效率更高，尤其是隨著誤差精度的進一步提高，高階時間格式的優勢將更加明顯。Vermeire 和 Nadarajah[148]發展了一種顯示-隱式混合時間推進方法。相比 IRK 方法，計算效率更高的高階時間離散格式相繼被提出[149-150]，這類方法在每個時間步只需進行一次雅可比矩陣計算。

2.2.4 高階格式的激波捕捉方法

對于強可壓縮流動，高階格式面臨的一大挑戰是：以極小的數值耗散來分辨小尺度脈動，在光滑區保持格式的高精度特性不被破壞，同時光滑地捕捉流場中出現的各種間斷。對于有限差分類格式、有限體積類格式以及有限元類格式，常用的激波捕捉策略是相似的，主要包括：限制器、WENO類插值/重構、人工黏性模型等。限制器類方法[151-152]以保證格式的TVD或TVB性質為目標，優點是穩定性好，常用于二階格式，但往往因為耗散過大，不利于實現高精度，因而在各類高階格式中應用較少。

在有限差分的框架下，WENO類方法成為最常用的可壓縮流動高階格式，獲得了巨大的成功。其核心思想就是將格式的設計模板分成若干子模板，在每個子模板上重構一個低階數值通量，通過某種光滑測試因子計算每個子模板的權，將所有子模板的數值通量非線性疊加，得到高階或者最優精度的數值通量。WENO格式具有如下特點，在光滑區，格式可以達到最優精度，在激波區，基本選擇了最光滑的子模板。而經典的WENO格式在面對各種不同的問題也逐漸暴露出一些不足之處，針對WENO格式的改進主要包括兩類：1)針對WENO格式在極值點附近降階的問題進行改進，如WENO-M[153]和WENO-Z[154]格式；2) WENO格式計算定常激波時，在激波下游會出現微弱的非物理波動，導致殘差不能收斂到滿意程度，對于激波噪聲，這種微弱的非物理波動則非常大，有時比激波噪聲的幅度高兩個數量級。張樹海等提出了兩種方法消除這種微弱的非物理波動，主要包括一種適合于五階WENO格式的光滑測試因子[155]和迎風型插值方法[156]。

而在有限體積的框架下，提高精度需要擴展模板，這使得WENO類方法的構造較為復雜，向高階格式推廣往往存在穩定性問題。當將此類方法推廣到有限元框架下時，為了保證計算效率，需要依賴額外的激波指示器[157]，先判斷那些需要進行修正的單元(即包含間斷的單元)，然后在這些單元上采用WENO重構，而在其他單元處采用原始的高階格式。因此，激波指示器的性能就很大程度上決定了基于WENO重構方法的計算精度和穩定性。在有限元的框架下引入WENO重構，初衷是為了利用WENO的高精度特性，有效保持有限元離散的精度，然而，實際應用發現，WENO重構很難推廣到三階精度以上。作者認為，這主要是因為WENO重構往往涉及多維高次多項式的內插或外插過程，容易導致計算結果的不穩定。Dumbser等[158-159]利用后驗思想MOOD(Multi-dimensional Optimal Order Detection)來對包含激波的單元進行修正，修正時將一個DG單元劃分成若干個子單元，在單元內部引入WENO等方法進行限制，所給出的算例最高精度達到10階，很好地保留了DG格式的高階特性。

捕捉激波的另外一種方法是人工黏性方法，即通過添加人工黏性項來抑制虛假振蕩。人工黏性方法最早可以追溯到1950年von Neumann 和 Richtmyer的工作[160]，隨后Baldwin和MacCormack[161]以及Jameson[162]沿著這一方向均取得了杰出的成果。這些早期的工作大多采用低階有限體積方法。而當格式精度提升后，有限差分和有限體積框架下的人工黏性方法則因耗散過大而較少采用。與之相反，人工黏性與高階有限元類方法的結合能夠在一個單元內分辨激波[163](即亞格子分辨率)，使得人工黏性這一經典概念煥發新生，引起了人們的極大興趣，也誕生了眾多方法[164-167]。人工黏性方法的關鍵是人工黏性系數的計算，需要設計能夠可靠衡量流場光滑性的指示器。目前常用的光滑性指示器大致有三種：1)流場關鍵參數(如密度)的導數；2)模態分解系數的衰減速率；3)熵增。研究表明[168]，這三種模型都有各自的缺點：1)基于導數的模型需要根據流動特性結合開關機制來減少對于小尺度脈動等的耗散；2)基于模態分解的模型往往存在高估光滑性的情況；3)基于熵增的模型則耗散偏大，且對于經驗參數的依賴性較強。而當格式階數較高(一般大于4階)的時候，各模型均能夠得到分辨率較好的結果，并且不同模型之間的差別也逐漸減小。

2.2.5 高階格式的計算代價分析

計算代價方面，目前，對不同格式計算代價的評估大多停留在以相同自由度個數為基準的前提下。例如，文獻[169]針對引擎內流的大渦模擬計算開展格式效率的對比，在誤差相當時，四階格式所需自由度數比二階格式低約一個數量級，也顯著低于基于二階格式的LES/RANS混合方法。然而需要指出的是，真正嚴格的對比則應是在達到相同誤差時衡量格式的計算效率。此外，評估計算代價還應考慮算法與硬件的結合。現代基于CPU的大規模并行計算集群雖然能夠有效壓縮計算周期，但計算成本仍然巨大。而DG、FR這類高階格式使用的網格單元數少，局部計算量集中，更利于充分利用現代計算架構的浮點運算能力，更適合GPU計算。基于CPU并行策略的有限體積求解器一般只能使用并行集群性能峰值的5%～10%，而使用GPU計算的FR格式計算速度最高可以達到性能峰值的58%[170]。

Vincent等[171]從誤差和計算代價的角度系統對比了FR和基于二階有限體積法的商用軟件STAR-CCM+在湍流計算中的計算成本，圖34給出了計算所用的成本(硬件購置成本×計算時間)和流場結果的對比。在計算成本相當的前提下，5階FR計算得到的結果更精細。

(a) 計算成本對比

(b) 密度等值面圖

需要指出的是，上述計算成本對比選取的算例屬于分離流動，更容易體現高階格式的優勢。而對于定常RANS計算，高階格式的計算代價是否具有競爭力仍需進一步研究。發展h-p自適應方法將是進一步提高有限元類高階格式在計算代價方面競爭力的有效途徑。

2.2.6 軟件實現

高階格式因為算法設計復雜，編程困難，在一定程度上也限制了其推廣和應用。鑒于此，一些研究團隊開發了高階格式的開源軟件，有助于幫助初學者快速完成基于高階格式的研究和計算任務。以下列舉部分開源軟件：

1) 有限差分類：中科院力學所的OpenCFD(WENO格式、緊致格式等)；

2) 有限元類：如英國帝國理工大學的Nektar++[172](DG格式)和PyFR[173](FR格式)，美國斯坦福大學的HiFiLES[174](FR格式)，瑞士洛桑聯邦理工大學的NUDG[175](DG格式)，美國德克薩斯A&M大學的deal.II[176](有限元方法一般框架)等。

3 發展趨勢及建議

本文針對航空航天CFD應用中的物理模型和計算方法等核心問題，評述了其發展現狀和面臨的挑戰。一方面，CFD的應用無論是深度和廣度都有了快速發展，其作用也越來越大；另一方面，CFD在關鍵的核心理論等方面都還面臨著諸多問題和困難，進展緩慢。這促使我們認真思考未來的CFD應該如何發展，這里將分別針對本文中的各個方向，給出一些思考、建議和意見，以期拋磚引玉，圖畫未來。

3.1 湍流模型

1) 對于使用者來說，建議重視如下問題：通過相關文獻等方式，了解湍流模型的基本性能和特點，注意模型中常數或參數的選擇原則、邊界條件正確的設置方法、對網格尤其是近壁網格和大流動梯度處等網格的限定和要求，盡量通過多種手段保證湍流模型的收斂性。

2) 繼續改進、完善和發展現有的渦黏性湍流模型，提高預測精度、擴展適應性和預測能力：根據近二十年來國外湍流模型的新發展，選擇評價較高、實用性較好、性能優越的新發展模型，如高魯棒性的橢圓松弛法、K-KL模型等，促進其實用化發展和推廣應用。對于現有經典兩方程模型，重點研究長度尺度和近壁模擬方法的評估和改進。

3) 支持RSM模型的長期研究，大力推廣其在航空航天領域的應用。

4) 在發展和驗證新的湍流模型中，重視DNS和LES的作用，使其能提供寶貴的數據庫、流動結構和機理模型，尤其是建模指導等。

5) 美國NASA等在2017年召開了Turbulence Modeling Symposium[177]，會議主要討論湍流模型的現狀、新思想以及應對未來的問題。與會者充分討論后，認為近期應該開展2個方面的研究： 1)不確定性量化(UQ)； 2)數據驅動建模。

3.2 轉捩模型

隨著對轉捩物理機理認知的逐漸深入和實驗數據的不斷豐富，轉捩模型已取得了豐碩的研究成果。目前已有數十種轉捩模型被提出，這些模型對特定流動特征下的轉捩預測具有較高的精度，特別是在二維簡單轉捩流動中。然而就目前的發展來看，這些轉捩模型距離在工程三維復雜外形中的應用尚有一定差距。主要困難如下：

1) 邊界層轉捩路徑的復雜化和影響因素的多樣化導致尚無關于轉捩系統完備的物理機制解釋；

2) 轉捩現象本身的高復雜性和強非線性，使得在基于物理的方程框架中對所有轉捩機制進行完整建模十分困難，特別是涉及線性和非線性效應及其相互作用機制時。因此，如何將轉捩內在機制合理準確地引入轉捩模型仍是個難題；

3) 轉捩模型的普適性嚴重不足，單類方法往往只適用于特定流動類型，對其他類型轉捩流動的預測性能顯著下降或完全失效；

4) 為使方程封閉，轉捩模型均不可避免地引入了大量模化參數，從而導致預測結果的不確定度。

建議在后續研究中應重點關注以下幾點：

1) 加強轉捩機制研究和理論建模，使得轉捩模型能夠反映更本質的物理機制；

2) 充分利用飛行/風洞實驗數據、DNS等高精度數值模擬手段、穩定性理論結論和機器學習等數據挖掘技術，校驗轉捩模型相關模型的正確性和準確性，實現轉捩模型的預測能力提升；

3) 普適性本身就是一個值得探討的問題，一方面可考慮針對性地發展特定流動類型轉捩模型，實現術業有專攻，另一方面引入混合模型，通過提取不同轉捩機制類型的特征參數，預先判定轉捩主導機制，實現分區轉捩預測；

4) 針對轉捩模型在工程應用中所面臨的普適性差、參數調整對經驗依賴性強等問題，開展轉捩模型的不確定度量化和敏感性分析。

3.3 通量求解

通量格式作為CFD領域的熱點問題在過去的幾十年間取得了豐碩的研究成果。但是，在針對流動分離、激波邊界層干擾等多維復雜流動開展數值模擬時，現有方法仍具有模擬精度不足、計算結果嚴重依賴于網格分布等缺陷。此外，尚無一種通量格式可以完全避免非物理的激波異常和“overheating”現象。為此，本文作者有如下建議：

1) 鑒于現有通量格式無法完全避免激波異常現象的出現，在模擬高速流動時，應盡量選取激波穩定性較佳的通量格式，并注意強激波附近網格長細比的合理設置，通過高質量的網格生成抑制求解過程中可能出現的非物理流場脈動。

2) 應深入挖掘膨脹流動模擬中非物理“overheating”現象出現的數學機制，從理論上找到改進這一缺陷的具體措施。

3) 推廣全速域通量格式在航空航天領域的應用，改善當前主流軟件在寬速域復雜流動模擬中的預測精度以及計算網格依賴性。

4) 針對當前通量格式在多維復雜流動模擬中精度不足的缺陷，大力發展真正多維全速域通量格式，進而提高CFD技術在航空航天領域的可靠性。

3.4 高階格式

1) 高階格式框架的選取取決于求解的具體問題：a)幾何外形簡單的大渦模擬或直接數值模擬建議選擇有限差分類格式；b)復雜構型的定常RANS計算，在能夠生成高質量結構網格的前提下可采用高階有限差分格式，如WCNS格式等。在幾何構型更為復雜的情況下，可嘗試選擇三階精度的DG格式，更高階精度的DG格式會面臨收斂性和魯棒性等問題；c)復雜構型的湍流精細模擬可選擇有限元類的高階格式。

2) 高階格式在向實際復雜問題推廣過程中遇到的困難是全方位的，涉及網格生成、激波捕捉、時間推進、計算穩定性、計算效率等各方面，尤其亟需發展高效率、低存儲、魯棒性強的隱式時間格式。高階格式面臨的這些困難在二階格式中已經存在，在高階格式的框架下則變得更加突出。

3) CFD其他關鍵技術比如通量格式、物理模型等在高階格式的框架下，性能也會發生不同程度的變化。因此，高階格式的發展還應注意結合這些CFD的關鍵環節綜合考慮。

4) 目前來看，大規模并行集群的性能并未得到充分的利用。在計算機硬件更新換代日新月異的今天，高階格式要提高競爭力，在方法發展過程中應當重視與計算機硬件的緊密結合，充分發掘硬件資源的計算潛能。