一种求解三维非稳态对流扩散反应方程的高精度有限差分格式*

魏剑英， 葛永斌

（宁夏大学 数学统计学院，银川 750021）

引 言

非稳态对流扩散反应方程是一类基本的发展方程，在生态环境、流体力学、生物数学等领域都有广泛的应用.由于解析解通常很难求出，因此研究各种有效实用的数值方法对此类方程进行求解就显得极为重要.

有限差分方法作为经典实用的数值方法之一，在计算流体力学中发挥着重要作用.目前关于非稳态对流扩散方程已经有非常多的研究报道，如文献[1-11].而关于非稳态对流扩散反应方程的报道相对还较少，并且大多是针对一维和二维问题的.如崔翔鹏等[12]针对一维问题，在单调迭代的差分方法上进行改进，构造了一种预估-校正差分方法，该方法的截断误差为O(τ2+h4)；Karaa[13]提出了求解二维问题的高精度局部一维格式，其截断误差为O(τ2+h4)，且是无条件稳定的；赵秉新[14]将一维方程利用指数变换消去反应项，采用高阶迎风紧致差分公式对对流项进行离散，采用高阶对称差分公式对扩散项进行离散，时间方向采用四阶Runge-Kutta 方法计算，所提格式的截断误差为O(τ4+h4)，且该格式适用于对流占优问题的计算；杨录峰等[15]对一维方程消去对流项后，利用四阶Padé公式，构造了一种无条件稳定的差分格式，其截断误差为O(τ2+h4)；杨晓佳等[16]针对一维非稳态方程，采用Taylor 级数展开式结合余项修正法，推导出一种截断误差为O(τ2+τh2+h4)的隐式差分格式，但该格式是条件稳定的.已有的文献中关于三维非稳态对流扩散反应方程的研究则更加少见.Kaya[17]利用添加一些特殊网格点的处理技巧，结合算子分裂技术，在时间方向用Crank-Nicolson(C-N)格式或向后差分公式离散，构造了求解三维对流扩散反应方程的二阶精度的隐格式，且该格式适合于对流和反应占优问题的计算；张亚刚[18]针对三维非稳态方程，空间利用四阶紧致差分公式，时间利用四阶向后差分公式，提出了一种截断误差为O(τ4+h4)的无条件稳定的紧致差分格式，但该格式为五层格式，所以除初始条件外，还需分别构造第一、第二、第三时间层上的格式，才能使时间向下推进.综上所述，关于三维非稳态对流扩散反应方程的高阶紧致有限差分格式的研究还很少见，有待于进一步深入研究.为此，本文将针对含有变对流项和反应项系数的三 维非稳态对流扩散反应方程，建立一种高精度紧致的有限差分格式.

1 差分格式的建立

考虑三维非稳态对流扩散反应方程：

初始条件为

边界条件为

然后，利用向前差分离散上式右端的时间导数项，利用式(5)计算空间方向的一阶导数项，利用式(6)离散二阶导数项，可以得到

再将式(10)代入式(8)，得到

最后把式(11)代入式(7)，舍去高阶项，得到

类似地，可以给出uy和uz在边界上的计算格式.

2 稳定性分析

下面利用Fourier 稳定性分析法对格式(12)的稳定性进行分析.在方程(1)中，假设对流项和反应项系数均为常数，分别为且f项精确成立，则格式(12)变为

在网格点(xi,yj,zk,tn)处，令

将式(16)、(17)分别代入式(15)，进行化简，得到

根据式(5)有

3 算 法

4 数值实验

本文将利用以下三个问题进行数值实验，验证上述高精度紧致差分格式(12)的稳定性和有效性，文中涉及到的误差及收敛阶定义如下:

右端项f(x,y,z,t)和 初边值条件均由精确解u(x,y,z,t)=etsin(πx)sin(πy)sin(πz)给出.

问题1 是含有非齐次项的三维非稳态问题，本文将对对流项系数取两种不同函数的情况进行计算：①p=tx，q=ty，r=tz； ②p=t(x+y+z)，q=t2xyz，r=sin[π(x+y+z+t)].对于上述两种情况，本文格式的计算结果均符合理论分析.表1 给出了当 τ=h2，t=0．25时 的L∞范 数误差、L2范数误差以及收敛阶的比较.可以看出，C-N 格式和BTCS 格式在空间上达到的精度是二阶，而本文格式在空间上达到的精度是四阶.表2 给出了当h=0．031 25,t=0．5时 ，不同时间步长τ 下 的L∞范 数误差、L2范数误差、收敛阶的比较及计算时间.从表中可以看出，本文方法的两种误差在时间上均达到了二阶精度，当时间步长较小时计算收敛慢，所用的计算时间相对较长.表3 给出了当t=1， 网格数取 3 2×32×32， 网格比 λ 分别取0.8，1.6，3.2，6.4 时，该问题的L∞范 数误差和L2范数误差.计算结果表明，对于给定的所有λ 值，本文格式同C-N 格式和BTCS 格式一样，仍保持收敛，即本文格式是无条件稳定的，这与理论分析的结果一致.

表1 问题1 当τ =h2,t=0．25 时 的 L∞误 差和 L2误差及收敛阶Table 1 L∞ errors, L2 errors and convergence rates for τ =h2, t=0．25 in problem 1

表2 问题1 当h =0．031 25,t=0．5时 ，本文格式的 L∞误 差、L 2误差、收敛阶及CPU 时间Table 2 L∞ errors, L2 errors, convergence rates and CPU time for h =0．031 25, t=0．5 in problem 1

表3 问题1 当h =1/32,t=1时 ，对不同网格比λ =τ/h2的 L∞误 差和 L2误差Table 3 L∞ errors and L2 errors for h =1/32,t=1 for different mesh ratio λ =τ/h2 values in problem 1

问题2考虑如下Gauss 脉冲问题：

问题2 是常系数问题，表4 给出了当 h=0．025时 对不同的 t,p,q,r,τ，本文格式与几种交替方向隐格式(ADI)的 L∞范 数误差和 L2范数误差的比较.其中 Douglas-Gunn ADI格式[1]为二阶精度无条件稳定的格式，HOC-ADI格 式[3]和 PHOC-ADI格式[4]均为时间二阶、空间四阶精度格式，HOC-ADI 格式是条件稳定的，而PHOC-ADI 格式是无条件稳定的，EHOC-ADI 格式[5]是时间二阶、空间四阶精度的无条件稳定的指数型格式，RHOC-ADI 格式[6]是时间二阶、空间四阶精度的无条件稳定的有理型格式.从表中可以看出，当p=q=r=0．8 时 ，本文格式的 L∞范 数误差和 L2范数误差比 D ouglas-Gunn ADI格式、EHOC-ADI 格式和RHOCADI 格式的都小，计算结果更精确；随着 p,q,r的增大，本文格式的计算误差比前四种ADI 格式的计算误差都小，但比RHOC-ADI 格式的计算误差稍大一些；当 p,q,r增大到800 和8000 时，HOC-ADI 格式的计算结果发散.表5 给出了当t =0．25,τ=h2,p=q=r=0．8时 ，本文格式的 L∞范 数误差和 L2范数误差及收敛阶，表中结果表明本文格式在空间上达到了四阶精度.图1和图2 给出了问题2 在计算区域1 ．2≤x,y≤1．8内 ，截面 z =1．5上对不同的 t ,p,q,r,τ，不同格式的计算解和解析解的等值线，其中虚线为解析解，实线为计算解.从图1 可以看出，当h=0．025,t=1．25, p=q=r=0．8,τ=6．25×10−3时 , P HOC-ADI 格 式、 E HOC-ADI格 式、 R HOC-ADI格式和本文格式的数值解均能与解析解吻合得很好，能准确捕获以(1.5,1.5)为中心的移动脉冲.从图2 可以看出，当h=0．025,t=1．25×10−4, p=q=r=8 000,τ=6．25×10−7时，RHOC-ADI 格式的数值解曲线基本与解析解曲线一致，本文格式的数值解与解析解之间的吻合度只产生了微弱的偏差，EHOC-ADI 格式的数值解与解析解之间产生的偏差稍微大一些，而PHOC-ADI 格式的数值解与解析解之间产生了相当大的偏差.图中结果表明，同为高精度紧致差分格式，当 h,α 固 定不变时，随着 p,q,r值的不断增大，相较于EHOC-ADI 格式和PHOC-ADI 格式，本文格式和RHOC-ADI 格式更适合求解此类波的传播问题.

图2 问题2 当h =0．025, t=1．25×10−4, p=q=r=8 000, τ=6．25×10−7 时 ，在区域1 ．2≤x,y≤1．8内，截面z=1.5 上的等值线图：(a) PHOC-ADI 格式和精确解；(b) EHOC-ADI 格式和精确解；(c) RHOC-ADI 格式和精确解；(d) 本文格式和精确解Fig.2 Numerical contour lines compared with exact solutions at z=1．5 in region 1．2≤x,y≤1．8 in problem 2 for h=0．025, t=1．25×10−4,p=q=r=8 000, τ=6．25×10−7: (a) PHOC-ADI scheme and exact solutions; (b) EHOC-ADI scheme and exact solutions; (c) RHOC-ADI scheme and exact solutions;(d) the present scheme and exact solutions

表4 问题2 当h =0．025 时，在不同参数下的 L∞误 差和 L2误差Table 4 L ∞ errors and L 2 errors for h =0．025 in problem 2

表5 问题2 当t =0．25,τ=h2,p=q=r=0．8时 ，本文格式的 L∞误 差和L 2误差及收敛阶Table 5 L∞ errors, L2 errors and convergence rates for t =0．25,τ=h2,p=q=r=0．8 in problem 2

图1 问题2 当h =0．025,t=1．25,p=q=r=0．8,τ=6．25×10−3 时，在区域 1 ．2≤x,y≤1．8内 ，截面 z =1．5上的等值线图：(a) PHOC-ADI 格式和精确解；(b) EHOC-ADI 格式和精确解；(c) RHOC-ADI 格式和精确解;(d) 本文格式和精确解Fig.1 Numerical contour lines compared with exact solutions at z=1．5 in region 1．2≤x,y≤1．8 in problem 2 for h=0．025, t=1．25,p=q=r=0．8,τ=6．25×10−3: (a) PHOC-ADI scheme and exact solutions; (b) EHOC-ADI scheme and exact solutions; (c) RHOC-ADI scheme and exact solutions; (d) the present scheme and exact solutions

问题3考虑如下非线性问题：

其中

精确解为u(x,y,z,t)=−2α[a2e−γαtnxπcos(nxπx)sin(nyπy)sin(nzπz)]d0， 初边值条件由其给出，其中，d0=1/[a1+a2e−γαtsin(nxπx)sin(nyπy)sin(nzπz)].

问题3 是三维Burgers 方程，在计算中取nx=ny=nz=3,a1=1,a2=0．1.表6 给出了该问题在不同空间网格数下，当τ =h2,α=0．1,t=0．5时 的L∞范数误差及收敛阶的比较.DHOC 格式[20]是条件稳定的，且时间上的精度为二阶、空间上的精度为四阶.从表中可以看出，DHOC 格式和本文格式在空间上的计算精度均达到了四阶.表7 给出了该问题当 τ =0．001,t=0．1时 ，对不同的 α本 文格式与DHOC 格式的L∞范数误差的比较，从表中可以看出，本文格式的计算误差与DHOC 格式的计算误差具有相同的数量级，但本文格式的计算误差更小一些.表8 给出了当h=0．062 5,α=0．1,t=2， 网格比 λ分别取1，2，4，6，8 时，针对该问题本文格式与DHOC 格式的L∞范 数误差和L2范数误差的比较.计算结果表明，随着λ 不断增大，DHOC 格式的计算结果发散，即格式是条件稳定的，而本文格式的计算结果仍保持收敛，即格式是无条件稳定的，这与两种格式各自的理论分析完全一致.充分说明了同为高精度紧致格式，相较于DHOC 格式，本文方法更适合求解此类非线性问题.

表6 问题3 当τ =h2, α=0．1, t=0．5时 的 L∞ 误 差和 L2误差及收敛阶Table 6 L∞ errors, L2 errors and convergence rates for τ =h2, α=0．1,t=0．5 in problem 3

表7 问题3 当τ =0．001, t=0．1时 ，对不同α 的 L∞误差Table 7 L∞ errors for τ =0．001, t=0．1 with different α values in problem 3

表8 问题3 当h =0．062 5, α=0．1, t=2 时 ，对不同网格比λ =τ/h2 的 L∞误 差和 L2误差Table 8 L ∞ errors and L 2 errors for h =0．062 5, α=0．1, t=2 with different mesh ratio λ =τ/h2 values in problem 3

5 结 论

本文针对含有变对流项和反应项系数的三维非稳态对流扩散反应方程，空间一阶导数采用四阶Padé公式计算，二阶导数利用一阶导数来逼近，时间采用Taylor 级数展开和余项修正，构造了一种高精度紧致有限差分格式，该格式被证明是无条件稳定的，其截断误差为即格式具有空间四阶精度，时间二阶精度.由于格式是两层的，因此可直接利用初始条件启动，不需要如文献[18]一样另外构造启动层的计算格式，从而使得计算过程更为简便.最后进行的数值实验验证了本文的理论结果.