一種求解一維理想磁流體方程組的保正拉氏方法

鄒世俊, 蔚喜軍, 戴自換

(1. 中國工程物理研究院研究生院，北京 100088; 2. 北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，北京 100088)

0 引言

可壓縮磁流體方程組(MHD)在天體物理、可控?zé)岷司圩儭⒔饘僖睙挼缺姸嗫蒲蓄I(lǐng)域中扮演著重要的角色。此方程組是將氣體動力學(xué)方程組與Maxwell 方程組結(jié)合起來，廣泛的應(yīng)用于描述導(dǎo)電流體在磁場中的動力學(xué)現(xiàn)象。當(dāng)忽略黏性，電阻及相對論效應(yīng)時，上述方程轉(zhuǎn)化為理想可壓縮磁流體方程組。這組方程組十分復(fù)雜以至于很難對其進(jìn)行數(shù)學(xué)分析和精確求解，因此，對于理想磁流體方程組設(shè)計(jì)精確魯棒的數(shù)值方法有著十分重要的意義。理想磁流體方程組可以被寫為如下形式

在物理上，上述系統(tǒng)所描述的物理現(xiàn)象中，密度ρ和熱力學(xué)壓力p總是非負(fù)的。但是，一個逼近上述系統(tǒng)的數(shù)值格式所得到的近似解并不能總是保持這種正性，這會引起系統(tǒng)失去雙曲性并導(dǎo)致數(shù)值不穩(wěn)定。為了解決這一問題，許多研究者提出了不同的數(shù)值技巧。然而，絕大多數(shù)此類方法都是在歐拉框架下提出的，即網(wǎng)格在空間中是固定的。在流體力學(xué)的數(shù)值模擬中，有另一種框架也被經(jīng)常使用，這就是拉格朗日框架。在拉氏框架下，計(jì)算網(wǎng)格隨流體運(yùn)動。因此，拉氏方法更加適合處理多介質(zhì)及自由界面問題。但是，在計(jì)算流體力學(xué)中保正拉氏格式很少被研究，更不用說是針對于MHD 系統(tǒng)設(shè)計(jì)的方法了。Bezard 和Desp′res[3]針對一維拉氏MHD 方程組設(shè)計(jì)了一種數(shù)值格式，這一格式滿足熵不等式并在適當(dāng)?shù)腃FL 條件下滿足保正性質(zhì)。Gallice[6]運(yùn)用簡單黎曼解的概念提出了一種熵保正Godunov 型格式，這一格式對拉氏和歐拉MHD 系統(tǒng)都可以保接觸間斷。近些年來，有一些保正拉氏格式[4,9]被提出。這些格式可以很好的運(yùn)作于歐拉方程并能得到高階精度。但是，所有這些格式都不能推廣到MHD 系統(tǒng)。

本文中，我們將關(guān)注的焦點(diǎn)放在針對于一維理想磁流體系統(tǒng)的保正拉氏算法設(shè)計(jì)上?；谖墨I(xiàn)[10]中在拉氏框架下針對二維可壓縮理想MHD 方程組設(shè)計(jì)的RKDG 方法，我們推導(dǎo)出針對一維理想MHD 方程組的拉氏格式。之后，針對拉氏格式推廣了HLLD 近似黎曼解[7]，這一黎曼解可以保持密度和熱力學(xué)壓力的正性。這種拉氏HLLD 數(shù)值通量的表達(dá)式與其在歐拉框架下的表達(dá)式并不相同。運(yùn)用這一特別的HLLD 數(shù)值通量，我們可以構(gòu)造出一種保正拉氏格式。最后，我們將會給出幾個一維中具有挑戰(zhàn)性的數(shù)值算例來證明我們方法的有效性。

本文結(jié)構(gòu)如下，在第1 節(jié)中，我們給出針對于一維理想MHD 方程組的拉氏格式。之后我們給出拉氏HLLD 數(shù)值通量的推導(dǎo)過程及其保正性質(zhì)。最后，通過運(yùn)用拉氏HLLD 數(shù)值通量，我們可以構(gòu)造出求解一維理想MHD 方程組的保正拉氏格式。在第2 節(jié)中，通過給出數(shù)值算例可以證明我們方法的保正性質(zhì)。最終，我們會在第3 節(jié)中給總結(jié)。

1 求解一維理想MHD 方程組的保正格式

1.1 一維理想MHD 方程組的拉氏格式

眾所周知，黎曼問題(8)真解的求解可能既困難又需要大量的計(jì)算，因此在構(gòu)造Godunov 型數(shù)值格式時經(jīng)常會使用近似黎曼解。基于HLLD 黎曼解[7]的良好性質(zhì)，我們在針對一維MHD 系統(tǒng)的方法中使用此黎曼解。

1.2 拉氏HLLD 通量及其保正性質(zhì)

圖1 拉氏框架下HLLD 通量的黎曼扇

引理1 拉氏HLLD 求解器具有保正性質(zhì)，如果SL和SR滿足

1.3 拉氏保正格式

我們將(7)中的Uh取為分片常數(shù)并運(yùn)用向前歐拉時間離散，可以得到拉氏框架下積分弱形式(6)的離散格式

圖2 應(yīng)用HLLD 求解器的拉氏格式在Ii×[tn,tn+1]中的波系演化

2 數(shù)值算例

在這一節(jié)中，我們給出幾個一維中具有挑戰(zhàn)性的數(shù)值算例來測試之前提出方法所具有的保正性質(zhì)。若不使用保正方法，則所有這些算例在計(jì)算時都會產(chǎn)生負(fù)壓。

2.1 超快稀疏波問題

在這一小節(jié)中，考慮一個存在超快稀疏波的激波管問題。其初始條件如下

在圖3 中，我們分別給出u0= 3 和u0= 3.1 時在t= 0.05 時刻密度ρ，速度分量ux和總壓的圖像。如文獻(xiàn)[7]中所述，在對眾多近似黎曼解進(jìn)行數(shù)值測試后可以發(fā)現(xiàn)，雖然在Mf= 3 時，所有測試過的近似黎曼解都不會產(chǎn)生負(fù)的密度和壓力，但是在Mf= 3.1 時有些卻無法計(jì)算下去。這是因?yàn)檫@些無法進(jìn)行計(jì)算的近似黎曼解并不保正。在圖3 中可以看到我們的格式對于Mf=3 和Mf=3.1 都運(yùn)算的很好，這一數(shù)值結(jié)果可以說明拉氏HLLD 黎曼解的保正性質(zhì)。注意到，圖3 中的圖像與文獻(xiàn)[7]中所給出的圖像有所不同，這是因?yàn)楸疚闹械慕Y(jié)果是運(yùn)用拉氏方法得到的，網(wǎng)格隨著流體移動到了計(jì)算區(qū)域兩端，在x=0 附近網(wǎng)格非常稀疏，x=0 附近密度趨于真空。

圖3 一維超快膨脹激波管問題的結(jié)果

2.2 一維真空激波管問題

在這一小節(jié)中，我們給出文獻(xiàn)[5]中提出的一維真空激波管問題的數(shù)值結(jié)果。問題的初始條件如下

計(jì)算區(qū)域取為[?0.5,0.5]。γ= 5/3 并應(yīng)用外流邊界條件。這一問題用來說明我們的拉氏保正格式可以處理擁有極低密度和壓力的問題。

在圖4 中，我們給出在t=0.1 時刻2 000 網(wǎng)格上計(jì)算得到的密度和壓力的數(shù)值結(jié)果。通過此圖可以觀察到，在x=?0.3 附近網(wǎng)格非常稀疏。與之前的數(shù)值算例類似，在拉氏方法中網(wǎng)格隨著流體移動，因此x=?0.3 附近趨于真空。通過和文獻(xiàn)[5]中的結(jié)果進(jìn)行比較，我們可以確信無論是低密度還是低壓力都捕捉得很好。注意到，如果不運(yùn)用保正方法來計(jì)算此問題，則計(jì)算程序?qū)驗(yàn)榉俏锢斫獾某霈F(xiàn)在幾個時間步內(nèi)崩潰。

圖4 運(yùn)用保正拉氏格式在2 000 網(wǎng)格上得到的一維真空激波管問題的密度和壓力

2.3 旋轉(zhuǎn)阿爾文波脈沖問題

本節(jié)，我們考慮文獻(xiàn)[1]中提出的旋轉(zhuǎn)阿爾文波脈沖問題。給出初始條件如下計(jì)算使用周期邊界條件，γ= 5/3。算例中，初始壓力比總能的萬分之一還要小。除此之外，該問題的解中還存在一個很強(qiáng)的阿爾文波間斷，大多數(shù)數(shù)值格式在計(jì)算此問題時都會產(chǎn)生負(fù)壓。

該問題的數(shù)值模擬采用了800 個計(jì)算網(wǎng)格。在圖5 中，我們給出拉氏保正格式在t=0.156 時刻得到的總能和壓力。注意到在解中存在兩個脈沖，因此成功的數(shù)值模擬應(yīng)該能很好的保持這一特點(diǎn)。與前面的數(shù)值算例類似，如果不使用保正方法計(jì)算此問題，則計(jì)算程序會在幾個時間步內(nèi)崩潰。

圖5 運(yùn)用保正拉氏格式在800 網(wǎng)格上，t=0.156 時，計(jì)算得到的旋轉(zhuǎn)阿爾文波脈沖問題的總能和熱壓

3 總結(jié)

本文中，我們給出了一維理想MHD 系統(tǒng)的一種拉氏形式，并提出了一種在適當(dāng)?shù)男盘査俣认驴梢员３置芏群蜔崃W(xué)壓力正性的拉氏HLLD 數(shù)值通量。運(yùn)用這一HLLD 通量，我們構(gòu)造出了一種針對一維理想MHD 系統(tǒng)的保正拉氏格式。最后，我們給出了幾個具有挑戰(zhàn)性的數(shù)值算例來證明本文提出的方法所具有的保正性質(zhì)。

在將來的工作中，我們將會探索如何對于多維MHD 系統(tǒng)構(gòu)造保正拉氏格式。在對拉氏MHD 系統(tǒng)特別是多維情形設(shè)計(jì)保證格式時存在數(shù)個困難。這些困難使得多維MHD 系統(tǒng)的保正拉氏格式構(gòu)造變?yōu)榱艘粋€十分復(fù)雜的工作。因此，這一工作仍然需要繼續(xù)研究。