变系数MHD-Boussinesq方程的二阶BDF时间离散格式的误差分析
Error Analysis of a Second-Order BDF Time Discretization Scheme for the Variable-Coefficient MHD-Boussinesq Equations
摘要: 基于外推线性化方法,本文研究了具有变系数的二维不可压缩磁流体动力学方程组和热对流扩散方程所耦合MHD-Boussinesq方程组的二阶BDF时间离散格式,理论上证明了该格式具有无条件稳定性。在合理的正则性假设下,给出了速度、压力、磁场和温度在离散范数意义下的最优时间收敛阶,证明了格式具有二阶时间精度。最后通过数值实验验证了理论分析的结果。
Abstract: Based on the extrapolated linearization method, this paper investigates a second-order BDF time-discrete scheme for the coupled MHD-Boussinesq system, which consists of the two-dimensional incompressible magnetohydrodynamic equations with variable coefficients and the heat convection-diffusion equation. The unconditional stability of the scheme is theoretically proven. Under reasonable regularity assumptions, optimal temporal convergence orders in the discrete norm are derived for the velocity, pressure, magnetic field, and temperature, demonstrating that the scheme achieves second-order temporal accuracy. Finally, numerical experiments are conducted to validate the theoretical findings.
文章引用:徐轩. 变系数MHD-Boussinesq方程的二阶BDF时间离散格式的误差分析[J]. 应用数学进展, 2026, 15(2): 65-79. https://doi.org/10.12677/aam.2026.152050

1. 引言

本文考虑具有粘性、电导率和热扩散系数依赖温度的变系数不可压缩磁流体动力学方程组和热对流扩散方程所耦合的MHD-Boussinesq方程组,该方程组在二维区域 Ω×( 0,T ] 上具有如下非线性耦合形式:

t u( υ( θ )u )+( u )u+p+B×curlB=θ e 2 +f (1.1)

t B+curl( κ( θ )curlB )curl( u×B )=g (1.2)

t θ( σ( θ )θ )+uθ=ε u 2 +ψ (1.3)

u=0,B=0, (1.4)

其中 Ω R 2 是一个有界凸区域。未知变量 u( x,t ),p( x,t ),B( x,t ) θ( x,t ) 分别表示速度场、压力、磁场和温度。粘性 υ( θ ) 、电导率 κ( θ ) 和热扩散率 σ( θ ) 是依赖于温度的正值函数。 εR 是一个常数, e 2 是由 e 2 = ( 0,1 ) T 给出的单位向量, θ e 2 代表浮力。给定的函数 f,g ψ 分别是磁感应的体积力、热源和外加电流,且在分布意义下满足 g=0 .方程(1.1)~(1.4)中的叉乘和旋度算子定义为:

u×B= u 1 B 2 u 2 B 1 ,             B×r=( r B 2 ,r B 1 ), curlB= B 2 x B 1 y ,             curlr= r y r x ,

其中 u=( u 1 , u 2 ),B=( B 1 , B 2 ) r 是一个标量函数。此外,我们赋予方程(1.1)~(1.4)如下初始条件:

u( x,0 )= u 0 ( x ),      B( x,0 )= B 0 ( x ),      θ( x,0 )= θ 0 ( x ) (1.5)

和边界条件:

u=0,     θ=0,Bn=0,     n×curlB=0,      Ω. (1.6)

在上述MHD-Boussinesq方程组中,方程(1.1)表示具有浮升力 θ e 2 、洛伦兹力 B×curlB 以及体积力 f 影响的动量守恒定律。方程(1.2)表明电磁场由麦克斯韦方程主导,而方程(1.3)描述了由速度和热源驱动的温度变化情况。无散度条件 u=0 表示流体的不可压缩特性, B=0 则意味着磁场中不存在磁单极子。显然,(1.1)~(1.4)式通过线性Boussinesq近似将不可压缩MHD方程与热方程耦合起来。该方程组在物理和工程领域具有广泛应用(参见文献[1]-[6])。

关于MHD-Boussinesq方程组及其相关问题的数值方法,已有一些研究结果。当 θ=0,ε=0 时,MHD-Boussinesq方程退化为不可压缩MHD方程,国内外学者有大量关于MHD方程组数值算法的研究,如Gao和Qiu在[7]中研究一阶Euler半隐式有限元格式,在低正则性假设下获得了半隐有限元格式的最优误差估计。Li和Luo在[8]中研究了一类MHD方程组的二阶半隐式Crank-Nicolson有限元格式,并给出了有限元误差估计。Zhang,Hou和Shan在[9]中对于三维MHD方程研究了二阶线性化Crank-Nicolson格式。当 B=0,ε=0 时,MHD-Boussinesq方程退化为Boussinesq方程。对于具有常系数粘性系数和热扩散率的Boussinesq方程,[10][11]分别研究了二阶Crank-Nicolson有限元格式和BDF2有限元格式。在上述研究中,所考虑的问题均基于常系数粘性、电导率和热扩散率的假设。针对变系数MHD-Boussinesq方程的数值算法,Li和Wang在[12]中提出了线性化且无条件稳定的一阶Euler半隐式有限元格式,其中用Mini有限元离散速度与压力,用分片线性有限元离散磁场和温度,在合理的正则性假设以及CFL条件下,得到了速度,磁场和温度的最优误差估计。

本文考虑具有变系数的二维不可压缩磁流体动力学方程组和热对流扩散方程所耦合MHD-Boussinesq方程组,给出一个基于外推线性化的二阶BDF时间离散格式,理论证明了该格式具有无条件稳定性和二阶时间精度,并通过数值实验验证了其有效性。

在本文中,符号 C 表示一个正常数,其与时间步长 τ 无关,且在不同位置取值可能不同。

2. 预备知识

我们使用粗体符号 H k ( Ω ) W k,p ( Ω ) L p ( Ω ) 表示向量值空间 H k ( Ω ) 2 W k,p ( Ω ) 2 L p ( Ω ) 2 ,其中Sobolev空间和Lebesgue空间 H k ( Ω ) 2 W k,p ( Ω ) L p ( Ω ) 按照经典方式定义(参见文献[1]),其范数分别记为 H k W k,p L p 。空间 H 0 1 ( Ω ) 中的函数在边界 Ω 上迹为0。

X= H 0 1 ( Ω ) X 0 ={ vX,v=0Ω }

W={ B H 1 ( Ω ),Bn=0Ω } W 0 ={ BW,B=0Ω }

M= L 0 2 ( Ω )={ q L 2 ( Ω ), Ω qdx =0 } Y= H 0 1 ( Ω )

其范数分别为:

v x = v L 2 q M = q L 2 θ Y = θ L 2

B W = curlB L 2 + B L 2 ,且 B=0 .

根据Poincare不等式和Sobolev嵌入定理,可得:

v L p C v L 2 vX 1p<+ (2.1)

φ L p C φ L 2 φY 1p<+ . (2.2)

对于 u,v,wX θ,φY 我们定义三线性项为:

b 1 ( u,v,w )= Ω ( u )vwdx ,     b 2 ( u,θ,φ )= Ω ( uθ )φdx .

u=0 ,根据分部积分,易知上述三重线性项满足反对称性质,即:

b 1 ( u,v,w )= b 1 ( u,w,v ) b 2 ( u,θ,φ )= b 2 ( u,φ,θ ) (2.3)

由此可得:

b 1 ( u,v,v )=0 b 2 ( u,θ,θ )=0 (2.4)

此外

b 1 ( u,v,w )= Ω ( u )vwdx b 2 ( u,θ,φ )= Ω ( uθ )φdx . (2.5)

且有以下不等式:

b 1 ( u,v,w )C u L 2 1/2 u L 2 1/2 v L 2 w L 2 (2.6)

b 1 ( u,θ,φ )C u L 2 1/2 u L 2 1/2 θ L 2 φ L 2 . (2.7)

为了得到时间收敛阶,我们将作出如下假设。

假设1:初始值满足

u 0 X 0 H 2 ( Ω ) B 0 W 0 H 2 ( Ω ) θ 0 Y H 2 ( Ω ) . (2.8)

假设2:粘性系数 υ( θ ) ,电导率 κ( θ ) 和热扩散率 σ( θ ) 是三个光滑函数且满足

α 1 υ( z ),κ( z ),σ( z )α,  | υ ( z ) |,| κ ( z ) |,| σ ( z ) |β,  zR . (2.9)

其中常数 α1 β0

假设3:MHD-Boussinesq方程组(1.1)~(1.4)存在唯一全局光滑解,且满足以下正则性条件:

( u,B ) L ( 0,T; H 2 ( Ω ) ) L 2 ( 0,T; H 3 ( Ω ) ) (2.10)

θ L ( 0,T; H 2 ( Ω ) ) L 2 ( 0,T; H 3 ( Ω ) ) (2.11)

( u t , B t ) L ( 0,T; L 2 ( Ω ) ) L 2 ( 0,T; H 1 ( Ω ) ) (2.12)

θ t L ( 0,T; L 2 ( Ω ) ) L 2 ( 0,T; H 0 1 ( Ω ) ) (2.13)

( u tt , B tt ) L 2 ( 0,T; L 2 ( Ω ) ) θ tt L 2 ( 0,T; L 2 ( Ω ) ) (2.14)

( u ttt , B ttt ) L 2 ( 0,T; L 2 ( Ω ) ) θ ttt L 2 ( 0,T; L 2 ) . (2.15)

注2.1:在条件(2.10)及 υ( z ),κ( z ),σ( z ) α 1 的前提下,[13]已证明当 f=g=0 ψ=0 时,MHD-Boussinesq方程组(1.1)~(1.4)满足正则性条件(2.12)~(2.15)。因此,我们要求 ( f,g ) L ( 0,T; L 2 ( Ω ) ) L 2 ( 0,T; H 1 ( Ω ) ) ψ L ( 0,T; L 2 ( Ω ) ) L 2 ( 0,T; H 1 ( Ω ) ) 以确保(2.16)~(2.15)成立。

3. 二阶BDF时间离散格式

本节给出用于求解方程组(1.1)~(1.4)的二阶BDF时间离散格式,并证明数值格式的无条件稳定性。

N Ν + 0= t 0 < t 1 << t N =T 为时间区间 [ 0,T ] 的部分,其固定时间步长为 τ=T/N ,且对 0nN t n =nτ . 对于定义在 [ 0,T ] 上的任意函数 g(t) ,我们定义

D τ g n+1 = 3 g n+1 4 g n + g n1 2τ g ^ n =2 g n g n1 1nN1 .

对任意 a,b,cR ,有以下恒等式成立:

2( 3a4b+c )a= | a | 2 | b | 2 + ( 2ab ) 2 ( 2bc ) 2 ( a2b+c ) 2

进而可得telescope公式:

( D τ g n+1 , g n+1 )= 1 4τ ( g n+1 L 2 2 g n L 2 2 + g ^ n+1 L 2 2 g ^ n L 2 2 )+ 1 4τ g n+1 2 g n + g n1 L 2 2 ,   1nN1. (3.1)

此外,由带积分余项的泰勒公式可得:

D τ g n+1 g t ( t n+1 )= 1 2τ t n1 t n+1 ( 2 ( t t n ) + 2 1 2 ( t t n1 ) 2 ) g ttt ( t )dt (3.2)

g ^ ( t n )g( t n+1 )= t n1 t n+1 ( 2 ( t t n ) + ( t t n1 ) ) g tt ( t )dt (3.3)

其中 ( t t n ) + =max{ t t n ,0 } . 则对任意范数 都成立

g ^ ( t n )g( t n+1 ) 2 C τ 3 t n1 t n+1 g tt ( t ) 2 dt (3.4)

D τ g n+1 g t ( t n+1 ) 2 C τ 3 t n1 t n+1 g ttt ( t ) 2 dt . (3.5)

基于上述记号并采用经典外推线性化方法,我们给出如下二阶BDF时间离散格式。

第一步:对于给定的 u 0 , B 0 θ 0 ,我们采用以下一阶欧拉时间离散格式求解 u 1 ,p, B 1 θ 1

1 τ ( u 1 u 0 )+( υ( θ 0 ) u 1 )+ u 0 u 1 + B 0 ×curl B 1 + p 1 = θ 0 e 2 + f 1 (3.6)

u 1 =0 (3.7)

1 τ ( B 1 B 0 )+curl( κ( θ 0 )curl B 1 )curl( u 1 × B 0 )= g 1 (3.8)

B 1 =0 (3.9)

1 τ ( θ 1 θ 0 )+( σ( θ 0 ) θ 1 )+ u 1 θ 1 =ε u 2 1 + ψ 1 . (3.10)

第二步:当 1nN1 时,对于给定的 u ^ n , B ^ n , θ ^ n ,我们通过以下方式求解 u n+1 ,p, B n+1 θ n+1

D τ u n+1 +( υ( θ ^ n ) u n+1 )+ u ^ n u n+1 + B ^ n ×curl B n+1 + p n+1 = θ ^ n e 2 + f n+1 (3.11)

u n+1 =0 (3.12)

D τ B n+1 +curl( κ( θ ^ n )curl B n+1 )curl( u n+1 × B ^ n )= g n+1 (3.13)

B n+1 =0 (3.14)

D τ θ n+1 +( σ( θ ^ n ) θ n+1 )+ u n+1 θ n+1 =ε u 2 n+1 + ψ n+1 . (3.15)

为证明上述二阶BDF时间离散格式的无条件稳定性和二阶时间收敛精度,我们需要用到[11]中所提供的离散Gronwall不等式。

引理3.1:设 a k , b k , γ k 为正实数且 B>0 ,使得对 n1

a n +τ k=1 n b k τ k=1 n γ k a k +B (3.16)

则当 τ 足够小满足 τ γ k <1 且令 σ k = ( 1τ γ k ) 1 时,对 a n +τ k=1 n b k exp( τ k=1 n γ k σ k )B 成立

a n +τ k=1 n b k exp( τ k=1 n γ k σ k )B . (3.17)

下述引理给出了所提数值格式的无条件稳定性。

引理3.2在假设2的条件下,对于充分小的 τ 1nN1 ,由(3.6)~(3.15)定义的数值解 u n+1 , B n+1 , θ n+1 满足:

u n+1 L 2 2 + B n+1 L 2 2 + θ n+1 L 2 2 + u ^ n+1 L 2 2 + B ^ n+1 L 2 2 + θ ^ n+1 L 2 2 + α 1 τ i=0 n ( u i+1 L 2 2 + B i+1 W 2 + θ i+1 L 2 2 ) C( u 0 L 2 2 + B 0 L 2 2 + θ 0 L 2 2 + u ^ 1 L 2 2 + B ^ 1 L 2 2 + θ ^ 1 L 2 2 )+Cτ i=0 n ( f i+1 L 2 2 + g i+1 L 2 2 + ψ i+1 L 2 2 ) (3.18)

证明:将(3.10)式两端乘以 2τ θ 1 并在 Ω 上积分,利用 ( ab,a )= 1 2 ( | a | 2 | b | 2 + | ab | 2 ) 可得:

θ 1 L 2 2 θ 0 L 2 2 +2 α 1 τ θ 1 L 2 2 2 α 1 τ( u 2 1 , θ 1 )+2τ( ψ 1 , θ 1 ) α 1 τ θ 1 L 2 2 +Cτ( ψ 1 L 2 2 + u 1 L 2 2 ) (3.19)

将(3.6)式,(3.8)式分别与 2τ u 1 2τ B 1 作内积并在 Ω 上积分,利用假设2,(3.1)式以及 ( a×curlb,c )=( c×a,curlb ) ,相加得:

( u 1 L 2 2 + B 1 L 2 2 )( u 0 L 2 2 + B 0 L 2 2 )+2 α 1 τ( u 1 L 2 2 + B 1 W 2 ) 2τ( θ 0 e 2 , u 1 )+2τ( f 1 , u 1 )+2τ( g 1 , B 1 ) α 1 τ( u 1 L 2 2 + B 1 W 2 )+Cτ( θ 0 L 2 2 + f 1 L 2 2 + g 1 L 2 2 ) (3.20)

将(3.19)和(3.20)相加,可得:

( u 1 L 2 2 + B 1 L 2 2 + θ 1 L 2 2 )+ α 1 τ( u 1 L 2 2 + B 1 W 2 + θ 1 L 2 2 ) u 0 L 2 2 + B 0 L 2 2 + θ 0 L 2 2 +Cτ u 1 L 2 2 + θ 1 L 2 2 +Cτ f 1 L 2 2 + g 1 L 2 2 + ψ 1 L 2 2 (3.21)

类似的,将(3.15)式两端乘以 4τ θ n+1 并在 Ω 上积分,利用(2.11)式和(3.1)式可得:

θ n+1 L 2 2 θ n L 2 2 + θ ^ n+1 L 2 2 θ ^ n L 2 2 +4 α 1 τ θ n+1 L 2 2 4 α 1 τ( u 2 n+1 , θ n+1 )+4τ( ψ n+1 , θ n+1 ) 3 α 1 τ θ n+1 L 2 2 +Cτ( ψ n+1 L 2 2 + u n+1 L 2 2 ) (3.22)

将(3.11)式,(3.13)式分别与 4τ u n+1 4τ B n+1 作内积并在 Ω 上积分,利用假设2,(3.1)式以及 ( a×curlb,c )=( c×a,curlb ) 可得:

( u n+1 L 2 2 + B n+1 L 2 2 )( u n L 2 2 + B n L 2 2 )+( u ^ n+1 L 2 2 + B ^ n+1 L 2 2 ) ( u ^ n L 2 2 + B ^ n L 2 2 )+4 α 1 τ( u n+1 L 2 2 + B n+1 W 2 ) 4τ( θ ^ n e 2 , u n+1 )+4τ( f n+1 , u n+1 )+4τ( g n+1 , B n+1 ) 3 α 1 τ( u n+1 L 2 2 + B n+1 W 2 )+Cτ( θ ^ n L 2 2 + f n+1 L 2 2 + g n+1 L 2 2 ) (3.23)

将(3.22)和(3.23)相加后求和,可得:

( u n+1 L 2 2 + B n+1 L 2 2 + θ n+1 L 2 2 )+( u ^ n+1 L 2 2 + B ^ n+1 L 2 2 + θ ^ n+1 L 2 2 ) + α 1 τ i=1 n ( u i+1 L 2 2 + B i+1 W 2 + θ i+1 L 2 2 ) u 1 L 2 2 + B 1 L 2 2 + θ 1 L 2 2 + u ^ 1 L 2 2 + B ^ 1 L 2 2 + θ ^ 1 L 2 2 +Cτ i=1 n ( u i+1 L 2 2 + θ i L 2 2 ) +Cτ i=1 n ( f i+1 L 2 2 + g i+1 L 2 2 + ψ i+1 L 2 2 ) (3.24)

注意到(3.21),可得:

( u n+1 L 2 2 + B n+1 L 2 2 + θ n+1 L 2 2 )+( u ^ n+1 L 2 2 + B ^ n+1 L 2 2 + θ ^ n+1 L 2 2 ) + α 1 τ i=0 n ( u i+1 L 2 2 + B i+1 W 2 + θ i+1 L 2 2 ) u 0 L 2 2 + B 0 L 2 2 + θ 0 L 2 2 + u ^ 1 L 2 2 + B ^ 1 L 2 2 + θ ^ 1 L 2 2 +Cτ i=0 n ( u i+1 L 2 2 + θ i L 2 2 ) +Cτ i=0 n ( f i+1 L 2 2 + g i+1 L 2 2 + ψ i+1 L 2 2 )

利用引理3.1中的离散Gronwall不等式,对于充分小的 τ( τ>0 ) ,我们得到离散能量不等式(3.18),从而完成引理3.2的证明。

4. 误差分析

n=0 时,在(1.1)~(1.4)中取 t= t 1 ,可以看出 ( u( t 1 ),p( t 1 ),B( t 1 ),θ( t 1 ) ) 满足:

u( t 1 )u( t 0 ) τ ( υ( θ( t 0 ) )u( t 1 ) )+u( t 0 )u( t 1 ) +B( t 0 )×curlB( t 1 )+p( t 1 )=θ( t 0 ) e 2 + f 1 + R u 1 ,          u 1 =0, (4.1)

B( t 1 )B( t 0 ) τ +curl( κ( θ( t 0 ) )curlB( t 1 ) )curl( u( t 1 )×B( t 0 ) ) = g 1 + R B 1 ,        B 1 =0,   (4.2)

θ( t 1 )θ( t 0 ) τ ( σ( θ( t 0 ) )θ( t 1 ) )+u( t 1 )θ( t 1 )=ε u 2 ( t 0 )+ ψ 1 + R θ 1 . (4.3)

其中,截断误差函数 R u 1 , R B 1 , R θ 1 由下式给出:

R u 1 = u( t 1 )u( t 0 ) τ u t ( t 1 )+( υ( θ( t 1 ) )u( t 1 ) )( υ( θ( t 0 ) )u( t 1 ) )+u( t 0 )u( t 1 ) u( t 1 )u( t 1 )+B( t 0 )×curl( B( t 1 ) )B( t 1 )×curl( B( t 1 ) )+θ( t 1 ) e 2 θ( t 0 ) e 2 ,

R B 1 = B( t 1 )B( t 0 ) τ B t ( t 1 )+curl( κ( θ( t 0 ) )curlB( t 1 ) )curl( κ( θ( t 1 ) )curlB( t 1 ) ) curl( u( t 1 )×B( t 0 ) )+curl( u( t 1 )×B( t 1 ) ),

R θ 1 = θ( t 1 )θ( t 0 ) τ θ t ( t 1 )+( σ( θ( t 1 )θ( t 0 ) )u( t 1 ) ).

1nN1 时,在(1.1)~(1.4)式中取 t= t n+1 ,可以看出 ( u( t n+1 ),p( t n+1 ),B( t n+1 ),θ( t n+1 ) )

满足

D τ u( t n+1 )( υ( θ ^ ( t n ) )u( t n+1 ) )+ u ^ ( t n )u( t n+1 ) + B ^ ( t n )×curlB( t n+1 )+p( t n+1 )=θ( t n ) e 2 + f n+1 + R u n+1 ,    u n+1 =0, (4.4)

D τ B( t n+1 )+curl( κ( θ ^ ( t n ) )curlB( t n+1 ) )curl( u( t n+1 × B ^ ( t n ) ) )= g n+1 + R B n+1 ,  B n+1 =0, (4.5)

D τ θ( t n+1 )( σ( θ ^ ( t n ) )θ( t n+1 ) )+u( t n+1 )θ( t n+1 )=ε u 2 ( t n+1 )+ ψ n+1 + R θ n+1 . (4.6)

其中,截断误差函数 R u n+1 , R B n+1 , R θ n+1 由下式给出:

R u n+1 = D τ u( t n+1 ) u t ( t n+1 )( υ( θ ^ ( t n )θ( t n+1 ) )u( t n+1 ) )+( u ^ ( t n )u( t n+1 ) )u( t n+1 )        +( B ^ ( t n )B( t n+1 ) )×curlB( t n+1 )+( θ( t n+1 ) θ ^ ( t n ) ) e 2 ,

R B n+1 = D τ B( t n+1 ) B t ( t n+1 )+curl( κ( θ ^ ( t n )θ( t n+1 ) )curlB( t n+1 ) )          curl( u( t n+1 )×( B ^ ( t n )B( t n+1 ) ) ),

R θ n+1 = D τ θ( t n+1 ) θ t ( t n+1 )( σ( θ ^ ( t n )θ( t n+1 ) )θ( t n+1 ) ).

基于taylor公式和解的正则性假设(2.10)~(2.15),并利用(3.4)~(3.5),可得截断误差的以下估计:

R u 1 L 2 2 + R B 1 L 2 2 + R θ 1 L 2 2 C τ 2 (4.7)

τ i=1 N1 ( R u i+1 L 2 2 + R B i+1 L 2 2 + R θ i+1 L 2 2 ) C τ 4 . (4.8)

e u n+1 =u( t n+1 ) u n+1 , e p n+1 =p( t n+1 ) p n+1 , e B n+1 =B( t n+1 ) B n+1 , e θ n+1 =θ( t n+1 ) θ n+1

e u 0 = e B 0 = e θ 0 =0

接下来,本文的主要结果由以下定理给出。

引理4.1:在假设1~3的条件下,当 τ 足够小时,时间离散系统(1.1)~(1.4)有如下的估计成立:

e u 1 L 2 2 + e B 1 L 2 2 + e θ 1 L 2 2 + α 1 τ( e u 1 L 2 2 + e B 1 W 2 + e θ 1 L 2 2 )C τ 4 (4.9)

证明:

n=0 时,将(4.1),(4.2),(4.3)式减去(3.6),(3.7),(3.8)式可得误差方程:

D τ e u 1 υ( θ( t 0 ) ) e u 1 = R u 1 (4.10)

D τ e B 1 +curl( κ( θ( t 0 ) ) )curl e B 1 = R B 1 (4.11)

D τ e θ 1 σ( θ( t 0 ) ) e θ 1 +( e u 1 )θ( t 1 )+ u 1 e θ 1 = R θ 1 +ε e u2 1 (4.12)

第一步:估计 e u 1 e B 1

将(4.10)式与 2τ e u 1 作内积并在 Ω 上积分,由 ( ab,a )= 1 2 ( | a | 2 | b | 2 + | ab | 2 ) 可以得到关于 e u 1 的不等式:

e u 1 L 2 2 +2 α 1 τ e u 1 L 2 2 2τ( R u 1 , e u 1 ) (4.13)

对上式右端各项进行精细估计。关键项的处理如下:

截断误差项,易证:

2τ( R u 1 , e u 1 ) 1 2 e u 1 L 2 2 +C τ 2 R u 1 L 2 2 (4.14)

将(4.14)代入(4.13)式中,我们可以得到:

e u 1 L 2 2 + α 1 τ e u 1 L 2 2 C τ 2 R u 1 L 2 2 (4.15)

将(4.11)式与 2τ e B 1 作内积并在 Ω 上积分,由 ( ab,a )= 1 2 ( | a | 2 | b | 2 + | ab | 2 ) 可以得到关于 e B 1 的不等式:

e B 1 L 2 2 +2 α 1 τ e B 1 W 2 2τ( R B 1 , e B 1 ) (4.16)

对上式右端各项进行精细估计。关键项的处理如下:

截断误差项,易证:

2τ( R u 1 , e u 1 ) 1 2 e B 1 L 2 2 +C τ 2 R B 1 L 2 2 (4.17)

将(4.17)式代入(4.16)式中,我们可以得到:

e B 1 L 2 2 + α 1 τ e B 1 W 2 C τ 2 R B 1 L 2 2 (4.18)

将(4.15)和(4.18)式相加,我们可以得到:

e u 1 L 2 2 +  e B 1 L 2 2 + α 1 τ e u 1 L 2 2 + α 1 τ e B 1 W 2 C τ 2 R u 1 L 2 2 +C τ 2 R B 1 L 2 2 (4.19)

第二步:估计 e θ 1

将(4.12)式与 2τ e θ 1 作内积并在 Ω 上积分,由 ( ab,a )= 1 2 ( | a | 2 | b | 2 + | ab | 2 ) 可以得到关于 e θ 1 的不等式:

e θ 1 L 2 2 +2 α 1 τ e θ 1 L 2 2 2τε( e u2 1 , e θ 1 )+2τ( R θ 1 , e θ 1 ) (4.20)

对上式右端各项进行精细估计。关键项的处理如下:

1) 热源项,由Young不等式可得:

2τε( e u2 1 , e θ 1 ) 1 2 α 1 τ e u 1 L 2 2 + 3 2 α 1 τ e θ 1 L 2 2 (4.21)

2) 截断误差项,易证:

2τ( R θ 1 , e θ 1 ) 1 2 e θ 1 L 2 2 +C τ 2 R θ 1 L 2 2 (4.22)

将(4.21)和(4.22)代入(4.20)中,我们可以将(4.20)式改写为:

e θ 1 L 2 2 + α 1 τ e θ 1 L 2 2 α 1 τ e u 1 L 2 2 +C τ 2 R θ 1 L 2 2 (4.23)

将(4.19)和(4.23)式相加,由Young不等式和(4.7)可以得到:

e u 1 L 2 2 + e B 1 L 2 2 + e θ 1 L 2 2 + α 1 τ e u 1 L 2 2 + α 1 τ e B 1 W 2 + α 1 τ e θ 1 L 2 2 C τ 2 ( R u 1 L 2 2 + R B 1 L 2 2 + R θ 1 L 2 2 )C τ 4

定理4.1:在假设1~3的条件下,当 τ 足够小时,时间离散系统(1.1)~(1.4)有如下的估计成立:

e u m+1 L 2 2 + e B m+1 L 2 2 + e θ m+1 L 2 2 +τ n=0 m ( e u n+1 L 2 2 + curl e B n+1 L 2 2 + e θ n+1 L 2 2 ) C 0 τ 4 (4.24)

0mN1 ,其中 C 是与时间步长 τ 无关的常数。

证明:

1nN1 时,将(4.4),(4.5),(4.6)式减去(3.6),(3.7),(3.8)式可得误差方程:

D τ e u n+1 ( ( υ( θ ^ ( t n ) )υ( θ ^ n ) )u( t n+1 ) )υ( θ ^ ( t n ) ) e u n+1 +( e ^ u n )u( t n+1 ) + u ^ n e u n+1 + e ^ B n ×curl( B( t n+1 ) )+ B ^ n ×curl e B n+1 = e ^ θ n + R u n+1 (4.25)

D τ e B n+1 +curl( κ( θ ^ ( t n ) )curl e B n+1 )+curl( κ( θ ^ ( t n )κ( θ ^ n ) )curlB( t n+1 ) ) curl( e u n+1 × B ^ ( t n ) )curl( u n+1 × e ^ B n )= R B n+1 (4.26)

D τ e θ n+1 ( ( σ( θ ^ ( t n ) )σ( θ ^ n ) )u( t n+1 ) )σ( θ ^ ( t n ) ) e θ n+1 +( e u n+1 )θ( t n+1 )+ u n+1 e θ n+1 = R θ n+1 +ε e u2 n+1 (4.27)

第一步:估计 e u n+1 e B n+1

将(4.25)式与 4τ e u n+1 作内积并在 Ω 上积分,由(2.6) (2.11)和(3.1)可以得到关于 e u n+1 的不等式:

( e u n+1 L 2 2 e u n L 2 2 + e ^ u n+1 L 2 2 e ^ u n L 2 2 )+4 α 1 τ e u n+1 L 2 2 4τ( ( υ( θ ^ ( t n ) )υ( θ ^ n ) )u( t n+1 ), e u n+1 )4τb( ( e ^ u n ),u( t n+1 ), e u n+1 )   4τ( e ^ B n ×curl( B( t n+1 ) ), e u n+1 )4τ( B ^ n ×curl e B n+1 , e u n+1 )   +4τ( e ^ θ n , e u n+1 )+4τ( R u n+1 , e u n+1 ) := i=1 5 I i n+1 ( e u n+1 ) 4τ( B ^ n ×curl e B n+1 , e u n+1 ) (4.28)

对上式右端各项进行精细估计。关键项的处理如下:

1) 非线性项,由拉格朗日中值定理、Holder不等式、Young不等式、(2.11)以及(2.14)可知:

I 1 n+1 ( e u n+1 )=4τ( ( υ( θ ^ ( t n ) )υ( θ ^ n ) )u( t n+1 ), e u n+1 ) 4βτ e ^ θ n L 2 u( t n+1 ) L e u n+1 L 2 3 5 α 1 τ e u n+1 L 2 2 +Cτ e ^ θ n L 2 2 u( t n+1 ) L 2 3 5 α 1 τ e u n+1 L 2 2 +Cτ e ^ θ n L 2 2 (4.29)

2) 对流项,由Holder不等式、Young不等式以及(2.14)可知:

I 2 n+1 ( e u n+1 )=4τb( ( e ^ u n ),u( t n+1 ), e u n+1 ) 4τ e ^ u n L 2 u( t n+1 ) L 3 e u n+1 L 6 3 5 α 1 τ e u n+1 L 2 2 +Cτ u( t n+1 ) L 3 2 e ^ u n L 2 2 3 5 α 1 τ e u n+1 L 2 2 +Cτ e ^ u n L 2 2 (4.30)

3) 洛伦兹力项,我们采用与估计 I 2 n+1 ( e u n+1 ) 的方法,可得:

I 3 n+1 ( e u n+1 )=4τ( e ^ B n ×curlB( t n+1 ), e u n+1 ) 4τ e ^ B n L 2 curlB( t n+1 ) L 3 e u n+1 L 6 3 5 α 1 τ e u n+1 L 2 2 +Cτ u( t n+1 ) L 3 2 e ^ B n L 2 2 3 5 α 1 τ e u n+1 L 2 2 +Cτ e ^ B n L 2 2 (4.31)

4) 易证:

I 4 n+1 ( e u n+1 )=4τ( e ^ θ n , e u n+1 ) 3 5 α 1 τ e u n+1 L 2 2 +Cτ e ^ θ n L 2 2 (4.32)

5) 截断误差项,易证:

I 5 n+1 ( e u n+1 )=4τ( R u n+1 , e u n+1 ) 3 5 α 1 τ e u n+1 L 2 2 +Cτ R u n+1 L 2 2 (4.33)

将(4.29)~(4.33)相加,我们将(4.28)式重新写为:

( e u n+1 L 2 2 e u n L 2 2 + e ^ u n+1 L 2 2 e ^ u n L 2 2 )+ α 1 τ e u n+1 L 2 2 Cτ e ^ θ n L 2 2 +Cτ e ^ u n L 2 2 +Cτ e ^ B n L 2 2 +Cτ R u n+1 L 2 2 4τ( B ^ n ×curl e B n+1 , e u n+1 ) (4.34)

将(4.26)式与 4τ e B n+1 作内积并在 Ω 上积分,利用(2.6) (2.11)和(3.1)可以得到关于 e B n+1 的不等式:

( e B n+1 L 2 2 e B n L 2 2 + e ^ B n+1 L 2 2 e ^ B n L 2 2 )+4 α 1 τ e B n+1 W 2 4τ( κ( θ ^ ( t n )κ( θ ^ n ) )curlB( t n+1 ),curl e B n+1 )   +4τ( curl( e u n+1 × B ^ ( t n ) ), e B n+1 )+4τ( curl( u n+1 × e ^ B n ), e B n+1 )+4τ( R B n+1 , e B n+1 ) := i=1 2 J i n+1 ( e B n+1 ) +4τ( curl( e u n+1 × B ^ ( t n ) ), e B n+1 )+4τ( curl( u n+1 × e ^ B n ), e B n+1 ) (4.35)

对上式右端各项进行精细估计。关键项的处理如下:

1) 非线性项,采用类似 I 1 n+1 ( e u n+1 ) 的证明方法,可以得到:

J 1 n+1 ( e B n+1 )=4τ( κ( θ ^ ( t n )κ( θ ^ n ) )curlB( t n+1 ),curl e B n+1 ) 4βτ e ^ θ n L 2 curlu( t n+1 ) L curl e B n+1 L 2 α 1 τ e B n+1 W 2 +Cτ e ^ θ n L 2 2 curlu( t n+1 ) L 2 α 1 τ e B n+1 W 2 +Cτ e ^ θ n L 2 2 (4.36)

2) 截断误差项,易证:

J 2 n+1 ( e B n+1 )=4τ( R B n+1 , e B n+1 ) α 1 τ e B n+1 W 2 +Cτ R B n+1 L 2 2 (4.37)

将(4.38)和(4.39)相加,我们将(4.35)式重新写为:

( e B n+1 L 2 2 e B n L 2 2 + e ^ B n+1 L 2 2 e ^ B n L 2 2 )+ α 1 τ e B n+1 W 2 Cτ e ^ θ n L 2 2 +Cτ R B n+1 L 2 2 +4τ( curl( e u n+1 × B ^ ( t n ) ), e B n+1 )+4τ( curl( u n+1 × e ^ B n ), e B n+1 ) (4.38)

将(4.34)与(4.38)式相加,我们可以得到:

( e u n+1 L 2 2 e u n L 2 2 + e ^ u n+1 L 2 2 e ^ u n L 2 2 )+( e B n+1 L 2 2 e B n L 2 2 + e ^ B n+1 L 2 2 e ^ B n L 2 2 ) + α 1 τ e u n+1 L 2 2 +2 α 1 τ e B n+1 W 2 Cτ e ^ θ n L 2 2 +Cτ e ^ u n L 2 2 +Cτ e ^ B n L 2 2 +Cτ R u n+1 L 2 2 +Cτ R B n+1 L 2 2 4τ( B ^ n ×curl e B n+1 , e u n+1 ) +4τ( curl( e u n+1 × B ^ ( t n ) ), e B n+1 )+4τ( curl( u n+1 × e ^ B n ), e B n+1 ) (4.39)

对于动量方程和磁场方程的耦合项,由 ( a×curlb,c )=( a×c,curlb ) 、Holder不等式、Young不等式以及 B=0 可知:

4τ( B ^ n ×curl e B n+1 , e u n+1 )+4τ( curl( e u n+1 × B ^ ( t n ) ), e B n+1 )+4τ( curl( u n+1 × e ^ B n ), e B n+1 ) =4τ( ( e u n+1 × e ^ B n ,curl e B n+1 )+( u n+1 × e ^ B n ,curl e B n+1 ) )=4τ( u( t n+1 )× e ^ B n ,curl e B n+1 ) 4τ u( t n+1 ) L e ^ B n L 2 curl e B n+1 L 2 α 1 τ e B n+1 W 2 +Cτ e ^ B n L 2 2 (4.40)

将耦合项的估计代入(4.39)式,可得:

( e u n+1 L 2 2 e u n L 2 2 + e ^ u n+1 L 2 2 e ^ u n L 2 2 )+( e B n+1 L 2 2 e B n L 2 2 + e ^ B n+1 L 2 2 e ^ B n L 2 2 ) + α 1 τ e u n+1 L 2 2 + α 1 τ e B n+1 W 2 Cτ e ^ θ n L 2 2 +Cτ e ^ u n L 2 2 +Cτ e ^ B n L 2 2 +Cτ R u n+1 L 2 2 +Cτ R B n+1 L 2 2 (4.41)

第二步:估计 e θ n+1

将(4.27)式与 4τ e θ n+1 作内积并在 Ω 上积分,利用(2.11)和(3.1)可以得到关于 e θ n+1 的不等式:

( e θ n+1 L 2 2 e θ n L 2 2 + e ^ θ n+1 L 2 2 e ^ θ n L 2 2 )+4 α 1 τ e θ n+1 L 2 2 4τ( ( ( σ( θ ^ ( t n ) )σ( θ ^ n ) )u( t n+1 ) ), e θ n+1 )   4τb( ( e u n+1 ),θ( t n+1 ), e θ n+1 )+4τ( R θ n+1 , e θ n+1 )+4τ( ε e u2 n+1 , e θ n+1 ) := i=1 4 L i n+1 ( e θ n+1 ) (4.42)

对上式右端各项进行精细估计。关键项的处理如下:

1) 非线性项,采用类似估计 I 1 n+1 ( e u n+1 ) 的方法,可以得到:

L 1 n+1 ( e θ n+1 )=4τ( ( υ( θ ^ ( t n ) )υ( θ ^ n ) )θ( t n+1 ), e θ n+1 ) 4βτ e ^ θ n L 2 θ( t n+1 ) L e θ n+1 L 2 3 4 α 1 τ e θ n+1 L 2 2 +Cτ e ^ θ n L 2 2 θ( t n+1 ) L 2 3 4 α 1 τ e θ n+1 L 2 2 +Cτ e ^ θ n L 2 2 (4.43)

2) 对流项,采用类似估计 I 2 n+1 ( e u n+1 ) 的方法,可以得到:

L 2 n+1 ( e θ n+1 )=4τb( ( e u n+1 ),θ( t n+1 ), e θ n+1 ) 4τ e u n+1 L 2 θ( t n+1 ) L 3 e θ n+1 L 6 3 4 α 1 τ e θ n+1 L 2 2 +Cτ u( t n+1 ) L 3 2 e u n+1 L 2 2 3 4 α 1 τ e θ n+1 L 2 2 +Cτ e u n+1 L 2 2 (4.44)

3) 洛伦兹力项,同理可得:

L 3 n+1 ( e θ n+1 )=4τb( ( e u n+1 ),θ( t n+1 ), e θ n+1 ) 4τ e u n+1 L 2 θ( t n+1 ) L 3 e θ n+1 L 6 3 4 α 1 τ e θ n+1 L 2 2 +Cτ u( t n+1 ) L 3 2 e u n+1 L 2 2 3 4 α 1 τ e θ n+1 L 2 2 +Cτ e u n+1 L 2 2 (4.45)

4) 截断误差项,易证:

L 4 n+1 ( e θ n+1 )=4τ( R θ n+1 , e θ n+1 ) 3 4 α 1 τ e θ n+1 L 2 2 +Cτ R θ n+1 L 2 2 (4.46)

将(4.43)~(4.46)式相加,我们将(4.42)式重新写为:

( e θ n+1 L 2 2 e θ n L 2 2 + e ^ θ n+1 L 2 2 e ^ θ n L 2 2 )+ α 1 τ e θ n+1 L 2 2 Cτ e ^ θ n L 2 2 +Cτ e u n+1 L 2 2 +Cτ R θ n+1 L 2 2 (4.47)

将(4.41)和(4.47)相加,我们可以得到:

( e u n+1 L 2 2 e u n L 2 2 + e ^ u n+1 L 2 2 e ^ u n L 2 2 )+( e B n+1 L 2 2 e B n L 2 2 + e ^ B n+1 L 2 2 e ^ B n L 2 2 ) +( e θ n+1 L 2 2 e θ n L 2 2 + e ^ θ n+1 L 2 2 e ^ θ n L 2 2 )+ α 1 τ e u n+1 L 2 2 + α 1 τ e B n+1 W 2 + α 1 τ e θ n+1 L 2 2 Cτ e u n+1 L 2 2 +Cτ e ^ θ n L 2 2 +Cτ e ^ u n L 2 2 +Cτ e ^ B n L 2 2 +Cτ R u n+1 L 2 2 +Cτ R B n+1 L 2 2 +Cτ R θ n+1 L 2 2 (4.48)

将(4.48)从 n=1 m( mN1 ) 求和,并由引理4.1和(4.8)式,我们有

e u m L 2 2 + e ^ u m L 2 2 + e B m L 2 2 + e ^ B m L 2 2 + e θ m L 2 2 + e ^ θ m L 2 2 + α 1 τ n=1 m ( e u n+1 L 2 2 + e B n+1 W 2 + e θ n+1 L 2 2 ) τ n=1 m ( e u n+1 L 2 2 + e ^ u n L 2 2 + e ^ B n L 2 2 + e ^ ( θ ) L 2 2 ) +C τ 4

利用离散的grownwall不等式,当 τ 充分小时,我们有:

e u m L 2 2 + e ^ u m L 2 2 + e B m L 2 2 + e ^ B m L 2 2 + e θ m L 2 2 + e ^ θ m L 2 2 + α 1 τ n=1 m ( e u n+1 L 2 2 + e B n+1 W 2 + e θ n+1 L 2 2 ) C τ 4

5. 数值结果

在本节中,我们给出数值结果以验证由定理4.1中推导的时间误差估计(4.9)。为简便起见,我们在单位正方形区域 Ω 中求MHD-Boussinesq方程(1.1)~(1.4)。此外,我们选取适当的 f g ψ 使得精确解 ( u,B,p,θ ) 具有以下形式:

u= ( y 2 e t , x 2 e t ) T , p=( 2x1 )( 2y1 ) e t , B= ( x 2 ( x1 ) 2 y( y1 )( 2y1 ) e t , y 2 ( y1 ) 2 x( x1 )( 2x1 ) e t ) T , θ=sin( πx )sin( πy )sin( t ).

粘性系数 υ( θ ) ,电导率 κ( θ ) 和热扩散率 σ( θ ) 分别由以下公式给出:

υ( θ )= e θ ,   κ( θ )= e θ ,   σ( θ )=1+θ

此外我们设定参数 ε=1 ,终止时间 T=1 。通过依次取逐渐缩小的网格尺寸 h=1/4 ,1/8 ,,1/ 64 ,并选取不同迭代次数 N=1/ h 2 ,使得 τ= h 2 。在此条件下,由公式(4.9)可得:

Table 1. Numerical errors and convergence rate of the velocity field

1. 速度场的数值误差与收敛速率

τ=h

u u n L 2

阶数

1/4

4.79E−03

-

1/8

1.24E−03

1.95

1/16

3.15E−04

1.98

1/32

7.93E−05

1.99

1/64

1.99E−05

1.99

Table 2. Magnetic field numerical error and convergence rate

2. 磁场数值误差与收敛速率

τ= h 2

B B n L 2

阶数

1/4

4.96E−04

-

1/8

1.47E−04

1.75

1/16

3.89E−05

1.92

1/32

9.86E−06

1.98

1/64

2.47E−06

1.99

Table 3. Numerical errors and convergence rate of the temperature field

3. 温度场的数值误差与收敛速率

τ= h 2

θ θ n L 2

阶数

1/4

5.57E−02

-

1/8

1.47E−02

1.93

1/16

3.72E−03

1.98

1/32

9.34E−04

1.99

1/64

2.34E−04

2.00

数值结果展示在表1~3中,这些表格给出了不同尺寸 h 对应的数值误差及其收敛速率。可以清楚地观察到,在 L 2 范数下呈现二阶收敛速率,这与理论预测的(4.23)式高度吻合。

参考文献

[1] Barleon, L., Casal, V. and Lenhart, L. (1991) MHD Flow in Liquid-Metal-Cooled Blankets. Fusion Engineering and Design, 14, 401-412. [Google Scholar] [CrossRef
[2] Hashizume, H. (2006) Numerical and Experimental Research to Solve MHD Problem in Liquid Blanket System. Fusion Engineering and Design, 81, 1431-1438. [Google Scholar] [CrossRef
[3] Smolentsev, S., Moreau, R., Bühler, L. and Mistrangelo, C. (2010) MHD Thermofluid Issues of Liquid-Metal Blankets: Phenomena and Advances. Fusion Engineering and Design, 85, 1196-1205. [Google Scholar] [CrossRef
[4] Davidson, P.A. (2001). An Introduction to Magnetohydrodynamics. Cambridge University Press.[CrossRef
[5] Titterton, D. and Weston, J.L. (2004) Strapdown Inertial Navigation Technology. Aerospace & Electronic Systems Magazine IEEE, 20, 33-34.
[6] Majda, A.J. and Bertozzi, A.L. (2001) Vorticity and Incompressible Flow. Cambridge University Press. [Google Scholar] [CrossRef
[7] Gao, H.D. and Qiu, W.F. (2019) A Linearized Energy Preserving Finite Element Method for the Dynamical Incompressible Magnetohydrodynamics Equations. Computer Methods in Applied Mechanics and Engineering, 346, 982-1001.
[8] Li, Y. and Luo, X. (2019) Second-Order Semi-Implicit Crank-Nicolson Scheme for a Coupled Magnetohydrodynamics System. Applied Numerical Mathematics, 145, 48-68. [Google Scholar] [CrossRef
[9] Zhang, Y.H., Hou, Y.R. and Shan, L. (2015) Numerical Analysis of the Crank-Nicolson Extrapolation Time Discrete Scheme for Magnetohydrodynamics Flows. Numerical Methods for Partial Differential Equations, 31, 2169-2208. [Google Scholar] [CrossRef
[10] Wang, S.H. and Li, Y. (2024) Optimal Convergent Analysis of a Linearized Euler Finite Element Scheme for the 2D Incompressible Temperature-Dependent MHD-Boussinesq Equations. Communications in Nonlinear Science and Numerical Simulation, 138, Article 108264. [Google Scholar] [CrossRef
[11] Heywood, J.G. and Rannacher, R. (1990) Finite-Element Approximation of the Nonstationary Navier-Stokes Problem. Part IV: Error Analysis for Second-Order Time Discretization. SIAM Journal on Numerical Analysis, 27, 353-384. [Google Scholar] [CrossRef
[12] Zhang, G.D., Yang, J. and Bi, C. (2018) Second Order Unconditionally Convergent and Energy Stable Linearized Scheme for MHD Equations. Advances in Computational Mathematics, 44, 505-540. [Google Scholar] [CrossRef
[13] Bian, D. and Liu, J. (2017) Initial-Boundary Value Problem to 2D Boussinesq Equations for MHD Convection with Stratification Effects. Journal of Differential Equations, 263, 8074-8101. [Google Scholar] [CrossRef