袁先旭,陳堅(jiān)強(qiáng),杜雁霞,郭啟龍,肖光明,傅亞陸,梁飛,涂國(guó)華,*

1. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000

2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000

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

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

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

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

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

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

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

1.1 轉(zhuǎn)捩預(yù)測(cè)模型

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

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

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

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

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

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

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

1.2 湍流模擬模型

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

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

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

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

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

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

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

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

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

2.1 氣/液兩相流模型

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

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

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

張阿漫等[30]基于Roe近似黎曼求解器建立了多相流SPH模型，引入耗散項(xiàng)限制因子模擬黏性流動(dòng)，保證了模擬大密度比復(fù)雜流動(dòng)界面時(shí)壓力場(chǎng)的光滑性和穩(wěn)定性。針對(duì)三維軸對(duì)稱問題，在柱坐標(biāo)系下對(duì)軸對(duì)稱SPH方程進(jìn)行了全新的推導(dǎo)，建立了軸對(duì)稱多相流模型[31]，并通過經(jīng)典的表面張力測(cè)試和氣泡上浮算例驗(yàn)證了模型的有效性。圖9[30]給出了Rayleigh-Taylor不穩(wěn)定性問題求解的壓力云圖變化，圖10[31]對(duì)比了不同Ar數(shù)和Bo數(shù)下的上浮氣泡最終形狀，與試驗(yàn)結(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方法求解流場(chǎng)信息，并發(fā)展了如圖11所示的虛粒子耦合算法來傳遞單元-粒子交界面的物理信息，模擬了彈性擋板液艙晃蕩減阻問題；同時(shí)將虛粒子耦合算法延伸應(yīng)用到了流體-結(jié)構(gòu)共軛傳熱問題中[33]，突破性地基于FEM-SPH方法對(duì)傳熱-流體-結(jié)構(gòu)耦合問題進(jìn)行了模擬分析，圖12給出了耦合過程中某一時(shí)刻的溫度場(chǎng)云圖。

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

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

2.2 飛機(jī)結(jié)冰模型

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

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

圖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ù)光滑表面、粗糙表面及不同潤(rùn)濕性表面上液滴飛濺的試驗(yàn)結(jié)果，提出了更具普適性的液滴發(fā)生飛濺的臨界條件預(yù)測(cè)公式，圖14給出了預(yù)測(cè)公式計(jì)算結(jié)果與已有文獻(xiàn)數(shù)據(jù)對(duì)比，驗(yàn)證了所提關(guān)系式的普適性。

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

2.3 流/熱/固耦合模型

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

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

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

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

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

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

圖17 基于熱格子Boltzmann方法的A320艙內(nèi)流場(chǎng)預(yù)測(cè)結(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ù)值預(yù)測(cè)模型及數(shù)值計(jì)算方法。圖18給出了不同溫度、不同壓力下高溫接觸熱阻的變化規(guī)律。

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

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

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

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

2.4 熱氣動(dòng)彈性預(yù)測(cè)模型

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

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

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

3 多物理場(chǎng)耦合計(jì)算模型與方法

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

3.1 高溫近壁流動(dòng)與材料表面催化燒蝕

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

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

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

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

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

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

3.2 稀薄過渡流計(jì)算模型與方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖28 間斷附近各點(diǎn)通量計(jì)算中最大最小光滑因子的比值[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)計(jì)算方法，朱華君等[70]利用與有限差分混合的方法，構(gòu)造了具有強(qiáng)激波捕捉能力的計(jì)算方法，結(jié)合了兩種方法的優(yōu)勢(shì)，顯著提高了這類高階精度計(jì)算方法在高超聲速流動(dòng)中的魯棒性和應(yīng)用能力，圖29給出了該方法在極高M(jìn)ach數(shù)問題計(jì)算中得到的密度云圖。

圖29 強(qiáng)激波捕捉高階精度CPR計(jì)算極高M(jìn)ach數(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方法黏性項(xiàng)離散的直接間斷有限元方法(DDG)，并將其推廣應(yīng)用于高階精度DG/FV混合方法中。研究表明DDG方法相比傳統(tǒng)DG黏性項(xiàng)離散的BR2方法具有單步計(jì)算量少，計(jì)算收斂所需步數(shù)少等優(yōu)點(diǎn)。典型旋成體低速層流繞流算例三階精度DDG-DG1/FV2相比BR2-DG1/FV2計(jì)算效率可以提高1.7倍，相比三階BR2-DGP2計(jì)算效率可以提高4.3倍(圖30)。

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

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

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

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

圖32 非結(jié)構(gòu)高階精度算法湍流圓柱繞流計(jì)算網(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)高精度算法理論框架，建立了不低于三階空間精度的緊致有限體積算法，其中涵蓋了保精度限制器、高效時(shí)間推進(jìn)方法、高階網(wǎng)格算法、高精度邊界處理方法等，保證了算法的分辨率、計(jì)算效率以及魯棒性。

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

5 總結(jié)與展望

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

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

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

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

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

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

致 謝

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