不可壓磁流體力學(xué)方程組的高精度數(shù)值解法

莊昕, 葛永斌, 袁冬芳

(1.河南工程學(xué)院理學(xué)院, 河南 鄭州451191; 2.寧夏大學(xué)數(shù)學(xué)統(tǒng)計(jì)學(xué)院, 寧夏 銀川750021;3.內(nèi)蒙古科技大學(xué)理學(xué)院, 內(nèi)蒙古 包頭014010))

1.引言

磁流體力學(xué)(Magnetohydrodynamics, 簡寫為MHD)是結(jié)合經(jīng)典流體力學(xué)和電動(dòng)力學(xué)的方法研究磁流體(導(dǎo)電流體)在磁場中的運(yùn)動(dòng)規(guī)律.磁流體是一種特殊的功能材料, 在天體物理、工業(yè)生產(chǎn)以及受控?zé)岷朔磻?yīng)等領(lǐng)域都具有十分廣泛的應(yīng)用.MHD是在非導(dǎo)電流體力學(xué)的基礎(chǔ)上, 研究導(dǎo)電流體中流場和磁場的相互作用, 進(jìn)行這種研究必須對經(jīng)典流體力學(xué)方程(Naiver-Stokes方程, 簡寫為N-S方程)加以修正, 以便得到MHD基本方程組.MHD方程組在計(jì)算流體力學(xué)領(lǐng)域扮演著非常重要的角色, 尋求其精確而穩(wěn)定的數(shù)值求解方法是眾多科研工作者夢寐以求的目標(biāo), 鑒于問題的復(fù)雜性, 其研究遠(yuǎn)不如單個(gè)方程和N-S方程組充分, 高精度的計(jì)算方法[1?7]就更為少見.

Armero和Simo[8]對不可壓MHD方程組進(jìn)行了研究, 文章結(jié)合有限元近似檢驗(yàn)了長期擴(kuò)散和無條件非線性穩(wěn)定問題的算法.Nesliturk和Tezer-Sezgin[9]用有限元方法求解不可壓MHD管道流動(dòng)問題, 流動(dòng)受到電動(dòng)力的驅(qū)動(dòng).目前在實(shí)際應(yīng)用中通常會遇到低雷諾數(shù)和高雷諾數(shù)問題, 很難同時(shí)對兩類問題取得數(shù)值解, 原因在于慣性項(xiàng)和粘性項(xiàng)的影響.Peaceman和Rachford[10]用ADI方法求解拋物型方程和橢圓型方程, 對于低雷諾數(shù)問題, 這種方法收斂于精確解的速度慢.為了克服PR-ADI格式的缺點(diǎn), DAI[11]基于ZHANG[12]的差分格式提出了一種新的ADI格式求解二維拋物型方程,DAI的兩層差分格式推廣了PR-ADI格式,并叫做推廣的ADI格式, 此格式克服了PR-ADI格式收斂速度慢的缺點(diǎn).Navarro和Cabezas-Gomez[13]應(yīng)用推廣的ADI方法求解低雷諾數(shù)下二維不可壓粘性MHD問題, 這類流動(dòng)問題的研究涉及生物、化學(xué)、物理學(xué)、電磁學(xué)和工程學(xué)等多個(gè)學(xué)科的知識, 屬于典型交叉學(xué)科研究范疇, 因此研究其精確穩(wěn)定的數(shù)值解法具有十分重要的科學(xué)意義和實(shí)際工程應(yīng)用價(jià)值.

本文提出了數(shù)值求解低雷諾數(shù)下二維不可壓粘性MHD方程組的高階緊致差分方法.將MHD方程組中密度-渦量-流函數(shù)形式的模型方程組經(jīng)過高階緊致差分離散后得到了空間四階精度, 時(shí)間二階精度的緊致差分格式.為了驗(yàn)證本文提出的高精度緊致差分方法的精確性和可靠性, 對有解析解的二維非定常不可壓MHD方程組的初邊值問題進(jìn)行數(shù)值實(shí)驗(yàn), 數(shù)值結(jié)果證明本文所建立的高階緊致格式精確有效并且無條件穩(wěn)定.高階緊致格式更適用于求解低磁雷諾數(shù)下不可壓粘性MHD問題, 本文正是基于這點(diǎn)進(jìn)行了深入的研究.

2.不可壓磁流體力學(xué)控制方程

在磁流體力學(xué)近似下, Maxwell方程組表示為:

上述Maxwell方程組中E為電場強(qiáng)度,B為磁感應(yīng)強(qiáng)度,D=ε0E+P為電感應(yīng)強(qiáng)度,H=B/μ0?M為磁場強(qiáng)度.J和ρc分別表示電流密度和電荷密度,ε0和μ0分別表示電解質(zhì)常數(shù)和磁導(dǎo)率,P,M分別表示極化強(qiáng)度和磁化強(qiáng)度.在電磁學(xué)中聯(lián)系電流密度J和電場強(qiáng)M之間關(guān)系的Ohm定律J=σ(E+V ×B)+ρcV, 其中σ為電導(dǎo)率,ρc為電荷密度,V是導(dǎo)電流體質(zhì)點(diǎn)的速度.可忽略電流電荷體密度ρc, 故電磁學(xué)Ohm定律在磁流體力學(xué)近似下可表示為

正因?yàn)樵诖帕黧w力學(xué)中研究的電磁現(xiàn)象是低頻的, 對于等離子體, ?？梢院雎云渲袉蝹€(gè)離子磁化性質(zhì)和極化性質(zhì)的影響, 即:P=0,M=0, 則

在方程(2.1)-(2.4)中, 方程(2.1)和(2.4)是一致的, (2.2)和(2.3)是一致的.故由(2.1)、(2.2)、(2.5)、(2.6)可得磁流體力學(xué)中控制磁場運(yùn)動(dòng)的方程―磁擴(kuò)散方程:

記η=, 稱為磁擴(kuò)散系數(shù), (2.8)式即磁擴(kuò)散方程.

另外, 在磁場作用下, 控制不可壓粘性流體流動(dòng)的動(dòng)力學(xué)方程為N-S方程:

其中,V是速度,p?為壓強(qiáng).控制磁流體流動(dòng)的方程組簡化為動(dòng)力學(xué)方程和磁擴(kuò)散方程的聯(lián)立方程組.引進(jìn)特征參數(shù)：特征長度L, 特征速度U, 特征磁場B0, 雷諾數(shù)和磁雷諾數(shù)Re及Rem,磁普朗特?cái)?shù)Prm, 將方程(2.4)、(2.8)-(2.10) 中所有物理量進(jìn)行無量綱化：

可得控制磁流體流動(dòng)的方程組, 記

針對二維不可壓磁流體力學(xué)方程組, 引入渦量ω, 流函數(shù)ψ和a, 電流密度j

則得磁流體力學(xué)方程組的電流密度-渦量-流函數(shù)形式的方程組：

其中,f1和f2為非齊次項(xiàng).

3.數(shù)值方法

首先在直角坐標(biāo)系中將二維計(jì)算區(qū)域用一簇平行于坐標(biāo)軸的直線族進(jìn)行網(wǎng)格剖分.設(shè)τ表示時(shí)間步長, 空間取等間距網(wǎng)格, 步長用h表示.則網(wǎng)格點(diǎn)為(xi,yj,tn),xi=ih,yj=jh,i,j=0,1··· ,m,h= 1/m,n ≥0 ;m表示x和y方向的網(wǎng)格等分?jǐn)?shù).為了方便起見, 用數(shù)字0, 1, 2, 3,4, 5, 6, 7, 8分別代表網(wǎng)格節(jié)點(diǎn)(xi,yj),(xi+1,yj),(xi?1,yj), (xi,yj?1),(xi+1,yj+1), (xi?1,yj+1),(xi?1,yj?1) , (xi+1,yj?1) , 二維空間離散子域圖參考文[4].

將方程(2.20)記為如下Poisson方程的形式：

其中g(shù)=g(x,y).當(dāng)Φ=ψ,g=?ω時(shí), 為流函數(shù)方程(2.20); 當(dāng)Φ=a,g=?jRem時(shí), 為流函數(shù)方程(2.22); 當(dāng)Φ=ω,g=Re[(ωt+ψyωx ?ψxωy)?N(ayjx ?axjy)] +f1為渦量方程(2.19);當(dāng)Φ=j,g=Rem[jt+(ψyjx ?ψxjy)]?(ayωx ?axωy)+f2, 為電流密度方程(2.21).

由文[2], 二階導(dǎo)數(shù)在點(diǎn)(xi,yj)可以利用下面四階緊致差分公式：

其中

將(3.2)和(3.3)代入方程(3.1)的左邊整理可

將(3.4)代入(3.5)式, 略去高階項(xiàng)整理可得方程(3.1)的四階緊致差分格式

注意到

則(3.6)式即

對于流函數(shù)方程(2.20), 取Φ=ψ,g=?ω結(jié)合(3.6)式得流函數(shù)方程的四階緊致差分格式：

同理可得方程(2.22)的四階緊致差分格式：

對于渦量方程(2.19), 取Φ=ω,g=Re[(ωt+ψyωx ?ψxωy)?N(ayjx ?axjy)] +f1,利用(3.8)式整理可得:

結(jié)合(2.19)-(2.22)式可得：

為保證格式具有四階精度, 避免網(wǎng)格節(jié)點(diǎn)超出給定的九點(diǎn)模板, 下面兩式可化為：

其中

將上述四式代入式(3.11), 整理后可得如下公式：

將上述離散公式代入(3.12)式可得渦量方程(2.19)的四階全隱緊致差分格式, 整理后可得下形式

上式中,n ≥1且系數(shù)An+1,Bn,Cn?1均為渦量在第n+1時(shí)間層上的值.注意到這是一個(gè)三層格式, 除了知道初始時(shí)刻的值外, 還須知道第一個(gè)時(shí)間步的值, 方可接著向下推進(jìn).所以上述推導(dǎo)過程的最后一步, 只需將時(shí)間項(xiàng)以向后差分來代替, 即時(shí)間離散用下式代替：且令即可得到渦量在第一個(gè)時(shí)間步上的離散格式：

對于電流密度方程(2.21)取Φ=j,g=Rem(jt+ψyjx ?ψxjy)?(ayωx ?axωy) +f2, 則由(3.8)得

格式的推導(dǎo)過程同渦量方程的格式推導(dǎo), 整理后可得如下公式：

將ψ,a,ω和j的各階偏導(dǎo)數(shù)在九點(diǎn)模板中的二階精度的離散公式差分公式代入(3.16)式可得電流密度方程(2.21)的四階全隱緊致差分格式, 整理后可得如下形式:

其中,n ≥1 且系數(shù)A′n+1,B′n,C′n?1均為電流密度在第n+1時(shí)間層上的值.注意到這是一個(gè)三層格式, 除了知道初始時(shí)刻的值外, 還須知道第一個(gè)時(shí)間步的值, 方可接著向下推進(jìn).所以上述推導(dǎo)過程的最后一步, 只需將時(shí)間項(xiàng)以向后差分來代替, 即時(shí)間離散用下式代替ut=un+1?un/τ+O(τ)：且令n=0即可得到電流密度在第一個(gè)時(shí)間步上的離散格式:

4.數(shù)值算例

為了驗(yàn)證本文的高精度緊致差分格式的精確性和可靠性, 針對有解析解的二維非定常不可壓磁流體力學(xué)方程組的初邊值問題進(jìn)行數(shù)值實(shí)驗(yàn).在方程(2.19)-(2.22)中, 令

ψ=sin(x+y ?t),

a=sin(x+y+t),

ω=2sin(x+y ?t),

j=2/Rem·sin(x+y+t).

將ψ,ω,a,j代入(2.20)和(2.22)式, 可得非齊次項(xiàng)函數(shù)

f1=2Re·cos(x+y ?t)?4sin(x+y ?t),

f2=?4/Rem·sin(x+y+t)?2cos(x+y+t).

假設(shè)問題具有Dirichlet邊界條件, 計(jì)算區(qū)域在(0,1)×(0,1)上, 且0≤T ≤1.

表4.1和4.2給出了當(dāng)時(shí)間步長取τ=h2,Re=Rem= 0.0001, 0.01,N= 1.0時(shí), 渦量ω, 流函數(shù)ψ,a, 及電流密度j的最大絕對誤差Eω,Eψ,Ea,Ej和收斂階Rate = ln(E1/E2)/ln 2, 其中E1和E2分別為粗網(wǎng)格(網(wǎng)格步長為h)和相鄰的細(xì)網(wǎng)格上(網(wǎng)格步長為h/2)的最大絕對誤差.從表4.1和4.2的數(shù)據(jù)結(jié)果可以看出, 采用本文方法計(jì)算渦量、流函數(shù)以及電流密度等的最大誤差已基本達(dá)到四階精度, 與文章理論分析結(jié)果一致.

表4.1 當(dāng)Re=Rem =0.0001,N =1.0時(shí), 渦量、流函數(shù)及電流密度的最大誤差和收斂階

表4.3-4.4給出了當(dāng)空間步長分別取1/16,N分別為時(shí)1.0,100,Re=Rem=0.01, 不同網(wǎng)格比r(r=τ/h2)下, 在t=1時(shí)刻的時(shí)間推進(jìn)步數(shù)以及渦量ω、流函數(shù)ψ,a, 及電流密度j的最大誤差Eω,Eψ,Ea,Ej.從表4.3和4.4可以看出, 在不同的網(wǎng)格比值下, 渦量、流函數(shù)以及電流密度等計(jì)算結(jié)果均收斂, 因此本文格式是無條件穩(wěn)定的全隱緊致差分格式.另外, 由于本文方法時(shí)間精度為二階, 為了空間上達(dá)到四階精度, 時(shí)間步長要和空間步長的平方保持同一數(shù)量級, 因此時(shí)間步長不能去太大.

表4.3 當(dāng)h =1/16,Re=Rem =0.01,N =1.0,t=1時(shí), 不同r, 渦量、流函數(shù)及電流密度的最大誤差和收斂階

表4.2 當(dāng)Re=Rem =0.01,N =1.0時(shí), 渦量、流函數(shù)及電流密度的最大誤差和收斂階

表4.4 當(dāng)h =1/16,Re=Rem =0.01,N =100,t=1時(shí), 不同r, 渦量、流函數(shù)及電流密度的最大誤差和收斂階

圖4.1和圖4.2分別給出當(dāng)τ=h2,t= 1,Re=Rem= 0.01,N= 1.0, 網(wǎng)格為16×16時(shí), 流函數(shù)ψ和a的精確解和數(shù)值解; 圖4.3給出渦量ω和電流密度j的誤差曲面.從流函數(shù)ψ和a的精確解圖和數(shù)值解圖對比可以看出這兩個(gè)量的數(shù)值解和解析解吻合的相當(dāng)好, 從渦量ω和電流密度j的誤差曲面圖可以看出渦量的誤差較小, 而電流密度的誤差相對較大.

圖4.3 渦量ω的誤差曲面(左)和電流密度j的誤差曲面(右)

圖4.2 流函數(shù)a的精確解(左)和數(shù)值解(右)

圖4.1 流函數(shù)ψ的精確解(左)和數(shù)值解(右)

5.結(jié)論

本文結(jié)合渦量-流函數(shù)方法提出了數(shù)值求解多變量、強(qiáng)耦合、非線性二維不可壓MHD問題的四階全隱式緊致差分方法, 并對有精確解的MHD問題進(jìn)行了有效的數(shù)值模擬.本文格式空間為四階精度, 時(shí)間為二階精度, 實(shí)驗(yàn)結(jié)果與精確解吻合的非常好, 結(jié)果證實(shí)了本文所提方法的精確性和穩(wěn)定性.由于高階緊致差分格式具有小的離散子域并且精度高, 邊界條件容易處理, 在計(jì)算靠近計(jì)算區(qū)域邊界的網(wǎng)格點(diǎn)時(shí)不需要考慮邊界外面的網(wǎng)格點(diǎn), 格式的穩(wěn)定性好.因此, 與傳統(tǒng)高精度差分方法相比, 本文全隱式高精度緊致格式更適用于求解低雷諾數(shù)下不可壓粘性MHD問題.本文研究工作為數(shù)值模擬不可壓粘性磁流體力學(xué)問題提供了一種切實(shí)可行的方法, 可直接推廣到三維磁流體力學(xué)方程組的數(shù)值求解中去, 我們將會另文報(bào)道.