國家數(shù)值風洞(NNW)工程中的CFD基礎科學問題研究進展

袁先旭,陳堅強,杜雁霞,郭啟龍,肖光明,傅亞陸,梁飛,涂國華,*

1. 中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室，綿陽 621000

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

隨著計算機硬件能力的飛速提升，計算流體力學(CFD)已逐漸成為支撐武器裝備、航空航天、交通運輸和節(jié)能環(huán)保等諸多軍民工業(yè)領域發(fā)展的核心手段。根據(jù)估算，波音飛機的設計過程中CFD已占據(jù)50%的份額，并預測隨著CFD技術(shù)的發(fā)展，未來CFD將達到70%的份額[1]。

“數(shù)值風洞”是實現(xiàn)CFD技術(shù)應用于解決實際工程問題的重要工具，其研發(fā)建設在美歐等西方國家已經(jīng)得到高度重視，目前國內(nèi)市場上常見的與CFD相關(guān)的系列軟件套裝(包括網(wǎng)格生成、流場解算、多學科耦合分析與優(yōu)化等)，基本被西方國家產(chǎn)品壟斷。2018年，中國啟動了由中國空氣動力研究與發(fā)展中心(CARDC)牽頭建設的國家數(shù)值風洞(National Numerical Windtunnel, NNW)工程[2]，旨在自主研發(fā)功能先進、種類齊全的CFD工業(yè)軟件系統(tǒng)，建成擁有中國自主知識產(chǎn)權(quán)、國內(nèi)開放共享的空氣動力數(shù)值模擬平臺。目前，已經(jīng)面向全國發(fā)布了自主知識產(chǎn)權(quán)通用計算流體力學軟件NNW-FlowStar、國內(nèi)首款可控開源流體工程軟件NNW-PHengLEI，在國內(nèi)多家單位開始應用。

NNW工程套裝軟件的開發(fā)、優(yōu)化與升級過程中催生了一批理論與關(guān)鍵技術(shù)創(chuàng)新，其中在軟件架構(gòu)、集成與測試、驗證與確認等方面的主要進展在文獻[3-6]中給出了較為全面的綜述。

為了加強NNW工程套裝軟件在國防工業(yè)與國民經(jīng)濟發(fā)展中的重要作用，CFD基礎理論方面的持續(xù)創(chuàng)新是必不可少的。以航空航天為例，總變差減小(Total Variation Diminishing, TVD)格式的提出及廣泛應用為高超聲速流動提供了高效魯棒的求解方法，增強了CFD對高速飛行器氣動力/熱的鑒定評估能力；而高精度數(shù)值方法和大渦模擬(LES)理論的成熟使涉及邊界層轉(zhuǎn)捩、大范圍分離湍流問題的CFD預測更加精確。因此，NNW工程在若干基礎科學問題研究方面投入了大量資源，通過聯(lián)合國內(nèi)眾多高校及研究機構(gòu)，在湍流與轉(zhuǎn)捩模型、多相多介質(zhì)模型、多物理場耦合模型、高階精度算法等方面取得了一系列階段性的原創(chuàng)研究成果，對NNW工程軟件系統(tǒng)的功能拓展和能力提升打下了堅實基礎。

本文簡要總結(jié)了NNW工程實施以來在若干基礎科學問題方面取得的主要研究進展。隨著工程進度進入后半程，滿足一定成熟度要求的基礎研究成果將會首先集成至NNW-PHengLEI軟件并面向全國開源共享，經(jīng)過充分的檢驗和改進后將進一步集成至面向流體工程應用的通用CFD軟件。

1 轉(zhuǎn)捩與湍流模型及計算方法

轉(zhuǎn)捩與湍流是流體力學領域的兩大“百年難題”，其復雜的流動機理至今尚未完全清楚[7]。為了滿足中國工業(yè)CFD軟件需求，NNW工程針對轉(zhuǎn)捩與湍流預測模型及計算方法做了很多亮點工作，相應預測能力與精度均得到提高。

1.1 轉(zhuǎn)捩預測模型

目前橫流轉(zhuǎn)捩模型預測速域主要集中在亞跨聲速范圍，多數(shù)模型不適用于超聲速/高超聲速橫流轉(zhuǎn)捩預測[8]。張毅鋒等針對超聲速/高超聲速情況下三維邊界層轉(zhuǎn)捩面臨的轉(zhuǎn)捩機制復雜、影響因素多、缺乏標定數(shù)據(jù)等難點問題，基于γ-Reθ轉(zhuǎn)捩模型，提出了合理的壓力梯度修正、普朗特數(shù)修正[9]，進行了橫流轉(zhuǎn)捩拓展，并開展了橫流轉(zhuǎn)捩模型參數(shù)不確定度分析[10]。陳堅強等[2,11-13]進一步拓展了高超聲速橫流轉(zhuǎn)捩數(shù)據(jù)，構(gòu)造了考慮表面粗糙度影響的高超聲速橫流轉(zhuǎn)捩判據(jù)，提出了C-γ-Reθ轉(zhuǎn)捩模型，并將該模型成功用于高超聲速大攻角橫流轉(zhuǎn)捩預測，實現(xiàn)對帶攻角錐的常規(guī)風洞試驗/飛行試驗、HIFiRE-5靜音風洞試驗的高超聲速三維邊界層轉(zhuǎn)捩預測，預測結(jié)果與試驗結(jié)果吻合良好(圖1)，并結(jié)合線性穩(wěn)定性理論對HyTRV升力體表面的高超聲速橫流轉(zhuǎn)捩開展了研究。

圖1 HIFiRE-5橫流轉(zhuǎn)捩預測Fig.1 Prediction of cross-flow transition on HIFiRE-5

在轉(zhuǎn)捩研究中，基于流動穩(wěn)定性理論的eN方法被認為具有較為堅實的理論基礎以及較好的精度。但傳統(tǒng)穩(wěn)定性分析與eN方法受限于積分路徑等，對三維復雜問題實用性較差。黃章峰等[14]開展了感受性、三維邊界層積分路線等方面的研究，基于廣義增長率的守恒特性和三維擾動傳播射線理論提出了改進的拋物化穩(wěn)定性方程(Parabolized Stability Equations, PSE)、eN方法，顯著提高了對一般三維復雜流動穩(wěn)定性分析的計算效率和轉(zhuǎn)捩預測可靠性。圖2[14]對比了原始PSE方法及改進方法的預測結(jié)果，其中x表示沿著積分路徑的流向坐標，A/Ax=110為橫流波波幅比值，ω為圓頻率，散點為直接數(shù)值模擬(Direct Numerical Simulation, DNS)計算數(shù)據(jù)，虛線為原始PSE計算結(jié)果，實線為改進方法(RTPSE)的計算結(jié)果。

圖2 新型PSE/eN方法預測結(jié)果[14]Fig.2 Prediction results of new PES/eN method[14]

畢衛(wèi)濤等[15]基于結(jié)構(gòu)系綜理論(SED)提出了一種構(gòu)建工程轉(zhuǎn)捩模型的新思路，即通過試驗和可靠的計算數(shù)據(jù)確定轉(zhuǎn)捩邊界層的結(jié)構(gòu)系綜,提煉反映轉(zhuǎn)捩邊界層物理狀態(tài)和相似性的多層結(jié)構(gòu)參數(shù),進而形成物理圖像清晰、定量描述精確的新型轉(zhuǎn)捩模型。其研究結(jié)果表明該思路能夠成功應用于刻畫由自由來流湍流誘發(fā)的平板邊界層強迫轉(zhuǎn)捩和有攻角的高超聲速尖錐轉(zhuǎn)捩兩類流動，圖3[15]給出了6°攻角尖錐迎風面熱流SED-SL模型預測結(jié)果與試驗測量的對比，其中圖3(a)為試驗測得的溫差ΔT云圖，圖3(b)為SED-SL模型計算得到的熱流系數(shù)Ch，二者對比說明預測的轉(zhuǎn)捩位置與試驗取得了初步的一致性。

圖3 6°攻角尖錐迎風面熱流SED-SL模型預測結(jié)果[15]Fig.3 Prediction results of SED-SL model on windward surface in 6° angle of attack cone flow[15]

由于實際流動轉(zhuǎn)捩對初值十分敏感，因此實際飛行條件和風洞試驗條件之間的來流差異需要在各種預測模型中加以考量。涂國華等[16]對風洞試驗和飛行試驗的轉(zhuǎn)捩N值進行了關(guān)聯(lián)，發(fā)展了一種適合高超鈍錐的轉(zhuǎn)捩N值關(guān)聯(lián)方法。此外，涂國華等[17]針對三維全局穩(wěn)定性分析的Jacobian系數(shù)矩陣巨大、特征值不易求解的難題，提出了一種全局穩(wěn)定性模態(tài)分解方法，可以有效實現(xiàn)復雜流動的全局穩(wěn)定性分析。

1.2 湍流模擬模型

陳十一團隊[18]開展了可壓縮各向同性湍流的多過程分析方法和高保真度大渦模擬模型研究，基于亥姆霍茲分解和Kavasznay分解方法實現(xiàn)了可壓縮各向同性湍流的多尺度性質(zhì)分析，誤差達到0.1%以內(nèi)。通過使用人工神經(jīng)網(wǎng)絡方法，并結(jié)合可壓縮湍流的物理機理，王建春等[19-20]發(fā)展了可壓縮各向同性湍流的高保真度大渦模擬模型，圖4給出了反卷積人工神經(jīng)網(wǎng)格(Deconvolution Artificial Neural Network, DANN)結(jié)構(gòu)示意圖。新的大渦模型在先驗驗證中的相關(guān)系數(shù)達到0.95以上，在后驗驗證中，可以精確地預測能譜、結(jié)構(gòu)函數(shù)等統(tǒng)計量，以及瞬態(tài)流場結(jié)構(gòu)(圖5)。

圖4 反卷積人工神經(jīng)網(wǎng)格結(jié)構(gòu)示意圖[19]Fig.4 Schematic of structure of deconvolutional artificial neural network[19]

圖5 DANN在各向同性湍流計算中的應用[19]Fig.5 Application of DANN in computation of isotropic turbulence[19]

針對高雷諾數(shù)壁湍流模擬網(wǎng)格需求過大的問題，吳霆等[21]實現(xiàn)了Werner-Wengle壁面應力模型的大渦模擬(LES)。

陳浩等[22]提出了一種基于分區(qū)混合和基于湍流尺度混合的雙重RANS/LES混合模型，合理引入了壁面影響修正、網(wǎng)格拉伸修正以及多重過渡函數(shù)。在具有強非定常特性的大范圍分離湍流模擬中，雙重RANS/LES混合模型的預測結(jié)果較傳統(tǒng)DES類方法有明顯改進。

張偉偉等[23-24]結(jié)合量綱分析、標度理論以及神經(jīng)網(wǎng)絡機器學習方法，發(fā)展了可以獨立于傳統(tǒng)湍流模型的渦黏封閉方法，實現(xiàn)了將數(shù)據(jù)驅(qū)動模型與Navier-Stokes(N-S)方程求解器的動態(tài)耦合，并驗證了方法對不同翼型/機翼/翼身組合體在不同馬赫數(shù)、雷諾數(shù)、攻角等流動狀態(tài)下的泛化能力。結(jié)果表明預測的氣動力系數(shù)符合源數(shù)據(jù)(圖6[24])，在一定流動狀態(tài)范圍內(nèi)可以重現(xiàn)傳統(tǒng)RANS模型的解，體現(xiàn)了機器學習方法在未來湍流模型化工作中的良好前景。

圖6 NACA0012翼型表面摩擦阻力系數(shù)分布對比(實線：SA模型，符號：機器學習)[24]Fig.6 Comparison of skin friction coefficients on NACA0012 airfoil (solid-lines: SA model, symbols: machine learning)[24]

黃生洪等[25-26]針對湍流邊界層入口條件生成問題，基于Kraichnan單波數(shù)脈動隨機算法模型提出了一種入口邊界湍流生成算法，即DSRFG(Discretizing and Synthesizing Random Flow Generation)，成功解決了人工合成湍流產(chǎn)生符合預定湍流能量譜特性和空間相關(guān)性的入口條件(圖7)，圖中f表示頻率，fS表示不同頻率Fourier模態(tài)的幅值，DSRFG從算法上精確滿足了入口邊界解的適定性(主要是湍流脈動的流量連續(xù)條件)。此外，DSRFG方法還具有良好的并行特性，目前已通過波數(shù)空間的log擴展，進一步使DSRFG方法能夠產(chǎn)生大跨度波數(shù)空間湍流譜范圍的脈動。

圖7 DSRFG算法結(jié)合高精度大渦模擬結(jié)果[25]Fig.7 High-order LES result of DSRFG method[25]

2 多相多介質(zhì)計算模型與方法

多相與多介質(zhì)問題研究是航空航天、能源動力等領域的重要理論和關(guān)鍵技術(shù)基礎[27]，相關(guān)問題的數(shù)模模擬也對CFD技術(shù)的發(fā)展提出了極大挑戰(zhàn)。為了拓展中國工業(yè)CFD軟件的應用領域，NNW工程重點針對跨介質(zhì)飛行器出/入水過程模擬、飛行器結(jié)冰與防除冰特性預測、流/熱/固耦合熱質(zhì)傳遞分析以及熱氣動彈性預測等難點問題，發(fā)展了多類多相多介質(zhì)計算模型及計算方法，有效提高了多相流數(shù)值預測的計算精度。

2.1 氣/液兩相流模型

準確捕捉氣液兩相界面是氣/液兩相流動過程模擬的核心問題。光滑粒子流體動力學方法(SPH)作為較流行的無網(wǎng)格粒子方法，具有拉氏計算準確描述物質(zhì)界面的優(yōu)勢，能自然滿足自由面邊界條件，非常適合于模擬介質(zhì)大變形等復雜現(xiàn)象。

施文奎等[28-29]發(fā)展了空間變光滑長度SPH方法和多自由度流固耦合模型，圖8[28]給出了基于發(fā)展的SPH方法模擬帶多自由度運動的三維返回艙水上回收問題的計算結(jié)果。將模擬結(jié)果與試驗結(jié)果進行了定量比對，驗證了方法和模型的有效性，分析了入水沖擊現(xiàn)象和規(guī)律，拓展了SPH方法在復雜入水問題中的應用能力。

圖8 返回艙水上回收模擬[28]Fig.8 Simulation of a space capsule water recovery[28]

張阿漫等[30]基于Roe近似黎曼求解器建立了多相流SPH模型，引入耗散項限制因子模擬黏性流動，保證了模擬大密度比復雜流動界面時壓力場的光滑性和穩(wěn)定性。針對三維軸對稱問題，在柱坐標系下對軸對稱SPH方程進行了全新的推導，建立了軸對稱多相流模型[31]，并通過經(jīng)典的表面張力測試和氣泡上浮算例驗證了模型的有效性。圖9[30]給出了Rayleigh-Taylor不穩(wěn)定性問題求解的壓力云圖變化，圖10[31]對比了不同Ar數(shù)和Bo數(shù)下的上浮氣泡最終形狀，與試驗結(jié)果吻合較好。

圖9 Rayleigh-Taylor不穩(wěn)定性問題壓力云圖[30]Fig.9 Pressure contours of Rayleigh-Taylor instability problem[30]

圖10 不同Ar數(shù)和Bo數(shù)下的最終氣泡形狀[31]Fig.10 Terminal bubble shape with different Ar number and Bo number[31]

劉謀斌等[32]發(fā)展了ES-FEM-SPH方法，利用ES-FEM方法處理結(jié)構(gòu)變形，SPH方法求解流場信息，并發(fā)展了如圖11所示的虛粒子耦合算法來傳遞單元-粒子交界面的物理信息，模擬了彈性擋板液艙晃蕩減阻問題；同時將虛粒子耦合算法延伸應用到了流體-結(jié)構(gòu)共軛傳熱問題中[33]，突破性地基于FEM-SPH方法對傳熱-流體-結(jié)構(gòu)耦合問題進行了模擬分析，圖12給出了耦合過程中某一時刻的溫度場云圖。

圖11 單元與對應虛粒子耦合示意圖[32]Fig.11 Schematic diagram of interaction between element segment and corresponding ghost particles[32]

圖12 溫度場云圖(ES-FEM-SPH方法)[33]Fig.12 Contour of temperature with ES-FEM-SPH method[33]

2.2 飛機結(jié)冰模型

飛機結(jié)冰是多個物理過程綜合作用下的復雜現(xiàn)象，既包含傳熱、傳質(zhì)和液體流動等宏觀過程，又涉及晶粒形核和生長等微觀過程[34]。

針對傳統(tǒng)方法難以對三維孔隙結(jié)構(gòu)進行定量表征的問題，李偉斌等[35]基于結(jié)冰孔隙呈球形、孔徑取值連續(xù)、孔徑服從特定分布等二維圖像的合理假設，提出了基于結(jié)冰二維圖像定量信息的三維微觀結(jié)構(gòu)建模方法。圖13給出了1 cm3結(jié)冰三維孔隙結(jié)構(gòu)的建模結(jié)果，d為孔隙直徑，與實驗結(jié)果的對比表明，該方法可行，為結(jié)冰精細化研究提供一種新途徑。

圖13 1 cm×1 cm×1 cm結(jié)冰的孔隙分布[35]Fig.13 Pores distribution of modeling ice with size of 1 cm×1 cm×1 cm[35]

易賢等[36]根據(jù)光滑表面、粗糙表面及不同潤濕性表面上液滴飛濺的試驗結(jié)果，提出了更具普適性的液滴發(fā)生飛濺的臨界條件預測公式，圖14給出了預測公式計算結(jié)果與已有文獻數(shù)據(jù)對比，驗證了所提關(guān)系式的普適性。

圖14 飛濺量與飛濺參數(shù)之間的關(guān)系[36]Fig.14 Relationship between splashing quantity and parameter[36]

2.3 流/熱/固耦合模型

未來高超聲速飛行器在長時間飛行過程中，不僅涉及傳統(tǒng)氣動熱、結(jié)構(gòu)傳熱及應力變形等多因素耦合的流/熱/固耦合問題，同時也需要進一步考慮新型復合防熱材料的跨尺度效應和艙內(nèi)熱效應等綜合熱效應問題[37]。

楊肖峰等[38]根據(jù)高焓氣流下表面反應傳熱過程的跨尺度動力學機理，建立了如圖15所示的表面效應的CFD/RMD耦合計算方法，獲得具有反應動力學物理意義的界面參量，提升了含多相催化效應的高超聲速飛行器氣動加熱預測能力。

圖15 求解表面跨尺度催化傳熱過程的CFD/RMD耦合計算[38]Fig.15 Transcale CFD/RMD coupling diagram of surface catalytic heat transfer process[38]

李明佳等[39]針對高溫氧氣氛下C/SiC復合材料內(nèi)碳纖維的氧化燒蝕現(xiàn)象，首先通過階梯逼近方法構(gòu)建了C/SiC復合材料的三維幾何模型，并基于格子Boltzmann方法建立了如圖16所示的碳纖維的擴散-燒蝕模型，實現(xiàn)了氧化燒蝕過程中氣體擴散和界面推移的耦合求解。

圖16 平紋編織C/SiC復合材料三維介觀燒蝕單元模型[39]Fig.16 3D ablation mesoscopic model of plain weave C/SiC composite[39]

張超等[40]發(fā)展了基于熱格子Boltzmann方法及CPU/GPU異構(gòu)并行算法，建立了一種適用于復雜構(gòu)型艙內(nèi)流動傳熱的數(shù)值模擬手段，實現(xiàn)了十億量級網(wǎng)格規(guī)模的流/固耦合問題高效求解。圖17給出了含人體模型的A320客艙內(nèi)流場與溫度場的精細化預測結(jié)果。

圖17 基于熱格子Boltzmann方法的A320艙內(nèi)流場預測結(jié)果[40]Fig.17 Prediction results of flow field in A320 cabin based on thermal lattice Boltzmann method[40]

Chen等[41]發(fā)展了高溫條件下固-固粗糙表面幾何形貌表征與構(gòu)建技術(shù)，并綜合考慮界面多尺度傳熱、間隙輻射換熱以及力/熱耦合特性等，建立了典型防隔熱材料固-固界面多尺度傳熱特性數(shù)值預測模型及數(shù)值計算方法。圖18給出了不同溫度、不同壓力下高溫接觸熱阻的變化規(guī)律。

圖18 HTA-C與ZrB2-SiC材料間接觸熱阻預測[41]Fig.18 Prediction of contact thermal resistance between HTA-C and ZrB2-SiC[41]

趙建寧等[42-43]建立了基于多點接觸的高溫界面間斷模型，研究獲取了高溫間斷效應的形成演化機制和影響規(guī)律；構(gòu)建了高效分析局部間斷傳熱和力/熱耦合問題的間斷有限元(DG)方法，突破了間斷問題分析中兼顧精度和效率等的方法瓶頸。圖19[42]給出了某一維間斷問題分析中DG與CG計算效率的對比結(jié)果。在此基礎上，進一步發(fā)展了分區(qū)耦合的間斷/連續(xù)伽遼金混合有限元(DG/CG)方法，形成了適應于復雜工程大規(guī)模計算的高效安全評估能力。圖20[43]展示了基于DG/CG的疏導式防熱結(jié)構(gòu)熱/力耦合分析結(jié)果。

圖19 某一維間斷問題分析中DG與CG計算效率對比[42]Fig.19 Comparison of calculation efficiency of DG and CG for 1D discontinuity problem[42]

圖20 基于DG/CG的疏導式防熱結(jié)構(gòu)熱/力耦合分析[43]Fig.20 DG/CG based thermal mechanical coupled analysis of dredging thermal protection structure[43]

2.4 熱氣動彈性預測模型

長時間氣動力/熱綜合作用下的熱氣動彈性問題給未來高超聲速飛行器的熱安全帶來了前所未有的挑戰(zhàn)[44]。為滿足工程化應用需求，飛行器熱氣動彈性計算不僅需要準確計算單個物理場，還要考慮各個物理場之間的耦合效應，給數(shù)值模擬手段的計算精度與效率均提出了巨大挑戰(zhàn)。

張偉偉等[45]發(fā)展了基于時空雙向耦合策略的高效熱氣動彈性計算方法，構(gòu)建了可適應于結(jié)構(gòu)受熱引發(fā)時變熱模態(tài)的高超聲速非定常氣動力模型。圖21[45]為適用于時變熱模態(tài)的氣動力降階模型示意圖及其與傳統(tǒng)CFD結(jié)果對比。結(jié)果表明，該方法突破了結(jié)構(gòu)模態(tài)隨時間變化，以及材料熱/力學特性隨溫度變化條件下非定常氣動力高效重建這一瓶頸問題，大幅度提升了全耦合熱氣動彈性預測方法的工程實用性。

圖21 適用于時變熱模態(tài)的氣動力降階模型及其與傳統(tǒng)CFD結(jié)果對比[45]Fig.21 Aerodynamic reduced order model for time-varying thermal modes and comparison between its results and CFD results[45]

3 多物理場耦合計算模型與方法

高超聲速飛行器周圍的氣體流動存在稀薄效應、湍流效應、催化效應、燒蝕效應、輻射效應、非平衡效應等復雜的物理現(xiàn)象與化學反應，這些效應相互關(guān)聯(lián)，耦合發(fā)生，在空氣動力學上統(tǒng)稱為多物理場。多物理場耦合問題的研究是模擬高超聲速飛行器真實飛行狀態(tài)的關(guān)鍵技術(shù)基礎。NNW工程根據(jù)各自關(guān)注問題的側(cè)重，開展了不同的物理模型和耦合方法研究。

3.1 高溫近壁流動與材料表面催化燒蝕

在高溫近壁流動方面，劉朋欣[46]、吳正園[47]等開展了高溫化學非平衡湍流邊界層直接數(shù)值模擬研究。劉朋欣等采用完全氣體模型和高溫空氣化學反應模型對比了高超聲速楔形體頭激波后的流動狀態(tài)，圖22[46]給出了平板邊界層瞬時流場的旋渦結(jié)構(gòu)圖，分析了化學非平衡效應對湍流統(tǒng)計量、湍流脈動量的影響趨勢，并分析了邊界層標度律，從圖23[46]可以看出，Van Driest變換后的速度分布趨于一致。吳正園等對高溫氣體效應與湍流的耦合作用機理進行了分析，研究表明,高溫氣體效應使邊界層內(nèi)平均溫度顯著降低,平均密度顯著升高，壁面平均壓強和脈動壓強均增加。

圖22 Multi算例瞬時流場旋渦結(jié)構(gòu)圖(Qcr=0.05，法向高度著色)[46]Fig.22 Instantaneous vortex structures in Multi case(Qcr=0.05, colored by normal height)[46]

圖23 Van Driest 變換后的速度分布[46]Fig.23 Profile of Van Driest transformed velocity[46]

李佳偉等[48]針對高超聲速進氣道前緣“Ⅳ型”激波干擾產(chǎn)生的氣動加熱與結(jié)構(gòu)傳熱多物理場耦合計算問題，發(fā)展了一種基于有限體積法的流-熱-固一體化計算方法，并提出了一種新的雙溫阻模型。該方法可以高效解決流場與結(jié)構(gòu)傳熱穩(wěn)態(tài)求解問題，較快計算出穩(wěn)態(tài)結(jié)構(gòu)與流場的溫度分布。計算結(jié)果表明，“Ⅳ型”激波干擾作用會大大增加壁面最大壓力系數(shù)和熱流系數(shù)，對高速飛行器的熱防護性能提出了更高的要求。

針對典型表面材料的催化/燒蝕特性以及與流場的耦合作用問題，崔志亮等[49]采用分子動力學模擬 MD 方法結(jié)合 ReaxFF 反應勢函數(shù)，研究了高焓氧原子與碳基材料碰撞過程中發(fā)生的表面催化燒蝕過程及機理，探究了氧原子能流密度、石墨烯堆疊層數(shù)以及不同石墨表面晶型結(jié)構(gòu)的影響。

磁流體控制與應用是空氣動力學與電磁學相結(jié)合的內(nèi)容，目前對耦合電磁場的高超聲速流動機理的認識并不透徹。丁明松等[50]研究了高溫氣體效應對高超聲速磁流體控制效率的影響，在低雷諾數(shù)下耦合求解了帶電磁源項的三維N-S流場控制方程和電場泊松方程，對比分析了不同氣體模型對磁流體控制的影響。研究表明，高溫氣體效應會極大地降低磁控增阻效果，會明顯地增強部分表面區(qū)域的磁控熱流減緩效果。

針對霍爾效應對高超聲速磁流體力學控制的影響問題，丁明松等[51]耦合求解各向異性霍爾電場泊松方程和帶電磁源項的高溫熱化學非平衡流動控制方程組，建立了高超聲速流動磁流體力學控制霍爾效應數(shù)值模擬方法。圖24[51]給出了環(huán)形電流密度的分布情況，計算結(jié)果表明，霍爾效應能改變洛倫茲力在等離子體中的分布，降低整體的磁控增阻特性；霍爾效應對高超聲速 MHD 控制的影響，與壁面導電性和壁面附近漏電層的“漏電”現(xiàn)象緊密相關(guān)。

圖24 環(huán)形電流密度大小[51]Fig.24 Density of annular electric current[51]

3.2 稀薄過渡流計算模型與方法

在跨流域飛行復雜流動方面，基于Navier-Stokes方程的CFD方法雖然有較高的計算效率，但對稀薄流動問題的物理描述不準確，亟需探索跨流域復雜流動問題的新型數(shù)值仿真模擬方法。

鐘誠文團隊[52-54]發(fā)展了計算高效、模型精確的確定論和統(tǒng)計論兩套跨流域多尺度高效數(shù)值方法。在確定論方面，提出了非結(jié)構(gòu)速度空間技術(shù)并發(fā)展了相應的積分補償策略，空腔流動三維速度空間網(wǎng)格如圖25所示[54]；發(fā)展了雙原子分子氣體的全流域守恒隱式算法和多重隱式算法，提高了宏觀量的預測精度和計算效率。

圖25 空腔流動三維速度空間網(wǎng)格(49 950個網(wǎng)格)[54]Fig.25 3D velocity space mesh for cavity flow (49 950 cells)[54]

在統(tǒng)計論方面，鐘誠文團隊[55]發(fā)展了二維非結(jié)構(gòu)網(wǎng)格下UGKWP (Unified Gas-Kinetic Wave-Particle) 算法，實現(xiàn)了非結(jié)構(gòu)網(wǎng)格下波粒機制的耦合及多尺度計算，并驗證了該方法在非結(jié)構(gòu)網(wǎng)格下對全流域流動的模擬能力，在連續(xù)區(qū)能夠還原為N-S方法。

費飛等[56-57]提出了基于BGK模型的多尺度隨機粒子方法，比較了DSMC、BGK、FPM 3種方法中黏性系數(shù)與時間步長的關(guān)系(圖26[57])，提高了隨機粒子方法在連續(xù)流區(qū)域的計算精度，進一步擴展了隨機粒子方法在跨流域流動中的計算能力。

圖26 DSMC、BGK、FPM 方法黏性系數(shù)與時間步長的關(guān)系[57]Fig.26 Relation of viscosity and time step for DSMC, BGK and FPM methods[57]

趙文文等[58]提出了一種改進的NCCR+ (Nonlinear Coupled Constitutive Relationsl) 模型，該模型具有模擬中等努森數(shù)非平衡高超聲速流動壓縮區(qū)和膨脹區(qū)的潛力。李廷偉等[59]針對統(tǒng)一氣體動理論方法計算效率低的問題，采用數(shù)據(jù)驅(qū)動的思想，發(fā)展了一種稀薄非平衡流非線性本構(gòu)關(guān)系求解方法。結(jié)果表明，該方法可同時提高計算效率和計算精度，應用潛力較高。

方明等[60]發(fā)展了一種采用組分追蹤加權(quán)因子方法的DSMC模型，并將其應用于考慮稀有氣體電離效應的三維復雜外形再入飛行器繞流計算，計算結(jié)果與飛行測量數(shù)據(jù)具有很好的一致性。

4 高階精度數(shù)值計算方法

高精度數(shù)值計算方法在精度、分辨率等方面與傳統(tǒng)方法相比有明顯的優(yōu)勢，在復雜非定常、多尺度問題大規(guī)模精細模擬方面展現(xiàn)出很強的發(fā)展?jié)摿Γ兄谔岣逳NW流場解算軟件的精確性。然而，高精度方法在魯棒性和高效計算方面仍然有很大的改進空間，有必要進行持續(xù)深入的研究。

4.1 結(jié)構(gòu)網(wǎng)格高精度算法

高精度有限差分(Finite Difference)是最為常用的一類結(jié)構(gòu)網(wǎng)格網(wǎng)精度算法，其中代表性的有WENO類格式等。隨著工程上對高速流動的關(guān)注度越來越高，構(gòu)造混合格式成為了有限差分方法的發(fā)展趨勢。郭啟龍等[61]在加權(quán)格式的算法框架上提出了一種基于特征變量的間斷識別因子，并通過引入歸一化函數(shù)極大地消除了間斷識別因子的問題相關(guān)性，基于這種間斷識別方法，可以構(gòu)造出既能穩(wěn)定捕捉激波又具有較低耗散的混合格式。相關(guān)算法已成功應用于轉(zhuǎn)捩控制、湍流邊界層DNS等[62-63]復雜流動研究中(圖27)。隨后，李辰等[64]將以上間斷識別算法擴展到更窄的四點模板格式，通過降低模板的寬度有效增加了格式的魯棒性。針對復雜外形模擬需求，李辰等[65-66]基于NND格式的模板，通過引入非線性加權(quán)和光滑度量因子改進，提出了WNND、WENO-L兩種面向?qū)嶋H工程應用的高精度有限差分格式。

圖27 高精度混合格式的應用Fig.27 Applications of high-order hybrid scheme

文獻[67-69]基于對WENO格式光滑因子的分析，提出了對正弦函數(shù)為常數(shù)的光滑因子設計準則，基于該準則給出了一個光滑因子通用公式，并構(gòu)造了七階WENO格式及更高階的WENO格式。新光滑因子在流場連續(xù)區(qū)域變化小而在間斷附近變化大，且形式更簡單，圖28[67]中顯示了一維間斷附近不同光滑因子得到的最大值與最小值之比分布，其中新光滑因子(圖中WENO7-S)的比值相比原始WENO7-JS結(jié)果高近2個量級，因此能為格式帶來更好的間斷穩(wěn)定性。該系列WENO格式與其線性格式具有同樣的近似色散關(guān)系，且同時具有計算量小、間斷穩(wěn)定性好和分辨率高的優(yōu)勢。

圖28 間斷附近各點通量計算中最大最小光滑因子的比值[67]Fig.28 Quotients of the largest smoothness indicator divided by the smallest one of every flux near discontinuities[67]

基于高階精度基于重構(gòu)的修正過程(Correction Procedure via Reconstruction, CPR)計算方法，朱華君等[70]利用與有限差分混合的方法，構(gòu)造了具有強激波捕捉能力的計算方法，結(jié)合了兩種方法的優(yōu)勢，顯著提高了這類高階精度計算方法在高超聲速流動中的魯棒性和應用能力，圖29給出了該方法在極高Mach數(shù)問題計算中得到的密度云圖。

圖29 強激波捕捉高階精度CPR計算極高Mach數(shù)問題密度云圖[70]Fig.29 Density distribution of extremely high Mach number problem using high-order CPR method[70]

4.2 非結(jié)構(gòu)網(wǎng)格高精度算法

邵帥[71]發(fā)展了用于非結(jié)構(gòu)DG方法黏性項離散的直接間斷有限元方法(DDG)，并將其推廣應用于高階精度DG/FV混合方法中。研究表明DDG方法相比傳統(tǒng)DG黏性項離散的BR2方法具有單步計算量少，計算收斂所需步數(shù)少等優(yōu)點。典型旋成體低速層流繞流算例三階精度DDG-DG1/FV2相比BR2-DG1/FV2計算效率可以提高1.7倍，相比三階BR2-DGP2計算效率可以提高4.3倍(圖30)。

圖30 旋成體標準算例密度殘值收斂曲線[71]Fig.30 Density residual convergence of laminar flow over a body of revolution[71]

龔小權(quán)等[72]研究了基于精確雅可比矩陣計算技術(shù)的GMRES隱式方法。研究表明基于左端項Jacobi矩陣精確求解的GMRES計算效率相比Jacobi矩陣近似求解的GMRES和LU-SGS具有明顯優(yōu)勢，典型亞聲速無黏ONEAR-M6算例(DG四階精度，馬赫數(shù)為0.4，0°攻角，300萬自由度)計算效率可提高一個數(shù)量級以上(圖31)。

圖31 ONERA-M6機翼DG(P3)的密度殘值收斂曲線[72]Fig.31 Density residual convergence history of DG(P3) for ONERA-M6 wing[72]

燕振國等[73]基于非結(jié)構(gòu)高階精度算法，發(fā)展了一種新的預處理策略，并由此構(gòu)建了高效隱式時間推進方法，顯著提高了大規(guī)模非定常模擬的計算效率，且提高了非結(jié)構(gòu)高階精度方法在湍流高保真模擬等領域的模擬能力，圖32展示了該方法在湍流圓柱繞流問題中計算得到的旋結(jié)構(gòu)。

圖32 非結(jié)構(gòu)高階精度算法湍流圓柱繞流計算網(wǎng)格和旋渦結(jié)構(gòu)[73]Fig.32 Computational mesh and vortex structures of turbulent circular cylinder simulations using unstructured high-order methods[73]

任玉新等[74-75]基于前期提出的結(jié)構(gòu)/非結(jié)構(gòu)高精度算法理論框架，建立了不低于三階空間精度的緊致有限體積算法，其中涵蓋了保精度限制器、高效時間推進方法、高階網(wǎng)格算法、高精度邊界處理方法等，保證了算法的分辨率、計算效率以及魯棒性。

文獻[76-80]針對DDG方法中的激波捕捉問題提出了一種基于強殘差的人工黏性技術(shù)，能夠捕獲激波，且隱式算法可以收斂到機器誤差。同時設計了一種基于本地矩陣存儲的隱式MPI并行算法，在二維、三維算例的測試中，計算效率提升了5%～10%，并且保證了千核并行效率在90%以上。

5 總結(jié)與展望

本文總結(jié)了NNW工程在CFD基礎科學問題研究工作上的主要進展：在轉(zhuǎn)捩與湍流模型及計算方法方面，發(fā)展了C-γ-Reθ轉(zhuǎn)捩模型，改進了eN方法的三維轉(zhuǎn)捩預測能力，建立了基于人工神經(jīng)網(wǎng)絡的大渦模擬模型等；在多相多介質(zhì)計算模型與方法方面，發(fā)展了空間變光滑長度SPH方法和多自由度流固耦合模型，發(fā)展了ES-FEM-SPH方法，建立了表面效應的CFD/RMD耦合計算方法等；在多物理場耦合計算模型與方法方面，發(fā)展了基于有限體積法的流-熱-固一體化計算方法，提出了基于BGK模型的多尺度隨機粒子方法等；在高階精度數(shù)值計算方法方面，構(gòu)造了高精度混合有限差分格式，發(fā)展了新型DDG/FV格式等。

展望未來，NNW工程基礎研究系統(tǒng)將在國內(nèi)各參研團隊的共同努力下，結(jié)合軟件推廣及實際應用中的基礎理論及關(guān)鍵技術(shù)瓶頸問題，持續(xù)開展深入研究，以期進一步擴展國產(chǎn)CFD軟件的能力邊界。具體來說，各個重點方向的發(fā)展思路為：

1) 面向轉(zhuǎn)捩與湍流，DNS能力邊界將隨著計算機性能的發(fā)展而大幅拓展，需要著力發(fā)展配套異構(gòu)設備的高效高分辨率時空算法等技術(shù)；(隱式)大渦模擬、RANS-LES混合模擬將成為工程應用的主流，并需突破“灰區(qū)”抑制、近壁模型等關(guān)鍵技術(shù)。

2) 面向多相多介質(zhì)流動，將更多地涉及微觀/介觀尺度的數(shù)值方法；同時需要進一步建立針對流/固界面大變形的非線性、不同相態(tài)轉(zhuǎn)變等非平衡問題的計算模型與方法。

3) 面向多物理場耦合流動，需要考慮耦合的物理過程逐步增多，需要建立能夠表征3種甚至4種物理過程的“耦合”的計算模型，重點涉及高溫非平衡、氣動噪聲、電磁流等。

4) 面向高精度算法，主要解決實用范圍擴展、魯棒性不滿足工程實用的難題，如基于非結(jié)構(gòu)網(wǎng)格滿足保正、保單調(diào)、能量/熵穩(wěn)定的高精度算法；此外，針對各種極端流動條件下的高精度算法也是需要重點研究的方向。

致 謝

感謝中國空氣動力研究與發(fā)展中心劉朋欣、李辰、向星皓、李明、燕振國、羅勇、李娜, 西北工業(yè)大學朱林陽, 中國科學技術(shù)大學黃生洪、薛正陽對本文的貢獻。