非結(jié)構(gòu)網(wǎng)格體心梯度求解方法的精度分析

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

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

0 引 言

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

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

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

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

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

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

1 梯度求解方法的理論精度推導(dǎo)

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

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

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

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

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

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

1.1 格林高斯法

對于在單元內(nèi)部和單元界面上一階連續(xù)的函數(shù)f,根據(jù)格林高斯公式可以得到:

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

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

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

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

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

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

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

1.2 節(jié)點格林高斯法

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

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

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

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

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

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

1.3 最小二乘法

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

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

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

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

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

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

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

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

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

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

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

從上述二維推導(dǎo)很容易擴展到三維。高斯和節(jié)點高斯方法中:

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

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

2 精度測試

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

2.1 梯度求解誤差分析

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

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

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

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

圖3 網(wǎng)格x向拉伸示意圖Fig.3 Illustration of grid stretching in x direction

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

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

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

2.2 收斂精度驗證

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

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

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

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

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

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

只有當(dāng)兩側(cè)體心和線中點在一條直線上并且使用式(7)的插值權(quán)函數(shù)時才有:

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

三維情況下有類似的推導(dǎo)。此時約束方程的個數(shù)為9個,分別為:

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

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

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

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

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

2.3 NACA0012翼型亞聲速繞流

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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

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