一維非定常流體力學(xué)數(shù)值模擬

付智豪，余 聰

(中山大學(xué) 物理與天文學(xué)院，廣東 珠海 519082)

數(shù)值模擬在天體流體力學(xué)中具有重要意義，傳統(tǒng)的計算流體力學(xué)(Computational Fluid Dynamics，CFD)包括有限差分法和有限體積法等，而后又發(fā)展了諸如格子玻耳茲曼 (lattice Boltzmann method，LBM)和離散玻耳茲曼(discrete Boltzmann method,DBM)等計算方法[1-10].Riemann問題是CFD中的基本問題，即解決初始間斷下物理量的演化過程.此類方法除了可以有效的處理間斷解，也能自洽的求解待求區(qū)域內(nèi)光滑變化的物理量，十分適合研究具有復(fù)雜空間結(jié)構(gòu)和時間結(jié)構(gòu)的流場.對于激波管內(nèi)的Riemann問題，可求解出其理論解，但缺點是在復(fù)雜的初始條件下，模擬常常需要耗費大量的時間.在實際研究中，人們常常無法求得理論解析解，這時可以采用近似解的方法代替理論解.Roe通量差分裂格式作為CFD中最廣泛使用的格式[2]，適用性較廣，且精度較高，可以此研究不同初始條件下問題解的情況.本文僅考慮理想情況下的流體，從一維Euler方程出發(fā)，推導(dǎo)了Riemann問題下速度、密度及壓力的理論解，然后介紹了Roe通量差分裂格式及算法，以Sod、Lax及Shu-Osher激波管為算例檢驗算法的有效性，并在最后介紹了拓展至高維流體運動的方法.

1 理論解

1.1 控制方程

流體的運動十分復(fù)雜，但此處僅考慮無黏性的理想流體運動，且為了簡化問題只考慮一維情況，于是可采用守恒型的一維Euler方程進(jìn)行描述[1]：

(1)

1.2 Riemann問題的解

在求解方程組(1)時，如果設(shè)置密度、速度和壓強的初值為間斷型條件，那么流體物理量的演化會出現(xiàn)激波、稀疏波和接觸間斷的現(xiàn)象，此類問題便稱為Riemann問題.如圓柱管內(nèi)有一薄膜，當(dāng)兩側(cè)壓強差大到使其破裂，會形成沖向一側(cè)的激波，同時在另一側(cè)形成膨脹的稀疏波，二者間存在一個接觸間斷區(qū).激波即密度、壓強和速度在波陣面上發(fā)生突變的壓縮波.

Riemann問題是CFD中的核心問題之一，求出其精確解對解決流體運動具有重大意義.有限體積法便是將計算區(qū)域劃分為若干個區(qū)域，將每個小塊的邊界當(dāng)作間斷面，進(jìn)行Riemann問題近似求解，最終合并區(qū)域得到整個流場的運動.對于不同的初始條件，其解可分為5種情況[4]：1) 左激波+接觸間斷+右激波；2) 左稀疏波+接觸間斷+右激波；3) 左激波+接觸間斷+右稀疏波；4) 左稀疏波+接觸間斷+右稀疏波；5) 左稀疏波+真空區(qū)域+右稀疏波.

以情況2為例，對應(yīng)于管內(nèi)薄膜破裂，圖1給出解的x～t圖.Ⅰ和Ⅱ區(qū)為未受擾動區(qū)域，分別對應(yīng)稀疏波波前區(qū)域和激波波前區(qū)域，滿足初始條件.Ⅲ區(qū)對應(yīng)于稀疏波過后的波后區(qū)域，Ⅳ區(qū)對應(yīng)于激波過后的波后區(qū)域，3條實線間的Ⅴ區(qū)對應(yīng)于稀疏波區(qū)域.其中在Ⅲ區(qū)和Ⅳ間的虛線表示接觸間斷，Ⅳ區(qū)和Ⅱ間的實線表示激波間斷.

圖1 左稀疏波+接觸間斷+右激波

1.3 數(shù)學(xué)解釋

若將式(1)改寫為擬線性形式[4]：

(2)

A(U)=L-1ΛL

(3)

(4)

1.4 等熵流動

(5)

若流體運動無黏性且滿足絕熱條件，即等熵流動[4]，能量方程可用

(6)

代替.γ為絕熱指數(shù)，于是由式(5)得

(7)

其中u為流體的運動速度，c為聲速.第1式為描述特征線的方程，第2式為特征線上的不變量，又稱Riemann不變量，上下標(biāo)分別對應(yīng)向右傳播和向左傳播.

1.5 激波間斷

實際上的激波具有厚度，但數(shù)值僅是氣體分子自由程的數(shù)倍，故數(shù)學(xué)上可近似處理成為間斷.考慮流體積分形式的守恒型方程，設(shè)激波從左往右運動，運動速度為Z，左側(cè)物理量為(u1,p1,ρ1),右側(cè)物理量為(u2,p2,ρ2)，則間斷處滿足Rankine-Hugoniot(激波)跳躍關(guān)系式，簡稱R-H關(guān)系式[4].

(8)

1.6 具體求解

(9)

(10)

于是激波、稀疏波前后速度-壓力關(guān)系可寫成統(tǒng)一形式：

(11)

其中下標(biāo)i對應(yīng)波前區(qū)域的物理量.式(9)、式(10)相加得

u1-u2=f(p0,p1,ρ1)+f(p0,p2,ρ2)=F(p0)

(12)

若p1>p2，當(dāng)p0>p1時，中間區(qū)域壓力大于兩側(cè)，說明兩側(cè)均為激波，對應(yīng)情況1；當(dāng)p1>p0>p2，中間壓力比右側(cè)大比左側(cè)小，說明右側(cè)形成激波，左側(cè)形成稀疏波，對應(yīng)情況2；當(dāng)p2>p0，中間壓力比兩側(cè)都小，說明兩側(cè)均為稀疏波，對應(yīng)情況4；當(dāng)p0≤0，由于實際上壓力不能小于零，故為真空區(qū)，對應(yīng)情況5.若p1

圖2畫出了F(p)隨p的變化關(guān)系，可見其為單調(diào)遞增函數(shù)，于是p之間的大小關(guān)系等價于F(p)之間的大小關(guān)系.故流體的演化情況，也可先求出F(p1)、F(p2)，再與F(p0)=u1-u2比較進(jìn)行判斷.

圖2 (1,0,1)和(0.125,0,0.1)初始條件下F(p)與p的關(guān)系

由式(12)求解出中心區(qū)壓力p0后，代回式(9)或式(10)可求解出中心區(qū)速度u0.于是激波后區(qū)域內(nèi)的密度為

(13)

激波間斷的速度為

(14)

稀疏波后區(qū)域內(nèi)的密度，由絕熱過程公式可得

(15)

波頭速度及波尾速度為

(16)

(17)

(18)

(19)

(20)

2 Roe通量差分裂格式

Riemann問題雖然具有解析解，但實際運用中計算量較大，處理問題效果不是很好.若采取近似解的方法，不僅能夠精確描述流體運動，編程效率也能大大提高，下文便介紹近似解的Roe方法.

2.1 方程差分化

編程算法上采用有限體積法，對式(1)空間導(dǎo)數(shù)部分采取一階差分[8]：

(21)

(22)

得

(23)

2.2 基本原理

2.3 Roe格式

Roe給出了常矩陣構(gòu)造方法[2]，先由左右點值求得平均值：

(24)

(25)

(26)

(27)

(28)

2.4 熵修正

(29)

ε為設(shè)置參數(shù)，在此取為10-7.

3 數(shù)值算例

為了檢驗Roe算法的有效性，對以下經(jīng)典算例進(jìn)行了數(shù)值模擬，計算與畫圖均采用Matlab語言[8].

3.1 Sod激波管

Sod激波管，由兩個不連續(xù)分隔的恒定狀態(tài)組成，屬于經(jīng)典的Riemann問題[3,9].圖3展示了Roe格式下的數(shù)值解，并以理論解檢驗算法精度.發(fā)現(xiàn)其數(shù)值解與理論解值相符，但在間斷區(qū)域內(nèi)過渡平滑，這是由于算法精度不夠高導(dǎo)致的.

圖3 Sod激波管

如圖3所示，位于x=0.18的物理量均存在著一個突變，這對應(yīng)于一個向右傳播的激波.在與激波運動相反的方向，由于右側(cè)分子的減少，流體向外膨脹，形成了一個向左傳播的稀疏波.-0.12

3.2 Lax激波管

同Sod激波管一樣，也是由兩個不連續(xù)分隔的恒定狀態(tài)組成,但初始條件不同.圖4展示了Roe格式下的數(shù)值解，并以理論解檢驗算法精度.其數(shù)值解與理論解相符，間斷處的過渡是由于算法精度所限.

圖4 Lax激波管

3.3 Shu-Osher激波管

該問題能夠很好地測試算法捕捉平滑流動下激波相互作用的能力.初始條件是一個強激波，初始位置在x=-0.4，傳播到密度呈正弦變化的背景介質(zhì)中.圖5展示了用Roe格式求解出的密度、速度和壓強的位置分布，其中密度的位置分布與文獻(xiàn)[9]中求解出的趨勢一致,說明Roe格式算法具有較強的適用范圍.

圖5 Shu-Osher激波管

4 高維拓展

對于高維情況，只需對Euler方程進(jìn)行改寫，然后采用相應(yīng)的Roe通量差分裂格式[9]及差分化方程，便可進(jìn)行高維流體模擬計算.如3維Euler方程可寫為

(30)

則擬線性形式為

(31)

對于x方向，有

(32)

當(dāng)w=0時，方程對應(yīng)于2維流體模擬；當(dāng)v=w=0時，方程對應(yīng)于1維流體模擬.

5 結(jié)束語

非線性的流體偏微分方程中間斷解是物理中十分具有挑戰(zhàn)的問題，能夠?qū)げň?xì)結(jié)構(gòu)進(jìn)行有效追蹤和捕捉的算法是當(dāng)前流體力學(xué)的熱點之一.本文從一維Euler方程推導(dǎo)出其Riemann理論解，揭示了流體內(nèi)的激波、稀疏波和接觸間斷等現(xiàn)象物理本質(zhì).為了進(jìn)一步描述不同情況下的流體運動，文章采用了Roe格式，再次對Riemann問題求解，以檢驗其有效性.通過比較，發(fā)現(xiàn)一維下Roe格式與理論解符合度較高，能很好地描述激波管內(nèi)物理量的演化過程.由此認(rèn)為Roe格式是很好的近似解，繼而將該格式拓展至2維和3維形式.而且Roe格式的簡單形式使得其不僅能解決非定常的Euler方程，在考慮磁場效應(yīng)后，通過對方程的改寫，也能輕易的拓展到磁流體動力學(xué)的計算[9].而且采用更高精度的激波捕捉格式，如TVD、MUSCL和WENO等[3,4]，Roe格式求解與理論解幾乎無異.正是因為Roe格式在求解流體動力學(xué)性質(zhì)上的優(yōu)越性，使得其成為當(dāng)前天體物理應(yīng)用數(shù)值方法的主流方法.

本文僅考慮了理想狀態(tài)下的流體運動，但實際中，由于微觀粒子的復(fù)雜相互作用，在宏觀上存在著不同相與組分的離子界面，比熱、黏性系數(shù)、熱傳導(dǎo)率等不再保持為常量，基本平衡方程不再滿足，出現(xiàn)了非平衡的動力學(xué)及熱力學(xué)效應(yīng).采用NS和Euler建模在解決上述問題時存在一定的局限性，但近來發(fā)展的DBM方法從多尺度建模，不僅提高了對激波精細(xì)結(jié)構(gòu)的描述，而且在深度和廣度上都超越了原有建模[11].