1. 引言
不可压缩磁流体力学(MHD)方程组是由流体力学的Navier-stokes方程与电磁学中的Maxwell方程耦合而成,描述了在磁场存在且没有表面张力的情况下,不可压缩导电流体的运动规律,在工业上有广泛的应用,如在天体物理学和地球物理学[1]、核反应堆中的液态金属冷却[2] [3]、金属的电磁铸件[4]等领域。设
为
中一个具有连通边界
的有界域,考虑以下随时间变化的不可压缩磁流体动力学(MHD)方程组:
(1)
(2)
(3)
(4)
(5)
其中
,
是一个正常数。u、p、B、J、E分别表示速度、压力、磁场、电流密度和电场,
表示雷诺数,
为电导率,
为磁导率,f为给定的外力。
从物理的角度来看,无散度条件div B = 0对于不可压缩MHD方程组是重要的,它保证了没有磁单极子的存在。但是,这一条件的小扰动可能会导致不可压缩MHD方程的数值模拟出现较大的误差。近年来,许多研究人员研究满足MHD系统的无散度方案。Hu等人在[5]中,基于不可压缩MHD系统的磁电公式,提出了一种混合有限元方法,用Nédélec棱元逼近电场E,用Raviart-Thomas元逼近磁场B,该全离散格式可以保持磁场的无散度条件。我们知道对于任何
,都有div(curl w)≡0。因此,其中一种构造满足无散度条件的策略是将磁场B改写为一个旋转场。通过这种方法,Mao等人在[6]中引入具有时间规范的磁矢量势A,使其满足
和
。
将上述关系带入(1)~(5)中,可推出如下的矢量势不可压缩MHD方程组[6]:
(6)
(7)
(8)
其中
为磁雷诺数。
对该模型,可考虑如下的初始条件和边界条件[6]:
(9)
对于上述矢量势MHD方程组的数值方法及其收敛性已有一些研究工作。Mao在[6]中构造了一种线性化的Crank-Nicolson有限元格式,通过离散的能量不等式证明了数值解的弱收敛性。Li等人在[7]中,证明了非光滑和非凸域二维问题的
范数的强收敛性。在[8]、[9]和[10]中分别研究了矢量势MHD方程组的一阶欧拉、二阶Crank-Nicolson和BDF有限元格式,并建立了速度场u和磁矢量势A的误差估计。因此,由于
,离散磁感应
可以精确地保持无散度条件。在MHD系统中,另一个约束条件是速度场的无散度条件,即
。基于Navier-Stokes方程的投影算法[11] [12],对于常用的MHD系统,Prohl在[13]中研究了一种解耦的投影格式,其中速度场u、压力p和磁感应B的计算都是完全解耦的,但该解耦格式,不能满足离散的能量不等式,使得[13]中的解耦投影方案不是无条件稳定的。在[14]和[15]中作者提出了无条件稳定的时间半离散和有限元全离散投影方案,分别得到了时间误差估计和时空误差估计。目前,还没有关于求解矢量势MHD方程组的投影算法研究。
本文基于[11] [12]中投影方法的思想,构造求解矢量势MHD方程组的一阶投影有限元格式,使得原MHD方程组中原速度场和磁场的数值逼近解在离散层面上均满足无散度条件。在格式构造上,通过半隐方法来处理非线性项,使得所构造的投影格式是无条件稳定的线性化格式,易于实现程序设计。理论上,在合理的正则性假设下,我们得到了速度和磁矢量势在L2-范数和H(curl)-范数下的一阶时间误差估计。最后通过数值算例验证了所得的收敛性结果。
2. 预备知识
对于
,
,记
为Sobolev空间,当p = 2时,我们用
表示
。
中的范数由
表示,用
表示
的内积,
中的范数用
表示。符号
来表示一个正常数,它与时间步长大小
无关。
引入函数空间:
其中W空间的范数可定义为:
定义如下双线性项:
和三线性项:
显然该三线性项满足:
令
且初值满足
,问题(6)~(9)的变分形式为:求解
和
使得
(10)
(11)
接下来我们用
代表从
到H的正交投影算子,那么可以定义Stokes算子
如下:
从Poincare’不等式,插值不等式和Sobolev嵌入定理可以得到以下不等式:
(12)
(13)
其中
和
。
3. 一阶投影时间离散格式和主要结果
在本节中,我们给出问题(10)~(11)的一阶投影时间离散格式,并证明数值格式的无条件稳定性和本文的主要结果。
3.1. 投影时间离散格式
设
为时间区间
的均匀划分,时间步长为
和
,
。对于任何函数序列
,定义
,其中
。
定义
和任意
。下面给出求解问题(10)~(11)的一阶投影时间离散格式:
对于
,分别求解
和
使得
(14)
(15)
其中边界条件为
和一个
,和
(16)
其中边界条件为
对于上述(16)式进行变形可得:
(17)
3.2. 稳定性分析
下面通过两个引理证明该投影格式的无条件稳定性。
引理1 投影格式(14)~(16)有唯一解
。
证明:对于任何
,问题(14)~(15)可以重写为
令,易得
是
上的一个范数。定义右端为
,可以得到
另一方面,通过柯西–施瓦茨不等式,
的连续性是明显的,如下:
从而根据Lax-Milgram定理,引理1得证。
引理2 对于
,设
和
为(14)~(16)的解。对于任意的
,有如下的能量不等式
证明:让(14)式和(15)式分别与
和
做内积,两式相加得
(18)
这里我们用到
和。
接着让(16)式分别与
和
做内积,两式相加得
(19)
将(18)和(19)两式相加得
(20)
根据
和
,可得
将上述方程带入(20)式并从
到
求和,我们可以得到引理2的证明。
3.3. 主要结果
为了方便误差估计,记:
为了给出本文的主要结果,我们对初始值和解的正则性做以下假设:
假设1 假设方程(14)~(16)存在唯一的局部强解,满足如下正则性:
(21)
(22)
本文的主要结果为:
定理1 假设矢量势MHD方程组(14)~(16)的解满足假设1的正则性条件,则存在正常数C > 0使得
(23)
其中
是不依赖于
的正常数。
4. 时间误差估计
在这一部分,我们要给出定理1的证明。首先引入离散的Gronwall不等式[16]:
引理3 当整k ≥ 0时,设
和B为非负数,使得
(24)
假设
和
,则有
(25)
注:若(24)式中右边的第一项和仅加到
,那么(25)式对所有的
及
成立。
4.1. 误差方程
为了完成定理1的证明,首先给出相关的误差方程。
对于
,在(10)~(11)中令
,可得
(26)
和
(27)
其中截断误差函数为
(26)式、(27)式分别减去(14)式和(15)式,得到误差方程为
(28)
和
(29)
根据(17)有
(30)
为了方便误差估计,方程(28)~(29)中的一些项可重写为
和
4.2. 定理1的证明
首先,通过Hölder不等式、正则性假设和泰勒公式,我们很容易得到截断误差函数的误差估计。
引理3 若假设1成立,那么我们有如下估计
(31)
接下来证明定理1。(28)和(29)分别与
和
做内积,两式相加,我们可以得到
(32)
然后对右端逐项进行估计:
通过Hölder不等式、(12)~(13)、Young不等式和(21)~(22),估计
项有
其中,
。
估计
项。
类似的,我们可以估计
、
、
和
。
接下来估计
、
和
。考虑方程左侧有,故将这三项中
与
捆绑在一起,令复杂的四线性项转化为更易分配的三线性项。
其中,
其中,
结合上述不等式并带入(32)式,有
(33)
下面对
进行相关估计,令(30)分别与
和
做内积,可得
其中我们用到
在
上和
在
,接着将两式相加得
(34)
将(33)式与(34)式相加可得
(35)
根据(30),我们可以得到
进一步有
将上述不等式带入(35)式,可得
(36)
其中,我们使用
将(36)式从
到
求和,根据离散的Gronwall不等式,存在一些正常数C > 0使得
(37)
其中,我们使用引理3和下面的估计:
现在,我们给出
的估计。通过
和
可得
其中,我们使用
最后,将投影算子
作用于(30)式,使用(13)和(37),完成定理1的证明。
5. 数值结果
下面将给出数值实验来验证理论推导的最优误差估计。为了简便,考虑在单位正方形[0, 1] × [0, 1]上的问题(6)~(9)。在计算过程中,通过取适当的
和
使得真解有如下形式:
在数值模拟中,我们使用有限元方法进行空间离散,引入有限元空间
以及
,并选取Mini元逼近速度和压力,用棱元近似磁矢量势A。此外,我们取
以及
。
Table 1. Numerical error and convergence order under L2
表1. L2范数下的数值误差和收敛阶数
N |
|
阶数 |
|
阶数 |
10 |
7.56E−02 |
− |
1.74E−01 |
− |
20 |
4.44E−02 |
0.77 |
8.41E−02 |
1.05 |
30 |
3.12E−02 |
0.87 |
5.53E−02 |
1.03 |
40 |
2.40E−02 |
0.91 |
4.12E−02 |
1.03 |
60 |
1.64E−02 |
0.94 |
2.72E−02 |
1.02 |
80 |
1.25E−02 |
0.95 |
2.03E−02 |
1.02 |
Table 2. Numerical error and convergence order under L2
表2. L2范数下的数值误差和收敛阶数
N |
|
阶数 |
|
阶数 |
10 |
5.92E−01 |
− |
5.33E−01 |
− |
20 |
3.25E−01 |
0.87 |
2.44E−01 |
1.13 |
30 |
2.26E−01 |
0.89 |
1.57E−01 |
1.08 |
40 |
1.69E−01 |
1.01 |
1.16E−01 |
1.06 |
60 |
1.05E−01 |
1.17 |
7.59E−02 |
1.04 |
80 |
7.25E−02 |
1.30 |
5.64E−02 |
1.03 |
其中,表2中
为了验证一阶时间收敛阶,我们固定网格尺寸
,确保空间误差相对于时间误差可以忽略不计。我们取时间步长
,其中N为时间迭代步数。对于定理1中给出的时间误差估计,分别取
,得到的数值误差及收敛阶如表1和表2所示。从表1和表2可以看出,当时间步长逐渐变小时,得到的一阶收敛阶与定理1中的理论相一致。