鄒世俊, 蔚喜軍, 戴自換

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

0 引言

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

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

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

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

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

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

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

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

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

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

1.3 拉氏保正格式

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

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

2 數(shù)值算例

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

2.1 超快稀疏波問(wèn)題

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

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

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

2.2 一維真空激波管問(wèn)題

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

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

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

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

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

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

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

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

3 總結(jié)

本文中，我們給出了一維理想MHD 系統(tǒng)的一種拉氏形式，并提出了一種在適當(dāng)?shù)男盘?hào)速度下可以保持密度和熱力學(xué)壓力正性的拉氏HLLD 數(shù)值通量。運(yùn)用這一HLLD 通量，我們構(gòu)造出了一種針對(duì)一維理想MHD 系統(tǒng)的保正拉氏格式。最后，我們給出了幾個(gè)具有挑戰(zhàn)性的數(shù)值算例來(lái)證明本文提出的方法所具有的保正性質(zhì)。

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