二维非稳态对流扩散方程的高阶紧致差分格式

魏剑英，王燕，葛永斌

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

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

有关非稳态对流扩散方程的有限差分法的文献报道已有很多，如文献[1]利用半离散法和Padé逼近，推导了一维方程的截断误差为O(τ5+h4)的隐格式，且是无条件稳定的； 文献[2]利用Taylor级数展开和待定系数法构造了一致四阶紧致格式并结合三阶TVD-Runge-Kutta方法求解了一维方程； 文献[3]针对二维方程，对空间导数项和时间导数项分别采用四阶离散和加权平均的方法，推导出一种截断误差为O(τ2+h4)的格式； 文献[4]采用半离散的方法，从一维定常问题出发，利用四阶紧致差分算子和Crank-Nicolson(C-N)格式，推导了一种求解二维常系数非稳态对流扩散方程的差分格式，截断误差为O(τ2+h4)，且是无条件稳定的； 文献[5]针对此类方程，构造了有理型紧致交替方向隐式(ADI)格式，截断误差为O(τ2+h4)，且是无条件稳定的； 文献[6]还提出了该方程无条件稳定的指数型紧致ADI差分格式，其截断误差依然为O(τ2+h4)； 文献[7]利用指数变换消除对流项后转化为扩散方程，再利用紧致公式和扩展的Simpson公式构造了一维方程的高阶格式，截断误差为O(τ4+h4)，且是无条件稳定的； 文献[8]针对二维非稳态常系数方程，空间方向直接采用六阶组合紧致差分公式进行计算，时间方向用C-N格式离散，所提格式的截断误差为O(τ2+h6)，尽管该格式空间达到了六阶精度，但是由于其时间只有二阶精度，因此为了保证空间精度达到六阶，其计算所需要采取的时间步长必须为O(h3)，即必须采用较小的时间步长，并且该方法仅适用于常系数问题。

本文针对二维非稳态变系数对流扩散方程，首先对空间二阶导数采用一阶导数的四阶逼近公式，而一阶导数采用四阶Padé逼近，时间项采用二阶向后差分公式(BDF)，得到一个时间二阶、空间四阶精度的紧致差分格式； 然后利用Taylor级数展开，对一阶和二阶空间导数项采用六阶紧致差分公式，时间项采用三阶BDF得到一个时间三阶、空间六阶精度的紧致差分格式．由于采用了向后差分，因此两种格式均是无条件稳定的。

1 差分格式的建立

考虑如下二维非稳态对流扩散方程：

ut-α(uxx+uyy)+p(x，y，t)ux+q(x，y，t)uy=f(x，y，t) (x，y)∈Ω，t≥0

(1)

初始条件：

u(x，y，0)=g(x，y) (x，y)∈Ω

(2)

边界条件：

u(x，y，t)=s(x，y，t) (x，y)∈∂Ω，t>0

(3)

其中：Ω∈[a1，a2]×[b1，b2]为R2上的矩形区域，∂Ω为Ω的边界，u(x，y，t)为待求未知量，p(x，y，t)，q(x，y，t)分别为x，y方向的对流项系数，α为扩散项系数(常数)，f(x，y，t)为源项，且p(x，y，t)，q(x，y，t)，f(x，y，t)，g(x，y)，s(x，y，t)均为已知函数。

(4)

(5)

1.1 四阶精度紧致格式

将方程(1)改写为

ut=f(x，y，t)+α(uxx+uyy)-p(x，y，t)ux-q(x，y，t)uy

(6)

对二阶导数uxx和uyy采用如下逼近：

(7)

(8)

将(7)，(8)式代入(6)式，考虑其在点(xi，yj)的值，得到

(9)

为了使(9)式空间保持四阶精度，对ux和uy的计算采用如下四阶Padé公式[9]

(10)

(11)

考虑(9)式在第n+1时间层的离散值，对ut采用如下二阶BDF[10-11]进行离散

(12)

舍去高阶项，可得

(13)

(14)

其中G为(9)式的右端项。

在文献[12]中已经证明当K≤6(K代表精度)时，BDF是无条件稳定的．另外，由式(13)可以看出，在每个时间步上，该格式为5点模板的全隐格式．为了得到u(x，y，t)的计算值，(14)式中出现的ux和uy除需知道其内点值，还需知道其在边界点处的值，为保持与内点差分格式同样的精度，对边界上的ux，uy采用一致四阶边界条件[2]计算。

1.2 六阶精度紧致格式

为了得到方程(1)的六阶精度的格式，需由Taylor级数展开得到一、二阶导数的六阶近似公式

(15)

(16)

(17)

(18)

将(15)-(18)式代入(6)式，有

(19)

(20)

(21)

(22)

(23)

(24)

(25)

将(20)-(25)式代入(19)式，考虑其在点(xi，yj)的值，即可得到如下半离散的紧致格式

(26)

为了使(26)式具有六阶精度，对ux和uy采用文献[9]中的六阶差分公式计算，uxx和uyy采用文献[13]中的六阶差分公式计算

(27)

(28)

(29)

(30)

考虑(26)式在第n+1时间层的值，对ut采用三阶BDF[10]进行离散

(31)

略去高阶项，可得

(32)

2 数值实验

本文将利用以下4个问题进行数值实验，验证(2，4)格式和(3，6)格式的稳定性和有效性，文中涉及到的误差及收敛阶定义如下：

1) 最大绝对误差：

2)L2范数误差：

3) 收敛阶：

问题1考虑如下非齐次对流扩散问题

其精确解为u(x，t)=e-txy(1-x)(1-y)，初边值条件及右端项f(x，y，t)均由精确解给出。

表1给出了问题1当τ=h2，t=0.5时的计算结果．从表中可以看出，本文(2，4)格式和文献[6]格式在空间上均达到了四阶精度，但本文(2，4)格式的计算结果比BTCS格式和文献[6]格式更好； 本文(3，6)格式在空间上达到了六阶精度(网格为64时的数值已经达到机器精度，所以收敛阶受到影响不准确)，且最大绝对误差比(2，4)格式的小2～3个数量级，计算结果更精确．当h=0.01，t=0.5时，不同时间步长τ下的L∞误差、收敛阶及CPU时间由表2给出．从表中可以看出，本文(2，4)格式时间上达到的精度为二阶，本文(3，6)格式时间上达到的精度为三阶，这与理论分析一致； 另外当τ较大时计算收敛慢，所用的CPU时间相对较长，计算过程中(3，6)格式需要计算在节点处的一、二阶导数值，因此所用的计算时间比(2，4)格式长。

表1 问题1当τ=h2，t=0.5时的最大绝对误差L∞及收敛阶

表2 问题1当h=0.01，t=0.5时的最大绝对误差L∞，收敛阶及CPU时间

问题2考虑如下变系数对流扩散问题

其精确解为u(x，t)=txy(1-x)(1-y)ex+y，初边值条件及右端项f(x，y，t)均由精确解给出。

表3给出了问题2当τ=h2，t=0.25时本文格式与C-N格式和BTCS格式的L∞误差及收敛阶的比较．从表中可以看出，C-N格式和BTCS格式在空间上达到的精度均为二阶，本文(2，4)格式在空间上达到的精度为四阶，而本文(3，6)格式空间上达到的精度为六阶，且(3，6)格式的计算结果误差更小，充分验证了本文两种格式的精确性和有效性．另外，上述两个例子的计算结果进一步说明低精度的启动步并不影响格式的整体精度。

表3 问题2当τ=h2，t=0.25时的L∞误差及收敛阶

问题3考虑如下高斯脉冲问题

图1 t=1，τ=2.5×10-3，Pe=2，1.2≤x，y≤1.8时的等值线图

图2 t=0.001，τ=2.5×10-6，Pe=2000，1.2≤x，y≤1.8时的等值线图

表4 问题3当h=0.02时在不同参数下的L∞误差和L2误差

问题4考虑如下非线性方程

表5给出了问题4当t=1，τ=h2时，本文(2，4)格式和(3，6)格式在α分别取1和0.1时的最大绝对误差与收敛阶．从表中可以看出，本文(2，4)格式在空间上达到的精度是四阶，(3，6)格式达到的精度是六阶(网格为64，α=1时最大绝对误差的值已经达到机器精度)，这与理论推导是一致的，从表中还可以看出，α=1时所用的CPU时间比α=0.1时的长，(3，6)格式的CPU时间比(2，4)格式的长．表6给出了当N=7，τ=0.01，α=1时，t=15和t=20的计算解与精确解的绝对误差．表中数据表明，本文(2，4)格式的计算误差比文献[14]格式小5～6个数量级，与文献[15]格式的计算结果大致相当； 本文(3，6)格式的计算结果比文献[14]格式小7～9个数量级，比文献[15]格式小2～3个数量级．因此，对于此类非齐次边界的Burgers方程，本文两种高精度紧致差分格式非常有效．图3给出了N=20，t=1，τ=0.01，α=0.1时本文两种格式的绝对误差和文献[15]格式的绝对误差比较，从图中可以看出，本文(2，4)格式的绝对误差与文献[15]格式的绝对误差大致相当，但本文(3，6)格式的计算结果更加精确，绝对误差更小。

表5 问题 4 当τ=h2，t=1时的最大绝对误差L∞，收敛阶和CPU时间

表6 问题4当N=7，τ=0.01，α=1时的绝对误差

续表6

图3 N=20，t=1，τ=0.01，α=0.1时的精确解与绝对误差

3 结论