高階精度有限差分方法幾何守恒律研究進(jìn)展

毛枚良, 姜 屹, 閔耀兵,*, 朱華君, 鄧小剛

(1. 中國空氣動(dòng)力研究與發(fā)展中心, 綿陽 621000；2.軍事科學(xué)院 系統(tǒng)工程研究院, 北京 100082; 3.軍事科學(xué)院, 北京 100091)

0 引 言

近年來，隨著計(jì)算機(jī)浮點(diǎn)運(yùn)算能力的逐步提升,對(duì)湍流等包含時(shí)空多尺度流動(dòng)結(jié)構(gòu)的精細(xì)模擬成為研究熱點(diǎn)，針對(duì)湍流流動(dòng)的分離渦模擬(DES)[1-2]、大渦模擬(LES)[3-4]以及直接數(shù)值模擬(DNS)[5-6]等較為精細(xì)的數(shù)值模擬方法得到了進(jìn)一步發(fā)展和完善。精細(xì)的湍流模擬方法需要嚴(yán)格控制所有離散過程中的數(shù)值誤差，相較于傳統(tǒng)的二階精度計(jì)算方法，高階精度的計(jì)算方法具有明顯優(yōu)勢(shì)。

高階精度應(yīng)用場(chǎng)景的現(xiàn)實(shí)強(qiáng)烈需求催生了一系列的高階精度計(jì)算方法。精度最高的譜方法[7]具有理想的指數(shù)階精度，卻難以應(yīng)用于復(fù)雜外形的流動(dòng)模擬中，目前譜方法多用于對(duì)其他高階精度數(shù)值計(jì)算方法的考核與驗(yàn)證。其他的數(shù)值方法按照其離散框架大致可以分為基于有限差分方法的、基于有限體積方法的和基于有限元方法的高階精度計(jì)算方法，以及近年來逐漸衍生出的DG方法[8-9]、譜差分方法[10-11]以及譜體積方法[12-13]。由于離散框架的不同，不同類型的高階精度數(shù)值方法提高計(jì)算精度的方式迥異，尤其是在多維問題中。得益于逐維離散特性，有限差分方法能以較小計(jì)算代價(jià)實(shí)現(xiàn)多維問題的高階精度計(jì)算[14-15]，但適應(yīng)復(fù)雜幾何外形的能力較弱[16]。當(dāng)網(wǎng)格不夠光滑時(shí)，高階精度有限差分方法還存在明顯的降階問題。因此，在相當(dāng)長的一段時(shí)間內(nèi)，基于有限差分方法的高階精度算法多用于繞簡單外形的類似于湍流等多尺度復(fù)雜流動(dòng)問題的模擬中[17-18]，鮮有在復(fù)雜外形流動(dòng)中成功應(yīng)用的例子。近年來的研究結(jié)果表明[19-22]，基于有限差分方法的高階精度計(jì)算格式對(duì)復(fù)雜網(wǎng)格的適應(yīng)能力較差的主要原因是有限差分方法離散框架下的幾何守恒律不易滿足。

有限差分方法中幾何守恒律(Geometry Conservation Law, GCL)的確切定義最早由Thomas和Lombard[23-24]于1978年提出，而由于幾何守恒未滿足而引起的幾何守恒誤差最早則由Steger于1977年發(fā)現(xiàn)[25-26]。幾何守恒律一般可以分為體積守恒律(Volume Conservation Law,VCL)和面積守恒律(Surface Conservation Law,SCL)[27]，體積守恒律表征了網(wǎng)格運(yùn)動(dòng)過程中網(wǎng)格變換雅克比和網(wǎng)格運(yùn)動(dòng)速度之間需滿足的約束關(guān)系，面積守恒律則反映了網(wǎng)格變換導(dǎo)數(shù)各分量之間需滿足的約束關(guān)系。

Thomas和Lombard[23-24]在首次明確幾何守恒律概念時(shí)著重關(guān)注的正是運(yùn)動(dòng)網(wǎng)格下的體積守恒律問題。在有限差分離散框架下，Thomas等[23-24]在已知網(wǎng)格運(yùn)動(dòng)速度的情況下通過求解微分型的體積守恒律方程得到網(wǎng)格變換雅克比，從而使得體積守恒律得以滿足。但網(wǎng)格變換雅克比在離散后的體積應(yīng)由當(dāng)前網(wǎng)格點(diǎn)的位置唯一確定，Thomas等通過求解體積守恒律方程得到的數(shù)值體積與其定義之間存在相容性的問題。針對(duì)有限差分方法中運(yùn)動(dòng)網(wǎng)格的體積守恒律問題，Abe等[28-29]通過改變與網(wǎng)格運(yùn)動(dòng)相關(guān)的導(dǎo)數(shù)(網(wǎng)格運(yùn)動(dòng)速度)的計(jì)算方式使得體積守恒律在有限差分離散下自動(dòng)滿足。在有限體積離散框架下，Wang[30]和Zhang[31]通過改變網(wǎng)格界面運(yùn)動(dòng)速度的計(jì)算方式來使得體積守恒律自動(dòng)滿足。最近Chang等[32]還對(duì)有限體積方法中運(yùn)動(dòng)網(wǎng)格的幾種體積守恒解決方案進(jìn)行了較為詳細(xì)的分析與比較。從解決體積守恒律問題的角度，Abe等[28-29]的解決方案與Wang[30]和Zhang[31]的方案基本思路一致，只不過分別是在有限差分和有限體積的框架下實(shí)現(xiàn)的。

面積守恒律反映的是各網(wǎng)格變換導(dǎo)數(shù)分量之間需滿足的約束關(guān)系。早在1977年Steger[25-26]就發(fā)現(xiàn)，在采用傳統(tǒng)網(wǎng)格變換導(dǎo)數(shù)算法利用二階中心有限差分格式離散計(jì)算坐標(biāo)系下的守恒律方程時(shí)會(huì)面臨自由流不守恒的問題，其守恒誤差是以二階精度趨于0的。傳統(tǒng)的網(wǎng)格變換導(dǎo)數(shù)算法會(huì)導(dǎo)致流場(chǎng)中存在著完全由計(jì)算網(wǎng)格引起的誤差源(Hindman[33-34]將其稱為幾何誘導(dǎo)誤差)，容易對(duì)計(jì)算流場(chǎng)產(chǎn)生非物理的干擾，嚴(yán)重時(shí)還會(huì)引起數(shù)值振蕩[35]導(dǎo)致計(jì)算過程不穩(wěn)定。為了有效解決這一問題，1978年P(guān)ulliam和Steger[36-37]給出了一種加權(quán)平均的網(wǎng)格變換導(dǎo)數(shù)算法，但美中不足的是Pulliam和Steger[36-37]的加權(quán)平均方法中同時(shí)包括二階中心差分和二階中心插值，其在網(wǎng)格塊邊界上的應(yīng)用存在困難，且不具有普適性，不能直接推廣到高階精度有限差分方法的應(yīng)用中。緊隨其后，1979年Thomas和Lombard[24]首次提出了網(wǎng)格變換導(dǎo)數(shù)的守恒計(jì)算形式，并指出如果網(wǎng)格變換導(dǎo)數(shù)和對(duì)流通量導(dǎo)數(shù)均采用二階中心差分離散，則能滿足面積守恒律。盡管Thomas和Lombard[24]的面積守恒律解決方案是針對(duì)二階精度有限差分格式提出的，但其在解決有限差分幾何守恒律問題的歷史上仍具有里程碑式的意義，其工作也被后續(xù)諸多研究人員所引用[27-29, 38-40]。關(guān)于高階精度有限差分方法中的幾何守恒律問題，直到二十多年后的2002年才由Visbal和Gaitonde[38]將Thomas和Lombard[24]的解決方案推廣到高階精度的有限差分格式計(jì)算中。Visbal和Gaitonde[38]再次重申了采用網(wǎng)格變換導(dǎo)數(shù)守恒計(jì)算形式的必要性，且明確指出網(wǎng)格變換導(dǎo)數(shù)和對(duì)流項(xiàng)必須采用相同的高階差分格式離散，并成功將其應(yīng)用到高階緊致有限差分格式的計(jì)算中，后續(xù)研究人員將該解決思路稱之為Visbal和Gaitonde的數(shù)值方法。2010年Nonomura等[39]將Visbal和Gaitonde的數(shù)值方法[38]首先應(yīng)用到WCNS[41-43]中，模擬了三角翼的繞流問題，數(shù)值計(jì)算結(jié)果表明WCNS能夠很好地滿足面積守恒律。同時(shí)Nonomura等[39]還研究了有限差分的WENO格式[44-45]的面積守恒特性，發(fā)現(xiàn)由于在計(jì)算數(shù)值通量時(shí)采用了包含網(wǎng)格變換導(dǎo)數(shù)的混合通量分裂方法，導(dǎo)致非線性的WENO格式難以滿足面積守恒律，其對(duì)復(fù)雜網(wǎng)格的適應(yīng)能力不如WCNS[39]。為了能夠適應(yīng)繞復(fù)雜構(gòu)型流動(dòng)的數(shù)值模擬，計(jì)算格式必須滿足面積守恒律[27]。為此Nonomura[46]和Jiang[47]均對(duì)WENO格式進(jìn)行了一些改進(jìn)，以便能夠滿足面積守恒律。

2011年Deng等[40]從理論上詳細(xì)分析了面積守恒律的滿足條件，并給出了滿足幾何守恒律的充分條件(Conservative Metric Method,CMM)，同時(shí)還指出Thomas和Lombard[24]的解決方案是CMM條件的一個(gè)二階精度特例，而Visbal和Gaitonde的數(shù)值方法[38]則是CMM條件的高階精度緊致差分格式的應(yīng)用實(shí)例，并首次將幾何守恒律的滿足范圍從內(nèi)點(diǎn)嚴(yán)格地?cái)U(kuò)展到邊界和臨近邊界的計(jì)算點(diǎn)上。隨后，Deng等[48]從幾何意義的角度給出了滿足幾何守恒律的對(duì)稱守恒網(wǎng)格導(dǎo)數(shù)算法(Symmetrical Conservative Metric Method,SCMM)，解決了網(wǎng)格變換導(dǎo)數(shù)和雅克比計(jì)算形式的唯一性問題。

本文主要介紹有限差分方法中幾何守恒律問題的研究進(jìn)展，重點(diǎn)將集中在有限差分框架下的網(wǎng)格變換導(dǎo)數(shù)和雅克比的具體離散形式和幾何守恒律對(duì)數(shù)值計(jì)算結(jié)果的影響。

1 幾何守恒律問題

笛卡爾坐標(biāo)系(t,x,y,z)下的一般守恒律方程的微分形式可以表述為：

其中Q為流場(chǎng)變量，E、F和G為關(guān)于Q的通量。為了便于在貼體結(jié)構(gòu)網(wǎng)格中進(jìn)行有限差分離散，需對(duì)守恒律方程(1)進(jìn)行坐標(biāo)變換。

經(jīng)坐標(biāo)變換(2)，計(jì)算坐標(biāo)系下的守恒律方程可表述為：

其中：

以及：

網(wǎng)格變換導(dǎo)數(shù)的具體表達(dá)式：

以及網(wǎng)格變換雅克比為：

網(wǎng)格運(yùn)動(dòng)速度的具體表達(dá)式為：

為了簡化公式的表達(dá)形式，將笛卡爾坐標(biāo)系下的坐標(biāo)表示為矢量形式：

r=xex+yey+zez

(9)

其中矢量ex、ey和ez分別為笛卡爾坐標(biāo)系下的三個(gè)不變的坐標(biāo)基，x、y和z分別為矢量r在笛卡爾坐標(biāo)系下三個(gè)坐標(biāo)基方向上的分量。將網(wǎng)格變換導(dǎo)數(shù)(6)表示為矢量的分量形式：

則容易得到坐標(biāo)變換導(dǎo)數(shù)(6)更為緊湊的矢量形式為：

ξ=rη×rζ,η=rζ×rξ,ζ=rξ×rη

(11)

以及網(wǎng)格變換雅克比(7)的緊湊表達(dá)形式：

J-1=rξ·ξ=rη·η=rζ·ζ

(12)

類似的，網(wǎng)格運(yùn)動(dòng)速度(8)也可以表述為：

對(duì)于流場(chǎng)變量φ和φ，在(偏)微分算子d的作用下滿足：

d(φφ)=φdφ+φdφ

(14)

在微分算子(14)作用下，將網(wǎng)格導(dǎo)數(shù)表達(dá)式(6)和(8)代入式(5)中，容易得到：

It=Ix=Iy=Iz=0

(15)

則可將計(jì)算坐標(biāo)系下的守恒律控制方程(3)進(jìn)一步表述為：

采用有限差分方法離散守恒律控制方程(16)需滿足的一個(gè)基本條件是均勻流保持：如果在計(jì)算域內(nèi)給定初值為均勻來流，且計(jì)算域的所有邊界也都為均勻來流邊界(則流場(chǎng)中任意一點(diǎn)的變量Q、E、F和G均為常數(shù))，則數(shù)值離散守恒律控制方程(16)應(yīng)該依然得到與初始值相同的均勻流。容易看出，均勻流問題本身是笛卡爾坐標(biāo)系下守恒律控制方程(1)的解，同時(shí)也是計(jì)算坐標(biāo)系下守恒律控制方程(16)的解。對(duì)于均勻流動(dòng)問題，計(jì)算坐標(biāo)系下守恒律控制方程(16)退化為：

均勻流保持要求變量Q恒為常數(shù)，即?Q/?τ=0對(duì)于任意的流場(chǎng)變量Q恒成立。方程(17)右端在微分意義下恒為0，均勻流保持要求方程(17)在有限差分離散下亦恒成立。守恒律控制方程(16)在有限差分算子離散下可以表述為：

其中δ1為有限差分求導(dǎo)算子，下標(biāo)表示求導(dǎo)方向。對(duì)于均勻流問題，式(18)退化為：

其中,

上標(biāo)“N”表示有限差分離散。式(20)對(duì)于任意變量Q恒成立意味著：

一般情況下，差分算子δ具有以下性質(zhì)：

δ(αφ)=αδφ

δ(φ+φ)=δφ+δφ

δξ(δηφ)=δη(δξφ)

(23)

其中α為常數(shù)。值得特別指出的是：差分算子δ一般還存在如下不等式：

δ(φφ)≠φδφ+φδφ

(24)

2 滿足幾何守恒律的條件

為了滿足式(21)，需要采用在數(shù)學(xué)上與網(wǎng)格導(dǎo)數(shù)(6)等價(jià)的守恒計(jì)算形式(以上標(biāo)“S1”表示)：

將網(wǎng)格導(dǎo)數(shù)的守恒計(jì)算形式(25)離散后代入式(20)有：

考慮到差分算子的性質(zhì)(23)，在滿足下列條件：

在數(shù)學(xué)上與網(wǎng)格導(dǎo)數(shù)(6)等價(jià)的守恒計(jì)算形式并不唯一，類似于式(25)，還存在另外一種守恒計(jì)算形式(以上標(biāo)“S2”表示)：

式(29)還可以表述為更為緊湊的矢量形式：

在網(wǎng)格變換導(dǎo)數(shù)(6)的幾種能滿足幾何守恒律的守恒計(jì)算形式中，本文推薦使用其對(duì)稱守恒計(jì)算形式S3(29)和(30)，原因在于S3能夠較好反映計(jì)算網(wǎng)格的特性(離散后的網(wǎng)格變換導(dǎo)數(shù)表現(xiàn)為計(jì)算網(wǎng)格單元的矢量面積)。

如果網(wǎng)格導(dǎo)數(shù)(6)的所有方向上的δ2算子均為二階中心差分時(shí)，在O點(diǎn)的矢量面積如圖1所示。

圖1 網(wǎng)格變換導(dǎo)數(shù)(6)在二階中心差分離散下的面積示意圖

在對(duì)網(wǎng)格導(dǎo)數(shù)的守恒計(jì)算形式(25)、(28)和(29)進(jìn)行離散時(shí)，記其內(nèi)層的差分算子為δ3，當(dāng)滿足如下條件時(shí)：

結(jié)合CMM條件(27)容易得到：

δ1=δ2=δ3

(32)

類似的，如果網(wǎng)格導(dǎo)數(shù)的守恒計(jì)算形式(25)、(28)和(29)的所有方向上的δ2和δ3算子均為二階中心差分時(shí)，在O點(diǎn)的矢量面積分別如圖2和圖3所示。對(duì)于質(zhì)量不佳的網(wǎng)格分布情況，圖3所示的矢量面積能更好地反映計(jì)算網(wǎng)格的離散特性。

圖2 網(wǎng)格變換導(dǎo)數(shù)(25)和(28)在二階中心差分離散下的面積示意圖

圖3 網(wǎng)格變換導(dǎo)數(shù)(29)在二階中心差分離散下的面積示意圖

以二階中心差分算子離散網(wǎng)格變換雅克比(7)，其在O點(diǎn)的體積如圖4所示。基于網(wǎng)格變換雅克比的矢量形式(12)給出其對(duì)稱計(jì)算形式：

圖4 網(wǎng)格變換雅克比(7)在二階中心差分離散下的體積示意圖

在二階中心差分算子離散下，網(wǎng)格變換雅克比的對(duì)稱計(jì)算形式(33)在O點(diǎn)的體積由六個(gè)部分體積組成，為敘述簡潔，本文僅給出其在ξ負(fù)方向上的部分體積，如圖5所示。

網(wǎng)格變換導(dǎo)數(shù)離散后表現(xiàn)為計(jì)算網(wǎng)格單元的矢量面積，網(wǎng)格變換雅克比離散后為網(wǎng)格單元的體積，且應(yīng)為由矢量面積封閉包圍的體積。由圖4和圖5可以看出，基于網(wǎng)格變換雅克比的定義式(7)和對(duì)稱計(jì)算形式(33)離散后的體積均不能被網(wǎng)格變換導(dǎo)數(shù)離散后的矢量面積封閉和包圍，據(jù)此本文給出網(wǎng)格變換雅克比的對(duì)稱守恒計(jì)算形式：

類似的，在二階中心差分算子離散下，網(wǎng)格變換雅克比的對(duì)稱守恒計(jì)算形式(34)在O點(diǎn)的體積也由六個(gè)部分體積組成，為敘述簡潔，本文僅給出其在ξ負(fù)方向上的部分體積，如圖6所示。不同于圖4和圖5所示的體積，基于網(wǎng)格變換雅克比的對(duì)稱守恒計(jì)算形式(34)離散后的體積可由圖3所示的網(wǎng)格變換導(dǎo)數(shù)的對(duì)稱守恒計(jì)算形式(29)離散后的矢量面積完全封閉和包圍。

圖5 網(wǎng)格雅克比(33)在二階中心差分離散下的部分體積示意圖

鑒于采用了網(wǎng)格變換導(dǎo)數(shù)和網(wǎng)格變換雅克比的對(duì)稱守恒計(jì)算形式，將上述滿足幾何守恒律的解決方案稱之為網(wǎng)格變換系數(shù)的對(duì)稱守恒計(jì)算方法[48](Symmetrical Conservative Metric Method,SCMM)，其詳細(xì)敘述如下：

(1) 網(wǎng)格變換導(dǎo)數(shù)采用其對(duì)稱守恒計(jì)算形式S3(29)或(30)；

(2) 網(wǎng)格變換雅克比采用其對(duì)稱守恒計(jì)算形式(34)；

(3) 方程(18)中通量差分算子和網(wǎng)格變換雅克比對(duì)稱守恒計(jì)算形式(34)的外層差分算子δ1、網(wǎng)格變換導(dǎo)數(shù)的對(duì)稱守恒計(jì)算形式(29)中外層差分算子δ2以及其內(nèi)層差分算子δ3在同一計(jì)算坐標(biāo)方向上完全相同，即同時(shí)滿足條件(27)和(31)。

對(duì)稱守恒網(wǎng)格導(dǎo)數(shù)算法(SCMM)是有限差分方法滿足幾何守恒律的充分條件，與網(wǎng)格質(zhì)量無關(guān)。還可以進(jìn)一步證明，在SCMM條件下采用高階精度有限差分格式離散的網(wǎng)格變換導(dǎo)數(shù)和雅克比是若干個(gè)二階精度格式離散值的線性組合[49-50]。

3 滿足幾何守恒律的特征對(duì)接邊界技術(shù)

將結(jié)構(gòu)網(wǎng)格技術(shù)應(yīng)用于繞復(fù)雜構(gòu)型的流動(dòng)模擬時(shí)，單塊網(wǎng)格的適應(yīng)能力極其有限，需要采用多塊對(duì)接結(jié)構(gòu)網(wǎng)格技術(shù)。高階精度有限差分方法對(duì)數(shù)值誤差的嚴(yán)格控制決定了其對(duì)網(wǎng)格對(duì)接邊界的處理方法的要求更為苛刻，因此發(fā)展了一種能夠滿足幾何守恒律的特征對(duì)接邊界處理技術(shù)[51]。

對(duì)于靜止網(wǎng)格而言，將計(jì)算坐標(biāo)系下的守恒律控制方程(16)中空間導(dǎo)數(shù)項(xiàng)歸類寫為：

經(jīng)過系列特征分析[51]，在網(wǎng)格對(duì)接邊界兩端(下標(biāo)L和R)的離散方程可以分別表述為：

其中：

λi為網(wǎng)格對(duì)接方向上的特征值，PQV為守恒變量對(duì)特征變量的導(dǎo)數(shù)矩陣，其在計(jì)算坐標(biāo)系下的詳細(xì)表達(dá)形式可參見Kim和Lee關(guān)于廣義特征邊界條件的文章[52]。

如果方程(35)中的右端項(xiàng)還存在黏性耗散項(xiàng)，則在對(duì)接邊界處的黏性項(xiàng)取其左右值的平均：

在上述離散下可以證明，QL和QR在數(shù)學(xué)上嚴(yán)格相等。在進(jìn)行時(shí)間推進(jìn)后將對(duì)邊界上的流場(chǎng)變量Q取為左右值的算術(shù)平均以消除數(shù)值計(jì)算過程中舍入誤差的影響：

基于上述特征對(duì)接邊界處理技術(shù)和滿足幾何守恒律的WCNS-E-5格式，完成了若干復(fù)雜拓?fù)浣Y(jié)構(gòu)的對(duì)接網(wǎng)格算例，其中30P-30N三段翼的流動(dòng)模擬結(jié)果如圖7所示。

圖7 30P-30N三段翼數(shù)值模擬結(jié)果(左：網(wǎng)格；右：壓力)

由圖7所示的流場(chǎng)等值線可以看出，即使網(wǎng)格對(duì)接不太光滑，由于采用了能夠滿足幾何守恒律的特征對(duì)接邊界處理技術(shù)，數(shù)值方法在邊界附近也能保持精度不降，流場(chǎng)在對(duì)接邊界附近也比較光滑(對(duì)接面兩側(cè)網(wǎng)格的粗細(xì)程度相差不大時(shí))。

4 幾何守恒律誤差的影響研究

4.1 幾何守恒律誤差對(duì)數(shù)值解影響的理論分析

眾所周知，高階精度有限差分格式可以成功應(yīng)用于簡單幾何外形流動(dòng)的精細(xì)模擬，但在復(fù)雜幾何外形上的應(yīng)用卻遇到了很大的困難。我們認(rèn)為其根本原因是復(fù)雜幾何外形的網(wǎng)格質(zhì)量難以滿足高階精度差分格式的要求。由此出發(fā)，我們開展了非充分光滑網(wǎng)格下，幾何守恒律誤差對(duì)數(shù)值計(jì)算結(jié)果精度影響的理論分析工作[53]。

以二維情況下流場(chǎng)變量u對(duì)笛卡爾坐標(biāo)x的導(dǎo)數(shù)?u/?x為例說明：

其在有限差分格式離散下的數(shù)值解?uN/?x可表述為：

其中網(wǎng)格變換雅克比在二維情況下的對(duì)稱守恒計(jì)算形式為：

如果計(jì)算網(wǎng)格為M階連續(xù)，不妨設(shè)網(wǎng)格坐標(biāo)y對(duì)計(jì)算坐標(biāo)系ξ方向的導(dǎo)數(shù)為M階連續(xù)，即：

其余網(wǎng)格導(dǎo)數(shù)均足夠光滑。通過Taylor展開，在N階精度有限差分格式離散下，當(dāng)幾何守恒律不滿足(δ1≠δ2)時(shí)，數(shù)值解的精度為：

(44)

(45)

4.2 幾何守恒律誤差對(duì)計(jì)算精度的影響

為了驗(yàn)證幾何守恒律誤差對(duì)計(jì)算精度的影響，采用三角函數(shù)生成不同連續(xù)程度的網(wǎng)格[53-54]，三角函數(shù)的指數(shù)n代表網(wǎng)格波動(dòng)的連續(xù)程度，如圖8所示。

圖8 網(wǎng)格波動(dòng)示意圖

給定解析初場(chǎng)及其導(dǎo)數(shù)：

(46)

采用四階中心差分格式計(jì)算得到導(dǎo)數(shù)?uN/?x，給出其相對(duì)于解析導(dǎo)數(shù)(46)的數(shù)值誤差，逐步加密網(wǎng)格，其統(tǒng)計(jì)誤差L∞范數(shù)的數(shù)值精度表現(xiàn)如圖9所示。

圖9 數(shù)值誤差

通過圖9所示的不同連續(xù)程度的網(wǎng)格下的數(shù)值精度可以看出，當(dāng)計(jì)算網(wǎng)格不是足夠光滑時(shí)，幾何守恒律的滿足能夠有效提升數(shù)值解的離散精度，其數(shù)值精度相較于幾何守恒律不滿足時(shí)要高出一階精度，同理論分析結(jié)論一致[53]。

4.3 幾何守恒律誤差對(duì)計(jì)算穩(wěn)定性的影響

采用高階精度有限差分格式(HDCS)模擬對(duì)數(shù)值誤差敏感的雙圓柱散射問題，在同一時(shí)刻的脈動(dòng)壓力均方根如圖10所示。

圖10 脈動(dòng)壓力均方根云圖

從圖10可以看出，如果幾何守恒律滿足，聲場(chǎng)光滑合理；如果幾何守恒律不滿足，聲場(chǎng)會(huì)被幾何守恒律誤差所污染。σ較小時(shí)(σ≤1×10-3)，污染主要集中在網(wǎng)格質(zhì)量不佳的奇點(diǎn)附近；σ越大，聲場(chǎng)被污染得越嚴(yán)重，σ>1×10-2時(shí)，整個(gè)聲場(chǎng)都被污染，σ=1時(shí)聲場(chǎng)已被嚴(yán)重污染而導(dǎo)致計(jì)算過程的不穩(wěn)定。

通過對(duì)雙圓柱散射問題的數(shù)值模擬，驗(yàn)證了幾何守恒律對(duì)高階精度有限差分格式計(jì)算的重要影響。當(dāng)網(wǎng)格質(zhì)量不佳時(shí)，由于幾何守恒律不滿足而引起的網(wǎng)格離散誤差可能對(duì)流場(chǎng)產(chǎn)生比較明顯的污染，嚴(yán)重時(shí)還可能引起計(jì)算過程的不穩(wěn)定和計(jì)算結(jié)果的發(fā)散。幾何守恒律的滿足能夠大幅提升高階精度有限差分格式在復(fù)雜網(wǎng)格計(jì)算中的適應(yīng)能力，計(jì)算結(jié)果也更為可信。

5 結(jié) 論

以實(shí)現(xiàn)高階精度有限差分?jǐn)?shù)值方法在繞復(fù)雜構(gòu)型流動(dòng)中的應(yīng)用為目標(biāo)，綜合運(yùn)用理論分析和數(shù)值驗(yàn)證手段，比較系統(tǒng)地研究了結(jié)構(gòu)網(wǎng)格中高階精度有限差分方法的幾何守恒律問題，得到以下認(rèn)識(shí)：

1) 滿足幾何守恒律的充分條件表明，不是所有的有限差分算子都能較為方便地滿足幾何守恒律，而鄧小剛等提出的系列WCNS格式對(duì)于滿足幾何守恒律具有明顯優(yōu)勢(shì)。

2) SCMM條件是滿足幾何守恒律的充分條件，同時(shí)網(wǎng)格變換導(dǎo)數(shù)和雅克比均采用其對(duì)稱守恒計(jì)算形式，具有唯一性?；赟CMM條件離散的矢量面積(由網(wǎng)格變換導(dǎo)數(shù)分量組成)完全封閉，網(wǎng)格變換雅克比離散后體積恰好由矢量面積所圍成，更符合計(jì)算網(wǎng)格的幾何意義。

3) 自由流保持只是滿足幾何守恒律的一種表現(xiàn)形式，是幾何守恒律滿足的必要條件，而非充分條件。嚴(yán)格來講，幾何守恒律的本質(zhì)含義應(yīng)當(dāng)包括矢量面積完全封閉、體積由矢量面積圍成和確定以及在網(wǎng)格運(yùn)動(dòng)過程中依然維持上述關(guān)系等。在靜止網(wǎng)格中SCMM條件能夠完全滿足幾何守恒律。

4) 幾何守恒律的滿足能夠有效提高高階精度有限差分格式在質(zhì)量欠佳網(wǎng)格中的計(jì)算精度，同時(shí)還能大幅提升高階精度格式在復(fù)雜網(wǎng)格中的計(jì)算穩(wěn)定性。

在下一步工作中，我們將繼續(xù)完善幾何守恒律的相關(guān)理論，同時(shí)重點(diǎn)開展幾何守恒律研究成果的應(yīng)用工作，如指導(dǎo)高階精度有限差分格式的構(gòu)造[55,57]、結(jié)構(gòu)網(wǎng)格質(zhì)量的檢測(cè)等，以繼續(xù)提升高階精度有限差分算法在繞復(fù)雜構(gòu)型流動(dòng)數(shù)值模擬的應(yīng)用能力。