張懷寶，王光學，王靖宇

(中山大學 物理學院， 廣東 廣州 510006)

傳統(tǒng)的可壓縮流體求解器通常采用以守恒變量為自變量的守恒形式的控制方程，其中為代表的國內(nèi)外比較著名的求解器有HOSTA[1]，TRIP[2]，OpenCFD[3]，F(xiàn)UN3D[4]，OpenFOAM[5]等。采用有限體積方法，積分形式下的控制方程一般可以表示為：

(1)

式中，Q為守恒變量，F(xiàn)為無黏通量，F(xiàn)d為黏性通量，S為源項，t為時間，n為單元控制面的法向矢量。通過鏈式法則，在時間項上引入雅克比矩陣J=?Q/?P，可以將上述守恒變量守恒形式的控制方程轉(zhuǎn)換為以原始變量P為自變量的守恒形式的控制方程[6-7]，即

(2)

對于一維Euler方程，上述各個變量的定義為：

(3)

在某些具體的數(shù)值運算中，都需要用到原始變量，例如半點處的變量重構多采用原始變量[8]；另外，對黏性流體的求解，動力黏度和傳熱系數(shù)都是溫度的函數(shù)[9]，而黏性應力的計算需要采用速度梯度；邊界條件，例如無滑移邊界，一般假設壁面壓強的法向梯度為零[10]。而采用守恒變量形式的控制方程時，需要在每一個時間步推進后，將更新的守恒變量轉(zhuǎn)換成原始變量，進行相關計算。特別地，當涉及熱化學非平衡問題時，溫度是總能的隱函數(shù)，直接求取溫度比較困難，實際一般采用迭代方法，時間耗費較大，而采用原始變量形式，無須該轉(zhuǎn)換過程[11]。另外，采用原始變量形式可以大大簡化時間隱格式中雅克比矩陣的構造，推導解析形式的雅克比矩陣相對容易，而實現(xiàn)數(shù)值形式的雅克比矩陣更為直接[11]。

本文采用的原始變量類型的控制方程基于壓力、速度和溫度，一方面該組變量的選取與Weiss等[6]在構造低速預處理矩陣時的選取保持一致，因此現(xiàn)有的計算平臺可以直接采用其發(fā)展的預處理方法，并將計算能力擴展到低馬赫數(shù)問題；另一方面，溫度的計算相對于其他變量，尤其針對熱化學問題，較為重要也較為困難[12-14]；最后采用速度作為原始變量，對應動量方程，較為直觀。在對低速預處理系統(tǒng)的相關研究中，Turkel[15]指出，預處理后的控制方程盡管對于定常問題滿足守恒性質(zhì)，但是對于非定常問題則不成立，原因是預處理矩陣的引入破壞了時間準確性。然而從數(shù)值分析發(fā)現(xiàn)，即使不對引入的雅克比矩陣進行低速預處理，基于原始變量的控制方程在離散之后，對于非定常問題也不能得到準確的數(shù)值解。數(shù)值試驗也發(fā)現(xiàn)，如果流場光滑，非定常的數(shù)值結果誤差或許并不明顯，但是在梯度較大或者存在強間斷時，流場解與準確解相比會出現(xiàn)較大差異。為了消除該數(shù)值方法帶來的誤差，獲得準確的非定常解，可以構造相應的雙時間步方法，虛擬時間項采用原始變量，而物理時間項采用守恒變量，在相鄰的兩個物理時間步內(nèi)作為定常問題求解。盡管虛擬時間推進過程中存在誤差，但是方程最終會收斂到物理時間上的守恒形式。計算結果表明，該雙時間步方法與守恒變量守恒形式下的雙時間步方法計算結果相同，兩種方法等價。

1 數(shù)值方法

1.1 守恒變量守恒形式控制方程的離散

積分公式(1)可以表示為：

(4)

式中：

R=-A(F-Fd)·ndA+?VSdV

(5)

R為右端項(或者殘差項)。 對某一個控制單元i，右端項在空間上離散，可以表示為：

(6)

式中，Aj為單元i的第j個控制面的面積，nj為該面的法向矢量，N為單元i的控制面總個數(shù)，Vcell為單元i的體積。 不考慮網(wǎng)格的運動，時間項的離散采用一階向后差分，即

(7)

(8)

第n+1時間步的流場解為：

(9)

1.2 原始變量守恒形式控制方程的離散

類似地，積分公式(2)可以表示為：

(10)

式中：

R=-A(F-Fd)·ndA+?VSdV

(11)

即右端項與式(5)相同，該項的空間離散可以采用式(6)表示。 不考慮網(wǎng)格的運動，時間項的離散采用一階向后差分，即

(12)

(13)

第n+1時間步的流場解為：

(14)

1.3 非定常解的時間準確性

不難看出，盡管J=?Q/?P，但是Q≠JP，而是Q=JP+M。其中：

(15)

因此原始變量方法得到的守恒變量值為：

(16)

(17)

(18)

1.4 基于原始變量的雙時間步方程

(19)

對式(19)進行一階顯式離散，得：

(20)

(21)

(22)

應用近似關系式

(23)

雙時間步方程的線性方程組最終表示為：

(24)

2 激波管問題與結果分析

非定常流動算例采用典型的一維激波管問題，其初值條件見表1。驅(qū)動段與非驅(qū)動段的分隔位置位于x=0，氣體常數(shù)R=287.06，比熱比γ=1.4，各個物理量均進行無量綱化處理。

表1 激波管問題初值條件

計算采用中山大學計算流體力學研究中心開發(fā)的非結構網(wǎng)格計算平臺MulPhy。針對激波管問題，該求解器只求解Euler方程，采用Median-Dual格點格式離散控制體，通量格式采用AUSM+-up[20]，控制方程可以采用原始變量或者守恒變量作為自變量，時間推進采用一階顯式格式，線性求解器使用GMRES[21]方法。計算網(wǎng)格為二維，網(wǎng)格數(shù)量為600×20，分別對應x與y方向，且數(shù)值實驗表明，該網(wǎng)格尺度下空間上的數(shù)值誤差與時間上的數(shù)值誤差相比，可以忽略不計。時間推進的步長為Δt=1×10-7，總共推進1000個時間步，對應物理時間t=1×10-4。雙時間步計算中，每兩物理時間步之間的收斂殘差降到初始殘差的10個量級以下。圖1(a)～(c)分別是壓力、速度、溫度的計算結果。在每張圖中，四組計算結果分別對應解析解，守恒變量方法(式(9))數(shù)值解，原始變量方法(式(14))數(shù)值解，原始變量+雙時間步方法(式(24))數(shù)值解。從結果對比中可以看出，采用原始變量方法與守恒變量方法得到的壓力和速度曲線略有差別，但差別不大，而溫度曲線差別較為明顯。特別地，在激波強間斷附近，原始變量方法會產(chǎn)生較大誤差。從圖1(c)可以看出，激波后的溫度出現(xiàn)明顯的上沖，而采用守恒變量方法得到的數(shù)值解與解析解保持一致。圖2為溫度對比放大圖，由圖2可以看出，激波的預測位置，也出現(xiàn)了偏差。進一步研究表明，隨著時間的繼續(xù)推進，激波位置的誤差會繼續(xù)增大。

(a) 壓力曲線 (a) Pressure curve

(b) 速度曲線 (b) Velocity curve

(c) 溫度曲線 (c) Temperature curve圖1 三種數(shù)值解與解析解的原始變量對比Fig.1 Comparison of the primitive variables using three numerical approaches against analytical solution

(a) 激波與接觸間斷附近 (a) In the vicinity of the shock and contact discontinuity

(b) 激波附近 (b) In the vicinity of the shock圖2 三種數(shù)值解與解析解的溫度對比放大圖Fig.2 Contrast magnification of the temperatures using three numerical approaches against analytical solution

將守恒變量方法的數(shù)值解作為準確解，計算原始變量方法數(shù)值解的誤差范數(shù)，結果見表2。其中相對誤差范數(shù)的參考量采用驅(qū)動段的初值條件，由于初速度為零，涉及速度的參考量均采用驅(qū)動段的初始聲速a=347.3。可以看出，采用原始變量方法，各個原始變量(p,u,T)與各個守恒變量(ρ,ρu,ρE)均有不同程度的誤差。

表2 激波管問題數(shù)值結果的誤差范數(shù)統(tǒng)計

(25)

的解。

因此將守恒變量方程(1)引入雙時間步方法并增加一組數(shù)值解，對兩組采用雙時間步方法的數(shù)值解進行比較，溫度曲線如圖3所示，可以看出兩曲線完全重合；再次對各個流場變量的誤差范數(shù)進行統(tǒng)計，L2范數(shù)都在1.0×10-10以下，表明兩種方法等價。

圖3 基于原始變量的雙時間步方法與基于守恒變量 的雙時間步方法的溫度計算結果對比Fig.3 Comparison of the temperatures predicted by primitive-variable based and conservative-variable based approaches both together with dual-time stepping technique

3 雙馬赫反射問題與結果分析

為進一步驗證基于原始變量的雙時間步方法能夠恢復準確的數(shù)值解，對Woodward等[22]首先提出的雙馬赫反射問題進行計算。雙馬赫反射問題的初始條件和邊界條件如圖4所示，初始時刻，Ma=10的激波與平板呈60°，并向平板方向運動，激波與平板壁面相交產(chǎn)生激波反射過程。運動激波前后的氣體狀態(tài)分別為：

其中：u為x方向速度，v為y方向速度。

二維計算域為[0,4]×[0,1]，采用均勻網(wǎng)格，網(wǎng)格量為1921×480，分別對應x與y方向，網(wǎng)格收斂性試驗表明，該網(wǎng)格尺度滿足空間精度要求。時間推進的物理時間步長為Δt=4×10-5，總共推進5000個時間步，對應物理時間t=0.2。雙時間步計算中采用定虛擬時間步Δτ=4×10-5，并發(fā)現(xiàn)每兩物理時間步之間迭代20個虛擬時間步后結果即已經(jīng)收斂，實際計算采用30個虛擬時間步。通量格式采用AUSMPW+[23]。

圖4 雙馬赫反射問題的初始條件和邊界條件Fig.4 Initial and boundary conditions of double Mach reflection

數(shù)值計算首先采用基于原始變量的雙時間步方法(式(24))，得到t=0.2時流場的密度云圖，如圖5所示(圖中30條等勢線表示的密度范圍為2～22)?？梢钥闯?，隨著時間推進，流場演化出豐富的流動現(xiàn)象，數(shù)值方法能夠清晰地捕捉到馬赫反射后形成的三叉點、入射激波、馬赫桿、反射激波和接觸間斷以及近壁面結構。

圖5 雙馬赫反射問題t=0.2時密度云圖Fig.5 Density contours at t=0.2 for the double Mach reflection

作為對比，采用基于守恒變量方程(1)的雙時間步方法重新計算，得到另一組數(shù)值解，然后對兩組數(shù)值解的密度變量進行誤差范數(shù)統(tǒng)計，L2范數(shù)值為8.14×10-12，表明盡管該問題物理現(xiàn)象非常復雜，兩種方法的計算結果仍然相等，證明兩種雙時間步方法等價。

4 結論

1)將守恒變量守恒形式控制方程的時間項引入原始變量作為自變量，非定常數(shù)值解無法保證時間準確性。

2)針對性地構造基于原始變量的雙時間步方法，能夠恢復準確的數(shù)值解，該方法等價于基于守恒變量的雙時間步方法。

3)通過一維激波管問題與雙馬赫反射的數(shù)值實驗，上述兩個結論均得到驗證?；谠甲兞康碾p時間步方法可以推廣到一般的非定常流動問題中。