李彩云，安 慰，劉學(xué)軍，*，呂宏強(qiáng)

(1.南京航空航天大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，模式分析與機(jī)器智能工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室，南京 211106；2.南京航空航天大學(xué) 航空學(xué)院，南京 210016)

0 引言

計(jì)算流體力學(xué)（computational fluid dynamics,CFD）的任務(wù)是數(shù)值求解流體力學(xué)控制方程，其已被廣泛應(yīng)用于航空工程、生物科學(xué)工程、水利工程等所有涉及流場(chǎng)流動(dòng)、熱交換等現(xiàn)象的領(lǐng)域。計(jì)算網(wǎng)格的質(zhì)量直接影響數(shù)值計(jì)算的收斂性和計(jì)算精度。對(duì)于簡(jiǎn)單流場(chǎng)，人工生成的網(wǎng)格通?？梢詽M足需求。然而，復(fù)雜的非定常流場(chǎng)具有多尺度結(jié)構(gòu)，大多數(shù)情況下在計(jì)算前無法給出合適的網(wǎng)格，僅憑人工經(jīng)驗(yàn)調(diào)整網(wǎng)格，準(zhǔn)確度低，因此有必要研究網(wǎng)格自適應(yīng)技術(shù)以提高求解精度。

網(wǎng)格自適應(yīng)技術(shù)主要有3 種：局部加密法（H 型）、局部升階法（P 型）和移動(dòng)網(wǎng)格法（R 型）。H 型方法是根據(jù)某些誤差指標(biāo)局部細(xì)化和粗化網(wǎng)格，較常見的方法有基于誤差估算的方法[1]、伴隨方程法[2]等，通常需要增加大量網(wǎng)格節(jié)點(diǎn)且改變網(wǎng)格拓?fù)浣Y(jié)構(gòu)。P 型方法是在數(shù)值解劇烈變化的區(qū)域使用更高階的數(shù)值格式，以提高計(jì)算精度；通常和H 型方法結(jié)合使用，如求解具有傳熱效應(yīng)的不可壓縮流的HP 型自適應(yīng)算法[3]。這兩類自適應(yīng)方法分別在增加網(wǎng)格節(jié)點(diǎn)和高階計(jì)算方面增大了計(jì)算復(fù)雜度。R 型方法通過某種策略移動(dòng)網(wǎng)格節(jié)點(diǎn)，使網(wǎng)格節(jié)點(diǎn)在流場(chǎng)解變化強(qiáng)烈的區(qū)域更密集，同時(shí)不改變網(wǎng)格量和網(wǎng)格拓?fù)浣Y(jié)構(gòu)，因此R 型方法在提高計(jì)算精度的同時(shí)，不會(huì)造成計(jì)算量的增加。

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

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

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

1 方法介紹

圖1 給出了基于機(jī)器學(xué)習(xí)的非定常流場(chǎng)網(wǎng)格自適應(yīng)框架（具體細(xì)節(jié)在2.2 節(jié)介紹）。首先在初始網(wǎng)格T0上使用DG 有限元方法統(tǒng)計(jì)得到間斷量Re0；根據(jù)初始網(wǎng)格和間斷量訓(xùn)練BPNN 回歸模型，得到任意節(jié)點(diǎn)位置和間斷量的映射關(guān)系；然后使用MMPDE 變分法移動(dòng)網(wǎng)格，每次移動(dòng)后的網(wǎng)格Ti通過回歸模型預(yù)測(cè)得到對(duì)應(yīng)網(wǎng)格節(jié)點(diǎn)間斷量Re,i，重復(fù)n次后得到網(wǎng)格Tn0；最后通過Laplacian 平滑法對(duì)Tn0進(jìn)行質(zhì)量?jī)?yōu)化，得到新的網(wǎng)格Tn，并返回至CFD 計(jì)算模塊中求解N-S 方程。本節(jié)對(duì)所用到的方法進(jìn)行分別介紹。

圖1 基于機(jī)器學(xué)習(xí)的非定常流場(chǎng)網(wǎng)格自適應(yīng)的框架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 方法是典型的高階數(shù)值方法，通過提高單元上數(shù)值解多項(xiàng)式的階數(shù)，增加相應(yīng)單元上解函數(shù)的自由度來提高空間精度。該方法易于處理任意網(wǎng)格和復(fù)雜幾何區(qū)域，適用于網(wǎng)格自適應(yīng)[20]。

可壓縮黏性流動(dòng)的N-S 方程可以表示為：

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

式中，uj(t)代表方程變量，表示每個(gè)單元內(nèi)數(shù)值解的自由度，?j(x)是對(duì)應(yīng)的基函數(shù)，N表示p階時(shí)基函數(shù)的個(gè)數(shù)?？刂品匠痰慕釻在空間和時(shí)間上分別進(jìn)行離散，其中uj(t)是關(guān)于時(shí)間的函數(shù)，而 ?j(x)為空間離散的函數(shù)。

對(duì)式（1）分部積分得到：

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

本文采用總體提升因子Rh表征單元的間斷量。由于Rh在單元內(nèi)是多項(xiàng)式分布，這里取Rh對(duì)應(yīng)ρ的量在單元內(nèi)積分得到Re，即網(wǎng)格間斷量。記Rh對(duì)應(yīng) ρ的量為Rh,1，則有：

其中，Se為單元面積，||·||2為L(zhǎng)2 范數(shù)。

為了說明間斷量與流場(chǎng)數(shù)值之間的關(guān)系，圖3 給出了二維圓柱繞流結(jié)構(gòu)網(wǎng)格在Ma∞=0.2、Re=200 下的計(jì)算結(jié)果?？梢姼唛g斷量區(qū)域與馬赫云圖中的數(shù)值變化劇烈區(qū)域相對(duì)應(yīng)，表明尾流區(qū)間斷量大的網(wǎng)格單元局部數(shù)值解精度不足，需要移動(dòng)網(wǎng)格節(jié)點(diǎn)實(shí)現(xiàn)局部加密。

圖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 的變分法移動(dòng)網(wǎng)格

1.2.1 模型介紹

物理域Ω的網(wǎng)格Th有N個(gè)單元，每個(gè)單元記為K，節(jié)點(diǎn)記為X，節(jié)點(diǎn)X上的值由網(wǎng)格函數(shù)u=u(x,t)得到。固定邊界節(jié)點(diǎn)、移動(dòng)內(nèi)部網(wǎng)格節(jié)點(diǎn)，使其符合函數(shù)值分布規(guī)律。對(duì)應(yīng)計(jì)算域 ?c中的網(wǎng)格Tc,h有相同單元數(shù)，每個(gè)單元Kc的節(jié)點(diǎn)記為 ξ，Th網(wǎng)格節(jié)點(diǎn)和Tc,h網(wǎng)格節(jié)點(diǎn)的對(duì)應(yīng)轉(zhuǎn)換關(guān)系為x=x(ξ)、ξ=ξ(x)。為防止出現(xiàn)網(wǎng)格交叉或折疊，后面討論使用轉(zhuǎn)換 ξ=ξ(x)，計(jì)算網(wǎng)格節(jié)點(diǎn)上的值為=u(x(ξ))。從計(jì)算網(wǎng)格單元Kc到物理網(wǎng)格單元K有可逆仿射映射FK:Kc→K，使得K=FK(Kc)。

1.2.2 度量張量及MMPDE 變分法

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

式中基函數(shù) ?i為線性拉格朗日多項(xiàng)式，滿足：

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

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

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

根據(jù)等分布和對(duì)齊條件[14]，滿足M均勻網(wǎng)格需要符合下列條件：

結(jié)合式（4）得到度量張量為：

要注意的是，這里的項(xiàng)HK,α是單元K處的海森矩陣。即，首先計(jì)算網(wǎng)格節(jié)點(diǎn)處函數(shù)解的梯度，由此得到的MK控制網(wǎng)格節(jié)點(diǎn)向梯度值大的位置移動(dòng)。根據(jù)均勻分布原理，本文目標(biāo)是得到M度量下的均勻網(wǎng)格。變分法移動(dòng)網(wǎng)格即利用求解泛函的極值來確定生成網(wǎng)格所需的坐標(biāo)變換。為避免其離散化過程丟失泛函的幾何結(jié)構(gòu)，用物理域和計(jì)算域?qū)?yīng)的網(wǎng)格單元映射替代傳統(tǒng)的節(jié)點(diǎn)坐標(biāo)變換，并使用MMPDE求解泛函極值，即基于MMPDE 的變分移動(dòng)網(wǎng)格方法[14]。

帶有M的網(wǎng)格泛函的通用形式為：

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

式中，Tc0,h為 ?c上的參考網(wǎng)格，Φh為當(dāng)前計(jì)算網(wǎng)格到物理網(wǎng)格的變換函數(shù)。

1.3 BPNN 回歸模型

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

神經(jīng)網(wǎng)絡(luò)的輸入層和輸出層神經(jīng)元個(gè)數(shù)分別由輸入和輸出數(shù)據(jù)維數(shù)確定，這里輸入數(shù)據(jù)為網(wǎng)格節(jié)點(diǎn)坐標(biāo)，輸出數(shù)據(jù)為網(wǎng)格節(jié)點(diǎn)對(duì)應(yīng)的間斷量；隱含層層數(shù)和每層神經(jīng)元個(gè)數(shù)根據(jù)經(jīng)驗(yàn)和實(shí)驗(yàn)測(cè)試確定。為保證神經(jīng)網(wǎng)絡(luò)實(shí)現(xiàn)非線性擬合，每個(gè)神經(jīng)元需選用一個(gè)非線性函數(shù)作為激活函數(shù)。關(guān)于不同激活函數(shù)的研究，詳見文獻(xiàn)[23]。常用的激活函數(shù)有Tanh 函數(shù)：

反向傳播（backpropagation,BP）神經(jīng)網(wǎng)絡(luò)是一種建立在梯度下降算法基礎(chǔ)上的多層前饋神經(jīng)網(wǎng)絡(luò)，用輸出層神經(jīng)元的預(yù)測(cè)誤差反向估計(jì)上一層隱含層神經(jīng)元的預(yù)測(cè)誤差，繼而逐層將誤差從輸出層反向傳播到隱含層，最后到輸入層，從而實(shí)現(xiàn)對(duì)連接權(quán)重的調(diào)整。權(quán)值調(diào)整差值為：

式中，負(fù)號(hào)表示梯度下降，η為學(xué)習(xí)率，?E/?wi為誤差關(guān)于當(dāng)前權(quán)值的偏導(dǎo)數(shù)。使用鏈?zhǔn)椒▌t可以求出誤差對(duì)每個(gè)神經(jīng)元權(quán)值的偏導(dǎo)，由式（14）得到新的權(quán)值，繼續(xù)正向傳播，進(jìn)行下一次權(quán)值調(diào)整，直到網(wǎng)絡(luò)的輸出誤差小于既定值，或者訓(xùn)練次數(shù)達(dá)到預(yù)設(shè)的學(xué)習(xí)次數(shù)為止。更多關(guān)于BPNN 的細(xì)節(jié)，如初始化權(quán)值、網(wǎng)絡(luò)結(jié)構(gòu)等，見文獻(xiàn)[24-25]。

2 基于BPNN 的非定常流場(chǎng)網(wǎng)格自適應(yīng)方法

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

2.1 Laplacian 網(wǎng)格平滑

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

圖4 附面層網(wǎng)格偏移Fig.4 Boundary layer grid offset

Laplacian 光滑將節(jié)點(diǎn)移動(dòng)到其鄰接單元中心，如圖5 所示，節(jié)點(diǎn)I的坐標(biāo)xi變換公式為：

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

式中，?K為xi鄰接單元集合，|?K|為鄰接單元個(gè)數(shù)。式（15）控制節(jié)點(diǎn)向密度大的方向光滑，因此可以避免光滑導(dǎo)致網(wǎng)格密度特征損失，且附面層網(wǎng)格基本維持單元特征。

對(duì)于狹長(zhǎng)形狀的網(wǎng)格，上述方法可能出現(xiàn)無效的移動(dòng)結(jié)果，這里引入單元質(zhì)量概念[9]度量網(wǎng)格質(zhì)量好壞。網(wǎng)格質(zhì)量通常選取為網(wǎng)格單元的最大角度差，復(fù)雜方法如網(wǎng)格外接圓半徑、變形矩陣等應(yīng)用也較為廣泛[11]。這里網(wǎng)格質(zhì)量根據(jù)單元角度定義，如四邊形平面角度αi(i=1,2,3,4)，定義 ηi=sinαi，則單元質(zhì)量?。?/p>

式（16）可有效識(shí)別單元中是否出現(xiàn)負(fù)角。若單元質(zhì)量提高則移動(dòng)網(wǎng)格節(jié)點(diǎn)，否則移回距離原節(jié)點(diǎn)的1/2 位置；移動(dòng)次數(shù)達(dá)到預(yù)設(shè)上限或者單元質(zhì)量差值達(dá)到閾值，則光滑停止。

2.2 基于BPNN 的非定常流場(chǎng)網(wǎng)格自適應(yīng)方法

為了使移動(dòng)網(wǎng)格方法獨(dú)立于網(wǎng)格函數(shù)，本文采用回歸模型代替誤差估計(jì)方法[13]，其訓(xùn)練過程如圖6所示（具體網(wǎng)絡(luò)參數(shù)見3.3 節(jié)）。T0為初始網(wǎng)格，Re0為DG 有限元方法在初始網(wǎng)格T0上離散N-S 方程求解一段時(shí)間所得的網(wǎng)格間斷量。由T0和Re0訓(xùn)練BPNN 回歸模型，網(wǎng)格節(jié)點(diǎn)坐標(biāo)為模型輸入數(shù)據(jù)，節(jié)點(diǎn)對(duì)應(yīng)統(tǒng)計(jì)間斷量為模型輸出數(shù)據(jù)。

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

根據(jù)間斷量進(jìn)行網(wǎng)格自適應(yīng)，需對(duì)變分法中的度量張量（式（10））進(jìn)行修改：

將上述訓(xùn)練所得回歸模型與變分法移動(dòng)網(wǎng)格相結(jié)合，每次網(wǎng)格移動(dòng)的節(jié)點(diǎn)坐標(biāo)通過回歸模型預(yù)測(cè)當(dāng)前位置的函數(shù)值。圖7 給出第n+1 次迭代過程，為當(dāng)前物理網(wǎng)格，為回歸模型對(duì)預(yù)測(cè)出的網(wǎng)格間斷量，為給定的參考網(wǎng)格。算法迭代過程如下：

圖7 基于BPNN 回歸的變分法移動(dòng)網(wǎng)格Fig.7 Variational Moving Grid Based on BPNN Regression

5）經(jīng)過N次迭代后，得到最終自適應(yīng)物理網(wǎng)格。

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

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

3 實(shí)驗(yàn)和分析

3.1 實(shí)驗(yàn)數(shù)據(jù)

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

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

圖10 采用結(jié)構(gòu)初始網(wǎng)格計(jì)算得到的非定常馬赫云圖Fig.10 Unsteady Mach cloud obtained by using structure initial grid computation

3.2 數(shù)據(jù)預(yù)處理

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

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

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

3.3 神經(jīng)網(wǎng)絡(luò)設(shè)置

BPNN 訓(xùn)練數(shù)據(jù)為初始網(wǎng)格上917 個(gè)節(jié)點(diǎn)坐標(biāo)及對(duì)應(yīng)統(tǒng)計(jì)間斷量，隱含層和輸出層的激活函數(shù)隱含層分別為Tanh 函數(shù)和Sigmoid 函數(shù)，訓(xùn)練誤差目標(biāo)值為Pgoal=0.000 01，學(xué)習(xí)率為lr=0.01，最大訓(xùn)練次數(shù)為Pepoch=1 000，驗(yàn)證次數(shù)為Pmaxfail=100。

一般情況下，BPNN 網(wǎng)絡(luò)層數(shù)加深容易導(dǎo)致梯度爆炸和梯度消失的問題，因此隱含層的數(shù)量通常設(shè)置為小于4。本文首先探究了適用于預(yù)測(cè)網(wǎng)格節(jié)點(diǎn)間斷量的網(wǎng)絡(luò)結(jié)構(gòu)，選擇預(yù)測(cè)結(jié)果的均方根誤差（RMSE）作為網(wǎng)絡(luò)結(jié)構(gòu)的評(píng)價(jià)指標(biāo)，其隱含層計(jì)算方式為：

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

圖13 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)選擇箱線圖Fig.13 Neural network structure selection boxplot

3.4 自適應(yīng)結(jié)果分析

傳統(tǒng)有限元網(wǎng)格的質(zhì)量通常以翹曲度、扭曲角以及雅可比比率等參數(shù)指標(biāo)來度量，而對(duì)于CFD 高精度方法的計(jì)算網(wǎng)格而言，網(wǎng)格質(zhì)量由計(jì)算收斂所得數(shù)值解的精度來評(píng)價(jià)，網(wǎng)格質(zhì)量過低會(huì)導(dǎo)致計(jì)算無法收斂。因此，本節(jié)以計(jì)算精度來分析評(píng)估網(wǎng)格質(zhì)量。

為了探究間斷量的數(shù)據(jù)預(yù)處理對(duì)網(wǎng)格移動(dòng)結(jié)果的影響，本文應(yīng)用基于BPNN 的非定常流場(chǎng)網(wǎng)格自適應(yīng)方法，分別對(duì)有、無數(shù)據(jù)預(yù)處理的間斷量進(jìn)行網(wǎng)格移動(dòng)。圖14（a、b）分別標(biāo)記了初始網(wǎng)格有、無預(yù)處理的不動(dòng)節(jié)點(diǎn)和高間斷量節(jié)點(diǎn)，圖14（c、d）分別為調(diào)整后的對(duì)應(yīng)網(wǎng)格，網(wǎng)格節(jié)點(diǎn)加密位置與高間斷量區(qū)域相對(duì)應(yīng)。

圖14 有無預(yù)處理自適應(yīng)網(wǎng)格結(jié)果對(duì)比Fig.14 Comparison of adaptive meshgrid results with and without preprocessing

將圖14（c、d）中的網(wǎng)格節(jié)點(diǎn)位置信息返回至CFD 計(jì)算模塊進(jìn)行N-S 方程求解，計(jì)算條件與實(shí)驗(yàn)數(shù)據(jù)的相同，升力系數(shù)、表面壓力系數(shù)等計(jì)算結(jié)果的輸出周期為10 個(gè)時(shí)間步長(zhǎng)。圖14（e～h）顯示了計(jì)算迭代5 000 步的對(duì)比結(jié)果。

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

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

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

表1 初始網(wǎng)格和有無數(shù)據(jù)預(yù)處理的自適應(yīng)網(wǎng)格計(jì)算耗時(shí)Table 1 CPU time of initial mesh and adaptive mesh with or without data preprocessing

綜上分析，本文所提基于BPNN 的網(wǎng)格自適應(yīng)方法適用于非定常流場(chǎng)網(wǎng)格調(diào)整。另外，基于DG 有限元方法統(tǒng)計(jì)獲得的間斷量指標(biāo)需經(jīng)過一定時(shí)間的非定常計(jì)算，該統(tǒng)計(jì)時(shí)長(zhǎng)難以嚴(yán)格控制為整數(shù)周期或準(zhǔn)周期，這一定程度上影響了網(wǎng)格優(yōu)化質(zhì)量。因此，本文提出的預(yù)處理過程旨在消除不同非定常計(jì)算時(shí)長(zhǎng)對(duì)網(wǎng)格優(yōu)化效果的影響，通過對(duì)數(shù)據(jù)進(jìn)行預(yù)處理并控制度量張量M，使得網(wǎng)格自適應(yīng)結(jié)果符合流場(chǎng)特征，從而提高計(jì)算精度。0

3.5 模型泛化性分析

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

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

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

BPNN 訓(xùn)練數(shù)據(jù)為初始網(wǎng)格節(jié)點(diǎn)及對(duì)應(yīng)間斷量，由于數(shù)據(jù)分布與前一算例相似，采用相同的網(wǎng)絡(luò)參數(shù)即可得到較好的訓(xùn)練結(jié)果（見圖17）?？紤]到附面層網(wǎng)格有相對(duì)成熟的網(wǎng)格生成規(guī)范，因此在該算例中，附面層的四邊形結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)保持固定（見圖15），自適應(yīng)后的網(wǎng)格如圖18 所示。初始網(wǎng)格和自適應(yīng)網(wǎng)格分別在相同初始條件下進(jìn)行30 000 步非定常計(jì)算，結(jié)果如圖19～圖21 所示。

圖17 BPNN 訓(xùn)練結(jié)果Fig.17 BPNN training results

圖18 混合網(wǎng)格圓柱繞流自適應(yīng)網(wǎng)格Fig.18 Adaptive hybrid grid of flow around the cylinder

圖19 采用初始網(wǎng)格計(jì)算得到的馬赫云圖Fig.19 Mach cloud obtained by initial grid computation

圖19 和圖20 分別為初始網(wǎng)格和自適應(yīng)網(wǎng)格計(jì)算30 000 步后的馬赫云圖，可以看到，網(wǎng)格自適應(yīng)方法在流場(chǎng)結(jié)構(gòu)更復(fù)雜的渦街區(qū)域?qū)崿F(xiàn)了網(wǎng)格的相對(duì)加密，計(jì)算所得的流場(chǎng)結(jié)構(gòu)更清晰。圖21 為升力系數(shù)曲線對(duì)比，可以發(fā)現(xiàn)，自適應(yīng)網(wǎng)格升力系數(shù)曲線（藍(lán)色實(shí)線）較初始網(wǎng)格升力系數(shù)曲線（紅色虛線），提前出現(xiàn)了顯著的非定常振蕩。

圖20 采用自適應(yīng)網(wǎng)格計(jì)算得到的馬赫云圖Fig.20 Mach cloud obtained by adaptive grid computation

圖21 初始網(wǎng)格和自適應(yīng)網(wǎng)格下的升力系數(shù)對(duì)比Fig.21 Comparison of CL between the initial grid and the adaptive grid

表2 給出了自適應(yīng)前后網(wǎng)格的非定常計(jì)算、BPNN訓(xùn)練與網(wǎng)格自適應(yīng)步驟的耗時(shí)對(duì)比。相比表1，神經(jīng)網(wǎng)絡(luò)訓(xùn)練平均耗時(shí)隨數(shù)據(jù)量增加呈線性增長(zhǎng)。在處理更大規(guī)模的網(wǎng)格時(shí)，網(wǎng)格自適應(yīng)所需的計(jì)算時(shí)間遠(yuǎn)小于非定常計(jì)算時(shí)間，因此本文所提網(wǎng)格自適應(yīng)方法在提高計(jì)算精度的同時(shí)，能夠保證計(jì)算效率。

表2 網(wǎng)格自適應(yīng)前后的各種計(jì)算耗時(shí)Table 2 Different CPU time before and after grid adaptive

4 結(jié)論

CFD 高精度數(shù)值模擬結(jié)果依賴流場(chǎng)方程求解時(shí)的數(shù)值格式和計(jì)算網(wǎng)格質(zhì)量。針對(duì)非定常流場(chǎng)計(jì)算網(wǎng)格調(diào)整困難的問題，本文研究了結(jié)合機(jī)器學(xué)習(xí)和變分法的移動(dòng)網(wǎng)格方法，在不增加網(wǎng)格節(jié)點(diǎn)數(shù)且不改變拓?fù)浣Y(jié)構(gòu)的前提下，移動(dòng)網(wǎng)格節(jié)點(diǎn)，具體方法為：基于DG 有限元方法，在原始網(wǎng)格上獲得間斷量指標(biāo)統(tǒng)計(jì)值，采用機(jī)器學(xué)習(xí)方法進(jìn)行一次性的R 型網(wǎng)格自適應(yīng)，同時(shí)采用Laplacian 平滑法保證網(wǎng)格質(zhì)量，并在調(diào)整后的網(wǎng)格上進(jìn)行最終的非定常計(jì)算。整個(gè)過程不涉及幾何守恒律問題，使得該網(wǎng)格自適應(yīng)方法獨(dú)立于CFD 計(jì)算，同時(shí)減少了非定常網(wǎng)格自適應(yīng)的時(shí)間耗費(fèi)。

本文通過圓柱繞流算例比較了調(diào)整前后二維圓柱繞流網(wǎng)格的非定常計(jì)算結(jié)果，渦街區(qū)域相對(duì)加密的自適應(yīng)網(wǎng)格更早出現(xiàn)非定?，F(xiàn)象，且流場(chǎng)結(jié)構(gòu)更為清晰。該實(shí)驗(yàn)結(jié)果表明，間斷量指標(biāo)統(tǒng)計(jì)值可作為非定常流場(chǎng)網(wǎng)格自適應(yīng)的參考依據(jù)，通過結(jié)合機(jī)器學(xué)習(xí)使得網(wǎng)格調(diào)整自動(dòng)化、智能化，調(diào)整后的網(wǎng)格在保證計(jì)算效率的同時(shí)，提高了CFD 計(jì)算精度。

本文網(wǎng)格自適應(yīng)方法以三角形（二維）為基礎(chǔ)單元，任意網(wǎng)格單元可以剖分成三角形單元進(jìn)行網(wǎng)格自適應(yīng)，且移動(dòng)過程中不產(chǎn)生負(fù)體積問題，能夠保持原始網(wǎng)格的拓?fù)浣Y(jié)構(gòu)。因此，所提方法對(duì)網(wǎng)格單元沒有限制，適用于結(jié)構(gòu)網(wǎng)格、非結(jié)構(gòu)網(wǎng)格以及混合網(wǎng)格。后續(xù)可以將本文方法拓展至三維，以四面體為基本單元進(jìn)行網(wǎng)格自適應(yīng)。