劉 君, 韓 芳, 夏 冰

(大連理工大學(xué) 航空航天學(xué)院, 遼寧 大連 116024)

0 引 言

美國在《2030年CFD遠景規(guī)劃》[1]中認為，自適應(yīng)網(wǎng)格技術(shù)和誤差估計是復(fù)雜流動和復(fù)雜幾何外形數(shù)值模擬的重大瓶頸技術(shù)。網(wǎng)格是計算流體力學(xué)(CFD)解決復(fù)雜外形工程應(yīng)用的關(guān)鍵，是該學(xué)科的重要內(nèi)容。在網(wǎng)格理論研究方面，除了傳統(tǒng)的生成算法外，近些年幾何守恒律(Geometric Conservation Law, GCL)問題也受到關(guān)注。它直接影響計算結(jié)果精度，是解決“誤差估計”的基礎(chǔ)。

1 幾何守恒律的研究現(xiàn)狀

幾何守恒律的概念最早是Thomas和Lombard于1978年在AIAA會議上提出[2]，隨后在AIAA期刊正式發(fā)表[3]。但對幾何守恒律問題的研究，最早可見于1961年Trulio和Trigger的工作報告中對一維坐標下“體積守恒(Conservation of Volume)”問題的討論[4]。

對于笛卡爾坐標系(t,x,y,z)的三維守恒型Euler方程：

(1)

應(yīng)用于復(fù)雜外形時，需要進行坐標變換，在計算空間中(τ,ξ,η,ζ)形式如下：

(2)

其中：

對于流動參數(shù)為常數(shù)的均勻流場，以上方程變成：

(3)

這是Thomas針對有限差分法給出的微分形式的GCL方程，但未對此進行深入討論，因為他們主要關(guān)注有限體積法在網(wǎng)格變形時涉及到的GCL問題。對于有限體積法，基于ALE(Arbitrary Lagrange-Euler)描述下的三維Euler方程：

(4)

其中，Ω為控制體，?Ω為控制邊界，n為控制體邊界外法向向量，dσ為體積微元，ds為面積微元，xc為網(wǎng)格面的運動速度。把均勻流場參數(shù)代入，得到積分形式的GCL方程：

(5)

根據(jù)式(5)，對于靜止網(wǎng)格xc=0，有限體積法不存在GCL問題，即使是運動網(wǎng)格，如果沒有變形，式(5)也自動滿足。國內(nèi)采用結(jié)構(gòu)網(wǎng)格與物體剛性固連運動開展動導(dǎo)數(shù)研究的文獻[5-8]和相關(guān)學(xué)術(shù)專著[9]采納了以上結(jié)論。但是，從本文推導(dǎo)來看，這種結(jié)論僅對梯度重構(gòu)界面通量的有限體積法成立，對于采用有限差分法重構(gòu)界面通量的高精度格式不一定正確。本文后面的例子表明，所有包含有四邊形的網(wǎng)格均存在幾何參數(shù)計算不相容引入誤差的風險。

下面考察有限差分法產(chǎn)生GCL問題的數(shù)學(xué)本質(zhì)。推導(dǎo)物理空間(t,x,y,z)的守恒控制方程變換到計算空間(τ,ξ,η,ζ)得到如下原始形式：

=UIt+FIx+GIy+HIz=0

(6)

只是在時間和空間導(dǎo)數(shù)連續(xù)條件下，使用了式(3)及以下恒等式：

(7)

從理論上講，帶有源項的方程(6)才是原始方程(1)的等價形式，但坐標變換系數(shù)大多采用網(wǎng)格位置的差分來近似求解，此時恒等式(3)和式(7)可能不成立，在這種情況下方程(2)與方程(1)不等價，把它作為曲線坐標系下的等價形式進行CFD計算，得到的數(shù)值解必然存在誤差。研究GCL算法的目的是消除由于幾何參數(shù)計算引起的誤差。

在GCL研究領(lǐng)域，把涉及網(wǎng)格變形和時間導(dǎo)數(shù)的式(3)稱為體積守恒律，把僅有空間導(dǎo)數(shù)的式(7)稱為面積守恒律[10]。動網(wǎng)格和靜網(wǎng)格均須考慮面積守恒律，這點已受到廣泛關(guān)注，亦是目前GCL研究的主要內(nèi)容。下面對相關(guān)研究文獻進行綜述。

有限差分法GCL研究可以追溯到CFD應(yīng)用于復(fù)雜外形流場模擬的初期。1977年Steger[11]采用二階中心格式來計算坐標變換系數(shù)，發(fā)現(xiàn)即使對均勻來流也存在質(zhì)量不守恒的現(xiàn)象，后來Hindman[12-13]等人進一步給出了這種誤差導(dǎo)致非物理解、引起數(shù)值振蕩的算例。此后，1978年Steger和Pulliam[14]采用一種加權(quán)平均算法計算坐標變換系數(shù)，解決了二階中心格式離散流動方程時遇到的GCL問題。隨著高精度格式進入實用階段，人們發(fā)現(xiàn)Steger的算法難以推廣，GCL問題重新引起重視。

2002年Visbal和Gaitonde[15]采用與流動方程相同的高階離散格式計算坐標變換系數(shù)，成功解決了高階緊致格式(CS)的GCL問題。2010年Nonomura等人把Visbal[16]的算法推廣到WCNS格式，但是在應(yīng)用于WENO格式時遇到難題，最后通過修改WENO格式來迎合GCL[17-18]。2011年，為解決緊致類格式的GCL問題，鄧小剛分析了構(gòu)建GCL算法的充分條件，提出坐標變換系數(shù)的守恒算法(Conservative Metric Method: CMM)[19]，并證明Steger等人的算法是CMM的二階精度特例，Visbal等人的算法是 CMM的高階精度應(yīng)用例證，并于2013年從計算幾何學(xué)角度提出坐標變換系數(shù)計算結(jié)果唯一的對稱守恒算法(Symmetrical Conservative Metric Method, SCMM)[20]。下面對這一算法簡單介紹。

(8)

理論上導(dǎo)數(shù)乘積滿足交換律：yηzζ=zζyη，差商離散以后交換律不再成立。為解決這個離散過程產(chǎn)生出來的問題，Steger等人和Thomas等人引入算子δ3構(gòu)造對稱算法，例如導(dǎo)數(shù)yηzζ的差商采用

(9)

(10)

(11)

同時還有對J-1的計算要求，這里不再詳述。由于緊致格式可寫成若干中心格式的組合，因此，按照上面離散坐標變換系數(shù)的高階格式可以表示為幾個矢量面積的線性組合，很容易將SCMM推廣到高階精度緊致格式。但是對于其它不滿足以上算子條件的高精度格式，很難應(yīng)用。

關(guān)于有限差分下的體積守恒律問題，Abe[21-22]等人借鑒Zhang[28]等人在有限體積框架下對GCL問題所做的工作，通過改變與網(wǎng)格運動相關(guān)的坐標變換系數(shù)的計算方式，使體積守恒律得到滿足。

綜上所述，有限差分法的GCL研究現(xiàn)狀是，鄧小剛等人成功地解決了面積守恒律問題[18-19]，Abe等人成功解決了體積守恒律問題，但以上工作都是針對緊致格式等中心格式的，對其它格式缺乏普適性。

2 構(gòu)造普適性GCL算法的難點

前面分析了GCL問題的數(shù)學(xué)本質(zhì)，下面通過具體實例討論在離散過程中幾何參數(shù)的誤差是如何產(chǎn)生的。以二維問題為例，采用中心差分格式計算圖1中(i,j)點坐標變換產(chǎn)生的導(dǎo)數(shù)，其中Jocobi行列式為：

(12)

其余系數(shù)ξx、ξy、ηx和ηy也可以采用類似的差分格式得到。在計算空間內(nèi)Δξ=Δη=1，內(nèi)層差分采用二階中心差分，外層考慮最簡單的一階迎風格式，對守恒律方程(6)進行數(shù)值離散得到：

(13)

在任意曲線條件下難以保證Ix=0和Iy=0，如果這些參數(shù)隨空間變化產(chǎn)生的誤差不同，即使是均勻流場也無法維持均勻和穩(wěn)定。

(14)

對于不變形網(wǎng)格，在二維情況下三角形網(wǎng)格、三維情況下四面體網(wǎng)格，采用單元梯度來構(gòu)造二階或者更高精度的，計算過程中不涉及相鄰網(wǎng)格的幾何參數(shù)，有限體積法的幾何參數(shù)計算表達式唯一且相容，不存在GCL問題。但是對于結(jié)構(gòu)網(wǎng)格基礎(chǔ)上的有限體積法，通過式(13)可以看出，如果沿用有限差分法，采用擴展模板的方法來構(gòu)造高精度格式，同樣存在不滿足GCL的風險，因此，前述文獻關(guān)于“不變形結(jié)構(gòu)網(wǎng)格不存在GCL問題”的結(jié)論值得進一步討論。

近兩年國內(nèi)學(xué)者在WENO格式的GCL研究方面取得進展[29-30]，給出的均勻流場算例結(jié)果誤差曲線不能保持直線，計算初期變化明顯，后期增長變緩。從研究思路看，和SCMM類似，采用調(diào)整坐標變換系數(shù)的算法來滿足GCL，同樣面臨把系數(shù)的解析值代入不一定滿足的問題。

3 基于等價方程的GCL補償算法

為了便于理解本文針對有限差分法建立基于修正方程的GCL補償算法，下面簡單介紹非結(jié)構(gòu)變形網(wǎng)格的GCL算法分類和建立GCL算法的思路。

劉君等人對非結(jié)構(gòu)變形網(wǎng)格的GCL進行了研究[23]，提出如下觀點：GCL方程(5)是流體力學(xué)在均勻流場條件下的退化方程，不是理論上必須滿足的新約束條件，而不滿足GCL產(chǎn)生非物理解的機理是計算過程中體積增量與面積運動形成的體積不等。采用圖3所示網(wǎng)格單元來說明。

通過以上簡單模型認識到計算過程中體積增量與面積運動形成體積不等是產(chǎn)生非物理解的機理，以此來考察國內(nèi)外文獻構(gòu)造的多種GCL算法，本質(zhì)都是修改右端項使得面積運動形成的體積等于體積增量ΔV。

2015年，Chang等人[24]對現(xiàn)有的有限體積的幾何守恒律算法進行了分類，按照誤差消除方式的不同，分為帶源項的“體積約束方法(Volume Constrained Method，VCM)”和不帶源項的“面積約束方法(Face Constrained Method，F(xiàn)CM)”。2016年，根據(jù)選擇的修正參數(shù)的不同，劉君[25]把以上算法進一步細分為四類：

(1) 1978年Thomas等人[1]提出的修正面積的GCL算法，引入平均面積：

(15)

這類算法對網(wǎng)格運動速度沒有限制，但僅適用于兩個時間層的格式。

(2) 1995年Farhat等人[26]構(gòu)造時間二階精度隱格式(BDF2)時遇到tn-1和tn+1時刻之間多套網(wǎng)格的 GCL問題，提出修正速度的GCL算法。為了滿足幾何相容，計算過程中每個時間層根據(jù)面積計算速度，推導(dǎo)過程復(fù)雜，涉及20個參數(shù)，增加了計算量。后來2006年Mavriplis等人[27]研究了Farhat算法，發(fā)現(xiàn)存在參數(shù)不唯一的問題，針對BDF2格式，有些情況下并非全都嚴格滿足幾何相容條件，他們提出如下網(wǎng)格運動速度算法：

(16)

(3) 2004年Zhang等人[28]采用時間兩層隱式格式(Crank-Nicolson格式)構(gòu)造隱式格式，提出同時修正面積和網(wǎng)格速度的算法：

(17)

(4) 2009年劉君[23]基于以上機理認識，提出修正體積的第四類GCL算法。例如對BDF2格式，采用如下體積代替tn+1時刻方程左端項的幾何體積Vn+1：

(18)

也可以選擇修正其它時間層的體積。

通過以上分類比較，可以看出前三類GCL算法(FCM方法)的思路是修改離散方程右端項中的幾何參數(shù)來滿足方程左端項的體積增量，劉君提出的第四類GCL算法(VCM方法)直接根據(jù)方程兩側(cè)的幾何參數(shù)不相容量來消除誤差。按照以上思路來分析有限差分法：

(1) 首先，微分形式的GCL方程同樣是流體力學(xué)控制方程的退化方程，不是理論上必須滿足的新約束條件，離散過程后幾何參數(shù)不相容是GCL的機理；

(2) SCMM算法思路是通過構(gòu)造坐標變換系數(shù)的離散格式滿足相容關(guān)系式：Ix=Iy=Iz=0，以此來消除曲線坐標系下方程(6)的右端項，類似于有限體積法前三類GCL算法；

守恒控制方程從物理空間(t,x,y)坐標系變換到計算空間(τ,ξ,η)坐標系的完整形式是：

(19)

(20)

如果不能保證以上GCL恒等式成立，離散以后和原始方程(1)不等價，稱為離散近似方程。

(21)

由于坐標變換系數(shù)和流動參數(shù)之間不相關(guān)，因此，在連續(xù)空間求導(dǎo)數(shù)，得到如下關(guān)系式：

(22)

實際上從式(1)到式(2)推導(dǎo)過程中也用到式(22)：

如果采用相同的算子離散，忽略高階小量后，式(22)依然成立：

(23)

下面按照這個左右兩側(cè)格式相同的原則來推導(dǎo)一階迎風格式的通量分裂：

(24)

按照式(23)展開計算方程(24)中的分裂格式，以其中一項為例：

整理后寫為：

(25)

這樣一來，補償源項以后的修正方程的一階迎風格式可以變成如下簡潔形式：

(26)

上式與原始方程的一階迎風格式比較，只是把所有離散點的坐標變換系數(shù)換成(i,j)點的參數(shù)。從前面對基于梯度重構(gòu)的有限體積法在靜止網(wǎng)格上沒有GCL誤差的機理分析，本文提出的基于修正方程的GCL補償算法，盡管初始離散中包含有相鄰點的坐標變換系數(shù)，但是最終形式和計算過程中不受相鄰點的坐標變換系數(shù)影響，二者的思路是一致的。

(27)

由于對一階顯式時間推進

(28)

式(27)可進一步簡化為：

(29)

以上是針對一階格式進行推導(dǎo)的，隨后，我們針對空間二階迎風格式及Runge-Kutta等高階時間精度格式也進行了推導(dǎo)，公式同樣成立。但尚未對隱式時間格式進行推導(dǎo)證明。

4 驗證算例

目前大家對GCL問題的研究大多基于高精度格式，實際上最早的GCL是為了消除均勻流場在曲線坐標系下計算時出現(xiàn)的因網(wǎng)格計算引起的誤差問題。因此，對GCL問題，有效的驗證算例應(yīng)該是均勻流場問題。對于均勻流場問題，流場參數(shù)導(dǎo)數(shù)為0，因此高精度格式和一階格式的表現(xiàn)是一樣的，如果一階格式出現(xiàn)GCL誤差，高精度格式只能減小誤差，不能消除，如果一階格式?jīng)]有出現(xiàn)GCL誤差，則證明網(wǎng)格是均勻網(wǎng)格，高精度格式也應(yīng)該沒有誤差。綜上，本文選擇采用一階格式計算均勻流場問題作為驗證算例。

上下邊界條件為固壁，左右邊界條件是基于Riemann問題求解的進、出口邊界。

控制方程分別采用離散近似方程式(20)與離散等價方程式(19)，并分別記為方法一與方法二。兩種計算方法皆分別采用包含F(xiàn)DS及FVS在內(nèi)的AUSM、HLLC、Roe、VanLeer四種格式計算，這四種格式都是CFD常用的基礎(chǔ)格式。

(30)

Δξ=1。時間離散采用一階顯式推進，CFL=0.9，輸出時間t=1.5。計算中皆采用無量綱的物理量。

對不同馬赫數(shù)Ma，計算非均勻網(wǎng)格圖4上的均勻流場問題，方法一計算結(jié)果密度誤差云圖如圖5至圖8所示，而方法二計算誤差全為0，故在此不再顯示。

由圖5至圖8可以看出，在馬赫數(shù)為0時，雖然四種格式計算的誤差各不相同，但都集中在網(wǎng)格中心圓周圍，因此處網(wǎng)格尺度變化最大，且正交特性最差。隨著馬赫數(shù)的增加誤差向下游傳播，直到馬赫數(shù)Ma>1，誤差沿著馬赫錐方向向下游傳播。

四種格式中，HLLC格式與Roe格式密度誤差云圖分布最為相似，這是因為這兩種格式都是基于近似Riemann解的思路提出的。但在進入超聲速后，兩者的區(qū)別越來越明顯，這是因為對Riemann問題，Roe格式認為間斷的左右狀態(tài)之間形成的是激波，通過構(gòu)造一個中間狀態(tài)來模擬計算；HLLC格式則認為間斷的左右狀態(tài)之間存在接觸間斷，通過構(gòu)造兩個中間狀態(tài)來進行計算，更能夠精確地模擬接觸間斷。

對運動網(wǎng)格，設(shè)計網(wǎng)格運動情況如圖9所示，網(wǎng)格以圖4為基準，如圖9所示進行每一時間步為5°的左右旋轉(zhuǎn)周期運動，正向最大旋轉(zhuǎn)角度為順時針旋轉(zhuǎn)30°，如圖9(a)所示，反向最大旋轉(zhuǎn)角度為逆時針旋轉(zhuǎn)30°，如圖9(b)所示。

同樣針對不同馬赫數(shù)、四種格式下方法一及方法二計算結(jié)果密度誤差絕對值的變化曲線如圖10所示。由于方法二計算密度誤差在四種格式下全為0，所以將四條曲線合為一條，以Method_2表示。

(a)Ma=0

(b)Ma=0.9

(c)Ma=1.1

(d)Ma=2.0

圖10動網(wǎng)格下密度誤差曲線
Fig.10Densityerrorsvarieswithtime

圖10所示的誤差曲線，每一時間點的誤差值是圖4及圖9所示區(qū)域所有點的誤差絕對值的平均值。由圖10可以看出，方法一在四種格式下計算動網(wǎng)格上的均勻流，均存在誤差，其中仍然是HLLC格式與Roe格式的計算結(jié)果最為接近。然后隨著馬赫數(shù)的逐漸增大，四種格式對誤差的計算結(jié)果也逐漸靠近。但對于方法二，不論采用何種格式，誤差總為0，因此可以用Method_2一條曲線表示。

5 結(jié) 論

通過理論推導(dǎo)和數(shù)值研究，得到如下結(jié)論：

(1)由于離散條件下坐標變換恒等式不再成立，導(dǎo)致目前CFD領(lǐng)域廣泛使用的曲線貼體坐標系下基本方程與直角坐標系下原始方程可能不等價。等價方程應(yīng)該帶有源項，由此提出離散等價方程概念。通過分析源項的產(chǎn)生機理，提出離散等價方程左右兩側(cè)格式需要滿足的準則。

(2)在準則指導(dǎo)下，根據(jù)坐標變換系數(shù)和流動參數(shù)之間呈線性關(guān)系的特點，提出新的耦合算法，并且采用算子理論進行證明。這種從經(jīng)典CFD角度看屬于凍結(jié)系數(shù)降低精度的算法，對于離散等價方程是保持格式精度并滿足幾何守恒律的。

(3)一階迎風格式數(shù)值算例表明，新的幾何守恒律算法對坐標變換系數(shù)沒有要求，可以直接使用準確的解析解，適合于FDS和FVS通量分裂格式。原則上算子理論推導(dǎo)出來的結(jié)論與具體格式無關(guān)，但用于高精度格式的效果還需要驗證。