非結構網格體心梯度求解方法的精度分析

趙 輝,張耀冰,陳江濤,鄧有奇

(中國空氣動力研究與發展中心,綿陽621000)

0 引 言

近些年來,使用非結構混合網格的計算流體力學(Computational Fluid Dynamics,CFD)在工程模擬中得到了廣泛的應用。非結構混合網格,沒有網格拓撲的限制,能夠高效地處理復雜外形問題。網格生成時需要較少的人工干預,網格生成的時間相對于多塊結構網格也大大減少,計算前期準備的工作量和所需時間顯著減少。非結構網格另一個比較吸引人的特點是便于開展基于流場解的網格自適應[1]。在兩個著名的航空標模計算會議、阻力會議和高升力會議上,超過一半的參與者使用基于非結構網格的流場計算器[2-16]。

現有的面向工程應用的非結構網格CFD求解器大多是基于空間二階精度的有限體積方法。為了達到空間的二階精度,使用分段線性重構多項式來近似流場變量的分布。梯度求解的精度在很大程度上決定了求解器的整體空間精度。常用的梯度方法主要有格林高斯法和最小二乘法兩種方法。

格林高斯法是利用格林高斯定理,將流場變量在體心的梯度轉換為變量在單元界面上的積分。格林高斯法最大的優點是它與計算通量的過程相似,因此不需要增加新的數據結構,可以提高計算效率。求解過程中只使用了與當前單元共面的單元信息,計算過程比較緊致,比較容易實現大規律并行計算。但是格林高斯法的缺點也很明顯。計算中只使用了共面的單元信息,這樣較窄的模板有可能導致較小的數值耗散,影響計算的魯棒性[17-18]。

為了改善傳統格林高斯法的求解精度和魯棒性問題,有學者提出了節點格林高斯法[17,19]。這種方法的思路是面心(二維情況下為線中點)的值用該面所有頂點的算術平均得到,而頂點的值又利用它周圍所有單元的體心值通過加權平均得到,這樣就充分利用了周圍控制體的所有已知信息,模板相較傳統的格林高斯法更寬。

最小二乘法[18,20-23]可以認為是一種k-exact重構(k＝1),是將流場變量在當前單元體心處做一階Taylor級數展開,隨后將該展開延拓到周圍鄰近單元(共面或者共點單元),并且要求滿足該單元的平均值約束條件。這樣得到一個未知量是體心梯度的最小二乘問題。

在現階段,關于梯度求解方法的文獻有很多,但是大多都是從數值測試的方面來評估各種方法的數值表現[17-27],比如Sozer[27]測試了幾種常用的梯度求解方法的數值表現,也得到了在任意網格拓撲下能保持一階精度的梯度求解方法。但是他們沒有過多涉及到方法本身的理論精度問題。國內的王年華等[30]在理論分析方面做了一定工作,側重點在數值測試方面,并對網格質量和網格類型的影響做了一定分析。本文結合前人的工作,重點在理論方面分析了幾種常見的梯度求解方法的精度,給出了精度推導的具體過程,并從數值方面對上述理論精度進行了驗證。在數值測試過程中,通過以當前單元體心為基準進行坐標局部縮放的做法,保證了在非結構網格上做精度測試時,網格拓撲能夠嚴格保持不變。這些研究內容對非結構網格梯度求解方法的理論分析和對比,具有一定的借鑒意義。本文的分析和數值驗證都是基于格心的有限體積法展開的。

1 梯度求解方法的理論精度推導

在本節中主要推導了幾種常用的梯度求解方法的理論精度,主要有格林高斯法、節點格林高斯法和最小二乘法。首先在二維情況下做了推導,最后將推導推廣到三維。

本文的研究考慮了直角坐標系下的Euler方程,控制方程為:

其中:U為守恒變量,F、G、H為無黏通量,定義如下:

在有限體積法中,將上述控制方程在每個網格單元Ωi內積分,得到:

在非結構網格中,為了達到空間二階精度,需要利用守恒變量的平均值信息,完成分段線性重構,即計算出守恒變量在單元體心的梯度。然后利用線性插值得到界面左右兩側的守恒變量,

其中:rL、rR分別表示從左右兩側體心到面心的矢量。最后可以使用成熟的Riemann求解器得到界面處的無黏通量Fn(UL,UR,n)。

1.1 格林高斯法

對于在單元內部和單元界面上一階連續的函數f,根據格林高斯公式可以得到:

其中n為控制體界面的外法線方向。定義梯度在單元上的平均值為表示單元內平均值,S表示控制體面積。

由于單元體心值是單元平均值的二階精度近似,式(4)可以改寫為:

在實際應用過程中,式(5)右端的積分∮fdl一般使用線中點的函數值代替,即:

因此下面的關鍵是求得fmid的近似。如果根據某種加權方法得到＝fmid＋O(Δ2),則根據格林高斯法計算得到的體心梯度為一階精度。否則計算將低于一階精度。

假定f1和f2分別為界面兩側函數的單元平均值或體心函數值(由于體心值是平均值的二階精度近似,兩者交換使用不影響下文的推導),r1和r2分別為從兩側體心到界面中點的矢徑。根據Taylor展開(只取一階項),若兩側體心和線中點在同一條直線上,即可以據此得到:

如果兩側體心和線中點不在同一條直線上,則無論距離權,還是使用面積權或者簡單算術平均,都達不到二階精度。

在很多情況下,特別是網格質量不是特別好的時候,“兩側體心和線中點不在同一條直線上”是經常發生的,這也給計算帶來很大的誤差。此時,即使對于線性分布的函數都不能準確還原fmid。這樣的話,單元體心梯度的求解也達不到一階精度。

1.2 節點格林高斯法

格林高斯法的關鍵是利用左右單元體心的值來求解線中點的值,這種做法雖然簡單,但只使用了與當前單元共面的單元的信息,如圖1左所示,周圍其他單元的信息并沒有用到。為了改善傳統格林高斯法的求解精度和魯棒性問題,有學者提出了節點格林高斯法。這種方法的思路是線中點的值用該線所有頂點的算術平均得到,而頂點的值又利用它周圍所有單元的體心值通過加權平均得到,這樣就充分利用了周圍控制體的所有已知信息,模板相較傳統的格林高斯法更寬,如圖1右所示。

圖1 格林高斯法(左)和節點格林高斯法(右)的模板點示意圖[28]Fig.1 Stencil illustration for Green-Gauss method(Left)and node-based Green-Gauss method(Right)[28]

在節點格林高斯法中,線中點的函數值由該線段端點的函數值算術平均得到:

如果f1和f2均有二階精度,則得到的至少有二階精度。

對于節點值f1和f2,可以根據周圍的節點值通過某種加權算法得到。f1＝,f2＝。加權函數的取法多種多樣,可以取1/rn,也可以通過單元面積加權。但是這些權重都不能達到二階精度,不能準確還原線性函數分布。Rausch[29]提出了一種滿足二階精度、能夠準確對線性函數完成插值(保線性)的權重函數,具體形式如下:

其中當前節點記為0,鄰近體心記為i,(i＝1,2,…,n)。不過需要指出的是,使用這種權函數需要做兩次基于節點的循環,計算量至少比使用距離權多了一倍。

1.3 最小二乘法

利用最小二乘法來估算體心梯度最早是由Barth[20-21]提出,可以理解為k-exact(k＝1)重構[24]。在當前單元及其領域內,假定線性函數分布,該函數分布不僅滿足當前單元的平均值約束,也滿足鄰近單元的平均值約束。二維情況下,當前單元記為0,鄰近單元記為i,(i＝1,2,…,N),滿足單元0上平均值約束的函數分布可以寫為:

將其在鄰近單元i上積分,滿足單元i的平均值約束:

將xi－x0記為Δxi,yi－y0記為Δyi,得到如下的最小二乘問題

其中:wi為權函數,可以取為單元i體心到單元0體心距離的函數。該最小二乘問題可以記為A x＝b。為了避免病態問題,上述最小二乘問題A矩陣使用Gram-Schmidt過程[23]分解為正交矩陣Q和上三角矩陣R的乘積,因此得到x＝R－1QTb。亦可以使用奇異值分解,通過矩陣A的廣義逆來估算體心梯度。

為了分析最小二乘法的精度,通過式(12)的法方程來分析求解精度[26]。

上述最小二乘問題的法方程為:

將fi－f0在(x0,y0)處做Taylor展開得到,經過簡單推導得到:

可見最小二乘法計算體心梯度達到一階精度,而且推導過程中不引入近似。對于任意的網格拓撲,都能準確計算線性函數的梯度。

在最小二乘法求解過程中,可以選擇與當前單元共面的單元,也可以選擇與當前單元共點的單元,這類似于格林高斯法和節點格林高斯法模板集的不同。只使用共面單元,導致模板集過窄,數值耗散較小,工程計算中可能遇到魯棒性問題。使用全部共點單元,計算量相對大一點,但是模板點更寬,計算的穩定性提高。

在實施過程中,權函數的選擇比較關鍵。在大拉伸比網格上,不加權的最小二乘法產生的左端項矩陣條件數過大。使用距離的倒數作為權重的話,可以有效改善矩陣的條件數問題,在各向異性網格上減輕計算的剛性問題,得到更好的梯度近似。

上面的分析給出了三種梯度求解方法的理論精度。格林高斯法中,只有當單元兩側體心和線中點在同一條直線上并且使用距離權插值到中點上時,計算得到的梯度才有一階精度。節點格林高斯法計算梯度的精度取決于節點上的函數值是通過何種加權方法得到。如果使用距離的倒數作為權重,計算得到的梯度不到一階精度。如果使用滿足二階精度的插值方法,計算得到的梯度為一階精度,能夠準確得到線性函數的梯度。但是這樣做付出的代價是必須循環兩次,計算量增大。對于最小二乘法,不管使用共面單元還是共點單元,不管在何種網格拓撲下,只要鄰近單元數目足夠,計算得到的梯度為一階精度,能夠準確得到線性函數的梯度。

從上述二維推導很容易擴展到三維。高斯和節點高斯方法中:

三維情況下,不管單元界面是三角形還是四邊形都有:

最小二乘法的精度分析可以更加直接地從二維推廣到三維,具體細節比較簡單,本文不再給出。在三維情況下,關于梯度求解方法的理論精度分析與二維情況一致。

2 精度測試

在上一節中,對三種梯度求解方法的理論精度做了推導,這一節將首先驗證不同梯度求解方法在網格變形情況下對線性函數梯度的求解精度。然后從數值方面在不同網格拓撲情況下對上節推導的理論精度進行驗證。最后通過低速無黏翼型繞流問題驗證在稍復雜問題下,各種梯度求解方法的表現。

2.1 梯度求解誤差分析

使用線性解析函數來驗證不同梯度求解方法在網格變形情況下預測梯度的能力?？紤]了5種梯度的求解方法,分別為格林高斯法(GreenGauss-CellBased),使用距離倒數權插值到節點的節點格林高斯法(GreenGauss-NodeBased-InverseDistance),使用保線性權插值到節點的節點格林高斯法(GreenGauss-NodeBased-Linear Preserving),使用共面單元的最小二乘法(LeastSquare-CellBased)和使用共點單元的最小二乘法(LeastSquare-NodeBased)。兩種最小二乘法都使用了距離倒數作為權函數。

在[0,1]×[0,1]的二維區域內,生成三角形的初始網格。網格變形方法參考了Cary[18]的處理方法,是通過代碼將初始網格點的位置由(xi,yi)平移到(xi＋(s－1)yi,yi),s為從1開始連續變化的正整數,稱為拉伸比例。在此過程中保持網格拓撲不變。驗證過程中考慮了四種不同的網格拓撲,分別記做Grid A、Grid B、Grid C、Grid D,如圖2所示。這里給出的是初始網格的示意圖。在網格還沒有變形的情況下,就已經有很多“兩側體心和線中點不在同一條直線上”的情況發生。

圖2 用來測試線性函數梯度求解的四套網格Fig.2 Grids used for gradient computation test of linearly varying function

圖3給出了s＝3時Grid C變形后的示意圖,網格在x方向拉伸嚴重。

圖3 網格x向拉伸示意圖Fig.3 Illustration of grid stretching in x direction

在測試中主要考慮不同梯度求解方法的影響,沒有關注邊界的處理。在格林高斯和節點格林高斯法中,在邊界附近單元梯度求解時,需要用到邊界上線段中點和節點的函數值。這里通過解析函數精確給出。在共面的最小二乘法中,有些邊界單元沒有足夠多的內部共面相鄰單元,無法給出預測梯度,在統計誤差時,這部分單元忽略。

測試使用的是線性函數f＝x＋sy。通過上面的理論分析可以預測,格林高斯法和使用距離倒數權插值到節點的節點格林高斯法不能準確預測該函數的梯度值,隨著網格變形加大,誤差更大。其他三種方法都能準確計算該函數梯度。圖4給出了4種網格拓撲下,計算得到的體心梯度誤差的1-范數和2-范數誤差隨著拉伸比例的變化情況,其中fy的誤差是除了拉伸比例s后的結果。格林高斯法和使用距離倒數權插值到節點的節點格林高斯法在四套網格上都不能準確預測體心梯度。特別是在Grid B,Grid C和Grid D上,隨著網格拉伸加大,誤差持續增加。后三種梯度求解方法在四種網格拓撲下,都準確預測了線性函數的體心梯度。隨著網格變形的加劇,計算誤差均在可以接受的范圍內。

圖4 不同梯度求解方法對線性函數梯度求解的誤差Fig.4 Gradient computation error for linear function with different methods

2.2 收斂精度驗證

非結構網格加密的過程中,很難保證局部網格的拓撲不變。對于計算域中的某一個參考點,網格加密后其周圍的網格拓撲很可能發生改變,這樣做數值精度測試是否嚴格是存疑的。Sozer[27]提出了一種新的思路。在測試過程中,保持網格不變,將測試函數在局部坐標系下進行縮放,以此來驗證方法的數值精度。但實施過程中,因為測試函數是定義在局部坐標系下的,過程稍微繁瑣。本文提出了一種新的數值精度驗證的方法。對于圖5所示的單元0(由紅色線段構成),在高斯法和共面最小二乘法中,單元0體心的梯度只與單元0以及與其共面的單元有關,在節點高斯法和共點最小二乘法中,單元0體心的梯度只與單元0以及與其共點的單元有關。將單元0和其相關單元組合在一起,以單元0的體心為基準點來進行局部的坐標縮放,通過單元0體心梯度的預測值和精確值的誤差來驗證方法的數值精度。在此過程中測試函數保持不變。如果只在計算域中的某一點上通過預測誤差來測試精度難免會受到函數分布的影響。在驗證過程中,將單元0及其相關單元在計算域內整體均勻平移,通過誤差的1-范數和2-范數來測試數值精度。驗證使用的測試函數是非線性函數:f＝(1＋x＋y＋xy)(sin2πx＋sin2πy)。

圖5 完成收斂精度測試的網格示意(Grid E)Fig.5 Illustration of grid used for the order of accuracy test(Grid E)

圖6給出了不同的梯度預測方法計算的梯度誤差隨著網格尺度的變化規律,其中h為當前單元0面積的平方根。使用保線性權插值到節點的節點格林高斯法,使用共面單元的最小二乘法和使用共點單元的最小二乘法這三種方法計算的梯度有著近似線性的網格收斂性,證明了梯度預測的精度為一階精度。格林高斯法和使用距離倒數權插值到節點的節點格林高斯法,梯度預測的誤差與預想的不一致。預計中這兩種方法的誤差隨著網格加密會逐漸減小,但收斂速度慢于其他后面三種方法,即梯度預測的精度低于一階。

圖6 不同梯度求解方法的數值精度驗證(Grid E)Fig.6 The order of accuracy for different gradient computation methods on Grid E

可以簡單地從理論上分析這種情況發生的原因。對于格林高斯法,預測的x方向體心梯度從式(3)出發,可以寫為:

其中:i＝1,…,N0表示單元0的邊界,j表示單元0的共面單元,兩者的公共面是i。將fi在當前單元體心0處做Taylor展開,可以得到:

只有當兩側體心和線中點在一條直線上并且使用式(7)的插值權函數時才有:

值得注意的是,式(21)也給出了格林高斯法中梯度求解為一階精度時權函數需要滿足的條件。在二維情況下,式(21)給出了四個約束方程,而未知數的個數為與當前單元共面的鄰近單元的個數。通過最小二乘法求得wi2(i＝1,…,N0)后,可以得到單元體心處的梯度為:

三維情況下有類似的推導。此時約束方程的個數為9個,分別為:

同樣通過最小二乘法求得wi2后,可以得到單元體心處的梯度為:

由于上述權函數是通過最小二乘法得到的,因此在一般情況下不能保證該方法得到的梯度精確滿足一階精度。但該權函數使用了當前單元和所有共面單元的幾何信息,并不像前面提到的距離加權或者面積加權等方法只使用當前單元和某一個共面單元信息。從式(20)可以推斷出,新權函數求解梯度的精度更高。這種在格林高斯法的框架下,使用新權函數求解體心梯度的方法將在我們后續的文章中詳細分析。

對于節點格林高斯法,通過類似的推導可以得到,當節點插值函數需要滿足式(25),體心的梯度預測才有一階精度。

其中:i＝1,…,N0表示單元0的節點編號,j＝1,…,Ni表示擁有節點i的單元編號,m＝1,…,Pi表示屬于單元0并且擁有節點i的邊界編號。當使用保線性權插值到節點上時,權函數滿足:

圖7 完成收斂精度測試的網格示意圖(Grid F)Fig.7 Illustration of grid used for the order of accuracy test(Grid F)

2.3 NACA0012翼型亞聲速繞流

最后一個算例是NACA0012翼型亞聲速無黏繞流。計算使用了逐漸加密的網格來驗證不同梯度求解方法對流場的模擬能力。計算采用了兩套網格,粗網格分別如圖8所示。其中Grid G是近似各向同性的三角形單元,Grid H是將各向異性的四邊形單元剖分成三角形單元。在不同的網格拓撲下進行計算可以更全面地驗證梯度求解方法的精度。計算的馬赫數為0.6,迎角為0°。理論上該算例的阻力系數應該為零,因此在分析中用阻力系數的誤差來分析數值精度。在計算中,使用共面單元的最小二乘法遇到了嚴重的魯棒性問題,無法得到收斂的結果,因此只給出了其他四種梯度求解方法的結果。

體心梯度的求解精度決定了求解器的整體精度。梯度求解為一階精度時,求解器整體有二階精度。梯度求解為零階精度時,求解器整體只有一階精度。圖9給出了不同梯度求解方法預測的阻力系數。坐標都是在對數坐標下的表示。橫坐標為網格的特征尺度,取為h＝1/N,其中N為總的單元數目。格林高斯法隨網格加密有明顯的非線性變化,說明此時還未完全進入網格收斂解區域,根據較密的兩套網格估算的數值精度也在一階左右。其他三種方法阻力系數隨著網格加密線性變化,且數值精度在二階左右。在一般的網格下,使用距離倒數權插值到節點的節點格林高斯法求解梯度時只有0階梯度,但是如果節點周圍的近似滿足:＝0,則該種方法計算得到的梯度也有一階精度,此時阻力系數的預測是二階精度。表2給出了在網格G上通過阻力系數計算得到的空間離散精度,進一步證實了上述的觀察。

圖8 NACA0012翼型繞流網格示意圖Fig.8 Illustration of grids used for simulation of flow around NACA0012 airfoil

表1 梯度求解方法的數值精度驗證(Grid F)Table 1 The order of accuracy for different methods on Grid F

為了進一步判別不同梯度求解方法對流動結構的刻畫能力,圖10和圖11分別給出了格林高斯法和使用共點單元的最小二乘法得到的馬赫數云圖和熵增云圖,其中熵增定義為:

在本算例中,流動是等熵的,理論上在計算域的任何位置均有ΔS＝0。這里只給出了在Grid H序列中細網格上的結果,在其他網格上的結果與此類似。雖然兩種梯度求解方法的收斂精度不同,但從圖10可以看到,兩者得到了基本一致的光滑的馬赫數分布。圖11的熵增云圖可以看出明顯的區別。格林高斯法得到的流場中在翼型附近和下游區域熵增明顯較大,而使用共點單元的最小二乘法得到的流場中熵增區域顯著較小。這與阻力系數的預測一起證明了格林高斯方法求解梯度的誤差較最小二乘法偏大。

圖9 NACA0012翼型繞流問題預測的阻力系數Fig.9 Drag coefficient for inviscid transonic flow around NACA0012 airfoil

表2 梯度求解方法的數值精度驗證(Grid G)Table 2 The order of accuracy for different methods(Grid G)

圖10 NACA0012翼型繞流的馬赫數云圖Fig.10 Computed Mach number contours for flow around NACA0012 airfoil

3 結 論

本文對二階精度非結構混合網格求解器中幾種常用的梯度求解方法進行了深入的理論分析。從理論精度方面,對幾種方法的截斷誤差精度進行了推導。使用保線性權插值到節點的節點格林高斯法、使用共面單元的最小二乘法和使用共點單元的最小二乘法,不管網格拓撲關系如何,都能保證梯度求解為一階精度。而格林高斯法和使用距離倒數權插值到節點的節點格林高斯法只有在特定的網格和權函數下,才能有一階精度。在一般的網格中,只有零階精度。

在得到上述理論之后,對各種梯度求解方法做了進一步的數值測試。首先驗證了不同方法在網格持續拉伸變形情況下,對線性函數梯度求解的誤差。格林高斯法和使用距離倒數權插值到節點的節點格林高斯法計算的梯度隨著網格拉伸,計算的梯度急劇惡化。而其他三種方法都能準確計算線性函數的梯度。

隨后從數值精度方面對上述的理論精度進行了驗證。這里提出了一種新的精度測試方法。以當前單元體心為基準點進行局部網格的縮放,可以精確地保持網格加密過程中的網格拓撲不變。數值測試驗證了上述推導的理論精度。對于格林高斯法和節點格林高斯法,推導得到了其滿足一階精度要求的權函數與網格尺度的約束關系。從這一約束出發,可以驗證某種權函數能否滿足一階精度要求。

最后通過NACA0012翼型亞聲速無黏繞流進一步分析了各種梯度求解方法的表現。使用保線性權插值到節點的節點格林高斯法和使用共點單元的最小二乘法求得的體心梯度均為一階精度,這樣解算器的整體精度為二階精度。使用距離倒數權插值到節點的節點格林高斯法只有在網格各向同性的情況下才有一階精度的體心梯度求解,否則低于一階。格林高斯法在一般的網格下體心梯度低于一階,解算器的整體精度也低于二階。使用共面單元的最小二乘法因為只使用共面單元,魯棒性較差,不能得到收斂的結果。

本文只是從梯度求解的精度方面對幾種方法進行了分析。在實際工程應用中,計算的魯棒性、計算量、并行傳輸的數據量等都是需要仔細考察的因素。后續文章將從這些方面對梯度求解方法進一步分析,并考察在邊界層等大拉伸比單元上的表現。文中也提出了一種在格林高斯法的框架下,盡可能提高梯度求解精度的權函數求解方法。在后續的文章中將對這一方法深入研究和分析。