基于機器學習的非定常流場網格自適應

李彩云，安 慰，劉學軍，*，呂宏強

(1.南京航空航天大學 計算機科學與技術學院，模式分析與機器智能工業和信息化部重點實驗室，南京 211106；2.南京航空航天大學 航空學院，南京 210016)

0 引言

計算流體力學（computational fluid dynamics,CFD）的任務是數值求解流體力學控制方程，其已被廣泛應用于航空工程、生物科學工程、水利工程等所有涉及流場流動、熱交換等現象的領域。計算網格的質量直接影響數值計算的收斂性和計算精度。對于簡單流場，人工生成的網格通常可以滿足需求。然而，復雜的非定常流場具有多尺度結構，大多數情況下在計算前無法給出合適的網格，僅憑人工經驗調整網格，準確度低，因此有必要研究網格自適應技術以提高求解精度。

網格自適應技術主要有3 種：局部加密法（H 型）、局部升階法（P 型）和移動網格法（R 型）。H 型方法是根據某些誤差指標局部細化和粗化網格，較常見的方法有基于誤差估算的方法[1]、伴隨方程法[2]等，通常需要增加大量網格節點且改變網格拓撲結構。P 型方法是在數值解劇烈變化的區域使用更高階的數值格式，以提高計算精度；通常和H 型方法結合使用，如求解具有傳熱效應的不可壓縮流的HP 型自適應算法[3]。這兩類自適應方法分別在增加網格節點和高階計算方面增大了計算復雜度。R 型方法通過某種策略移動網格節點，使網格節點在流場解變化強烈的區域更密集，同時不改變網格量和網格拓撲結構，因此R 型方法在提高計算精度的同時，不會造成計算量的增加。

Cao 和Huang[4]總結和比較了各種R 型自適應方法，將其分為兩類：基于速度的方法和基于位置的方法。前者建立關于網格速度的網格方程，通過對速度場積分得到網格節點的位置。如Miller 等[5]提出的移動有限元法（moving finite element,MFE），通過最小化網格速度xt和網格函數的L2 范數生成移動網格；Thomas 和Lombard[6]提出集合守恒定律（gather conservation law,GCL）方法，根據網格速度xt的雅可比矩陣Jt所滿足的恒等式得到移動網格。雖然這類網格自適應方法的網格密度和控制函數之間的關系更直接，但因為自適應網格是網格速度場時間集成的結果，后續時間步的網格移動受前期網格移動結果的影響，最終可能產生傾斜網格[4]。基于位置的方法直接控制網格節點的坐標，一個典型方法是變分法，通過求解泛函的極小值得到網格計算域到物理域的映射，從而得到移動后網格節點的位置。這類方法中，Winslow[7]提出了一種基于變量擴散的等勢方法求解泛函的極小值；在此基礎上，Thompson 等[8]使用泊松方程控制網格密度與方向改進泛函；Brackbill 和Saltzman[9]提出了一種結合了網格密度、平滑度和正交性的網格泛函方法；Dvinsky[10]使用諧波映射量作為網格函數；Brackbill[11]將Winslow 的方法與調和映射法結合，提出了一個同時考慮網格密度和對齊條件的泛函。Huang[12]使用各向同性和等分布條件構建網格泛函[13]，用于研究一種新的離散化和求解策略[14]，保留了連續函數的基本幾何結構并取得了較好結果。針對CFD 流場網格沒有準確的函數擬合流場值的問題，Wu 等[15]對Huang[14]的方法進行改進，用反向傳播神經網絡（backpropagation neural network,BPNN）回歸模型替代擬合函數法，在二維定常流場自適應網格上取得較好結果。

現有針對非定常流場數值模擬的網格自適應方法通常每隔一段時間步就進行一次網格調整，網格函數利用每個周期初始時間步的流場解進行擬合。因此所得網格會滯后于流場解，需要不斷地自適應處理以減少時間差。Tang[16]在每段時間步上進行多次自適應處理，使用插值方法獲得網格解，從而減小時間差，但是存在數值解插值導致誤差增大的問題。韓志熔等[17]用流場速度的紊亂度作為指示變量來構造網格自適應指示函數，采用H 型網格自適應方法提高局部網格密度；這種方法提高了計算精度，但同時也存在H 型方法增大計算量的問題。Jacobson 等[18]提出了跨聲速顫振的網格自適應方法，基于表示標量場插值誤差的多尺度度量進行網格調整。

本文提出采用間斷伽遼金（discontinuity Galerkin,DG）有限元法[19]對可壓縮Navier-Stokes（N-S）方程進行數值離散，實現高階精度；將變分法和機器學習相結合，實現非定常流場網格的一次性網格自適應，以滿足所有時間步的CFD 計算；同時避免在非定常計算過程中因調整網格而涉及的幾何守恒律問題。具體而言，在初始網格上使用DG 方法試算一段時間，直至出現周期性非定常現象，然后統計非定常流場中局部數值解精度的間斷量；由網格信息和間斷量訓練BPNN 回歸模型；使用MMPDE（moving mesh partial differential equation）變分法進行網格自適應，并優化網格質量，得到最終的自適應網格。該方法有望緩解傳統網格自適應存在的計算量巨大的問題，在不改變CFD 計算復雜度的情況下提升計算精度。

1 方法介紹

圖1 給出了基于機器學習的非定常流場網格自適應框架（具體細節在2.2 節介紹）。首先在初始網格T0上使用DG 有限元方法統計得到間斷量Re0；根據初始網格和間斷量訓練BPNN 回歸模型，得到任意節點位置和間斷量的映射關系；然后使用MMPDE 變分法移動網格，每次移動后的網格Ti通過回歸模型預測得到對應網格節點間斷量Re,i，重復n次后得到網格Tn0；最后通過Laplacian 平滑法對Tn0進行質量優化，得到新的網格Tn，并返回至CFD 計算模塊中求解N-S 方程。本節對所用到的方法進行分別介紹。

圖1 基于機器學習的非定常流場網格自適應的框架Fig.1 The framework for grid adaptation of unsteady flow fields based on machine learning

圖2 交界面的左右單元變量Fig.2 Left and right cell variables of the interface

1.1 DG 有限元方法

DG 方法是典型的高階數值方法，通過提高單元上數值解多項式的階數，增加相應單元上解函數的自由度來提高空間精度。該方法易于處理任意網格和復雜幾何區域，適用于網格自適應[20]。

可壓縮黏性流動的N-S 方程可以表示為：

為了使上述方程組封閉，補充理想氣體狀態方程p=(γ ?1)ρ[E?(u2+v2)/2]，常態下空氣的比熱為γ=cp/cv=1.4。本論文采用Bassi 和Rebay 提出的 BR2 格式，引入輔助變量 Θ=?U。對式（1）采用DG 有限元方法進行空間離散，得到守恒變量的高階格式：

式中，uj(t)代表方程變量，表示每個單元內數值解的自由度，?j(x)是對應的基函數，N表示p階時基函數的個數。控制方程的解U在空間和時間上分別進行離散，其中uj(t)是關于時間的函數，而 ?j(x)為空間離散的函數。

對式（1）分部積分得到：

式中，?K為單元K的邊界集合，e為單元K的任意一條邊界，為單元K上的基函數，k為該基函數的最高階數。

本文采用總體提升因子Rh表征單元的間斷量。由于Rh在單元內是多項式分布，這里取Rh對應ρ的量在單元內積分得到Re，即網格間斷量。記Rh對應 ρ的量為Rh,1，則有：

其中，Se為單元面積，||·||2為L2 范數。

為了說明間斷量與流場數值之間的關系，圖3 給出了二維圓柱繞流結構網格在Ma∞=0.2、Re=200 下的計算結果。可見高間斷量區域與馬赫云圖中的數值變化劇烈區域相對應，表明尾流區間斷量大的網格單元局部數值解精度不足，需要移動網格節點實現局部加密。

圖3 圓柱繞流間斷量分布及馬赫云圖（Ma∞=0.2,Re=200）Fig.3 Discontinuity distribution and Mach cloud diagram of flow around the cylinder (Ma∞=0.2,Re=200)

1.2 基于MMPDE 的變分法移動網格

1.2.1 模型介紹

物理域Ω的網格Th有N個單元，每個單元記為K，節點記為X，節點X上的值由網格函數u=u(x,t)得到。固定邊界節點、移動內部網格節點，使其符合函數值分布規律。對應計算域 ?c中的網格Tc,h有相同單元數，每個單元Kc的節點記為 ξ，Th網格節點和Tc,h網格節點的對應轉換關系為x=x(ξ)、ξ=ξ(x)。為防止出現網格交叉或折疊，后面討論使用轉換 ξ=ξ(x)，計算網格節點上的值為=u(x(ξ))。從計算網格單元Kc到物理網格單元K有可逆仿射映射FK:Kc→K，使得K=FK(Kc)。

1.2.2 度量張量及MMPDE 變分法

利用插值誤差估計可以獲得最佳度量張量M[12]。首先考慮在簡單情況下，對于函數u=u(x)，任意單元K的頂點ai,i=1,···,n的線性拉格朗日插值為：

式中基函數 ?i為線性拉格朗日多項式，滿足：

記單元K的中點為XK，根據泰勒定理，將u(x)展開并根據特征分解性質和范數性質，對任意q∈[1,∞]，插值誤差可最終表示為：

式中，H(u,xK)為xK處的海森值，為FK的雅可比矩陣，h.o.t.表示高階項，常數C=O(1)。

由于各向異性方法可以比各向同性方法提供更清晰的誤差界限，特別是當物理解顯示各向異性特征時，例如在一個方向上變化比其他方向更快，所以各向異性誤差界限提供了許多自適應算法分析設計所需的信息。采用索伯列夫空間中的插值理論，結合各向異性方法，將式（6）插值誤差范圍改為：

根據等分布和對齊條件[14]，滿足M均勻網格需要符合下列條件：

結合式（4）得到度量張量為：

要注意的是，這里的項HK,α是單元K處的海森矩陣。即，首先計算網格節點處函數解的梯度，由此得到的MK控制網格節點向梯度值大的位置移動。根據均勻分布原理，本文目標是得到M度量下的均勻網格。變分法移動網格即利用求解泛函的極值來確定生成網格所需的坐標變換。為避免其離散化過程丟失泛函的幾何結構，用物理域和計算域對應的網格單元映射替代傳統的節點坐標變換，并使用MMPDE求解泛函極值，即基于MMPDE 的變分移動網格方法[14]。

帶有M的網格泛函的通用形式為：

式中，J=?ξ/?x是 ξ=ξ(x)的雅可比矩陣，G是給定的光滑函數（對每一個參數）。由計算網格單元Kc到物理網格單元K的映射關系，式（11）可以直接近似為：

式中，Tc0,h為 ?c上的參考網格，Φh為當前計算網格到物理網格的變換函數。

1.3 BPNN 回歸模型

由于前文介紹的MMPDE 變分法需使用網格函數計算每次移動后網格節點上的函數值，無法直接應用于CFD 中，因此大部分自適應網格方法采用特定函數擬合流場，如Cao 等[22]用函數u(x,t)=tanh[50(3x?|y|?t)]表示翼型激波位置。而在CFD 數值仿真中，非定常流場數值變化劇烈、區域復雜，難以用某種特定形式的函數進行擬合，因而MMPDE 變分法無法直接用于CFD 流場網格自適應。結合1.1 節介紹的DG有限元法，本文提出采用衡量網格單元計算精度的間斷量數值作為網格自適應的依據；考慮到神經網絡具有較強的逼近非線性函數的能力，并具有自適應學習、較強的魯棒性和容錯性的特點，本文使用神經網絡擬合非線性的間斷量，為CFD 流場網格自適應提供了一種新的有效途徑。tanh(x)=(ex?e?x)/(ex+e?x)，Sigmoid 函數：f(x)=1/(1+e?x)，Relu函數：R elu=max(0,x)等。

神經網絡的輸入層和輸出層神經元個數分別由輸入和輸出數據維數確定，這里輸入數據為網格節點坐標，輸出數據為網格節點對應的間斷量；隱含層層數和每層神經元個數根據經驗和實驗測試確定。為保證神經網絡實現非線性擬合，每個神經元需選用一個非線性函數作為激活函數。關于不同激活函數的研究，詳見文獻[23]。常用的激活函數有Tanh 函數：

反向傳播（backpropagation,BP）神經網絡是一種建立在梯度下降算法基礎上的多層前饋神經網絡，用輸出層神經元的預測誤差反向估計上一層隱含層神經元的預測誤差，繼而逐層將誤差從輸出層反向傳播到隱含層，最后到輸入層，從而實現對連接權重的調整。權值調整差值為：

式中，負號表示梯度下降，η為學習率，?E/?wi為誤差關于當前權值的偏導數。使用鏈式法則可以求出誤差對每個神經元權值的偏導，由式（14）得到新的權值，繼續正向傳播，進行下一次權值調整，直到網絡的輸出誤差小于既定值，或者訓練次數達到預設的學習次數為止。更多關于BPNN 的細節，如初始化權值、網絡結構等，見文獻[24-25]。

2 基于BPNN 的非定常流場網格自適應方法

網格劃分是數值模擬至關重要的一步，它直接影響著后續數值計算結果的精確性。一方面由于流動的連續性被離散化，數值模擬的精度取決于網格的節點密度和分布；另一方面，網格單元的形狀（包括偏斜率、長寬比等）對數值模擬精度也有重要影響，如偏斜率（skewness）為實際節點形狀與同等體積等邊形節點的差別，高的偏斜率會降低解的精確性，并降低收斂性。本節基于邊界網格正交性原則和Laplacian網格光滑方法提高網格質量，提出基于BPNN 的非定常流場網格自適應方法。

2.1 Laplacian 網格平滑

在黏性流場中，邊界附近流動參數變化劇烈，且通過數值求解橢圓型方程生成的曲線坐標系的坐標線與邊界耦合，所以邊界網格要求滿足正交性條件，使用細密的貼體網格劃分。考慮到初始網格的附面層網格有相對成熟的網格生成規范，在自適應調整過程中需固定附面層網格節點Xfix，否則可能出現附面層網格偏斜率過大導致離散誤差增大的情況（如圖4）。此外，畸形網格不僅影響精度，而且可能導致計算不收斂，所以在移動過程中應盡量提高網格質量。本文采用Laplacian 光滑法優化網格形狀。

圖4 附面層網格偏移Fig.4 Boundary layer grid offset

Laplacian 光滑將節點移動到其鄰接單元中心，如圖5 所示，節點I的坐標xi變換公式為：

圖5 Laplacian 光滑鄰接單元示意圖Fig.5 Schematic of Laplacian smooth adjacency unit

式中，?K為xi鄰接單元集合，|?K|為鄰接單元個數。式（15）控制節點向密度大的方向光滑，因此可以避免光滑導致網格密度特征損失，且附面層網格基本維持單元特征。

對于狹長形狀的網格，上述方法可能出現無效的移動結果，這里引入單元質量概念[9]度量網格質量好壞。網格質量通常選取為網格單元的最大角度差，復雜方法如網格外接圓半徑、變形矩陣等應用也較為廣泛[11]。這里網格質量根據單元角度定義，如四邊形平面角度αi(i=1,2,3,4)，定義 ηi=sinαi，則單元質量取：

式（16）可有效識別單元中是否出現負角。若單元質量提高則移動網格節點，否則移回距離原節點的1/2 位置；移動次數達到預設上限或者單元質量差值達到閾值，則光滑停止。

2.2 基于BPNN 的非定常流場網格自適應方法

為了使移動網格方法獨立于網格函數，本文采用回歸模型代替誤差估計方法[13]，其訓練過程如圖6所示（具體網絡參數見3.3 節）。T0為初始網格，Re0為DG 有限元方法在初始網格T0上離散N-S 方程求解一段時間所得的網格間斷量。由T0和Re0訓練BPNN 回歸模型，網格節點坐標為模型輸入數據，節點對應統計間斷量為模型輸出數據。

圖6 回歸模型的訓練過程Fig.6 The training process of the regression model

根據間斷量進行網格自適應，需對變分法中的度量張量（式（10））進行修改：

將上述訓練所得回歸模型與變分法移動網格相結合，每次網格移動的節點坐標通過回歸模型預測當前位置的函數值。圖7 給出第n+1 次迭代過程，為當前物理網格，為回歸模型對預測出的網格間斷量，為給定的參考網格。算法迭代過程如下：

圖7 基于BPNN 回歸的變分法移動網格Fig.7 Variational Moving Grid Based on BPNN Regression

5）經過N次迭代后，得到最終自適應物理網格。

在迭代過程中，網格節點數和單元數保持不變。網格自適應后，使用Laplacian 平滑方法提高網格質量，整體流程如圖8 所示。由圖8 可見，本文提出的基于機器學習的非定常流場自適應方法完全獨立于CFD 計算。在初始網格上采用DG 有限元法離散N-S 方程，統計得到網格間斷量，輸入到基于機器學習的網格自適應模塊中；自適應過程中，采用回歸模型代替CFD 計算獲得每次移動后的網格節點值；自適應模塊所得網格節點位置信息返回到CFD 計算中，重新進行非定常流場計算，完成一次性網格調整。

圖8 算法迭代的整體流程圖Fig.8 Overall flow chart of algorithm iteration

3 實驗和分析

3.1 實驗數據

本節基于DG 有限元法對二維圓柱繞流進行數值模擬，采用2.2 節提到的基于BPNN 的非定常流場網格自適應方法測試。來流條件為Ma∞=0.2、Re=200。壁面采用無滑移絕熱邊界條件，遠場采用特征邊界條件，本算例采用結構網格進行計算，單元總數為868，節點數為917，第一層網格節點的物面距為圓柱直徑的0.01。采用四核并行計算，如圖9 所示，分區網格單元數分別為213、211、211 和223，負載較為均衡。時間步長取 Δt=0.01，時間離散采用一階隱式向后差分格式（BDF），空間離散選取三階DG 方法。試算一段時間直到出現如圖10 所示的1 個周期的非定常流場，統計所得網格間斷量如圖11（a）所示。

圖9 二維圓柱繞流初始網格Fig.9 Two-dimension initial grid of flow around the cylinder

圖10 采用結構初始網格計算得到的非定常馬赫云圖Fig.10 Unsteady Mach cloud obtained by using structure initial grid computation

3.2 數據預處理

由于初始值的影響，圖11（a）中遠場間斷量增大，需要對數據進行處理并做平滑，結果如圖11（b）所示。由于遠場附近間斷量均接近0，所以對應遠場附近度量張量M大小相近（如圖12（a）度量張量二維矢量圖），根據均勻分布準則，對間斷量進行變換。由圖11（b）可知，圓柱尾部的間斷量數值在y軸方向符合高斯概率密度函數（PDF）曲線，整體在x方向符合高斯累積分布函數（CDF）曲線，因此對間斷量Re進行預處理，對應的控制度量張量M大小為：

圖12 預處理前后度量張量示例Fig.12 Example of metric tensors before and after preprocessing

式（18）中，用權重系數λ1、λ2保持間斷量Re特征，取期望 μPDF1、μCDF1為 0，標準差 σPDF1、σCDF1分別為5 和1；預處理后的間斷量如圖11（c）所示，其空間分布更能合理反映流場的高梯度變化特性，且適用于MMPDE 網格自適應方法。式（19）中，用權重系數λM控制放縮M的幅度，取期望 μPDF2、μCDF2為0，標準差 σPDF2、σCDF2分別為8 和1，對應度量張量二維矢量如圖12（b）所示。另外，由于物面附近度量張量較大，且呈線性變化，所以固定附面層網格節點時，外層單元在網格自適應過程中不會產生畸變。

3.3 神經網絡設置

BPNN 訓練數據為初始網格上917 個節點坐標及對應統計間斷量，隱含層和輸出層的激活函數隱含層分別為Tanh 函數和Sigmoid 函數，訓練誤差目標值為Pgoal=0.000 01，學習率為lr=0.01，最大訓練次數為Pepoch=1 000，驗證次數為Pmaxfail=100。

一般情況下，BPNN 網絡層數加深容易導致梯度爆炸和梯度消失的問題，因此隱含層的數量通常設置為小于4。本文首先探究了適用于預測網格節點間斷量的網絡結構，選擇預測結果的均方根誤差（RMSE）作為網絡結構的評價指標，其隱含層計算方式為：

式中，yi和分別為真實值和回歸模型的預測值，m為樣本個數。圖13 顯示了在相同數據集上對不同隱含層設置的神經網絡進行十折交叉驗證的RMSE 箱線圖。由圖13 可知，具有一個隱含層的神經網絡預測誤差較大，這表明出現了欠擬合現象。隨著該隱含層神經元數量的增加，預測誤差逐漸降低；而增加隱含層也能夠顯著減小預測誤差，并保持穩定隱含層；但當隱含層過多或每個隱含層的神經元過多時，會出現過擬合，從而導致預測精度降低。因此，本文選擇三個隱含層[10,10,10]預測任意位置的網格節點間斷量。

圖13 神經網絡結構選擇箱線圖Fig.13 Neural network structure selection boxplot

3.4 自適應結果分析

傳統有限元網格的質量通常以翹曲度、扭曲角以及雅可比比率等參數指標來度量，而對于CFD 高精度方法的計算網格而言，網格質量由計算收斂所得數值解的精度來評價，網格質量過低會導致計算無法收斂。因此，本節以計算精度來分析評估網格質量。

為了探究間斷量的數據預處理對網格移動結果的影響，本文應用基于BPNN 的非定常流場網格自適應方法，分別對有、無數據預處理的間斷量進行網格移動。圖14（a、b）分別標記了初始網格有、無預處理的不動節點和高間斷量節點，圖14（c、d）分別為調整后的對應網格，網格節點加密位置與高間斷量區域相對應。

圖14 有無預處理自適應網格結果對比Fig.14 Comparison of adaptive meshgrid results with and without preprocessing

將圖14（c、d）中的網格節點位置信息返回至CFD 計算模塊進行N-S 方程求解，計算條件與實驗數據的相同，升力系數、表面壓力系數等計算結果的輸出周期為10 個時間步長。圖14（e～h）顯示了計算迭代5 000 步的對比結果。

圖14（e、f）分別顯示了有、無數據預處理的自適應網格升力系數曲線，相較于初始網格升力系數曲線（紅色虛線），經兩種方法調整后的網格升力系數曲線均提前出現非定常振蕩。從統計意義上講，圓柱繞流的合理網格分布應為上下對稱，帶有預處理流程的網格優化結果生成的最終網格基本符合上下對稱特征。另外對于高階DG 方法而言，高階基函數在不同網格內的分布跟單個網格單元的節點編號相關，即對稱網格單元中基函數的分布往往并不對稱，導致上下兩部分區域的數值解存在理論上的差別。因此，不同于有限體積法，采用高階DG 方法進行非定常流場計算不會因網格對稱而延遲或減弱非定常現象。從該圓柱算例的結果可見，經預處理的網格計算結果在1 000 步左右已經出現非定常效應，而未經預處理的網格計算結果在1 500 步左右才出現，說明數據預處理后的網格分布更合理，因此提高了計算精度。

圖14（g、h）顯示了分別采用兩種網格計算5 000步得到的馬赫云圖，可見圖14（g）尾流區清晰度高于圖14（h）。根據間斷量化指標，間斷量較大的區域需要較密的網格分布；間斷量較小的區域數值解相對光滑，網格可以稀疏分布。在圓柱繞流中，渦街區域需要更密的網格。對比結果顯示，數據預處理后的自適應網格獲得了更為精確的卡門渦街，可見預處理后的間斷量可以更合理地反映流場結構的物理規律。

表1 給出了相同條件下，初始網格和有無數據預處理的自適應網格進行5 000 步計算的耗時情況。其中，用于CFD 計算的硬件設備為Intel Core i7-9700 CPU 3.00 GHz，采用 8 核并行計算；用于神經網絡預測及網格自適應的硬件設備為Intel Core i9-10900K CPU 3.70 GHz，采用串行計算。從表1 可以看出，3 種網格的非定常流場計算耗時接近，自適應后的網格在提高計算精度的同時，沒有明顯降低計算效率，且數據預處理后的網格比未經預處理的網格計算耗時更少。

表1 初始網格和有無數據預處理的自適應網格計算耗時Table 1 CPU time of initial mesh and adaptive mesh with or without data preprocessing

綜上分析，本文所提基于BPNN 的網格自適應方法適用于非定常流場網格調整。另外，基于DG 有限元方法統計獲得的間斷量指標需經過一定時間的非定常計算，該統計時長難以嚴格控制為整數周期或準周期，這一定程度上影響了網格優化質量。因此，本文提出的預處理過程旨在消除不同非定常計算時長對網格優化效果的影響，通過對數據進行預處理并控制度量張量M，使得網格自適應結果符合流場特征，從而提高計算精度。0

3.5 模型泛化性分析

本節采用更復雜的非定常流場算例進一步驗證本文所提模型的泛化能力。該算例為相同來流條件（Ma∞=0.2、Re=200）下的圓柱繞流，采用四邊形/三角形混合網格，由貼體的四邊形結構網格和外場的三角形非結構網格構成，節點數為4 053，網格單元數為7 496，如圖15 所示。采用一階DG 方法，統計所得網格間斷量如圖16 所示。

圖15 非定常圓柱繞流混合初始網格Fig.15 Hybrid initial grid of unsteady flow around the cylinder

圖16 混合網格在Ma=0.2、Re=200 來流條件下的網格間斷量Fig.16 Mesh discontinuity of hybrid grid at Ma=0.2,Re=200

BPNN 訓練數據為初始網格節點及對應間斷量，由于數據分布與前一算例相似，采用相同的網絡參數即可得到較好的訓練結果（見圖17）。考慮到附面層網格有相對成熟的網格生成規范，因此在該算例中，附面層的四邊形結構網格節點保持固定（見圖15），自適應后的網格如圖18 所示。初始網格和自適應網格分別在相同初始條件下進行30 000 步非定常計算，結果如圖19～圖21 所示。

圖17 BPNN 訓練結果Fig.17 BPNN training results

圖18 混合網格圓柱繞流自適應網格Fig.18 Adaptive hybrid grid of flow around the cylinder

圖19 采用初始網格計算得到的馬赫云圖Fig.19 Mach cloud obtained by initial grid computation

圖19 和圖20 分別為初始網格和自適應網格計算30 000 步后的馬赫云圖，可以看到，網格自適應方法在流場結構更復雜的渦街區域實現了網格的相對加密，計算所得的流場結構更清晰。圖21 為升力系數曲線對比，可以發現，自適應網格升力系數曲線（藍色實線）較初始網格升力系數曲線（紅色虛線），提前出現了顯著的非定常振蕩。

圖20 采用自適應網格計算得到的馬赫云圖Fig.20 Mach cloud obtained by adaptive grid computation

圖21 初始網格和自適應網格下的升力系數對比Fig.21 Comparison of CL between the initial grid and the adaptive grid

表2 給出了自適應前后網格的非定常計算、BPNN訓練與網格自適應步驟的耗時對比。相比表1，神經網絡訓練平均耗時隨數據量增加呈線性增長。在處理更大規模的網格時，網格自適應所需的計算時間遠小于非定常計算時間，因此本文所提網格自適應方法在提高計算精度的同時，能夠保證計算效率。

表2 網格自適應前后的各種計算耗時Table 2 Different CPU time before and after grid adaptive

4 結論

CFD 高精度數值模擬結果依賴流場方程求解時的數值格式和計算網格質量。針對非定常流場計算網格調整困難的問題，本文研究了結合機器學習和變分法的移動網格方法，在不增加網格節點數且不改變拓撲結構的前提下，移動網格節點，具體方法為：基于DG 有限元方法，在原始網格上獲得間斷量指標統計值，采用機器學習方法進行一次性的R 型網格自適應，同時采用Laplacian 平滑法保證網格質量，并在調整后的網格上進行最終的非定常計算。整個過程不涉及幾何守恒律問題，使得該網格自適應方法獨立于CFD 計算，同時減少了非定常網格自適應的時間耗費。

本文通過圓柱繞流算例比較了調整前后二維圓柱繞流網格的非定常計算結果，渦街區域相對加密的自適應網格更早出現非定常現象，且流場結構更為清晰。該實驗結果表明，間斷量指標統計值可作為非定常流場網格自適應的參考依據，通過結合機器學習使得網格調整自動化、智能化，調整后的網格在保證計算效率的同時，提高了CFD 計算精度。

本文網格自適應方法以三角形（二維）為基礎單元，任意網格單元可以剖分成三角形單元進行網格自適應，且移動過程中不產生負體積問題，能夠保持原始網格的拓撲結構。因此，所提方法對網格單元沒有限制，適用于結構網格、非結構網格以及混合網格。后續可以將本文方法拓展至三維，以四面體為基本單元進行網格自適應。