付智豪，余 聰

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

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

1 理論解

1.1 控制方程

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

(1)

1.2 Riemann問題的解

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

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

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

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

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

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

(2)

A(U)=L-1ΛL

(3)

(4)

1.4 等熵流動(dòng)

(5)

若流體運(yùn)動(dòng)無黏性且滿足絕熱條件，即等熵流動(dòng)[4]，能量方程可用

(6)

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

(7)

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

1.5 激波間斷

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

(8)

1.6 具體求解

(9)

(10)

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

(11)

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

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

(12)

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

圖2畫出了F(p)隨p的變化關(guān)系，可見其為單調(diào)遞增函數(shù)，于是p之間的大小關(guān)系等價(jià)于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問題雖然具有解析解，但實(shí)際運(yùn)用中計(jì)算量較大，處理問題效果不是很好.若采取近似解的方法，不僅能夠精確描述流體運(yùn)動(dòng)，編程效率也能大大提高，下文便介紹近似解的Roe方法.

2.1 方程差分化

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

(21)

(22)

得

(23)

2.2 基本原理

2.3 Roe格式

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

(24)

(25)

(26)

(27)

(28)

2.4 熵修正

(29)

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

3 數(shù)值算例

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

3.1 Sod激波管

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

圖3 Sod激波管

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

3.2 Lax激波管

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

圖4 Lax激波管

3.3 Shu-Osher激波管

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

圖5 Shu-Osher激波管

4 高維拓展

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

(30)

則擬線性形式為

(31)

對(duì)于x方向，有

(32)

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

5 結(jié)束語

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

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