不可压磁流体力学方程组的高精度数值解法

庄昕, 葛永斌, 袁冬芳

(1.河南工程学院理学院, 河南 郑州451191; 2.宁夏大学数学统计学院, 宁夏 银川750021;3.内蒙古科技大学理学院, 内蒙古 包头014010))

1.引言

磁流体力学(Magnetohydrodynamics, 简写为MHD)是结合经典流体力学和电动力学的方法研究磁流体(导电流体)在磁场中的运动规律.磁流体是一种特殊的功能材料, 在天体物理、工业生产以及受控热核反应等领域都具有十分广泛的应用.MHD是在非导电流体力学的基础上, 研究导电流体中流场和磁场的相互作用, 进行这种研究必须对经典流体力学方程(Naiver-Stokes方程, 简写为N-S方程)加以修正, 以便得到MHD基本方程组.MHD方程组在计算流体力学领域扮演着非常重要的角色, 寻求其精确而稳定的数值求解方法是众多科研工作者梦寐以求的目标, 鉴于问题的复杂性, 其研究远不如单个方程和N-S方程组充分, 高精度的计算方法[1−7]就更为少见.

Armero和Simo[8]对不可压MHD方程组进行了研究, 文章结合有限元近似检验了长期扩散和无条件非线性稳定问题的算法.Nesliturk和Tezer-Sezgin[9]用有限元方法求解不可压MHD管道流动问题, 流动受到电动力的驱动.目前在实际应用中通常会遇到低雷诺数和高雷诺数问题, 很难同时对两类问题取得数值解, 原因在于惯性项和粘性项的影响.Peaceman和Rachford[10]用ADI方法求解抛物型方程和椭圆型方程, 对于低雷诺数问题, 这种方法收敛于精确解的速度慢.为了克服PR-ADI格式的缺点, DAI[11]基于ZHANG[12]的差分格式提出了一种新的ADI格式求解二维抛物型方程,DAI的两层差分格式推广了PR-ADI格式,并叫做推广的ADI格式, 此格式克服了PR-ADI格式收敛速度慢的缺点.Navarro和Cabezas-Gomez[13]应用推广的ADI方法求解低雷诺数下二维不可压粘性MHD问题, 这类流动问题的研究涉及生物、化学、物理学、电磁学和工程学等多个学科的知识, 属于典型交叉学科研究范畴, 因此研究其精确稳定的数值解法具有十分重要的科学意义和实际工程应用价值.

本文提出了数值求解低雷诺数下二维不可压粘性MHD方程组的高阶紧致差分方法.将MHD方程组中密度-涡量-流函数形式的模型方程组经过高阶紧致差分离散后得到了空间四阶精度, 时间二阶精度的紧致差分格式.为了验证本文提出的高精度紧致差分方法的精确性和可靠性, 对有解析解的二维非定常不可压MHD方程组的初边值问题进行数值实验, 数值结果证明本文所建立的高阶紧致格式精确有效并且无条件稳定.高阶紧致格式更适用于求解低磁雷诺数下不可压粘性MHD问题, 本文正是基于这点进行了深入的研究.

2.不可压磁流体力学控制方程

在磁流体力学近似下, Maxwell方程组表示为:

上述Maxwell方程组中E为电场强度,B为磁感应强度,D=ε0E+P为电感应强度,H=B/µ0−M为磁场强度.J和ρc分别表示电流密度和电荷密度,ε0和µ0分别表示电解质常数和磁导率,P,M分别表示极化强度和磁化强度.在电磁学中联系电流密度J和电场强M之间关系的Ohm定律J=σ(E+V ×B)+ρcV, 其中σ为电导率,ρc为电荷密度,V是导电流体质点的速度.可忽略电流电荷体密度ρc, 故电磁学Ohm定律在磁流体力学近似下可表示为

正因为在磁流体力学中研究的电磁现象是低频的, 对于等离子体, 常可以忽略其中单个离子磁化性质和极化性质的影响, 即:P=0,M=0, 则

在方程(2.1)-(2.4)中, 方程(2.1)和(2.4)是一致的, (2.2)和(2.3)是一致的.故由(2.1)、(2.2)、(2.5)、(2.6)可得磁流体力学中控制磁场运动的方程―磁扩散方程:

记η=, 称为磁扩散系数, (2.8)式即磁扩散方程.

另外, 在磁场作用下, 控制不可压粘性流体流动的动力学方程为N-S方程:

其中,V是速度,p∗为压强.控制磁流体流动的方程组简化为动力学方程和磁扩散方程的联立方程组.引进特征参数：特征长度L, 特征速度U, 特征磁场B0, 雷诺数和磁雷诺数Re及Rem,磁普朗特数Prm, 将方程(2.4)、(2.8)-(2.10) 中所有物理量进行无量纲化：

可得控制磁流体流动的方程组, 记

针对二维不可压磁流体力学方程组, 引入涡量ω, 流函数ψ和a, 电流密度j

则得磁流体力学方程组的电流密度-涡量-流函数形式的方程组：

其中,f1和f2为非齐次项.

3.数值方法

首先在直角坐标系中将二维计算区域用一簇平行于坐标轴的直线族进行网格剖分.设τ表示时间步长, 空间取等间距网格, 步长用h表示.则网格点为(xi,yj,tn),xi=ih,yj=jh,i,j=0,1··· ,m,h= 1/m,n ≥0 ;m表示x和y方向的网格等分数.为了方便起见, 用数字0, 1, 2, 3,4, 5, 6, 7, 8分别代表网格节点(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=g(x,y).当Φ=ψ,g=−ω时, 为流函数方程(2.20); 当Φ=a,g=−jRem时, 为流函数方程(2.22); 当Φ=ω,g=Re[(ωt+ψyωx −ψxωy)−N(ayjx −axjy)] +f1为涡量方程(2.19);当Φ=j,g=Rem[jt+(ψyjx −ψxjy)]−(ayωx −axωy)+f2, 为电流密度方程(2.21).

由文[2], 二阶导数在点(xi,yj)可以利用下面四阶紧致差分公式：

其中

将(3.2)和(3.3)代入方程(3.1)的左边整理可

将(3.4)代入(3.5)式, 略去高阶项整理可得方程(3.1)的四阶紧致差分格式

注意到

则(3.6)式即

对于流函数方程(2.20), 取Φ=ψ,g=−ω结合(3.6)式得流函数方程的四阶紧致差分格式：

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

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

结合(2.19)-(2.22)式可得：

为保证格式具有四阶精度, 避免网格节点超出给定的九点模板, 下面两式可化为：

其中

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

将上述离散公式代入(3.12)式可得涡量方程(2.19)的四阶全隐紧致差分格式, 整理后可得下形式

上式中,n ≥1且系数An+1,Bn,Cn−1均为涡量在第n+1时间层上的值.注意到这是一个三层格式, 除了知道初始时刻的值外, 还须知道第一个时间步的值, 方可接着向下推进.所以上述推导过程的最后一步, 只需将时间项以向后差分来代替, 即时间离散用下式代替：且令即可得到涡量在第一个时间步上的离散格式：

对于电流密度方程(2.21)取Φ=j,g=Rem(jt+ψyjx −ψxjy)−(ayωx −axωy) +f2, 则由(3.8)得

格式的推导过程同涡量方程的格式推导, 整理后可得如下公式：

将ψ,a,ω和j的各阶偏导数在九点模板中的二阶精度的离散公式差分公式代入(3.16)式可得电流密度方程(2.21)的四阶全隐紧致差分格式, 整理后可得如下形式:

其中,n ≥1 且系数A′n+1,B′n,C′n−1均为电流密度在第n+1时间层上的值.注意到这是一个三层格式, 除了知道初始时刻的值外, 还须知道第一个时间步的值, 方可接着向下推进.所以上述推导过程的最后一步, 只需将时间项以向后差分来代替, 即时间离散用下式代替ut=un+1−un/τ+O(τ)：且令n=0即可得到电流密度在第一个时间步上的离散格式:

4.数值算例

为了验证本文的高精度紧致差分格式的精确性和可靠性, 针对有解析解的二维非定常不可压磁流体力学方程组的初边值问题进行数值实验.在方程(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)式, 可得非齐次项函数

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

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

假设问题具有Dirichlet边界条件, 计算区域在(0,1)×(0,1)上, 且0≤T ≤1.

表4.1和4.2给出了当时间步长取τ=h2,Re=Rem= 0.0001, 0.01,N= 1.0时, 涡量ω, 流函数ψ,a, 及电流密度j的最大绝对误差Eω,Eψ,Ea,Ej和收敛阶Rate = ln(E1/E2)/ln 2, 其中E1和E2分别为粗网格(网格步长为h)和相邻的细网格上(网格步长为h/2)的最大绝对误差.从表4.1和4.2的数据结果可以看出, 采用本文方法计算涡量、流函数以及电流密度等的最大误差已基本达到四阶精度, 与文章理论分析结果一致.

表4.1 当Re=Rem =0.0001,N =1.0时, 涡量、流函数及电流密度的最大误差和收敛阶

表4.3-4.4给出了当空间步长分别取1/16,N分别为时1.0,100,Re=Rem=0.01, 不同网格比r(r=τ/h2)下, 在t=1时刻的时间推进步数以及涡量ω、流函数ψ,a, 及电流密度j的最大误差Eω,Eψ,Ea,Ej.从表4.3和4.4可以看出, 在不同的网格比值下, 涡量、流函数以及电流密度等计算结果均收敛, 因此本文格式是无条件稳定的全隐紧致差分格式.另外, 由于本文方法时间精度为二阶, 为了空间上达到四阶精度, 时间步长要和空间步长的平方保持同一数量级, 因此时间步长不能去太大.

表4.3 当h =1/16,Re=Rem =0.01,N =1.0,t=1时, 不同r, 涡量、流函数及电流密度的最大误差和收敛阶

表4.2 当Re=Rem =0.01,N =1.0时, 涡量、流函数及电流密度的最大误差和收敛阶

表4.4 当h =1/16,Re=Rem =0.01,N =100,t=1时, 不同r, 涡量、流函数及电流密度的最大误差和收敛阶

图4.1和图4.2分别给出当τ=h2,t= 1,Re=Rem= 0.01,N= 1.0, 网格为16×16时, 流函数ψ和a的精确解和数值解; 图4.3给出涡量ω和电流密度j的误差曲面.从流函数ψ和a的精确解图和数值解图对比可以看出这两个量的数值解和解析解吻合的相当好, 从涡量ω和电流密度j的误差曲面图可以看出涡量的误差较小, 而电流密度的误差相对较大.

图4.3 涡量ω的误差曲面(左)和电流密度j的误差曲面(右)

图4.2 流函数a的精确解(左)和数值解(右)

图4.1 流函数ψ的精确解(左)和数值解(右)

5.结论

本文结合涡量-流函数方法提出了数值求解多变量、强耦合、非线性二维不可压MHD问题的四阶全隐式紧致差分方法, 并对有精确解的MHD问题进行了有效的数值模拟.本文格式空间为四阶精度, 时间为二阶精度, 实验结果与精确解吻合的非常好, 结果证实了本文所提方法的精确性和稳定性.由于高阶紧致差分格式具有小的离散子域并且精度高, 边界条件容易处理, 在计算靠近计算区域边界的网格点时不需要考虑边界外面的网格点, 格式的稳定性好.因此, 与传统高精度差分方法相比, 本文全隐式高精度紧致格式更适用于求解低雷诺数下不可压粘性MHD问题.本文研究工作为数值模拟不可压粘性磁流体力学问题提供了一种切实可行的方法, 可直接推广到三维磁流体力学方程组的数值求解中去, 我们将会另文报道.