第1页 / 共10页
第2页 / 共10页
第3页 / 共10页
第4页 / 共10页
第5页 / 共10页
第6页 / 共10页
第7页 / 共10页
第8页 / 共10页
二维不可压类流函数-速度形式\磁流体力学方程组的紧致差分算法.pdf
http://www.paper.edu.cn
二维不可压类流函数 -速度形式
磁流体力学方程组的紧致差分算法
虞培祥1,2，田振夫1
1 复旦大学力学与工程科学系，上海 200433
2 复旦大学数学科学学院，上海 200433
摘要：基于二维不可压类流函数 -速度形式的磁流体力学方程组，本文分别建立了一种 2 阶 5
点常系数紧致差分格式以及一种高阶 9 点紧致差分格式。通过与具有解析解的算例进行比较，
本文提出的方法均能符合理论上的预期精度。
关键词：流体力学；紧致差分方法；完整磁流体力学模型；流函数 -速度形式；不可压
中图分类号： O361.3
Compact finite difference methods for the
quasi-streamfunction-velocity formulation of the
two-dimensional incompressible full
magnetohydrodynamic equations
Yu Pei-Xiang1,2, Tian Zhen-Fu1
1 Department of Mechanics and Engineering Science, Fudan University, Shanghai 200433
2 School of Mathematical Sciences, Fudan University, Shanghai 200433
Abstract: Based on the quasi-streamfunction-velocity formulation, a five-point constant
coefficient second-orderand a nine-point high-order compact finite difference methods are
developed for the two-dimensional steady-state incompressible full magnetohydrodynamics
equations. Both of the approaches are verified by solving the test problem with the analytical
solutions.
Key words: Fluid mechanics; compact finite difference method; full magnetohydrodynamics
equations; streamfunction-velocity formulation; incompressible
基金项目： 教育部高等学校博士学科点专项科研基金 (20100071110017)
作者简介： 虞培祥（1986-），男，博士后，主要研究方向：计算流体力学，磁流体力学。通信作者：田振夫（1963-），男，教
授，主要研究方向：计算流体力学，磁流体力学。
- 1 -
http://www.paper.edu.cn
0 引言
在我们的前期工作中，已经推导了二维不可压磁流体力学 (MHD) 的类流函数 -速度形式
控制方程组。在本文中，我们将为此方程组建立紧致差分格式。由于该方程组是从二维不可压
流函数 -速度形式 Navier-Stokes(N-S) 方程推广而来，因此可以借鉴这方面的最新研究成果。
对于流函数 -速度形式 N-S 方程，其特点是 4 阶非线性偏微分方程，因此采用传统的离散
方法，难以获得紧致差分格式，因而给计算带来了不便和难度。近年来，一些学者提出不同的
紧致差分格式来尝试解决这一问题。对于定常流动问题，Gupta 和 Kalita[1] 基于双调和方程
的离散方法 (即将流函数及其一阶导数同时作为未知函数求解)，提出了流函数 -速度方程的一
种常系数 9 点 2 阶紧致格式，但该格式系数矩阵不是对角占优，因此求解效率较低。为此，我
们 [2] 提出了一种流函数变量的 5 点 2 阶常系数紧致格式，该格式的显著特点是其关于流函数
的系数矩阵严格对角占优，且是对称正定，数值验证结果表明，该格式的收敛速度约为 Gupta
等 9 点常系数紧致差分格式的 2 倍，且构造方法简单、实用；此外，我们 [3] 还将这一算法与
四阶对流扩散型方程型紧致格式相组合，提出了适用于求解自然对流问题的新型算法。对于非
定常问题，Ben-Artzi 等 [4] 基于 Stephenson 格式建立了时间、空间均为 2 阶的纯流函数形式
N-S 方程的 9 点紧致差分格式，该方法的缺点是时间方向采用了预测 -校正方法，是条件稳定
格式，计算时间有严格限制。Kalita 和 Gupta[5] 将 Gupta 的 9 点 2 阶紧致差分格式推广到非
定常问题，得到了一种时间、空间均为 2 阶的流函数 -速度形式 N-S 方程的 9 点紧致格式，此
类格式仍存在流函数变量代数方程组系数矩阵不能对角占优缺陷，且对格式稳定性分析明显存
在缺陷。对于一般曲线坐标系，Pandit[6] 基于流函数方程提出了一种 2 阶紧致差分格式。然
而，在该格式中，流函数的一阶导数的差分格式采用了 2 阶中心差分格式，虽然最终的格式也
是 2 阶，但在整体精度和迭代效率上都难及 Gupta 在直角坐标系下的格式。最近, 我们 [7] 基
于一阶导数 4 阶中心差分格式，借鉴 5 点 2 阶常系数紧致格式的构造方法，在极坐标下构造
了一种 5 点的紧致格式，且格式系数仅依赖于坐标架，与流动参量 Re 无关，因此适用于求解
中等 Re 数的流动问题。对于高精度格式的研究，由于纯流函数或流函数 -速度形式 N-S 方程
是 4 阶偏微分方程的缘故，构造 4 阶 9 点紧致格式需要很强的技巧，2010 年，Ben-Artzi 等
[8] 首先提出了一种 4 阶 9 点紧致格式。
在本文中，我们首先将我们之前提出的 5 点 2 阶常系数格式推广到不可压 MHD 方程组
中。在此基础上，我们还将重点关注 Ben-Artzi 等的 4 阶 9 点紧致格式，并将其也应用到不可
压 MHD 问题中。
- 2 -
1 控制方程组
http://www.paper.edu.cn
在前期工作中，我们已经导出了如下的类流函数 -速度形式的无量纲定常不可压 MHD 方
)
(
[
(
u
[
[
)
(
∂4ψ
∂y4 = Re
Bx
∂y4 = Rem
∂4ϕ
∂2By
∂x2 +
∂2
∂x2
∂2v
∂y2
∂2v
∂x2 +
(
∂2By
+ By
∂y2
uBy vBx
+ v
(
)
∂2Bx
∂x2 +
∂2
∂y2
+
∂2u
∂x2 +
∂2u
∂y2
(
∂2Bx
∂y2
uBy vBx
+ f
)]
)]
)]
程组：
∂4ψ
∂x4 + 2
∂4ψ
∂x2∂y2 +
Ha2
Rem
∂4ϕ
∂x2∂y2 +
∂4ϕ
∂x4 + 2
∂ψ
u =
∂y
v = ∂ψ
∂x
∂ϕ
∂y
Bx =
By = ∂ϕ
∂x
(1)
(2)
(3)
(4)
(5)
(6)
(7)
其中，ψ 代表流函数，u 和 v 为速度分量，ϕ 为磁流函数，Bx 和 By 分别是磁场强度的分量。
Re 为雷诺数，Rem 为磁雷诺数，Ha 为哈特曼数。f 为额外的外力项。
需要特别指出，如果在边界上始终满足电流密度 J = Rem(uByvBx)，即 ( ∂2ϕ
∂x2 + ∂2ϕ
∂y2 ) =
Rem(uBy vBx)，则方程 2可简化为:
∂2ϕ
∂y2 = Rem
∂2ϕ
∂x2 +
(
∂ψ
∂x
∂ϕ
∂y
∂ψ
∂y
∂ϕ
∂x
)
这一形式与文献 [9] 中曾提到的控制方程完全相同，由于在该文献中磁流函数用 A 表示，因此
又称其为 ψ A 方法。
在本文中，我们将研究一般的类流函数 -速度形式的不可压 MHD 方程组 (1-6)。
2 2 阶紧致差分方法
在本节中，我们将推导类流函数 -速度形式控制方程组的 2 阶紧致差分格式。方便起见，
我们在 x 和 y 方向采用均匀网格尺寸 h。
在 [2] 中, 我们已经为不可压流函数 -速度形式建立了一种 2 阶 5 点常系数的紧致差分格
)
(
12
h2
式，即对于方程 (1)，当 Ha 0 时，有如下的紧致差分格式：
yv δ2
δ2
xψ + δ2
yψ
=
xδyu
(
12
h2 (δxv δyu) +
Re
xv + δ2
u(δ2
(
δxδ2
yv) v(δ2
1
2
)
)
xu + δ2
yu)
+ f + O(h2)
(8)
- 3 -
其中，所有算子（如 δx）均为标准的二阶中心差分算子，速度分量的离散采用如下的四阶紧致
差分格式：
http://www.paper.edu.cn
(1 +
(1 +
h2
6
h2
6
δ2
y)u = δyψ
x)v = δxψ
δ2
对于不可压 MHD 问题, 新增的右端项可表示成 ^f
[
(
)
(
+ By
)]
∂2Bx
∂y2
Bx
∂2By
∂x2 +
∂2By
∂y2
∂2Bx
∂x2 +
注意到该部分与流体的对流项非常相似，因此可采用如下的二阶精度的离散方法：
^f = Ha2
Rem
(
^f =
Ha2
Rem
)
)
)
Bx(δ2
xBy + δ2
xBx + δ2
yBx)
+ O(h2)
(
)
将式 (12) 代入式 (8), 即可得方程 (1) 的二阶精度格式：
12
yv δ2
h2
δ2
xψ + δ2
yψ
=
xδyu
(
yBy) By(δ2
(
δxδ2
yv) v(δ2
1
2
12
h2 (δxv δyu) +
Re
xv + δ2
u(δ2
(
对于磁流函数 ϕ 方程 (2)，类似可得
+
Ha2
Rem
)
δ2
xϕ + δ2
yϕ
=
(
12
h2
Bx(δ2
xBy + δ2
xu + δ2
yu)
yBy) By(δ2
(
(
h2 (δxBy δyBx) +
δxδ2
uBy vBx
Rem(δ2
x + δ2
y)
1
2
12
xBx + δ2
yBx)
)
yBy δ2
xδyBx
+ O(h2)
)
+ f + O(h2) (13)
)
(9)
(10)
(11)
(12)
(14)
(15)
(16)
其中，磁场强度 Bx 和 By 也采用与速度相同的离散方法：
(1 +
(1 +
h2
6
h2
6
δ2
y)Bx = δyϕ
x)By = δxϕ
δ2
综合式 (9)、(10) 和 (13)-(16)，我们就得到了类流函数 -速度形式定常不可压 MHD 方程的二
阶紧致差分格式.
3 高阶紧致差分方法
在本节中，我们将参考 Ben-Artzi 等 [8] 提出的 4 阶紧致差分格式，继续推导类流函数 -速
度形式控制方程组的 4 阶紧致差分格式。方便起见，我们同样在 x 和 y 方向采用均匀网格尺
寸 h。
- 4 -
http://www.paper.edu.cn
在 [8] 中, 流函数方程的双调和项可离散为：
( δ2
)
∂4ψ
∂x4 + 2
∂4ψ
∂x2∂y2 +
∂4ψ
∂y4 =
12
h2
+6δ2
xψ δ2
yψ + 2δxδ2
yψ + δyu δxv
yv 2δ2
xδ2
xδyu + O(h4)
(17)
其中，速度分量同样采用式 (9) 和 (10) 来计算.
当完全没有外力项（包括磁力项）时，方程 (1) 的右端项代表流体流动的对流项，在 [8]
中，作者提出了如下高阶差分格式：
(
(
[
[
(
u
∂3ψ
∂x3 +
3
δ2
x^v +
2
δ2
xδyψ + δ2
u
)
(
)
(
x ^u + δxδy ^v 3
∂3ψ
∂x∂y2
15
h2
δxψ + ^v
+ v
Re
= Re
+v
)]
∂3ψ
∂x2∂y
+ δxδ2
+
∂3ψ
∂y3
)
))]
(
y ^v δxδy ^u
yψ δ2
δyψ ^u
δ2
y ^u +
2
15
h2
+ O(h4)
(18)
与之前的速度近似值不同，此处的速度近似值 ^u 和 ^v 需要采用高于 4 阶的差分格式. 在 [8] 中,
作者给出了一种五点六阶的 Padé 内点格式，而在次边界上，采用单边的五阶格式。对于速度
项 ^v（2 i N 2, 1 j N 1），其具体格式为
10ψ0,j + 9ψ1,j 18ψ2,j ψ3,j
30h
ψN3,j + 18ψN2,j 9ψN1,j 10ψN,j
^v1,j +
3
10
^v2,j =
^vN1,j +
^vN,j =
3
10
^v0,j +
6
10
^vN2,j +
1
10
3
6
10
10
1
^vi1,j + ^vi,j +
3
ψi+1,j ψi1,j
对于速度项 ^u（1 i N 1, 2 j N 2），其格式为
^vi+1,j = 14
1
3
2h
9
30h
ψi+2,j ψi2,j
4h
1
9
^ui,2 = 10ψi,0 + 9ψi,1 18ψi,2 ψi,3
3
10
^ui,1 +
^ui,0 +
6
10
^ui,N2 +
1
10
3
6
10
10
1
^ui,j1 + ^ui,j +
3
^ui,N1 +
3
10
^ui,j+1 =
1
3
[
(
30h
+
2h
1
9
30h
14
9
ψi,j+2 ψi,j2
^ui,N = ψi,N3 + 18ψi,N2 9ψi,N1 10ψi,N
ψi,j+1 ψi,j1
)
(
15
h2
^By 3
^By δxδy
))]
+ O(h4)
^Bx +
)
^Bx
4h
对于 MHD 问题，类似地，我们可以得到
(
^Bx
(
yϕ δ2
δyϕ ^Bx
15
h2
其中， ^By（2 i N 2, 1 j N 1）采用如下格式：
^f = Ha2
Rem
xδyϕ + δ2
δ2
x
3
2
^Bx + δxδy
δxϕ + ^By
+ δxδ2
^By +
+ ^By
δ2
x
δ2
y
2
y
10ϕ0,j + 9ϕ1,j 18ϕ2,j ϕ3,j
3
10
^By 2,j =
^By 0,j +
^By 1,j +
6
10
6
10
^By N2,j +
1
10
3
10
1
^By i1,j + ^By i,j +
3
1
3
^By N1,j +
^By N,j =
3
10
^By i+1,j = 14
9
30h
ϕN3,j + 18ϕN2,j 9ϕN1,j 10ϕN,j
ϕi+1,j ϕi1,j
2h
30h
ϕi+2,j ϕi2,j
4h
1
9
(19)
(20)
(21)
(22)
- 5 -
http://www.paper.edu.cn
^Bx（1 i N 1, 2 j N 2）采用如下格式：
^Bxi,2 = 10ϕi,0 + 9ϕi,1 18ϕi,2 ϕi,3
3
10
^Bxi,N1 +
3
10
^Bxi,j+1 =
14
9
^Bxi,N = ϕi,N3 + 18ϕi,N2 9ϕi,N1 10ϕi,N
30h
ϕi,j+1 ϕi,j1
30h
ϕi,j+2 ϕi,j2
+
1
9
2h
4h
(23)
将式 (17), (18) 和 (23) 代入 Eq. (1), 我们就得到了该方程的四阶紧致差分格式。
1
3
^Bxi,0 +
^Bxi,1 +
6
10
6
10
^Bxi,N2 +
1
10
3
10
1
^Bxi,j1 + ^Bxi,j +
3
( δ2
(
xψ δ2
12
yψ
h2
[
(
= 12
δyu δxv
h2
(
[
δ2
xδyψ + δ2
+^v
Ha2
Rem
δ2
x^v +
(
+Re
^Bx
3
2
3
2
δ2
x
^u
(
yψ
xδ2
+ 6δ2
)
) 2δxδ2
(
)
δxψ + ^v
(
x ^u + δxδy ^v 3
15
h2
2
^By +
15
h2
+ ^By
δ2
xδyA + δ2
x
^Bx + δxδy
(25)
(26)
(27)
(28)
(
yψ δ2
)
15
h2
)
))]
y ^v δxδy ^u
δyψ ^u
(
yA δ2
y
+ δxδ2
δ2
y
^Bx +
15
h2
δyA ^Bx
^By δxδy
))]
yv + 2δ2
xδyu
+ δxδ2
δ2
y ^u +
δxA + ^By
^By 3
2
( δ2
)
^Bx
)
+ f + O(h4)
(24)
对于磁流函数方程 (2)，同样可被视作特殊的双调和方程，其左端项可采用如下离散：
∂4ϕ
∂x4 + 2
∂4ϕ
∂y4 =
∂4ϕ
∂x2∂y2 +
yϕ + δyBx δxBy
yBy 2δ2
)
(
对于右端项, 我们将项 (R = uBy vBx) 作为一个整体进行考虑，即：
xϕ δ2
yϕ + 2δxδ2
12
h2
+6δ2
(
)
xδ2
Rem
∂2
∂x2 +
∂2
∂y2
(uBy vBx) = Rem
∂2R
∂x2 +
∂2R
∂y2
xδyBx + O(h4)
为了得到右端项的高精度格式，我们需要获取 R 的二阶偏导数的高阶格式。参考 [10]，我
们可获得如下关于隐式格式：
^S0,j + 11 ^S1,j =
13R0,j 27R1,j + 15R2,j R3,j
13RN,j 27RN1,j + 15RN2,j RN3,j
h2
11 ^SN1,j + ^SN,j =
^Si1,j + 10 ^Si,j + ^Si+1,j =
h2
12
h2 (Ri1,j 2Ri,j + Ri+1,j)
以及
Si,0 + 11 Si,1 =
13Ri,0 27Ri,1 + 15Ri,2 Ri,3
13Ri,N 27Ri,N1 + 15Ri,N2 Ri,N3
h2
11 Si,N1 + Si,N =
Si,j1 + 10 Si,j + ^Si,j+1 =
h2
12
h2 (Ri,j1 2Ri,j + Ri,j+1)
- 6 -
http://www.paper.edu.cn
(29)
(30)
其中， ^S 和 S 分别代表 ∂2R
∂y2 的近似值. 于是，式 (26) 可表示成：by
(uBy vBx) = Rem( ^S + S) + O(h4)
∂2
∂y2
将式 (25) 和 (29) 代入 (2) 中, 就可以获得该方程的四阶紧致差分格式：
)
∂2
∂x2 +
∂x2 和 ∂2R
(
( δ2
)
(
xA δ2
yA
δyBx δxBy
xδ2
+ 6δ2
) 2δxδ2
yA
Rem
12
h2
= 12
h2
yBy + 2δ2
xδyBx Rem( ^S + S) + O(h4)
至此，我们已经建立了类流函数 -速度形式不可压 MHD 方程的四阶紧致差分格式，包括
式 (9)、(10)、(15)、(16)、(19)、(20)、(22)、(23)、(24)、(27)、(28) 和 (30)。尽管在上述离
散过程中，我们需要计算速度分量和磁场强度分量的两种不同精度的近似值，但由于这些近似
值的差分格式构成的代数方程组均为三对角方程，因此我们可以采用求解三对角方程组的专门
方法快速计算这些结果，由此并不会增加明显的计算量。
4 算法描述
在本文中，我们采用内外迭代的方法 [11] 来获取文中提到的算法的计算结果。为简便起
见，我们只描述四阶算法的具体求解流程。实际上，我们注意到在所有的离散方程组中，式
(24) 与 (30) 最复杂，其所花费的计算时间也最多，因此我们将采用多重网格方法 [12, 13] 来
加速这两个代数方程组的计算。
x, ^Bn
y 均为已知，则式 (9)、(10)、(15)、(16)、(19)、
假定 ψn, ϕn, ^un, ^vn, ^Bn
y , un, vn, Bn
x Bn
(20)、(22)、(23)、(24)、(27)、(28) 和 (30) 的数值解可以通过以下流程来获得：
1. 利用多重网格法求解式 (24) 来获取 ψn+1，其中右端项均采用前一迭代步的结果 (如 un、
vn、^un 和 ^vn 等)。
2. 利用多重网格法求解式 (30) 来获得 ϕn+1，其右端项的计算采用如下方法。
(a) 计算 R = uBy vBx，其中速度与磁场强度均采用高精度近似值的结果；
(b) 根据式 (27) 和 (28)，求解 R 的二阶偏导数，从而得到式 (30) 中的 Rem( ^S + S)；
(c) 式 (30) 的其余项采用前一迭代步的结果。
3. 根据式 (9) 和 (10)，更新速度的四阶近似值 un+1 and vn+1，并根据式 (19) 和 (20)，更
新速度的更高阶近似值 ^un+1 和 ^vn+1。
4. 类似地，更新磁场分量的四阶近似值和更高阶近似值。
5. 从 n = 0、1、2、... 开始，重复上述迭代过程 (1) 到 (4) 直到相邻两次迭代步的最大误差
小于收敛判据 ε. 在本文中, 选取 ε = 2 10
14.
- 7 -
http://www.paper.edu.cn
5 算法验证
为了验证我们提出的算法的准确性，我们考察一个含解析解的定常问题. 对于定常不可压
类流函数 -速度形式 MHD 方程组 (1)-(6), 可构造如下的解析解:
ψ = ex+yRe,
A = (x + y)2, Bx = 2(x + y), By = 2(x + y)
u = ex+yRe,
v = ex+yRe
于是可得到外力项 f 为
f = 4ex+y
(31)
(32)
65 65
表 1: 各个变量在不同计算网格下的 L2 误差及收敛阶
17 17
收敛阶
ψ
u
v
A
Bx
By
ψ
u
v
A
Bx
By
2.687E-06
9.038E-06
9.091E-06
1.592E-06
3.646E-06
3.646E-06
2.125E-08
2.498E-07
2.491E-07
1.805E-08
4.770E-08
4.770E-08
1.97
1.97
1.97
1.96
1.98
1.98
3.91
3.93
3.92
4.16
4.06
4.06
33 33
2 阶格式
6.838E-07
2.311E-06
2.324E-06
4.077E-07
9.215E-07
9.215E-07
4 阶格式
1.413E-09
1.642E-08
1.641E-08
1.009E-09
2.864E-09
2.864E-09
收敛阶
1.98
1.98
1.98
1.98
1.98
1.98
4.12
3.96
3.96
3.67
3.84
3.84
1.732E-07
5.863E-07
5.893E-07
1.036E-07
2.333E-07
2.333E-07
8.139E-11
1.058E-09
1.057E-09
7.933E-11
1.997E-10
1.997E-10
在 [0, 1] [0, 1] 的计算域内，我们对第一类边界条件下的上述问题进行了数值求解，其中
我们选取 Re = 1、Ha = 5 和 Rem = 100。表 1给出了各个变量在不同计算网格下的 L2 误差，
并由此可以得到计算得到的收敛阶。数据表明，二阶格式的计算收敛阶基本在 2 附近，而四阶
格式的收敛阶在 4 左右，这些结果都符合我们的理论预期，由此表明本算法是准确和可靠的。
6 结论
在本文中，我们分别推广了文献 [2] 和 [8] 中的关于流函数 -速度形式的不可压 Navier-
Stokes 方程的 2 阶 5 点常系数紧致格式和 9 点高精度紧致差分格式，从而分别提出了基于不
可压流函数 -速度形式定常 MHD 方程的二阶紧致差分格式和四阶紧致差分格式。数值实验表
- 8 -
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
