李廷偉，張莽，趙文文，陳偉芳，蔣勵劍

1. 浙江大學 航空航天學院，杭州 310027 2.中國電子科技集團公司第五十四研究所 航天系統(tǒng)與應(yīng)用專業(yè)部，石家莊 050081 3.中國運載火箭技術(shù)研究院 研究與發(fā)展中心，北京 100076

錢學森根據(jù)努森數(shù)Kn定義流體的稀薄程度，將流動區(qū)域劃分為連續(xù)流域(Kn<0.01)、滑移流域(0.0110)，稀薄非平衡流域包含除連續(xù)流域以外的3大流域[1]。

連續(xù)流域符合流體力學連續(xù)介質(zhì)假設(shè)條件，連續(xù)介質(zhì)流體力學是其理論基礎(chǔ)，其中最具有代表性的控制方程為Navier-Stokes(N-S)方程，形成了以廣義牛頓定律和傅里葉熱傳導定律為基礎(chǔ)的基于NSF(Navier-Stokes-Fourier)關(guān)系的連續(xù)流數(shù)值模擬計算方法，此方法的發(fā)展為其他各種方法的建立提供了寶貴經(jīng)驗。稀薄非平衡流域中氣體稀薄屬性逐漸凸顯，例如物面出現(xiàn)了十分顯著的速度滑移與溫度跳躍現(xiàn)象。此時在物面附近的努森層內(nèi)，連續(xù)介質(zhì)假設(shè)已經(jīng)失效，準確地說是基于3大守恒方程(質(zhì)量守恒、動量守恒和能量守恒)推導出來的黏性應(yīng)力項與熱流項已經(jīng)不能再簡單使用低階宏觀物理量(速度、溫度)的梯度線性表征，即線性本構(gòu)關(guān)系已經(jīng)不再適用于精準描述稀薄非平衡流動問題。

對稀薄非平衡流問題的研究主要圍繞稀薄氣體流動控制方程——Boltzmann方程開展，它是分子氣體動力學的基本方程，可以不受努森數(shù)的限制對連續(xù)流到自由分子流整個流域進行統(tǒng)一描述。Boltzmann方程是一個復(fù)雜的七維積分/微分方程，大部分研究均是對其直接或者間接求解，發(fā)展形成了多種數(shù)值計算方法與理論，統(tǒng)一氣體動理論格式(Unified Gas-Kinetic Scheme， UGKS)[2]是其中一種代表性方法。UGKS方法使用Boltzman方程的Bhatnagar-Gross-Krook(BGK)[2]類模型方程作為控制方程，用有限個離散速度替代整個速度空間，在一定條件下均能給出較為準確的流動特征，具備求解各流域多尺度問題的能力[3]。UGKS方法采用BGK類模型方程基于當?shù)胤e分解對通量進行構(gòu)造，將分子運動與碰撞過程自然地耦合在一起，物理意義更加明確，突破了DSMC方法[4]解耦處理所帶來的計算時間步長小于分子的平均碰撞時間、計算單元網(wǎng)格尺度小于分子平均自由程的限制。但由于UGKS方法在求解過程中依賴速度分布函數(shù)對速度空間進行離散，這對計算與存儲的要求非常高，多維與高速條件下計算效率較低。

隨著計算機理論與技術(shù)的迅猛發(fā)展，大數(shù)據(jù)時代啟發(fā)了人類思考問題新的思維，基于數(shù)據(jù)驅(qū)動和機器學習的研究方式也應(yīng)運而生。研究人員開始通過機器學習的方法解決流體力學領(lǐng)域中難以解決的問題，目前在氣動優(yōu)化設(shè)計[5-9]、湍流問題[10-15]、非定常氣動力建模[16-20]等領(lǐng)域上已經(jīng)開展了很多卓有成效的研究工作。

以N-S方程和UGKS方法為理論基礎(chǔ)，提出了一種適用于稀薄非平衡流數(shù)值模擬的基于數(shù)據(jù)驅(qū)動非線性本構(gòu)關(guān)系(Data-driven Nonlinear Constitutive Relations, DNCR)的數(shù)值計算方法?；诹鲌鎏卣鲄?shù)采用機器學習方法對N-S方程線性的黏性應(yīng)力項與熱流項進行非線性離散重構(gòu)，理論適用范圍可以覆蓋較大來流努森數(shù)條件，通過耦合非線性本構(gòu)關(guān)系求解宏觀守恒方程得到待預(yù)測狀態(tài)稀薄非平衡流動數(shù)值解。與傳統(tǒng)機器學習建模方法摒棄基本物理規(guī)律直接對數(shù)據(jù)進行訓練與預(yù)測的思想相比，DNCR方法未拋棄物理規(guī)律而是對流體本構(gòu)關(guān)系進行修正后耦合守恒方程迭代求解，同時保留了數(shù)據(jù)建模與物理建模的優(yōu)點，方法物理意義更加清晰明確。

在DNCR方法中，采用的具體機器學習算法為極端隨機樹[21-22](Extremely Randomized Trees)算法。在機器學習算法中，針對所研究的具體問題選取低冗余、高代表性的特征參數(shù)對于機器學習算法的泛化性能和運行效率具有重要影響。本文擬通過二維頂蓋驅(qū)動方腔流算例對高維非線性建模涉及的特征參數(shù)選取、參數(shù)調(diào)優(yōu)開展相關(guān)驗證與研究工作，選取若干典型狀態(tài)對極端隨機樹模型的泛化性能開展研究，并評估相關(guān)模型與方法的計算精度與計算效率。

1 數(shù)據(jù)驅(qū)動非線性本構(gòu)關(guān)系計算方法

DNCR基于N-S方程和UGKS方法的數(shù)值模擬計算結(jié)果作為流場樣本訓練數(shù)據(jù)，采用機器學習方法構(gòu)建熱流/應(yīng)力項與流場特征參數(shù)的高維復(fù)雜非線性回歸關(guān)系模型，對N-S方程本構(gòu)關(guān)系進行非線性修正。通過耦合數(shù)據(jù)驅(qū)動的非線性本構(gòu)關(guān)系求解宏觀量守恒方程得到稀薄非平衡流數(shù)值解。

DNCR方法計算流程如圖1所示，ψ代表在N-S流場數(shù)值模擬結(jié)果中提取的特征參數(shù)，作為機器學習模型的輸入量，ΔΘ為標記參數(shù)，作為機器學習模型的輸出量，即機器學習模型建立了ψ與ΔΘ之間的復(fù)雜高維回歸關(guān)系，與非線性本構(gòu)物理建模函數(shù)的區(qū)別在于其本構(gòu)關(guān)系沒有具體數(shù)學表達式，而是全流域存在離散當?shù)赜成洇住う?，所形成的離散本構(gòu)映射關(guān)系在全流場計算域適用。而模型泛化性能取決于訓練集特征代表性與特征參數(shù)選取，理論上若特征參數(shù)選擇未涉及到任何與空間網(wǎng)格與當?shù)刈鴺酥抵苯酉嚓P(guān)的參數(shù)且訓練集包含相似的非平衡流動特征，機器學習模型就具備一定的遷移預(yù)測能力。

圖1 DNCR示意圖Fig.1 Schematic of DNCR

DNCR方法的一個重要特點是訓練與預(yù)測過程相對獨立，圖1中紅線與綠線分別表示了訓練與預(yù)示流程。模型訓練過程通過對包含不同典型流動特征的基礎(chǔ)流場數(shù)據(jù)集開展訓練，獲得ψ與ΔΘ之間的復(fù)雜回歸關(guān)系。而預(yù)測過程則首先采用N-S求解器對待預(yù)測狀態(tài)開展初步計算，獲得待預(yù)測當?shù)靥卣髦郸祝缓蟛捎靡延柧毶傻幕貧w關(guān)系ψ→ΔΘ得到待預(yù)測流場當?shù)貥擞浿郸う?，最終通過時間推進方式求解質(zhì)量、動量與能量守恒方程，耦合非線性本構(gòu)關(guān)系完成計算，守恒方程為

(1)

式中：Q為守恒量；E、F為無黏通量；Ev、Fv為黏性通量。控制方程中黏性應(yīng)力張量與熱流項表達式為

(2)

式中：τ為黏性應(yīng)力項；μ為黏性系數(shù)；T為溫度；u、v為速度；κ為熱傳導系數(shù)；q為熱流項；下標Tag表示待預(yù)測流場當?shù)貥擞浿?。?1)與式(2)隨著時間推進，相關(guān)熱流、應(yīng)力張量與流場梯度量向標記值趨近，最終完成計算收斂，獲得待預(yù)測流場定常解。

當物體或擾動源的速度大于擾動信息的傳播時，擾動無法向上游傳播所形成的強壓縮間斷稱為激波。一維激波結(jié)構(gòu)是最簡單、最基本的非平衡流動現(xiàn)象之一，是研究本構(gòu)關(guān)系與非平衡流動的典型算例，可以用來對模型進行驗證[23]。

N-S、UGKS、DNCR相同網(wǎng)格下，采用單原子氣體Ar作為模擬氣體，其動力黏性系數(shù)μ采用逆冪律公式[24]計算：

(3)

單原子氣體Ar物性參數(shù)取μ0=2.272×10-5Pa·s，T0=300 K，γ=1.667，Pr=0.667，R=208.16 J/(kg·K)，其中μ0為Ar在溫度T0下的動力黏性系數(shù)，T0為Sutherland溫度常量，γ為比熱比，Pr為普朗特數(shù)，R為氣體常數(shù)，逆冪律冪次取ε=0.72。計算狀態(tài)如表1所示。

表1 一維Ar激波結(jié)構(gòu)計算狀態(tài)

圖2 激波結(jié)構(gòu)密度分布Fig.2 Density distributions in shock wave

圖3 激波結(jié)構(gòu)速度分布Fig.3 Velocity distributions in shock wave

圖4 激波結(jié)構(gòu)壓強分布Fig.4 Pressure distributions in shock wave

圖5 激波結(jié)構(gòu)溫度分布Fig.5 Temperature distributions in shock wave

2 訓練數(shù)據(jù)集生成

本文選取圖6所示頂蓋驅(qū)動方腔流動作為測試算例，計算網(wǎng)格如圖7所示，計算網(wǎng)格為61×61均勻網(wǎng)格，速度離散為28×28。選取5組稀薄非平衡流動狀態(tài)生成訓練數(shù)據(jù)集，計算狀態(tài)如表2所示。

圖6 頂蓋驅(qū)動方腔流示意圖Fig.6 Schematic of lid-driven cavity flow

圖7 二維頂蓋驅(qū)動方腔流計算網(wǎng)格Fig.7 Grid of two-dimensional lid-driven cavity flow

表2 二維方腔計算狀態(tài)Table 2 Calculation state of two-dimensional cavity

3 機器學習模型

3.1 CART與極端隨機樹算法

X與Y分別為特征輸入和標記輸出變量，并且Y為連續(xù)變量，給定訓練數(shù)據(jù)集：D={(x1,y1),(x2,y2),…,(xN,yN)}。

回歸樹是用于研究回歸問題的決策樹模型，決策樹是一種樹形結(jié)構(gòu)的基本分類與回歸算法。在機器學習中決策樹是一個預(yù)測模型，代表的是對象屬性(即特征參數(shù))與對象值(即標記值)之間的一種映射關(guān)系。一個回歸樹對應(yīng)著輸入空間(即特征空間)的一個劃分以及在劃分的單元上的輸出值。假設(shè)將訓練數(shù)據(jù)集所在的輸入空間劃分為M個單元：R1,R2,…,RM，這里的單元可以理解為從最頂端往下依次劃分構(gòu)建的M個分支，每個分支只對應(yīng)一個輸出值。遞歸地將每個區(qū)域劃分為兩個子區(qū)域并決定每個子區(qū)域的輸出值c，構(gòu)建二叉決策樹。

1) 選擇最優(yōu)切分變量j與切分點s：特征空間包含多個特征參數(shù)即特征變量，決策樹的構(gòu)建就是每層進行二叉樹的劃分，劃分時采用某特征變量能夠使誤差最小，此變量即為最優(yōu)切分變量；切分點指的是切分變量進行左右劃分的值，小于此值劃分為左子樹，大于此值劃分為右子樹，切分時某切分點能夠使誤差最小即為最優(yōu)切分點。通過求解:

(4)

遍歷變量j，對固定切分變量j掃描切分點s，選擇使式(4)最小的(j,s)。

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

R1(j,s)={x|x(j)≤s}

R2(j,s)={x|x(j)>s}

(5)

3) 繼續(xù)對兩個子區(qū)域調(diào)用步驟1)、2)，直至滿足條件。

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

(6)

極端隨機樹在生成CART過程中選擇劃分屬性不再是選擇最優(yōu)屬性而是完全隨機選擇。采用訓練數(shù)據(jù)集進行并行訓練得到t個CART，對本文回歸問題模型最終采用平均法：

(7)

3.2 特征參數(shù)選擇

在DNCR方法中，流場特征參數(shù)構(gòu)成特征空間{ψ}，熱流、應(yīng)力修正所需各項構(gòu)成標記空間{ΔΘ}。選取低冗余、高代表性的特征參數(shù)對機器學習算法的泛化性能和運行效率具有重要影響。首先從物理機理上考慮，選擇表征流場稀薄非平衡特征的參數(shù)，例如基于當?shù)亓鲌鎏荻鹊木植颗瓟?shù)KnGLL(ρ,P,T)作為主要的流場特征參數(shù)，然后嘗試加入其他參數(shù)并最終選取適當參數(shù)集作為流場特征參數(shù)建立特征空間。

本文采用方差過濾準則進行特征選擇，方差越小，表示該特征的數(shù)據(jù)差異越小，可以認為這個特征對區(qū)分樣本貢獻不大，導致預(yù)測的效果越差，因此可以剔除此特征或降低特征權(quán)重。另外考慮到真實物理問題中不同物理量之間數(shù)據(jù)的量級可能存在較大差異，量級較大導致方差較大，為了避免這一不利因素，使用標準差與極差的商值消除量級影響進行特征選擇。

(8)

式中:X={x1,x2,…,xn}代表某一特征參數(shù);n代表樣本點的數(shù)量;μ為X={x1,x2,…,xn}的算數(shù)平均值。標準差與極差做商的目的在于消除數(shù)據(jù)量綱的影響。

圖8 特征參數(shù)標準差對比Fig.8 Comparison of standard deviation of characteristic parameters

(9)

式中：Ypredict代表預(yù)測的標記值，Yactual代表真實的標記值，均可看作形如Y={y1,y2,…,yn}的向量，yi代表某網(wǎng)格點標記值，n代表計算網(wǎng)格點的數(shù)量。

圖9 余弦相似度對比Fig.9 Comparison of cosine similarity

3.3 極端隨機樹參數(shù)調(diào)優(yōu)

在研究具體問題時，不同機器學習算法中都有許多的參數(shù)需要人為設(shè)定，即使是同一種算法，面向不同回歸問題的參數(shù)配置也有所區(qū)別，最終得到的模型性能表現(xiàn)往往也存在十分明顯的差別。

對于極端隨機樹算法來說，主要參數(shù)包括框架參數(shù)：基學習器的數(shù)量n_estimators、度量分裂的標準criterion、是否采用袋外樣本來評估模型的好壞oob_score；決策樹參數(shù)：決策樹的每個節(jié)點隨機選擇的最大特征數(shù)max_features、葉子結(jié)點上應(yīng)有的最少樣例數(shù)min_samples_leaf、決策樹最大深度max_depth等參數(shù)。本文重點對n_estimators與max_features兩個參數(shù)進行研究，其余參數(shù)均采取默認設(shè)置，研究方法為單一變量法。

本文采用可決系數(shù)[27]評價參數(shù)調(diào)節(jié)后的模型優(yōu)劣情況?？蓻Q系數(shù)表示一個隨機變量與多個隨機變量關(guān)系的數(shù)字特征，即可決系數(shù)值越大越接近于1時說明特征參數(shù)對標記參數(shù)的解釋程度越高。

(10)

采用Kn=0.7與Kn=1.3兩組狀態(tài)對機器學習模型進行訓練，對Kn=1.0狀態(tài)進行預(yù)測。選取n_estimators在[1,1 000]范圍內(nèi)，R2隨n_estimators變化趨勢如圖10所示?？梢钥闯鯮2最終趨近于一個比較穩(wěn)定的值，擬合較好且隨著基學習器的數(shù)量增加未出現(xiàn)嚴重過擬合現(xiàn)象。如圖11所示，模型的訓練時間與基學習器數(shù)目基本滿足線性關(guān)系，考慮到模型訓練的時間成本，本文中選取n_estimators=300。

圖10 R2隨基學習器數(shù)目變化曲線Fig.10 R2 variation curve with n_estimators

圖11 訓練時間隨基學習器數(shù)目變化曲線Fig.11 Train time variation curve with n_estimators

根據(jù)前文中確定使用的12個特征值，max_features在[1,12]范圍中取值。為避免偶然誤差，訓練1 000次后取平均值，得到R2值隨max_features變化趨勢如圖12所示?？梢钥闯鲈贒NCR方法中極端隨機樹的預(yù)測準確性隨著max_features變化大致呈現(xiàn)遞增的趨勢，max_features的大小對模型訓練預(yù)測時間的影響可以忽略不計，在本文中選擇max_features=12。其余模型參數(shù)及其他待預(yù)測物理量的模型參數(shù)調(diào)優(yōu)過程與之類似。

圖12 R2隨最大特征數(shù)變化曲線Fig.12 R2 variation curve with max_features

4 DNCR計算精度與極端隨機樹泛化性能分析

圖13 預(yù)測情況Fig.13 Prediction situation

為更好地反映預(yù)測值與真實值的對比情況，采用均方誤差(Mean Square Error，MSE)和均方根誤差(Root Mean Square Error，RMSE)可以衡量預(yù)測值同真值之間的偏離程度。

(11)

(12)

標記值S1～S11的MSE、RMSE、R2情況如表3所示。從表中可以看出MSE與RMSE與標記值本身數(shù)據(jù)量級有關(guān)，并且與其量級相比，MSE與RMSE的值處于較低水平即預(yù)測效果較好，通過觀察R2值也以很好佐證，除S11預(yù)測最差外，其余均能達到0.98以上水平。

表3 標記值誤差特性Table 3 Error characteristics of predicted values

將每個網(wǎng)格點預(yù)測得到的標記值讀入DNCR計算程序迭代求解守恒方程，最終得到收斂待預(yù)測Kn=1.0時的流場。為了直觀對比計算精度，如圖14所示，選擇y=0.5 m截線上DNCR與NS、UGKS方法的物理量分布進行對比。

圖14 Kn=1.0計算結(jié)果對比 (y=0.5 m)Fig.14 Comparison of Kn=1.0 calculation results(y=0.5 m)

采用式(13)計算評估精度提升情況：

(13)

式中：若ξ趨于1，則表明DNCR方法預(yù)測值與UGKS基本一致；若ξ趨于0，則表示DNCR預(yù)測值與N-S方程計算結(jié)果趨于一致。

通過計算得到y(tǒng)=0.5 m截線上U、P、T的結(jié)果，DNCR相比N-S精度提升分別為92.23%、97.77%、90.80%。

為進一步評估極端隨機樹的泛化能力，選取Kn=0.5與Kn=1.5兩個超出訓練集范圍的努森數(shù)條件開展DNCR計算，重點考察DNCR方法的訓練集外延數(shù)據(jù)泛化性能。其中y=0.5 m截線DNCR與N-S、UGKS方法的y方向速度、溫度與壓力物理量分布對比分別如圖15、圖16所示。DNCR較N-S方程預(yù)測精度提升情況依次為73.59%、 95.71%、71.20%和89.12%、94.36%、87.69%，精度計算公式見式(13)，可以看出在Kn=0.5與Kn=1.5條件下，DNCR較N-S方程預(yù)測精度提升百分比與Kn=1.0狀態(tài)時相比略有下降，但仍保持較高水平，即DNCR方法中采用極端隨機樹模型具有較好的泛化性能。

圖15 Kn=0.5計算結(jié)果對比 (y=0.5 m)Fig.15 Comparison of Kn=0.5 calculation results (y=0.5 m)

圖16 Kn=1.5計算結(jié)果對比 (y=0.5 m)Fig.16 Comparison of Kn=1.5 calculation results (y=0.5 m)

計算效率方面，由于DNCR方法與N-S求解方法的區(qū)別僅在于對線性應(yīng)力與熱流項進行了非線性修正，因此其迭代求解速度與N-S方程基本相當。但DNCR方法需要首先采用N-S方程對待預(yù)測狀態(tài)進行預(yù)計算并提取流場特征值作為機器學習模型的輸入數(shù)據(jù)，最終才能在N-S計算結(jié)果上耦合數(shù)據(jù)驅(qū)動的非線性本構(gòu)關(guān)系求解宏觀量守恒方程得到稀薄非平衡流數(shù)值解，因此DNCR計算時間較N-S方程會有所增加，若不考慮機器學習模型訓練時間，DNCR與N-S方程計算效率能夠保持同一量級。

5 結(jié)論與展望

機器學習中的極端隨機樹模型對基于數(shù)據(jù)驅(qū)動非線性本構(gòu)關(guān)系(DNCR)的復(fù)雜高維非線性回歸問題具有較好的建模能力，其中特征參數(shù)的選擇與模型參數(shù)的調(diào)優(yōu)是機器學習建模中十分重要的關(guān)鍵問題，直接影響模型性能與精度。本文在DNCR方法基礎(chǔ)上通過不同努森數(shù)條件下二維頂蓋驅(qū)動方腔流算例對高維非線性建模涉及的特征參數(shù)選取、參數(shù)調(diào)優(yōu)開展相關(guān)驗證與研究工作，選取若干典型狀態(tài)(包括訓練集外延狀態(tài))對極端隨機樹模型的泛化性能開展研究，評估了相關(guān)模型與方法的計算精度與計算效率。計算結(jié)果初步表明通過特征參數(shù)篩選與模型調(diào)參，DNCR中的極端隨機樹模型不僅對訓練集努森數(shù)上下限范圍內(nèi)中間狀態(tài)具有較好的預(yù)測能力，對訓練狀態(tài)范圍以外的外延狀態(tài)預(yù)測效果也較好，體現(xiàn)出了一定的泛化能力。同時，DNCR方法在模型訓練完成后計算效率與N-S方程基本保持同一量級，能夠同時提高現(xiàn)有非平衡流數(shù)值方法的計算精度與計算效率，具有較高應(yīng)用潛力。

然而，作為機器學習方法在物理建模過程中的應(yīng)用案例，這類方法仍普遍面臨以下問題：① 與已知物理定律關(guān)聯(lián)困難；② 公開高可信度數(shù)據(jù)集匱乏。因此，針對現(xiàn)有DNCR方法在機器學習模型泛化性能、訓練集數(shù)據(jù)來源與數(shù)據(jù)成本及氣固邊界條件物理適定性方面的研究尚存在較大空白，需要有針對性地開展相關(guān)理論與應(yīng)用創(chuàng)新性研究，進一步提升DNCR方法非平衡流動描述能力。