面向稀薄流非線性本構預測的機器學習方法

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

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

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

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

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

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

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

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

1 數據驅動非線性本構關系計算方法

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

DNCR方法計算流程如圖1所示，ψ代表在N-S流場數值模擬結果中提取的特征參數，作為機器學習模型的輸入量，ΔΘ為標記參數，作為機器學習模型的輸出量，即機器學習模型建立了ψ與ΔΘ之間的復雜高維回歸關系，與非線性本構物理建模函數的區別在于其本構關系沒有具體數學表達式，而是全流域存在離散當地映射ψ→ΔΘ，所形成的離散本構映射關系在全流場計算域適用。而模型泛化性能取決于訓練集特征代表性與特征參數選取，理論上若特征參數選擇未涉及到任何與空間網格與當地坐標值直接相關的參數且訓練集包含相似的非平衡流動特征，機器學習模型就具備一定的遷移預測能力。

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

DNCR方法的一個重要特點是訓練與預測過程相對獨立，圖1中紅線與綠線分別表示了訓練與預示流程。模型訓練過程通過對包含不同典型流動特征的基礎流場數據集開展訓練，獲得ψ與ΔΘ之間的復雜回歸關系。而預測過程則首先采用N-S求解器對待預測狀態開展初步計算，獲得待預測當地特征值ψ，然后采用已訓練生成的回歸關系ψ→ΔΘ得到待預測流場當地標記值ΔΘ，最終通過時間推進方式求解質量、動量與能量守恒方程，耦合非線性本構關系完成計算，守恒方程為

(1)

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

(2)

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

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

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

(3)

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

表1 一維Ar激波結構計算狀態

圖2 激波結構密度分布Fig.2 Density distributions in shock wave

圖3 激波結構速度分布Fig.3 Velocity distributions in shock wave

圖4 激波結構壓強分布Fig.4 Pressure distributions in shock wave

圖5 激波結構溫度分布Fig.5 Temperature distributions in shock wave

2 訓練數據集生成

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

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

圖7 二維頂蓋驅動方腔流計算網格Fig.7 Grid of two-dimensional lid-driven cavity flow

表2 二維方腔計算狀態Table 2 Calculation state of two-dimensional cavity

3 機器學習模型

3.1 CART與極端隨機樹算法

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

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

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

(4)

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

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

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

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

(5)

3) 繼續對兩個子區域調用步驟1)、2)，直至滿足條件。

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

(6)

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

(7)

3.2 特征參數選擇

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

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

(8)

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

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

(9)

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

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

3.3 極端隨機樹參數調優

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

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

本文采用可決系數[27]評價參數調節后的模型優劣情況。可決系數表示一個隨機變量與多個隨機變量關系的數字特征，即可決系數值越大越接近于1時說明特征參數對標記參數的解釋程度越高。

(10)

采用Kn=0.7與Kn=1.3兩組狀態對機器學習模型進行訓練，對Kn=1.0狀態進行預測。選取n_estimators在[1,1 000]范圍內，R2隨n_estimators變化趨勢如圖10所示。可以看出R2最終趨近于一個比較穩定的值，擬合較好且隨著基學習器的數量增加未出現嚴重過擬合現象。如圖11所示，模型的訓練時間與基學習器數目基本滿足線性關系，考慮到模型訓練的時間成本，本文中選取n_estimators=300。

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

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

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

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

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

圖13 預測情況Fig.13 Prediction situation

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

(11)

(12)

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

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

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

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

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

(13)

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

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

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

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

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

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

5 結論與展望

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

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