不可壓縮磁流體方程的全解耦和無條件能量穩(wěn)定格式

周 帥,陳建華,王 琨,張國棟

(煙臺大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院,山東 煙臺 264005)

磁流體動(dòng)力學(xué)是研究等離子體等導(dǎo)電流體與電磁場的相互作用的物理學(xué)分支,在天體物理、熱核反應(yīng)、工業(yè)應(yīng)用[1-3]等眾多物理和工程分支中有著廣泛的應(yīng)用。導(dǎo)電流體可以在磁場的作用下產(chǎn)生電流,同時(shí)感應(yīng)電流與磁場相互作用,進(jìn)而形成洛倫茲力,并能改變流體的行為。因此,磁流體方程(MHD)將不可壓縮的納維-斯托克斯(Navier-Stokes)方程與麥克斯韋(Maxwell)方程耦合起來。近年來有許多關(guān)于MHD方程數(shù)值方法的研究[4-12]。對于穩(wěn)態(tài)的MHD方程,GUNZBURGER等[13]研究了標(biāo)準(zhǔn)的Galerkin有限元離散。SCH?TZAU等[14]提出了一種求解非凸多面體MHD方程的混合有限元方法。WACKER等[15]考慮了一種局部投影的方法來穩(wěn)定線性化MHD方程的兩個(gè)不足條件。

用數(shù)值方法求解非線性耦合的磁流體方程,一個(gè)具有挑戰(zhàn)性的問題是設(shè)計(jì)一種解耦線性的算法,同時(shí)保持能量穩(wěn)定性,即能量耗散定律能夠保持在離散的水平。全隱格式和耦合的半隱格式是無條件能量穩(wěn)定的,但是它們在每一個(gè)時(shí)間層上都需要解決一個(gè)耦合系統(tǒng),這會導(dǎo)致昂貴的時(shí)間消耗。所以,解耦是首先要解決的問題,在保持能量穩(wěn)定的同時(shí)也節(jié)省了時(shí)間成本。對于非定常納維-斯托克斯方程的求解,投影方法是最流行的解耦方法之一[16-18]。AN等[19]提出了一階線性耦合的投影方法求解單鞍點(diǎn)型磁流體方程,并證明了該方法的誤差估計(jì),由于方程是耦合的,所以在計(jì)算過程中時(shí)間消耗會比較大。CHOI等[5]發(fā)展了兩種基于投影方法求解單鞍點(diǎn)型磁流體方程的無條件穩(wěn)定的耦合格式。文獻(xiàn)[20-21]給出了單鞍點(diǎn)型磁流體方程的一階全離散投影格式的穩(wěn)定性證明和誤差分析。

本文針對雙鞍點(diǎn)型磁流體方程提出了一個(gè)線性的、無條件能量穩(wěn)定的、全解耦格式。對于模型的耦合困難,采用投影方法將壓力與速度解耦和拉格朗日乘子與磁場解耦,并且引入一個(gè)一階精度的穩(wěn)定項(xiàng)來穩(wěn)定磁場和速度場的解耦計(jì)算。對于模型的非線性困難,采用隱式-顯式方法來處理非線性項(xiàng),并且在非線性項(xiàng)的線性化處理中保持反對稱結(jié)構(gòu)。

1 磁流體方程和能量定律

設(shè)Ω是d(d=2,3)中的有界且單連通區(qū)域,考慮下面的磁流體方程:

(1)

(2)

·u=0,

(3)

·B=0,

(4)

u|?Ω=0,B×n|?Ω=0,r|?Ω=0,

(5)

u|t=0=u0(x),B|t=0=B0(x)。

(6)

其中,(x,t)∈Ω×(0,T],u為速度場,p為壓力場,B為磁場,r為拉格朗日乘子,物理參數(shù)Re,Rm和s分別為流體雷諾數(shù)、磁雷諾數(shù)和耦合系數(shù),n表示?Ω的單位外法向。

H(curl,Ω)={ω∈L2(Ω)d:×ω∈L2(Ω)d},

注1r是與約束·B=0相關(guān)的拉格朗日乘子,對式(2)取散度,得到·Bt+Δr=0,由于得到r≡0。

系統(tǒng)(1)～(6)遵循能量耗散定律:令式(1)與u作L2內(nèi)積,式(2)與sB作L2內(nèi)積,利用式(3)～(4)和分部積分公式,得到

對這兩個(gè)等式求和得

2 時(shí)間半離散格式

本節(jié)提出時(shí)間半離散格式,并證明其無條件能量穩(wěn)定性。對系統(tǒng)(1)～(6)構(gòu)造如下格式,給定(un,Bn,pn,rn),通過以下步驟計(jì)算(un+1,Bn+1,pn+1,rn+1):

(7)

(8)

求解un+1和pn+1滿足:

(9)

·un+1=0,

(10)

un+1·n|?Ω=0,

(11)

求解Bn+1和rn+1滿足:

(12)

·Bn+1=0,

(13)

rn+1|?Ω=0。

(14)

將上述三個(gè)式子聯(lián)立可以得到數(shù)值格式中的式(7)、(8)。

注2對式(9)取散度,可得pn+1滿足

(15)

(16)

注3對式(12)取散度,可得rn+1滿足

(17)

利用邊界條件rn+1|?Ω=0求解rn+1。然后計(jì)算Bn+1滿足

(18)

下面證明格式(7)～(14)的穩(wěn)定性。

定理1格式(7)～(14)具有如下能量穩(wěn)定性:

‖un+1‖2+s‖Bn+1‖2+ δt2‖pn+1‖2+sδt2‖rn+1‖2+

‖un‖2+s‖Bn‖2+ δt2‖pn‖2+sδt2‖rn‖2。

(19)

(20)

這里用到了分部積分

和

(21)

這里用到了

格式中(u·)u保持了反對稱結(jié)構(gòu),因此取檢驗(yàn)函數(shù)并與之做內(nèi)積,并利用((u·)v,v)=0的結(jié)論可以在穩(wěn)定性分析中將此非線性項(xiàng)消去,從而能夠簡化證明過程。

接下來,將式(9)重寫為

等式兩邊與自身做內(nèi)積,得

(22)

同理可得

(23)

將式(20)和(21)組合,利用式(22)和(23)得到

‖un+1‖2-‖un‖2+s‖Bn+1‖2-s‖Bn‖2+δt2‖pn+1‖2-δt2‖pn‖2+

(24)

其中

(25)

把式(25)代入式(24),可得式(19)。

3 全離散格式

定義協(xié)調(diào)有限元空間為

Ch:={B∈H0(curl,Ω):B|K∈Nk(K),K∈Th},

此外,有限元空間Xh和Mh必須滿足inf-sup條件

(26)

其中常數(shù)α與h無關(guān)。由于Sh?Ch,有限元空間Ch和Sh自然滿足inf-sup條件

(27)

基于格式(7)～(14)的全離散有限元格式如下:

(28)

(29)

(30)

(31)

(32)

(33)

注4對于式(29)中的對流項(xiàng),使用以下三線性形式:

因此可得

b(u,v,v)=0,u∈L2(Ω)d,v∈H1(Ω)d。

(34)

上式結(jié)合式(30)可得

(35)

同理可得

(36)

可清楚地看到格式(28)～(33)是線性的,全解耦的。

定理2 格式(28)～(33)具有如下的無條件能量穩(wěn)定性:

(37)

(38)

(39)

(40)

同理可得

(41)

最后將式(38)～(41)結(jié)合,可得

(42)

其中

(43)

將式(43)代入式(42)得

(44)

定理得證。

4 數(shù)值實(shí)驗(yàn)

在數(shù)值實(shí)驗(yàn)中,對于速度u,壓力p,以及拉格朗日乘子r,使用標(biāo)準(zhǔn)的P1b-P1-P1元,對于磁場B使用最低階的Nédélec′s元。這里,速度和壓力滿足inf-sup條件(26),磁場和拉格朗日乘子滿足inf-sup條件(27)。

4.1 時(shí)間收斂階測試

使用Ω=[0,1]×[0,1]上的解析解來計(jì)算時(shí)間收斂階,選擇精確解:

u=(ycos(t),xexp(t)),B=(yexp(t),xcos(t)),p=0,r=0,

表1 解耦格式的誤差和收斂階

4.2 空間收斂階測試

用Ω=[0,1]×[0,1]上的解析解來計(jì)算空間收斂階,選擇精確解:

u=(sin(2πy)sin(πx)sin(πx)exp(t),-sin(2πx)sin(πy)sin(πy)exp(t)),

B=(sin(πx)cos(πy)exp(t),-sin(πy)cos(πx)exp(t)),

p=(2x-1)(2y-1)cos(t),

r=0。

取s=Re=Rm=1,T=1,表2給出了解耦格式中速度的H1誤差‖eu‖H1,磁場的L2、H(curl)誤差(‖eb‖L2,‖eb‖curl)和壓力的L2誤差‖ep‖L2,它們具有一階精度,速度的L2誤差‖eu‖L2具有二階精度,拉格朗日乘子r接近于0。

表2 解耦格式當(dāng)δ t=h時(shí)的誤差和收斂階

同理,在Ω=[0,1]×[0,1]上使用相同的解析解來計(jì)算已有的耦合算法,方程如下:

取s=Re=Rm=1,T=1,表3給出了耦合算法中速度的H1誤差‖eu‖H1,磁場的L2、H(curl)誤差(‖eb‖L2,‖eb‖curl)和壓力的L2誤差‖ep‖L2,它們具有一階精度,速度的L2誤差‖eu‖L2具有二階精度,拉格朗日乘子r接近于0。

表3 耦合算法當(dāng)δ t=h時(shí)的誤差和收斂階

表2與表3相比較能夠得出本文提出的全解耦算法與已有的耦合算法在精度上大致相同,但在運(yùn)行速度方面,全解耦算法要比已有的耦合算法運(yùn)行速度快,所用的時(shí)間相對較短。

4.3 穩(wěn)定性測試

u0=(x2(x-1)2y(y-1)(2y-1),-y2(y-1)2x(x-1)(2x-1)),

B0=(sin(πx)cos(πy),-sin(πy)cos(πx)),

p0=0,

r0=0。

圖1 解耦格式的能量曲線

4.4 Hartmann流

Hartmann流描述了不可壓縮粘性導(dǎo)電流體在外部磁場Bd作用下沿均勻矩形截面管道的流動(dòng)?？紤]邊界條件

u=0,y=±1,

B×n=Bd×n, ?Ω。

該模型具有精確解

取Lx=4和G=1,考慮以下四種情形:

Ha=0.1,Re=Rm=0.1,s=1;

Ha=1.0,Re=Rm=1.0,s=1;

Ha=10,Re=Rm=10,s=1;

Ha=90,Re=Rm=30,s=9。

圖2和圖3分別給出了該格式在δt=1/100,h=1/128四種情形下的數(shù)值解u1和B1,可以發(fā)現(xiàn)它們與精確解很好地吻合。

圖2 數(shù)值解u1與精確解

圖3 數(shù)值解B1與精確解

5 結(jié)束語

對雙鞍點(diǎn)型磁流體方程提出了一個(gè)全解耦的線性高效格式,并嚴(yán)格證明了它的無條件穩(wěn)定性。該格式的主要特點(diǎn)是將耦合的非線性雙鞍點(diǎn)磁流體系統(tǒng)轉(zhuǎn)化為幾個(gè)線性橢圓問題,并在非線性項(xiàng)的線性化過程中保持了反對稱結(jié)構(gòu)。通過一系列的數(shù)值模擬,包括收斂性試驗(yàn)、能量穩(wěn)定性試驗(yàn)和Hartmann流試驗(yàn),驗(yàn)證了格式的穩(wěn)定性和收斂性。該格式的誤差估計(jì)將是后續(xù)工作的主要目標(biāo)。