張懷寶，王光學(xué)，王靖宇

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

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

(1)

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

(2)

對(duì)于一維Euler方程，上述各個(gè)變量的定義為：

(3)

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

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

1 數(shù)值方法

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

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

(4)

式中：

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

(5)

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

(6)

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

(7)

(8)

第n+1時(shí)間步的流場(chǎng)解為：

(9)

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

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

(10)

式中：

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

(11)

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

(12)

(13)

第n+1時(shí)間步的流場(chǎng)解為：

(14)

1.3 非定常解的時(shí)間準(zhǔn)確性

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

(15)

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

(16)

(17)

(18)

1.4 基于原始變量的雙時(shí)間步方程

(19)

對(duì)式(19)進(jìn)行一階顯式離散，得：

(20)

(21)

(22)

應(yīng)用近似關(guān)系式

(23)

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

(24)

2 激波管問(wèn)題與結(jié)果分析

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

表1 激波管問(wèn)題初值條件

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

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

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

(c) 溫度曲線 (c) Temperature curve圖1 三種數(shù)值解與解析解的原始變量對(duì)比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ù)值解與解析解的溫度對(duì)比放大圖Fig.2 Contrast magnification of the temperatures using three numerical approaches against analytical solution

將守恒變量方法的數(shù)值解作為準(zhǔn)確解，計(jì)算原始變量方法數(shù)值解的誤差范數(shù)，結(jié)果見(jiàn)表2。其中相對(duì)誤差范數(shù)的參考量采用驅(qū)動(dòng)段的初值條件，由于初速度為零，涉及速度的參考量均采用驅(qū)動(dòng)段的初始聲速a=347.3?？梢钥闯?，采用原始變量方法，各個(gè)原始變量(p,u,T)與各個(gè)守恒變量(ρ,ρu,ρE)均有不同程度的誤差。

表2 激波管問(wèn)題數(shù)值結(jié)果的誤差范數(shù)統(tǒng)計(jì)

(25)

的解。

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

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

3 雙馬赫反射問(wèn)題與結(jié)果分析

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

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

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

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

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

圖5 雙馬赫反射問(wèn)題t=0.2時(shí)密度云圖Fig.5 Density contours at t=0.2 for the double Mach reflection

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

4 結(jié)論

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

2)針對(duì)性地構(gòu)造基于原始變量的雙時(shí)間步方法，能夠恢復(fù)準(zhǔn)確的數(shù)值解，該方法等價(jià)于基于守恒變量的雙時(shí)間步方法。

3)通過(guò)一維激波管問(wèn)題與雙馬赫反射的數(shù)值實(shí)驗(yàn)，上述兩個(gè)結(jié)論均得到驗(yàn)證。基于原始變量的雙時(shí)間步方法可以推廣到一般的非定常流動(dòng)問(wèn)題中。