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

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

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

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

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

1 磁流體方程和能量定律

設Ω是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相關的拉格朗日乘子,對式(2)取散度,得到·Bt+Δr=0,由于得到r≡0。

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

對這兩個等式求和得

2 時間半離散格式

本節(jié)提出時間半離散格式,并證明其無條件能量穩(wěn)定性。對系統(tǒng)(1)～(6)構造如下格式,給定(un,Bn,pn,rn),通過以下步驟計算(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)

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

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

(15)

(16)

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

(17)

利用邊界條件rn+1|?Ω=0求解rn+1。然后計算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保持了反對稱結構,因此取檢驗函數(shù)并與之做內積,并利用((u·)v,v)=0的結論可以在穩(wěn)定性分析中將此非線性項消去,從而能夠簡化證明過程。

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

等式兩邊與自身做內積,得

(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é)調有限元空間為

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

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

(26)

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

(27)

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

(28)

(29)

(30)

(31)

(32)

(33)

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

因此可得

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

(34)

上式結合式(30)可得

(35)

同理可得

(36)

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

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

(37)

(38)

(39)

(40)

同理可得

(41)

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

(42)

其中

(43)

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

(44)

定理得證。

4 數(shù)值實驗

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

4.1 時間收斂階測試

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

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

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

4.2 空間收斂階測試

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

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 解耦格式當δ t=h時的誤差和收斂階

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

取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 耦合算法當δ t=h時的誤差和收斂階

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

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流描述了不可壓縮粘性導電流體在外部磁場Bd作用下沿均勻矩形截面管道的流動?？紤]邊界條件

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 結束語

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