蔣勵(lì)劍，趙文文，陳偉芳，堯少波

浙江大學(xué) 航空航天學(xué)院，杭州 310027

錢學(xué)森[1]根據(jù)來流氣體稀薄程度將氣體流域分為連續(xù)流域、滑移流域、過渡流域和自由分子流域，除連續(xù)流域外的其他3大流域又被合稱為稀薄非平衡流域。在連續(xù)流域，氣體滿足連續(xù)介質(zhì)假設(shè)，因而在宏觀連續(xù)介質(zhì)流體力學(xué)的基礎(chǔ)下，形成了以Navier-Stokes(N-S)方程為控制方程的連續(xù)流數(shù)值模擬計(jì)算方法，結(jié)合物面滑移邊界條件，可將N-S方程的適用領(lǐng)域擴(kuò)展到滑移流域。隨著氣體稀薄程度的增加，連續(xù)介質(zhì)假設(shè)在過渡流域和自由分子流域失效。當(dāng)特征尺度逼近分子平均自由程時(shí)，在微觀氣體動(dòng)理論的基礎(chǔ)上形成了直接模擬蒙特卡洛(Direct Simulation of Monto Carlo, DSMC)方法[2]，以微觀粒子的運(yùn)動(dòng)和碰撞直接模擬宏觀流動(dòng)。DSMC能夠在稀薄流域給出與真實(shí)飛行試驗(yàn)較為吻合的流場(chǎng)與氣動(dòng)數(shù)據(jù)，然而其統(tǒng)計(jì)特點(diǎn)導(dǎo)致計(jì)算需要消耗大量的內(nèi)存和機(jī)時(shí)，限制了該方法在工程領(lǐng)域近連續(xù)和多尺度流動(dòng)中的應(yīng)用。以Boltzmann方程為控制方程的統(tǒng)一氣體動(dòng)理論格式(Unified Gas-Kinetic Scheme, UGKS)方法是一套介于宏觀和微觀之間的介觀方法[3]。然而UGKS在求解過程中依賴速度分布函數(shù)在速度空間的離散，同樣給存儲(chǔ)和計(jì)算負(fù)荷帶來了很大的壓力。針對(duì)稀薄非平衡流域計(jì)算效率與計(jì)算精度的沖突，構(gòu)建同時(shí)具備高效率與高精度的數(shù)值方法在工程領(lǐng)域顯得尤為重要。

近年來，隨著大數(shù)據(jù)科學(xué)的興起、大規(guī)模并行計(jì)算及GPU設(shè)備的普及，機(jī)器學(xué)習(xí)方法在計(jì)算機(jī)視覺、自然語言處理等各個(gè)領(lǐng)域迅速得到應(yīng)用與發(fā)展。人工智能也被認(rèn)為是21世紀(jì)3大尖端技術(shù)(基因工程、納米科學(xué)、人工智能)之一[4]。在流體力學(xué)領(lǐng)域，機(jī)器學(xué)習(xí)也表現(xiàn)出巨大的潛力。在流動(dòng)控制領(lǐng)域，誕生了面向流動(dòng)控制基于機(jī)器學(xué)習(xí)的系統(tǒng)辨識(shí)與降階模型、基于遺傳規(guī)劃的主動(dòng)流動(dòng)控制、基于人工神經(jīng)網(wǎng)絡(luò)與深度強(qiáng)化學(xué)習(xí)的主動(dòng)流動(dòng)控制等方法[5]。在改善RANS模型方面，機(jī)器學(xué)習(xí)方法的應(yīng)用思路可以分為兩類[6]：一是采用增加源項(xiàng)的方式改變模型的控制方程形式，如Tracey等[7]利用機(jī)器學(xué)習(xí)重建RANS模型方程的源項(xiàng)，提出了開發(fā)湍流模型函數(shù)形式的新方法；二是構(gòu)造偏差函數(shù)，再與原RANS結(jié)果相加作為雷諾應(yīng)力預(yù)示結(jié)果，如Wang等[8]用隨機(jī)森林方法學(xué)習(xí)RANS模型偏差量，提高了RANS模型的準(zhǔn)確性；Xiao等[9]針對(duì)RANS和DNS結(jié)果的差異，利用機(jī)器學(xué)習(xí)方法并引進(jìn)了物理約束，對(duì)RANS結(jié)果進(jìn)行了修正。另一方面，也有學(xué)者選擇不再修正現(xiàn)有模型，而是直接尋找流場(chǎng)物理量與湍流之間的關(guān)系，如Zhu等[10]構(gòu)建黑箱模型，采用徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)，建立了基于機(jī)器學(xué)習(xí)的高維數(shù)據(jù)驅(qū)動(dòng)網(wǎng)絡(luò)模型；Ling等[11]利用張量不變基將伽利略不變性嵌入到預(yù)測(cè)的各向異性張量中，證明滿足伽利略不變性的神經(jīng)網(wǎng)絡(luò)在部分條件下具有更好的預(yù)測(cè)精度。而在稀薄流域，李廷偉等[12]采用極端隨機(jī)樹模型對(duì)非線性本構(gòu)關(guān)系的復(fù)雜高維非線性回歸問題進(jìn)行了建模，對(duì)基于數(shù)據(jù)驅(qū)動(dòng)非線性本構(gòu)方程數(shù)值計(jì)算方法(DNCR)進(jìn)行了初步驗(yàn)證。Zhang等[13]基于DSMC數(shù)據(jù)利用稀疏回歸方法推導(dǎo)宏觀控制方程，證明數(shù)據(jù)驅(qū)動(dòng)與分子模擬結(jié)合具有廣闊的發(fā)展前景。由于研究剛剛起步，這些適用于稀薄非平衡流動(dòng)問題的計(jì)算方法目前僅應(yīng)用于簡單的流動(dòng)問題，尚未擴(kuò)展到較全面和復(fù)雜的多尺度稀薄非平衡流動(dòng)。

在各種機(jī)器學(xué)習(xí)模型中，決策樹是一種基本的分類與回歸方法[14]。決策樹算法的概念最初可以追溯到1966年提出的CLS學(xué)習(xí)系統(tǒng)，而可用于回歸預(yù)測(cè)的CART算法在1984年由Breiman等提出。2001年，Breiman[15]繼續(xù)在此基礎(chǔ)上提出隨機(jī)森林算法，其本質(zhì)相當(dāng)于分類樹的組合，通過幾個(gè)模型降低泛化誤差，針對(duì)數(shù)值屬性的特征，將不同特征的樣本分配到對(duì)應(yīng)的分支。由于隨機(jī)森林的計(jì)算效率、精度等綜合性能較好，在如環(huán)境保護(hù)[16]、推薦系統(tǒng)[17]、道路交通情況預(yù)測(cè)[18]、結(jié)構(gòu)裂紋監(jiān)測(cè)[19]、動(dòng)力裝置可靠性分析[20]等許多領(lǐng)域都有應(yīng)用。對(duì)比隨機(jī)森林，極端隨機(jī)樹[21]則在隨機(jī)提取樣本的基礎(chǔ)上隨機(jī)提取特征，因而隨機(jī)性與抗干擾能力更強(qiáng)。

綜上所述，現(xiàn)有DNCR方法采用了極端隨機(jī)樹方法，直接提取流場(chǎng)物理量與物理量梯度作為特征值與標(biāo)記值，然而這些物理量大部分天然不具備伽利略不變性，使得所訓(xùn)練模型的預(yù)測(cè)范圍嚴(yán)重依賴于訓(xùn)練集坐標(biāo)系的固有方向，限制了模型的泛化性能。針對(duì)上述問題，本文在原始DNCR方法的基礎(chǔ)上提出了運(yùn)用旋轉(zhuǎn)不變量進(jìn)行訓(xùn)練與學(xué)習(xí)的改進(jìn)DNCR方法。改進(jìn)DNCR方法仍以N-S方程與UGKS方法作為理論基礎(chǔ)與數(shù)據(jù)來源，采用旋轉(zhuǎn)不變量代替特征值與標(biāo)記值中包含空間坐標(biāo)方向信息的流場(chǎng)物理量，使得所得訓(xùn)練模型能夠適用于任意旋轉(zhuǎn)后的坐標(biāo)系框架，增強(qiáng)了模型的泛化能力和遷移學(xué)習(xí)能力。為了驗(yàn)證改進(jìn)方法的有效性，本文以不同稀薄來流條件和不同幾何模型下的頂蓋方腔驅(qū)動(dòng)流和二維圓柱/橢圓柱繞流為例，開展了特征參數(shù)選取、模型參數(shù)調(diào)優(yōu)、泛化性能評(píng)估等工作，并評(píng)估改進(jìn)方法與原有DNCR方法的計(jì)算精度提升情況。計(jì)算結(jié)果表明：采用旋轉(zhuǎn)不變量進(jìn)行訓(xùn)練和學(xué)習(xí)后，改進(jìn)DNCR方法在模型的泛化性能和流場(chǎng)物理量預(yù)示精度上明顯優(yōu)于采用極端隨機(jī)樹直接學(xué)習(xí)流場(chǎng)物理量的原始DNCR方法。

1 DNCR計(jì)算方法與旋轉(zhuǎn)不變性修正

數(shù)據(jù)驅(qū)動(dòng)非線性本構(gòu)方程數(shù)值計(jì)算方法的基本原理是將N-S方程和UGKS方法的計(jì)算結(jié)果作為流場(chǎng)樣本訓(xùn)練數(shù)據(jù)，采用機(jī)器學(xué)習(xí)方法構(gòu)建熱流/應(yīng)力項(xiàng)與流場(chǎng)特征參數(shù)的高維復(fù)雜非線性回歸關(guān)系模型，對(duì)N-S方程本構(gòu)關(guān)系進(jìn)行非線性修正，最后通過求解耦合了數(shù)據(jù)驅(qū)動(dòng)非線性回歸關(guān)系模型的宏觀守恒方程得到待預(yù)測(cè)狀態(tài)的稀薄非平衡流數(shù)值解[12]。三維N-S流動(dòng)控制方程在笛卡爾直角坐標(biāo)系下的具體表達(dá)式為

(1)

式中：Q為守恒量；E、F、G為無黏通量；Ev、Fv、Gv為黏性通量。其黏性應(yīng)力張量和熱流項(xiàng)的表達(dá)式為

(2)

隨著氣體稀薄程度變高，努森數(shù)增大，基于連續(xù)介質(zhì)假設(shè)的N-S流動(dòng)控制方程逐漸失效。在稀薄非平衡流域，流動(dòng)問題主要通過Boltzmann方程來描述，該方程運(yùn)用概率速度分布函數(shù)來描述稀薄氣體分子的運(yùn)動(dòng)狀態(tài)，是氣體動(dòng)理學(xué)的基礎(chǔ)控制方程[22]。UGKS采用由BGK[23]簡化的Boltzmann模型方程計(jì)算氣體分子的遷移和碰撞過程：

e-(t-t′)/τdt′+e-t/τf(x-ut,u,0)

(3)

式中：分布函數(shù)f為坐標(biāo)x、分子速度u和時(shí)間t的函數(shù)，用來表示t時(shí)刻以速度u到達(dá)坐標(biāo)x的分子數(shù)密度；g為平衡態(tài)分布函數(shù)，依賴于宏觀物理量；τ為松弛時(shí)間，量級(jí)和分子平均碰撞時(shí)間相同。然而，UGKS將分子運(yùn)動(dòng)和碰撞過程耦合在一起的同時(shí)，也使用了額外的速度離散空間，使得在面對(duì)復(fù)雜高速流動(dòng)問題時(shí)收斂較慢，計(jì)算效率較低。

以二維問題為例，對(duì)熱流/應(yīng)力項(xiàng)進(jìn)行如下非線性修正：

(4)

將式(4)代入質(zhì)量、動(dòng)量與能量守恒方程中進(jìn)行迭代計(jì)算至收斂，最終獲得DNCR方法的流場(chǎng)預(yù)示結(jié)果。由于非線性修正式(4)的約束，最終收斂得到的流場(chǎng)解將趨近于訓(xùn)練模型給出的應(yīng)力、熱流與流場(chǎng)梯度。

在原始DNCR方法中，由于特征值和標(biāo)記值均直接選取流場(chǎng)的物理量與物理量梯度，不具備天然的伽利略不變性。在面對(duì)不同的待預(yù)測(cè)狀態(tài)與幾何模型時(shí)，只有當(dāng)待預(yù)測(cè)狀態(tài)的坐標(biāo)系方向與訓(xùn)練集坐標(biāo)系方向相同時(shí)，才能得到較為理想的預(yù)測(cè)結(jié)果。而當(dāng)坐標(biāo)系發(fā)生諸如旋轉(zhuǎn)變化時(shí)，預(yù)測(cè)結(jié)果可能無法收斂，導(dǎo)致訓(xùn)練模型適用范圍受限，泛化性能亟待提升。

為了使訓(xùn)練得到的機(jī)器學(xué)習(xí)模型能夠在任意旋轉(zhuǎn)坐標(biāo)系下適用，提高模型的泛化性能，拓展模型的適用范圍。本文對(duì)原方法中的特征值和標(biāo)記值進(jìn)行了優(yōu)化與旋轉(zhuǎn)不變性處理，其中DNCR方法的計(jì)算流程如圖1所示。

圖1 DNCR方法計(jì)算流程圖

(5)

由于應(yīng)力Δτij為對(duì)稱張量，可以將其表示為Δτij=VΛVT形式，其中Λ為特征值λ1、λ2、λ3構(gòu)成的對(duì)角陣，代表應(yīng)力Δτij的主應(yīng)力方向，不隨坐標(biāo)變換而變化，由此得到3個(gè)標(biāo)記值λ1、λ2、λ3；V為特征值對(duì)應(yīng)的特征向量v1、v2、v3構(gòu)成的矩陣，特征值對(duì)應(yīng)的特征向量兩兩正交，并且在此問題下表現(xiàn)為旋轉(zhuǎn)矩陣，即|V|=-1，V記錄了當(dāng)前坐標(biāo)系的方位信息。由于矩陣無法直接作為標(biāo)記值參與訓(xùn)練和預(yù)測(cè)，將旋轉(zhuǎn)矩陣變換為四元數(shù)(具體變換過程見附錄A)，因而得到4個(gè)標(biāo)記值q0、q1、q2、q3，代表應(yīng)力Δτij的坐標(biāo)方向信息。速度梯度的合成量采取與應(yīng)力同樣的處理方式。

在預(yù)測(cè)過程中，對(duì)N-S計(jì)算收斂得到的KnP、Knρ、udp/dx等23個(gè)旋轉(zhuǎn)不變量作為機(jī)器學(xué)習(xí)模型輸入；將UGKS與N-S的應(yīng)力、熱流、速度梯度、溫度梯度的張量特征值差量λi和由N-S特征向量變換為UGKS特征向量的旋轉(zhuǎn)矩陣Vq(Vq由四元數(shù)qj表示，具體變換過程詳見附錄A)作為模型輸出。以應(yīng)力張量為例，對(duì)DNCR續(xù)算所需的Δτij，有

(6)

2 機(jī)器學(xué)習(xí)模型選擇與參數(shù)調(diào)優(yōu)

改進(jìn)DNCR方法使用極端隨機(jī)樹作為學(xué)習(xí)模型，極端隨機(jī)樹本質(zhì)上是多個(gè)決策樹的集合，決策樹是一種基本的分類與回歸方法。決策樹的生成就是遞歸地構(gòu)建二叉決策樹的過程，回歸樹采用平方誤差最小化準(zhǔn)則。每一棵回歸樹對(duì)應(yīng)著輸入空間(即特征空間)的一個(gè)劃分以及在劃分的單元上的輸出值。假設(shè)已將輸入空間劃分為M個(gè)單元，并且在每個(gè)單元 上有一個(gè)固定的輸出值，于是回歸樹模型可表示為[14]

(7)

在輸入空間的劃分確定后，一般采用平方誤差來表示回歸樹對(duì)訓(xùn)練集的預(yù)測(cè)誤差，用使平方誤差最小的準(zhǔn)則來求解每個(gè)單元的最優(yōu)輸出。通常使用的最小二乘回歸樹算法，在訓(xùn)練數(shù)據(jù)集的輸入空間中，遞歸的將每個(gè)區(qū)域分為兩個(gè)子區(qū)域并決定其輸出值，來構(gòu)建二叉決策樹：

1) 選擇最優(yōu)切分變量j與切分點(diǎn)s，求解

(8)

遍歷變量，對(duì)固定的切分變量 掃描切分點(diǎn)，選擇使式(8)達(dá)到最小值的對(duì)。

2) 用選定的劃分區(qū)域并決定相應(yīng)的輸出值：

(9)

(10)

3) 繼續(xù)對(duì)兩個(gè)字區(qū)域調(diào)用步驟1)和步驟2)，直至滿足停止條件。

4) 將輸入空間劃分為M個(gè)區(qū)域R1,R2,…,RM，生成決策樹：

(11)

極端隨機(jī)樹作為二叉決策樹的一種集成算法，對(duì)所有訓(xùn)練樣本直接隨機(jī)選擇分叉特征得到結(jié)果，實(shí)現(xiàn)簡單，不容易產(chǎn)生過擬合，易找到對(duì)結(jié)果影響較大的特征，并且可并行化、速度快。

本文對(duì)基學(xué)習(xí)器數(shù)目以及最大特征數(shù)進(jìn)行調(diào)優(yōu)，通過可決系數(shù)的大小對(duì)兩個(gè)超參數(shù)進(jìn)行尋優(yōu)，可決系數(shù)公式為

(12)

圖2 基學(xué)習(xí)器數(shù)目調(diào)優(yōu)

圖3 最大特征數(shù)調(diào)優(yōu)

3 算例與分析

UGKS方法的數(shù)據(jù)作為訓(xùn)練集中最重要數(shù)據(jù)來源之一，需要確保能夠?qū)ο”》瞧胶饬鲌?chǎng)進(jìn)行較為準(zhǔn)確的描述，因此首先與文獻(xiàn)中DSMC[25]數(shù)據(jù)進(jìn)行對(duì)比，驗(yàn)證UGKS方法的準(zhǔn)確性。計(jì)算采用單原子氬氣圓柱繞流，來流條件Ma=5，Kn=0.1，來流分子數(shù)密度為n=1.294 4×1021/m3，黏性系數(shù)μ∞=2.117×10-5Pa·s，壁面溫度為Tw=273 K，壁面與DSMC保持一致均設(shè)置為漫反射邊界條件。UGKS方法物理空間網(wǎng)格剖分為50×64，速度空間離散網(wǎng)格剖分為93×93。

圖4和圖5給出了UGKS和DSMC兩種方法沿圓柱表面從駐點(diǎn)到后緣的物面應(yīng)力θ、熱流q結(jié)果對(duì)比曲線。其中橫坐標(biāo)為圓柱表面位置與圓心連線與y軸夾角角度，縱坐標(biāo)分別為表面熱流與黏性應(yīng)力值。從計(jì)算結(jié)果可以看出，UGKS結(jié)果與DSMC結(jié)果吻合較好，驗(yàn)證了UGKS計(jì)算模塊的可靠性。

圖4 UGKS與DSMC熱流對(duì)比

圖5 UGKS與DSMC應(yīng)力對(duì)比

本文所發(fā)展的改進(jìn)DNCR方法采用旋轉(zhuǎn)不變量作為特征值和標(biāo)記值，在任意方向的坐標(biāo)系下均能夠得到相同的流場(chǎng)結(jié)果，與實(shí)際物理現(xiàn)象吻合。原有DNCR方法直接采用流場(chǎng)物理量作為機(jī)器學(xué)習(xí)模型特征值和標(biāo)記值，對(duì)訓(xùn)練集坐標(biāo)系發(fā)生旋轉(zhuǎn)、平移后的待預(yù)測(cè)狀態(tài)無法開展準(zhǔn)確預(yù)測(cè)。為了便于方法之間對(duì)比，算例僅展示固定坐標(biāo)系下的對(duì)比結(jié)果。需要說明的是，改進(jìn)DNCR方法針對(duì)任意坐標(biāo)系(旋轉(zhuǎn)、平移)后的待預(yù)測(cè)狀態(tài)都能得到相同的流場(chǎng)計(jì)算結(jié)果。

為了反映新的特征參數(shù)構(gòu)造所帶來的模型泛化性能提升，選取二維圓柱和二維方腔兩種外形進(jìn)行驗(yàn)證和對(duì)比，并對(duì)部分待預(yù)測(cè)狀態(tài)進(jìn)行適當(dāng)?shù)膸缀瓮庑巫兓?。二維圓柱采用直徑0.304 8 m(12 in)下的數(shù)據(jù)為訓(xùn)練集，分別預(yù)測(cè)直徑 0.304 8 m下的二維圓柱外形，長軸0.609 6 m、短軸0.304 8 m的二維橢圓柱外形兩種幾何形狀在不同Kn數(shù)下的流場(chǎng)結(jié)果。二維方腔以外形 1 m×1 m的數(shù)據(jù)為訓(xùn)練集，分別預(yù)測(cè)1 m×1 m，2 m×1 m兩種幾何外形方腔在不同Kn數(shù)下的流場(chǎng)結(jié)果。最終所得結(jié)果的精度提升采用式(14)進(jìn)行評(píng)估：

(13)

3.1 二維圓柱氬氣高超聲速繞流

首先驗(yàn)證二維圓柱外形下的預(yù)測(cè)結(jié)果，計(jì)算采用單原子氣體氬氣，來流馬赫數(shù)與來流溫度分別為Ma=10,T0=198.439 K，氣體的物性參數(shù)取μ0=5.071 158×10-5Pa·s，γ=1.667，Pr=0.667，R=208.16 J/(kg·K)，逆冪律冪次ε=0.734，UGKS計(jì)算模型采用簡單硬球模型。模型的訓(xùn)練數(shù)據(jù)來自Kn=0.73與Kn=1.09這2組狀態(tài)，網(wǎng)格大小為33×91(周向×徑向)，圓柱直徑為0.304 8 m，第1層網(wǎng)格間距為1×-3m。

第1組待預(yù)測(cè)稀薄非平衡狀態(tài)為Kn=0.91同尺寸二維圓柱繞流，除來流密努森數(shù)變化外，其余氣體屬性相同，計(jì)算網(wǎng)格如圖6所示。圖7與圖8給出了不同計(jì)算方法溫度云圖對(duì)比,DNCR-RI(DNCR-Rotarional Invariants)為當(dāng)前采用旋轉(zhuǎn)不變量的改進(jìn)DNCR方法。為了更好的表現(xiàn)N-S、UGKS、改進(jìn)DNCR流場(chǎng)結(jié)果的差異，截取圖中網(wǎng)格對(duì)稱軸y=-x截線(45°截線)，對(duì)比截線上各項(xiàng)流場(chǎng)物理量的結(jié)果，如圖9～圖12所示。

圖6 直徑0.304 8 m二維圓柱網(wǎng)格示意圖

圖7 N-S與UGKS圓柱溫度對(duì)比

圖8 DNCR-RI與UGKS圓柱溫度對(duì)比

由計(jì)算結(jié)果可以看出，在來流Kn=0.91條件下，圓柱繞流頭部弓形激波附近區(qū)域表現(xiàn)出明顯稀薄非平衡特征，激波厚度由于平均分子自由程增大明顯變大。與UGKS與DNCR方法相比，N-S方程預(yù)示得到的激波位置更加靠近物面，激波脫體距離更小，表明N-S方程的連續(xù)性假設(shè)在這個(gè)來流條件下部分失效，無法對(duì)稀薄非平衡流場(chǎng)進(jìn)行準(zhǔn)確描述，而UGKS與DNCR流場(chǎng)結(jié)構(gòu)較為一致。對(duì)比N-S與UGKS計(jì)算結(jié)果，采用精度提升式(13)分別考察DNCR計(jì)算精度，可以得到壓力精度提升49.09%，溫度精度提升61.02%，x軸方向速度U精度提升43.85%，y軸方向速度V精度提升48.26%，表明現(xiàn)有DNCR方法在相同外形、不同Kn數(shù)下有較好的預(yù)測(cè)表現(xiàn)。

圖9 DNCR與UGKS圓柱截線y=-x上壓力對(duì)比

圖10 DNCR與UGKS圓柱截線y=-x上溫度對(duì)比

圖11 DNCR與UGKS圓柱截線y=-x上速度U對(duì)比

圖12 DNCR與UGKS圓柱截線y=-x上速度V對(duì)比

圖13和圖14反映了N-S、UGKS、DNCR 3種計(jì)算方式得到的摩擦阻力系數(shù)Cf和表面熱流對(duì)比，可以看到DNCR方法得到的結(jié)果與UGKS十分吻合，表明DNCR方法在氣動(dòng)熱和摩阻預(yù)測(cè)方面也呈現(xiàn)了較好的預(yù)測(cè)能力。

為了測(cè)試改進(jìn)DNCR方法模型的泛化性能，針對(duì)外形改變的情況進(jìn)行了預(yù)測(cè)和結(jié)果對(duì)比。第2組待預(yù)測(cè)狀態(tài)采用Kn=0.91，外形為長軸0.609 6 m、短軸0.304 8 m的二維橢圓柱，其余各項(xiàng)設(shè)置相同，計(jì)算網(wǎng)格如圖15所示，重點(diǎn)考察進(jìn)DNCR方法模型在變外形條件下泛化性能提升情況。圖16和圖17給出了不同計(jì)算方法溫度云圖對(duì)比。

圖13 圓柱表面沿x方向摩擦阻力系數(shù)

圖18～圖21中DNCR-RV(DNCR-Rotational Variants)為原有采用旋轉(zhuǎn)可變物理量的DNCR方法，可以明顯看出改進(jìn)DNCR方法得到的結(jié)果趨勢(shì)與UGKS結(jié)果更吻合。

若采用精度提升式(13)進(jìn)行量化分析，改進(jìn)DNCR方法對(duì)壓力的精度提升為-10.01%，對(duì)溫度的精度提升為71.19%，速度U的精度提升為63.80%，速度V的精度提升為53.36%。而原始DNCR方法的壓力精度提升為-90.52%，溫度精度提升為71.14%，速度U精度提升為57.07%，速度V精度提升為-42.76%。通過精度提升對(duì)比可以發(fā)現(xiàn)，當(dāng)待預(yù)示幾何外形發(fā)生變化時(shí)，原始DNCR方法以含有空間坐標(biāo)系方向信息的數(shù)據(jù)作為訓(xùn)練集，得到的預(yù)測(cè)結(jié)果精度提升不如采用基于旋轉(zhuǎn)不變量改進(jìn)DNCR方法，尤其是對(duì)諸如壓力與y方向速度V的預(yù)示。計(jì)算結(jié)果初步表明，改進(jìn)DNCR方法擁有更強(qiáng)的模型泛化性能，針對(duì)變幾何構(gòu)型的流場(chǎng)預(yù)示能力明顯增強(qiáng)。

圖14 圓柱表面沿x方向熱流分布

圖15 長軸0.609 6 m、短軸0.304 8 m二維橢圓柱網(wǎng)格示意圖

圖16 N-S與UGKS橢圓柱溫度對(duì)比示意圖

圖17 DNCR-RI與UGKS橢圓柱溫度對(duì)比示意圖

圖18 橢圓柱流場(chǎng)截線y=-x上壓力對(duì)比

圖19 橢圓柱流場(chǎng)y=-x截線溫度對(duì)比

圖20 橢圓柱流場(chǎng)截線y=-x速度對(duì)比

圖21 橢圓柱流場(chǎng)y=-x截線速度對(duì)比

3.2 頂蓋方腔驅(qū)動(dòng)流

為進(jìn)一步驗(yàn)證改進(jìn)DNCR方法的計(jì)算能力與泛化性能，對(duì)二維頂蓋方腔驅(qū)動(dòng)流進(jìn)行了計(jì)算與分析。計(jì)算采用單原子氣體氬氣，初始流場(chǎng)溫度T0=300 K，氣體物性參數(shù)取μ0=2.272×10-5Pa·s，γ=1.667，Pr=0.667，R=208.16 J/(kg·K)，逆冪律冪次取ε=0.81，UGKS計(jì)算模型采用簡單硬球模型。模型的訓(xùn)練數(shù)據(jù)來自Kn=0.003 2、Re=44.0與Kn=0.002 375、Re=59.3這2組狀態(tài)，物理計(jì)算網(wǎng)格為61×61的均勻網(wǎng)格，速度離散空間網(wǎng)格為28×28均勻分布，計(jì)算域尺寸為1 m×1 m。

第1組待預(yù)測(cè)狀態(tài)Kn=0.002 85、Re=49.4，幾何尺寸保持不變，同為1 m×1 m的方腔，其余各項(xiàng)設(shè)置相同，計(jì)算網(wǎng)格如圖22所示。圖23～圖26分別展示了x=0.5 m截線處壓力、溫度、速度U、速度V的對(duì)比。

圖22 1 m×1 m方腔示意圖 Fig.22 Schematic diagram of a 1 m×1 m squarecavity

圖23 1 m×1 m方腔截線x=0.5 m上壓力

圖24 1 m×1 m方腔截線x=0.5 m上溫度

圖25 1 m×1 m方腔截線x=0.5 m上速度U

圖26 1 m×1 m方腔截線x=0.5 m上速度V

為量化DNCR計(jì)算精度提升情況，借助精度提升計(jì)算式(13)，得到x=0.5 m處壓力精度提升68.15%，溫度精度提升55.89%，速度U精度提升31.74%，速度V精度提升-20.38%，速度V截線處絕對(duì)值較小，由于特征值和標(biāo)記值不涉及方向，所以預(yù)示結(jié)果較差，而絕對(duì)值較大的速度U結(jié)果明顯較好，當(dāng)?shù)伛R赫數(shù)精度提升了31.66%，受速度V影響較小。圖27～圖30給出了方腔流場(chǎng)截線y=0.5 m上流場(chǎng)物理量分布曲線。

圖27 1 m×1 m方腔截線y=0.5 m上壓力

圖28 1 m×1 m方腔截線y=0.5 m上溫度

圖29 1 m×1 m方腔截線y=0.5 m上速度U

如圖27～圖30所示，y=0.5 m截線上壓力精度提升77.47%，溫度精度提升58.66%，速度U精度提升69.43%，速度V精度提升44.49%，整體預(yù)示結(jié)果較好。

圖30 1 m×1 m方腔截線y=0.5 m上速度V

為了考察方法的泛化性能，對(duì)變外形情形進(jìn)行了預(yù)測(cè)和對(duì)比。第2組待預(yù)測(cè)狀態(tài)為Kn=0.002 85，Re=49.4，計(jì)算幾何外形為2 m×1 m的方腔，其余各項(xiàng)設(shè)置相同，網(wǎng)格如圖31所示。截取相對(duì)位置與第1組算例相同，分別給出x=1 m、y=0.5 m處的壓力、溫度、速度U、速度V結(jié)果，其中x=1 m截線計(jì)算結(jié)果如圖32～圖35所示。

圖31 2 m×1 m方腔示意圖

圖32 2 m×1 m方腔截線x=1 m上壓力

圖33 2 m×1 m方腔截線x=1 m上溫度

圖34 2 m×1 m方腔截線x=1 m上速度U

在x=1 m處壓力精度提升80.32%，溫度精度提升45.12%，速度U精度提升36.57%，速度V精度提升39.65%。y=0.5 m的截線結(jié)果如圖36～圖39所示。

在y=0.5 m處壓力精度提升54.63%，溫度精度提升4.49%，速度U精度提升42.49%，速度V精度提升3.72%。

圖35 2 m×1 m方腔截線x=1 m上速度V

圖36 2 m×1 m方腔截線y=0.5 m上壓力

圖37 2 m×1 m方腔截線y=0.5 m上溫度

圖38 2 m×1 m方腔截線y=0.5 m上速度U

圖39 2 m×1 m方腔截線y=0.5 m上速度V

通過對(duì)固定坐標(biāo)系不同狀態(tài)下二維圓柱與方腔算例的驗(yàn)證，從流場(chǎng)關(guān)鍵位置截線對(duì)比壓力、溫度、速度分布可以看出：改進(jìn)DNCR方法對(duì)不同來流狀態(tài)、訓(xùn)練集相同外形、變幾何外形的流場(chǎng)參數(shù)均能作出較好的預(yù)測(cè)，且當(dāng)外形發(fā)生改變時(shí)，預(yù)測(cè)效果優(yōu)于原有DNCR方法，表明改進(jìn)模型擁有更好的泛化性能。更為重要的是，若坐標(biāo)系發(fā)生旋轉(zhuǎn)時(shí)，改進(jìn)DNCR方法的流場(chǎng)結(jié)果與提升效果保持不變，原始DNCR方法則囿于特征值和標(biāo)記值不具備旋轉(zhuǎn)不變性而無法進(jìn)行預(yù)測(cè)，因此未在文中進(jìn)行量化比較，這也從另一方面反映了改進(jìn)DNCR特征量具備旋轉(zhuǎn)不變性的重要意義。

4 結(jié) 論

本文在原始DNCR方法的基礎(chǔ)上發(fā)展了基于旋轉(zhuǎn)不變量的改進(jìn)DNCR方法，對(duì)機(jī)器學(xué)習(xí)模型的特征值和標(biāo)記值進(jìn)行了重新篩選與處理，使其滿足旋轉(zhuǎn)不變性約束，最終使得訓(xùn)練模型的預(yù)示能力與范圍得到增強(qiáng)，解決了原始DNCR方法不適用于旋轉(zhuǎn)坐標(biāo)后的流場(chǎng)精確預(yù)示這一難題。改進(jìn)DNCR方法的計(jì)算流程與原始方法一致，即將修正應(yīng)力與熱流預(yù)測(cè)結(jié)果代入守恒方程進(jìn)行迭代并計(jì)算至收斂，在模型訓(xùn)練完成后，DNCR總計(jì)算時(shí)長等于N-S方程計(jì)算時(shí)長+模型預(yù)測(cè)時(shí)長+DNCR續(xù)算時(shí)長，預(yù)測(cè)時(shí)長相較N-S與DNCR計(jì)算時(shí)長可忽略不計(jì)，因此總計(jì)算時(shí)長與N-S方程處于同一量級(jí)，但計(jì)算精度較N-S方程得到大幅提升。以固定坐標(biāo)系下的氬氣二維高超聲速圓柱繞流和頂蓋方腔驅(qū)動(dòng)流為例，針對(duì)原有外形和變幾何外形的待預(yù)示狀態(tài)，將改進(jìn)的DNCR方法結(jié)果與N-S方程、UGKS方法及原有DNCR方法進(jìn)行了對(duì)比，結(jié)果表明使用旋轉(zhuǎn)不變量能夠顯著提升訓(xùn)練模型對(duì)外形變化的適應(yīng)能力，使其表現(xiàn)出更好的泛化性能，進(jìn)一步提高了DNCR方法的適用場(chǎng)景與應(yīng)用范圍。