CFD中統(tǒng)計(jì)誤差的數(shù)值精度分析

閔耀兵，馬燕凱，李松

中國(guó)空氣動(dòng)力研究與發(fā)展中心，綿陽(yáng) 621000

計(jì)算流體力學(xué)(CFD)是一門利用數(shù)值方法模擬流動(dòng)問(wèn)題的科學(xué)，其所依賴的數(shù)值方法在理論設(shè)計(jì)時(shí)需滿足一定的精度要求，而理論設(shè)計(jì)的精度在CFD應(yīng)用中不一定能夠達(dá)到，因此經(jīng)常需要對(duì)算法的理論設(shè)計(jì)精度進(jìn)行數(shù)值驗(yàn)證[1-6]。算法精度的數(shù)值驗(yàn)證過(guò)程多以加密網(wǎng)格的方式進(jìn)行，即通過(guò)逐步加密網(wǎng)格，使得數(shù)值誤差逐漸減小，由網(wǎng)格加密前后數(shù)值誤差的比值和網(wǎng)格尺度的比值可以計(jì)算出算法的數(shù)值精度[1-6]。一般情況下，相比于選取某固定點(diǎn)的數(shù)值誤差而言，對(duì)整個(gè)計(jì)算域內(nèi)的數(shù)值誤差進(jìn)行統(tǒng)計(jì)賦范更具有代表性和普適性。通常情況下，數(shù)值誤差的統(tǒng)計(jì)方式一般可取其L1范數(shù)、L2范數(shù)和L∞范數(shù)，并認(rèn)為統(tǒng)計(jì)誤差的各范數(shù)在描述算法的數(shù)值精度上具有等價(jià)性。因此，在CFD中對(duì)算法進(jìn)行數(shù)值精度驗(yàn)證時(shí)，有些研究人員僅考察統(tǒng)計(jì)誤差的L1范數(shù)[6-7]或者L2范數(shù)[8-10]，也有同時(shí)考察L1范數(shù)和L∞范數(shù)的[4-5]以及同時(shí)考察L2范數(shù)和L∞范數(shù)的[11]，當(dāng)然也有更為嚴(yán)謹(jǐn)?shù)难芯空咄瑫r(shí)基于L1范數(shù)、L2范數(shù)和L∞范數(shù)等3種常用范數(shù)進(jìn)行數(shù)值精度分析的[1-2,12]，卻鮮有文獻(xiàn)僅參考L∞范數(shù)。

然而，在對(duì)某具體的CFD算法的理論設(shè)計(jì)精度進(jìn)行數(shù)值驗(yàn)證時(shí)，由于CFD算法的理論設(shè)計(jì)多基于光滑流場(chǎng)假設(shè)，在遭遇流場(chǎng)間斷(如激波和接觸間斷等)[8]、網(wǎng)格不連續(xù)(如拐折和突然拉伸等)[12-13]等問(wèn)題時(shí)，算法的數(shù)值精度可能會(huì)存在不同程度的降階問(wèn)題[8,11-13]。即便是在光滑流場(chǎng)中，在極值點(diǎn)附近采用非線性加權(quán)插值也可能會(huì)產(chǎn)生降階問(wèn)題[14]。諸如此類的降階問(wèn)題通常只會(huì)出現(xiàn)于計(jì)算區(qū)域中的某個(gè)局部，不會(huì)導(dǎo)致整個(gè)計(jì)算區(qū)域的數(shù)值精度全部降階[13]。此時(shí)統(tǒng)計(jì)誤差的各范數(shù)在數(shù)值精度上一般不再具有等價(jià)性，各范數(shù)的數(shù)值精度與其賦范方式有關(guān)，不同賦范方式可能會(huì)給出不同的數(shù)值精度[13]，且一般情況下L1范數(shù)的數(shù)值精度最高，而L∞范數(shù)的數(shù)值精度最低，其中的原因和關(guān)系本文將給出詳細(xì)分析。

針對(duì)統(tǒng)計(jì)誤差各范數(shù)的數(shù)值精度不一致的問(wèn)題[13]，本文通過(guò)對(duì)統(tǒng)計(jì)誤差的具體賦范方式進(jìn)行理論分析，指出：當(dāng)且僅當(dāng)整個(gè)計(jì)算區(qū)域內(nèi)各點(diǎn)的數(shù)值精度完全一致時(shí)，統(tǒng)計(jì)誤差各范數(shù)的數(shù)值精度才相等，此時(shí)統(tǒng)計(jì)誤差的各范數(shù)在數(shù)值精度上具有等價(jià)性；而當(dāng)全場(chǎng)的數(shù)值精度不完全一致時(shí)，統(tǒng)計(jì)誤差的各范數(shù)的數(shù)值精度并不一致，其具體精度與賦范方式有關(guān)，且一般情況下表現(xiàn)為L(zhǎng)1范數(shù)的數(shù)值精度最高，而L∞范數(shù)的數(shù)值精度最低。本文的理論分析結(jié)果很好地解釋了統(tǒng)計(jì)誤差各范數(shù)的數(shù)值精度之間的關(guān)系，同時(shí)也為CFD算法精度驗(yàn)證時(shí)的誤差統(tǒng)計(jì)方式提供了理論參考。

在理論分析的基礎(chǔ)上，本文分別針對(duì)全場(chǎng)各點(diǎn)精度不完全一致和完全一致的情形進(jìn)行了驗(yàn)證，驗(yàn)證的結(jié)果與理論分析結(jié)論相符。然后對(duì)極值點(diǎn)附近的非線性加權(quán)插值引起的降階問(wèn)題進(jìn)行了具體分析，利用本文的理論分析結(jié)論，從數(shù)值精度上較好地解釋了非線性加權(quán)插值在極值點(diǎn)附近的降階問(wèn)題。

1 統(tǒng)計(jì)誤差的數(shù)值精度分析

當(dāng)計(jì)算區(qū)域中離散的網(wǎng)格尺度為h時(shí)，記整個(gè)離散區(qū)域中誤差統(tǒng)計(jì)的總點(diǎn)數(shù)為Nh，則隨著計(jì)算網(wǎng)格的加密，網(wǎng)格尺度h減小而總點(diǎn)數(shù)Nh增大，即Nh與1/h成正變關(guān)系。將每個(gè)點(diǎn)的數(shù)值誤差ei簡(jiǎn)單表示為

ei=cihrii=1,2,…,Nh

(1)

式中：ri為第i點(diǎn)的數(shù)值精度，一般情況下ri均為整數(shù)(?i)；ci為與網(wǎng)格尺度h無(wú)關(guān)的誤差常數(shù)，一般可以表示為求解變量及其導(dǎo)數(shù)的函數(shù)。

(2)

(3)

(4)

1.1 各點(diǎn)精度不完全一致的情形

在實(shí)際CFD計(jì)算中，由于離散邊界、流場(chǎng)間斷以及網(wǎng)格不連續(xù)等問(wèn)題的存在，經(jīng)常會(huì)出現(xiàn)各點(diǎn)數(shù)值精度ri并不完全相同的情況，即有

min(ri)

(5)

為便于分析，不妨設(shè)整個(gè)計(jì)算區(qū)域中總共存在n個(gè)數(shù)值精度(當(dāng)各點(diǎn)數(shù)值精度不完全相同時(shí)，n≥2)，其數(shù)值由小到大分別記為R1，R2，…，Rn，即

R1=min(ri)

(6)

一般情況下，隨著計(jì)算網(wǎng)格的加密，精度值Rj以及其數(shù)量n均保持不變。不妨記數(shù)值精度為Ri的數(shù)值誤差的點(diǎn)數(shù)為Ni，則有

(7)

(8)

不妨記

(9)

由于ci為與網(wǎng)格尺度h無(wú)關(guān)的常系數(shù)，故隨著計(jì)算網(wǎng)格的加密，Cj也基本保持不變。由式(9)有

(10)

在網(wǎng)格加密過(guò)程中，不妨設(shè)網(wǎng)格尺度h的加密倍數(shù)rfc為

(11)

式中：下標(biāo)c和f分別對(duì)應(yīng)粗網(wǎng)格和細(xì)網(wǎng)格，網(wǎng)格的加密過(guò)程使得rfc>1。當(dāng)網(wǎng)格尺度加密倍數(shù)為rfc時(shí)，整個(gè)計(jì)算區(qū)域中誤差統(tǒng)計(jì)的總點(diǎn)數(shù)Nh的變化關(guān)系為

(12)

式中：dN為計(jì)算網(wǎng)格的加密維數(shù)，dN=1,2,3分別對(duì)應(yīng)一維、二維和三維網(wǎng)格。同時(shí)記數(shù)值精度為Rj階的誤差點(diǎn)數(shù)Nj隨網(wǎng)格的加密滿足：

(13)

由于式(7)及Nj

dj≤dNj=1,2,…,n

(14)

隨著計(jì)算網(wǎng)格的加密，誤差統(tǒng)計(jì)的總點(diǎn)數(shù)Nh要比各精度的點(diǎn)數(shù)Nj增加得更快。

根據(jù)定義(12)及(13)，容易得到

(15)

(16)

為便于分析，不妨記

(17)

進(jìn)一步容易得到

(18)

(19)

(20)

將式(18)代入(20)中，并整理得到

(21)

由上述分析以及網(wǎng)格尺度關(guān)系式(11)，容易得到

(22)

由于dN≤3以及式(14)，容易知道dj≤3(?j)，特別地，對(duì)于本文數(shù)值算例中考慮的二維問(wèn)題(dN=2)有dj≤2，又m≥1，則式(22)中的數(shù)值精度Om在絕大多數(shù)情況下可以表示為

(23)

當(dāng)m趨于正無(wú)窮大時(shí)，有

(24)

考慮到一般情況下式(14)及式(23)成立，故有

O1≥O2≥…≥O∞

(25)

即當(dāng)各點(diǎn)計(jì)算精度不完全相同時(shí)，其統(tǒng)計(jì)誤差L1范數(shù)的數(shù)值精度O1最高，L2范數(shù)的數(shù)值精度O2次之，以此類推，L∞范數(shù)的數(shù)值精度O∞最低。

1.2 各點(diǎn)精度一致的情形

當(dāng)全場(chǎng)精度均一致時(shí)，有

min(ri)=max(ri)=R

(26)

類似于式(8)中的分析，容易得到

(27)

式中：

(28)

類似于Cj、Ch也為與網(wǎng)格尺度h無(wú)關(guān)的常數(shù)。此時(shí)統(tǒng)計(jì)誤差Lm范數(shù)的數(shù)值精度Om可以表示為

(29)

考慮到網(wǎng)格間距關(guān)系式(11)容易得到

Om=R

(30)

與各點(diǎn)計(jì)算精度不完全相同時(shí)的情形不一樣，當(dāng)各點(diǎn)計(jì)算精度均相同時(shí)，其統(tǒng)計(jì)誤差的各范數(shù)均為R階精度，此時(shí)統(tǒng)計(jì)誤差的各范數(shù)之間在數(shù)值精度上具有等價(jià)性。

2 數(shù)值算例

針對(duì)上述關(guān)于統(tǒng)計(jì)誤差數(shù)值精度的理論分析，設(shè)計(jì)相應(yīng)的數(shù)值試驗(yàn)進(jìn)行驗(yàn)證。類似于Mao等[13]的做法，設(shè)計(jì)不同連續(xù)程度的網(wǎng)格，計(jì)算各點(diǎn)的幾何守恒律誤差并進(jìn)行誤差統(tǒng)計(jì)，通過(guò)逐步加密網(wǎng)格以得到統(tǒng)計(jì)誤差的各范數(shù)的數(shù)值精度表現(xiàn)。

2.1 網(wǎng)格生成

本文中采用的周期性二維計(jì)算網(wǎng)格基于如下算法解析生成：

doj=1,Nj

doi=1,Ni

if(|x0|<0.4) then

y2=sinκ(x2)

else

y2=0

end if

xi,j=x0

yi,j=y0+A·y1·y2

end do

end do

(31)

式中：Ni和Nj分別為計(jì)算坐標(biāo)下ξ方向和η方向的網(wǎng)格分布點(diǎn)數(shù)，實(shí)際網(wǎng)格生成過(guò)程中取Ni=Nj=N+1；A為波動(dòng)幅值(A=0時(shí)為均勻直角網(wǎng)格)，在本文中均取A=0.2；函數(shù)y2的指數(shù)κ取不同值(在本文中指數(shù)κ取1,2,3)時(shí)可以得到在x=±0.4處不同連續(xù)程度的網(wǎng)格，函數(shù)y2描述的網(wǎng)格線如圖1所示。

圖1 函數(shù)y2分布

2.2 數(shù)值離散

二維網(wǎng)格中幾何守恒律誤差可以表示為(以Ix為例)

(32)

式中：

(33)

(34)

將式(34)的差分算子δ均取為如下的二階精度中心差分格式：

(35)

(36)

2.3 結(jié)果分析

(37)

(38)

同理有

(39)

而當(dāng)κ=3時(shí)，由于R1=Rn，則由式(26)對(duì)應(yīng)的分析易知

(40)

按照算法(31)生成的網(wǎng)格均滿足d1=1，下面本文考慮另外一種網(wǎng)格中統(tǒng)計(jì)誤差的數(shù)值精度問(wèn)題。取κ=3，在算法(31)的基礎(chǔ)上，按照如下方法修正x=±0.4、y=0處的網(wǎng)格分布：

表1 幾何守恒誤差的統(tǒng)計(jì)精度(κ=1)

表2 幾何守恒誤差的統(tǒng)計(jì)精度(κ=2)

表3 幾何守恒誤差的統(tǒng)計(jì)精度(κ=3)

if(|x0|=0.4&y0=0) then

end if

(41)

式中：γ為網(wǎng)格間距指數(shù)；B為與網(wǎng)格間距無(wú)關(guān)的常數(shù)，其取值方式為

(42)

式中：參數(shù)B的選取使得在x=±0.4、y=0處的網(wǎng)格波動(dòng)在加密過(guò)程中不至于被機(jī)器字長(zhǎng)所湮沒(méi)(雙精度浮點(diǎn)運(yùn)算)。按照算法(41)生成的計(jì)算網(wǎng)格在間斷附近的分布如圖2所示。

圖2 局部網(wǎng)格示意圖(γ=1)

表4 幾何守恒誤差的統(tǒng)計(jì)精度(γ=1)

(43)

(44)

(45)

(46)

基于算法(41)生成的網(wǎng)格的數(shù)值試驗(yàn)結(jié)果見(jiàn)表4～表6。當(dāng)γ=1和γ=2時(shí)，表4和表5中所示的統(tǒng)計(jì)誤差的各范數(shù)對(duì)應(yīng)的數(shù)值精度分別與表達(dá)式(44)和式(45)中的完全相同。當(dāng)γ=3時(shí)，表6中給出的統(tǒng)計(jì)誤差的各范數(shù)的數(shù)值精度與理論分析(46)相同，值得注意的是表6中L1范數(shù)的數(shù)值精度為二階，與理論分析(22)相同，而式(23)此時(shí)給出的數(shù)值精度應(yīng)為三階，此時(shí)式(22)和式(23)給出的數(shù)值精度并不相同。

2.4 算例應(yīng)用

針對(duì)WENO格式的非線性加權(quán)插值[4]在極值點(diǎn)附近可能會(huì)降階的問(wèn)題，Henrick等[14]設(shè)計(jì)了相應(yīng)的數(shù)值算例進(jìn)行驗(yàn)證，并指出WENO格式[4]中為避免分母為零而引入的小量ε的取值對(duì)數(shù)值精度的影響：ε取值相對(duì)較大時(shí)(如ε=10-6)，非線性權(quán)近似于其線性最優(yōu)權(quán)，WENO格式在極值點(diǎn)附近一般不會(huì)降階；而ε取值較小時(shí)(如ε=10-40)，非線性權(quán)遠(yuǎn)遠(yuǎn)偏離其線性最優(yōu)權(quán)，WENO格式在極值點(diǎn)附近會(huì)存在降階問(wèn)題。針對(duì)目前已得到廣泛應(yīng)用且同樣采用非線性加權(quán)思想構(gòu)造的WCNS-E-5格式[3,18-19]，有必要測(cè)試其在極值點(diǎn)附近的數(shù)值精度情況[18]，為了準(zhǔn)確反映WCNS-E-5格式[3]的非線性權(quán)在極值點(diǎn)附近的特性，取ε=10-300(四倍精度浮點(diǎn)運(yùn)算)。

依據(jù)算法(31)取A=0生成均勻網(wǎng)格，參照Henrick等[14]的測(cè)試方法給定周期性初值：

(47)

表5 幾何守恒誤差的統(tǒng)計(jì)精度(γ=2)

表6 幾何守恒誤差的統(tǒng)計(jì)精度(γ=3)

采用WCNS-E-5格式[3]離散ux并將其相對(duì)于解析值的數(shù)值誤差統(tǒng)計(jì)于表7中,表中的統(tǒng)計(jì)誤差的數(shù)值精度滿足：

(48)

值得指出的是，Henrick等[14]的數(shù)值實(shí)驗(yàn)中得出的統(tǒng)計(jì)誤差的數(shù)值精度與式(48)所述一致，只是Henrick等[14]只給出了m=1,2,∞這3種情況下的精度，且沒(méi)有對(duì)統(tǒng)計(jì)誤差的L1、L2和L∞范數(shù)所表現(xiàn)的不同的數(shù)值精度給出解釋與說(shuō)明。

由統(tǒng)計(jì)誤差的數(shù)值精度式(22)和式(23)易知R1=3，即統(tǒng)計(jì)區(qū)域中的最低數(shù)值精度為三階。對(duì)于周期性初值式(47)而言，在極值點(diǎn)附近u′=0、u″≠0、u?≠0，Henrick等[14]的分析指出此時(shí)非線性加權(quán)的WENO格式只有三階精度，則采用類似非線性加權(quán)方式構(gòu)造的WCNS格式在此類極值點(diǎn)附近也應(yīng)僅具有三階精度[18]，與表7中的最低數(shù)值精度為三階相吻合。

由于本算例基于二維網(wǎng)格，故dN=2，由式(22)易知d1=1，即在極值點(diǎn)附近WCNS格式[3]的降階以條帶狀形式出現(xiàn)，條帶的寬度與格式的計(jì)算模板有關(guān)，但條帶寬度不影響統(tǒng)計(jì)誤差的數(shù)值精度。

表7 統(tǒng)計(jì)誤差的數(shù)值精度

3 結(jié) 論

針對(duì)數(shù)值誤差統(tǒng)計(jì)時(shí)采用的不同賦范方式可能表現(xiàn)出不同數(shù)值精度的問(wèn)題，從理論上詳細(xì)分析了統(tǒng)計(jì)誤差賦范方式對(duì)其數(shù)值精度的影響，并給出了一般規(guī)律：當(dāng)且僅當(dāng)整個(gè)計(jì)算區(qū)域內(nèi)的數(shù)值精度完全一致時(shí)，統(tǒng)計(jì)誤差各范數(shù)的數(shù)值精度才相等，此時(shí)統(tǒng)計(jì)誤差的各范數(shù)在數(shù)值精度上才具有等價(jià)性；而當(dāng)全場(chǎng)的數(shù)值精度不完全一致時(shí)，統(tǒng)計(jì)誤差各范數(shù)的數(shù)值精度并不一致，其具體精度與賦范方式有關(guān)，且一般情況下表現(xiàn)為L(zhǎng)1范數(shù)的數(shù)值精度最高，L2范數(shù)的數(shù)值精度次之，以此類推，L∞范數(shù)的數(shù)值精度最低。

基于統(tǒng)計(jì)誤差數(shù)值精度的理論分析結(jié)論，本文有針對(duì)性地設(shè)計(jì)了相應(yīng)的數(shù)值試驗(yàn)予以驗(yàn)證，對(duì)于全場(chǎng)精度完全一致和不完全一致的情況，數(shù)值試驗(yàn)均給出了與理論分析完全吻合的結(jié)果，表明本文的研究結(jié)論能夠?yàn)镃FD中算法的精度驗(yàn)證提供理論依據(jù)。