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

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

0 引 言

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

高階精度應用場景的現(xiàn)實強烈需求催生了一系列的高階精度計算方法。精度最高的譜方法[7]具有理想的指數(shù)階精度，卻難以應用于復雜外形的流動模擬中，目前譜方法多用于對其他高階精度數(shù)值計算方法的考核與驗證。其他的數(shù)值方法按照其離散框架大致可以分為基于有限差分方法的、基于有限體積方法的和基于有限元方法的高階精度計算方法，以及近年來逐漸衍生出的DG方法[8-9]、譜差分方法[10-11]以及譜體積方法[12-13]。由于離散框架的不同，不同類型的高階精度數(shù)值方法提高計算精度的方式迥異，尤其是在多維問題中。得益于逐維離散特性，有限差分方法能以較小計算代價實現(xiàn)多維問題的高階精度計算[14-15]，但適應復雜幾何外形的能力較弱[16]。當網(wǎng)格不夠光滑時，高階精度有限差分方法還存在明顯的降階問題。因此，在相當長的一段時間內，基于有限差分方法的高階精度算法多用于繞簡單外形的類似于湍流等多尺度復雜流動問題的模擬中[17-18]，鮮有在復雜外形流動中成功應用的例子。近年來的研究結果表明[19-22]，基于有限差分方法的高階精度計算格式對復雜網(wǎ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)格運動過程中網(wǎng)格變換雅克比和網(wǎng)格運動速度之間需滿足的約束關系，面積守恒律則反映了網(wǎng)格變換導數(shù)各分量之間需滿足的約束關系。

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

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

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

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

1 幾何守恒律問題

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

其中Q為流場變量，E、F和G為關于Q的通量。為了便于在貼體結構網(wǎng)格中進行有限差分離散，需對守恒律方程(1)進行坐標變換。

經坐標變換(2)，計算坐標系下的守恒律方程可表述為：

其中：

以及：

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

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

網(wǎng)格運動速度的具體表達式為：

為了簡化公式的表達形式，將笛卡爾坐標系下的坐標表示為矢量形式：

r=xex+yey+zez

(9)

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

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

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

(11)

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

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

(12)

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

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

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

(14)

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

It=Ix=Iy=Iz=0

(15)

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

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

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

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

其中,

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

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

δ(αφ)=αδφ

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

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

(23)

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

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

(24)

2 滿足幾何守恒律的條件

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

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

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

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

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

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

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

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

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

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

δ1=δ2=δ3

(32)

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

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

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

以二階中心差分算子離散網(wǎng)格變換雅克比(7)，其在O點的體積如圖4所示?；诰W(wǎng)格變換雅克比的矢量形式(12)給出其對稱計算形式：

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

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

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

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

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

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

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

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

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

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

3 滿足幾何守恒律的特征對接邊界技術

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

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

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

其中：

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

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

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

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

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

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

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

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

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

以二維情況下流場變量u對笛卡爾坐標x的導數(shù)?u/?x為例說明：

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

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

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

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

(44)

(45)

4.2 幾何守恒律誤差對計算精度的影響

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

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

給定解析初場及其導數(shù)：

(46)

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

圖9 數(shù)值誤差

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

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

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

圖10 脈動壓力均方根云圖

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

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

5 結 論

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

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

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

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

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

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