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

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

1.1 格林高斯法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中當(dāng)前節(jié)點(diǎn)記為0,鄰近體心記為i,(i＝1,2,…,n)。不過需要指出的是,使用這種權(quán)函數(shù)需要做兩次基于節(jié)點(diǎn)的循環(huán),計(jì)算量至少比使用距離權(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)得到:

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

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

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

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

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

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

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

2 精度測試

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

2.1 梯度求解誤差分析

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

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

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

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

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

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

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

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

2.2 收斂精度驗(yàn)證

非結(jié)構(gòu)網(wǎng)格加密的過程中,很難保證局部網(wǎng)格的拓?fù)洳蛔?。對于?jì)算域中的某一個(gè)參考點(diǎn),網(wǎng)格加密后其周圍的網(wǎng)格拓?fù)浜芸赡馨l(fā)生改變,這樣做數(shù)值精度測試是否嚴(yán)格是存疑的。Sozer[27]提出了一種新的思路。在測試過程中,保持網(wǎng)格不變,將測試函數(shù)在局部坐標(biāo)系下進(jìn)行縮放,以此來驗(yàn)證方法的數(shù)值精度。但實(shí)施過程中,因?yàn)闇y試函數(shù)是定義在局部坐標(biāo)系下的,過程稍微繁瑣。本文提出了一種新的數(shù)值精度驗(yàn)證的方法。對于圖5所示的單元0(由紅色線段構(gòu)成),在高斯法和共面最小二乘法中,單元0體心的梯度只與單元0以及與其共面的單元有關(guān),在節(jié)點(diǎn)高斯法和共點(diǎn)最小二乘法中,單元0體心的梯度只與單元0以及與其共點(diǎn)的單元有關(guān)。將單元0和其相關(guān)單元組合在一起,以單元0的體心為基準(zhǔn)點(diǎn)來進(jìn)行局部的坐標(biāo)縮放,通過單元0體心梯度的預(yù)測值和精確值的誤差來驗(yàn)證方法的數(shù)值精度。在此過程中測試函數(shù)保持不變。如果只在計(jì)算域中的某一點(diǎn)上通過預(yù)測誤差來測試精度難免會(huì)受到函數(shù)分布的影響。在驗(yàn)證過程中,將單元0及其相關(guān)單元在計(jì)算域內(nèi)整體均勻平移,通過誤差的1-范數(shù)和2-范數(shù)來測試數(shù)值精度。驗(yàn)證使用的測試函數(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ù)測方法計(jì)算的梯度誤差隨著網(wǎng)格尺度的變化規(guī)律,其中h為當(dāng)前單元0面積的平方根。使用保線性權(quán)插值到節(jié)點(diǎn)的節(jié)點(diǎn)格林高斯法,使用共面單元的最小二乘法和使用共點(diǎn)單元的最小二乘法這三種方法計(jì)算的梯度有著近似線性的網(wǎng)格收斂性,證明了梯度預(yù)測的精度為一階精度。格林高斯法和使用距離倒數(shù)權(quán)插值到節(jié)點(diǎn)的節(jié)點(diǎn)格林高斯法,梯度預(yù)測的誤差與預(yù)想的不一致。預(yù)計(jì)中這兩種方法的誤差隨著網(wǎng)格加密會(huì)逐漸減小,但收斂速度慢于其他后面三種方法,即梯度預(yù)測的精度低于一階。

圖6 不同梯度求解方法的數(shù)值精度驗(yàn)證(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è)體心和線中點(diǎn)在一條直線上并且使用式(7)的插值權(quán)函數(shù)時(shí)才有:

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

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

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

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

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

其中:i＝1,…,N0表示單元0的節(jié)點(diǎn)編號,j＝1,…,Ni表示擁有節(jié)點(diǎn)i的單元編號,m＝1,…,Pi表示屬于單元0并且擁有節(jié)點(diǎn)i的邊界編號。當(dāng)使用保線性權(quán)插值到節(jié)點(diǎn)上時(shí),權(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翼型亞聲速繞流

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

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

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

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

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

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

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

表2 梯度求解方法的數(shù)值精度驗(yàn)證(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)格求解器中幾種常用的梯度求解方法進(jìn)行了深入的理論分析。從理論精度方面,對幾種方法的截?cái)嗾`差精度進(jìn)行了推導(dǎo)。使用保線性權(quán)插值到節(jié)點(diǎn)的節(jié)點(diǎn)格林高斯法、使用共面單元的最小二乘法和使用共點(diǎn)單元的最小二乘法,不管網(wǎng)格拓?fù)潢P(guān)系如何,都能保證梯度求解為一階精度。而格林高斯法和使用距離倒數(shù)權(quán)插值到節(jié)點(diǎn)的節(jié)點(diǎn)格林高斯法只有在特定的網(wǎng)格和權(quán)函數(shù)下,才能有一階精度。在一般的網(wǎng)格中,只有零階精度。

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

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

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

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