求解一维对流扩散反应方程的高阶紧致格式

赵秉新

(宁夏大学数学计算机学院，银川 750021)

对流扩散及反应现象广泛存在于自然界和工程应用中，如在大气、河流污染中污染物质的扩散，受温盐扩散影响的大洋环流，室内空气的调节，核工业中核反应堆的冷却及工业生产中的化学气相沉积中均存在复杂的对流扩散反应现象。描述这些问题的典型方程为流体动力学方程，而对流扩散反应方程(CDR，convection-diffusion-reaction)是其中的基本方程。由于无法求得对流扩散反应方程的解析解，因此，对此类方程的数值求解具有十分重要的理论和实际应用价值。

目前，求解CDR方程的常用数值方法为有限差分法、有限元法及有限体积法等，其中有限差分法作为一种发展成熟且易于应用的方法受到了广泛的关注。数值求解对流扩散反应型方程不同程度地存在着产生数值扩散和振荡解的可能，如经典二阶中心格式会产生非物理振荡，导致计算结果不可靠［1-2］等。减小数值振荡的有效方法是构造稳定的、且能有效反映对流项的迎风效应的数值格式，要求差分格式的设计要保持流体流动的传输特征。常用的一阶迎风格式是一种无条件稳定的格式，但却只有一阶精度，且存在无法有效捕捉微小物理量等缺陷。因此，要想对问题进行较成功的数值模拟，同时满足实际问题中高精度的要求，就必须对计算格式的求解效率及精度提出更高要求。高阶紧致迎风差分格式通过使用较少的网格基架点便能达到提高精度的要求，保持了迎风格式的优势，具有边界无需特殊处理以及能达到与谱方法相近的分辨率等优点，是数值计算方法研究的一个重要方向。Noye等［3］利用加权修正技术，提出了一种三阶半隐格式，但其是条件稳定的。陈国谦等［4］通过对对流系数和源项作2阶修正，消除截断误差中的低阶项，给出了对流扩散方程的4阶指数型格式。王彩华［5］通过在低精度格式基础上进行简单修正得到了针对一维无源对流扩散方程的4种高精度差分格式。通过在低阶格式中的源项中引入紧致修正项，杨茉等［6］提出了一种紧致修正方法。Ding等［7］提出了一种O(h4+τ4)阶无条件稳定的格式，但适用于对流占优问题的求解［8］，且存在求解效率不高的问题。

本文通过将原方程作简单转换，对得到的对流扩散方程中的对流项采用四阶紧致迎风格式离散，扩散项采用四阶Padé格式离散;之后，对于空间半离散格式，时间方向采用四阶龙格库塔方法计算，从而整体达到了O(h4+τ4)阶精度。经数值验证，该格式具有良好性能。该格式构造方法简单，容易理解，且易于推广到高维情形。

1 差分格式的构造

本文考虑一维非定常对流扩散反应方程:

其中:Ω =［0，l］×［0，T］;φ(x)为充分光滑的函数;常扩散系数a＞0;c(x，t)表示相速度;r为反应速度;g1，g2为常数。

以τ=Δt为时间步长，空间取均匀等间距网格，步长为 h= Δx=1/m，网格点为(xi，tn)，其中xi=ih，i=0，1，…，m，tn=nτ，n≥0。对方程(1)做如下变换:令u=exp(-rt)φ，代入方程(1)中消去反应项，对定解条件(2)和(3)作相应变换，得

方程(4)为未知量φ的对流扩散方程，记Q(x，t)=e-rtf(x，t)为源项。以下给出式(4)～(6)的差分格式，回代求解方程(1)～(3)。

1.1 空间离散

将方程(4)中的扩散项φxx采用四阶精度Padé格式离散:

其中Si为二阶偏导数φxx在i点的逼近值。该格式的离散矩阵为对角占优三对角矩阵，可直接采用追赶法求解。

对对流项进行离散时，考虑迎风效应而采用四阶紧致迎风格式［13］离散。例如对流项c(x，t)φx可表示成

其中Si利用式(7)计算。该格式基于3个网格基架点，具有四阶精度。

1.2 时间方向离散

Runge-Kutta方法最初用于常微分方程的数值求解。当把时间处理成独立变量时，可将式(4)进行空间离散后的半离散方程视为只依赖于时间的常微分方程，此时，可采用Runge-Kutta方法进行时间推进求解。本研究将式(4)的半离散格式写为

其中Lh为空间半离散算子，且

Qi为源项Q(x，t)在i点的离散值。

对半离散格式(8)，采用四步四阶Runge-Kutta方法求解，具体为:

其中Lh(φn)为φ在n时间层时式(9)的值。

2 数值验证

为验证格式的性能，考察3个有精确解的算例，并与已有格式的结果进行对比。

2.1 算例1

在x∈［0，1］上考虑如下对流扩散反应问题:

其中源项 f(x，t)=e-t[(r-1)x2+2(cx-a)]，其解析解为 u(x，t)=x2e-t。

本文对以下3种工况进行计算:

1)a=10-4，c=102，r=1，对流占优情形;

2)a=10-4，c=1，r=102，反应占优情形;

3)a=10-4，c=102，r=102，对流及反应占优情形。

对于上述3种工况，数值计算的结果均与精确解吻合得很好。图1(a)、(b)分别给出在工况1)的对流占优情形下，步长h=1/20和1/40时，问题(11)分别在t=1.5和3.0时刻数值结果与解析解的对比。可以看出，在粗细2种网格下，数值解均与精确解均吻合得很好。

图1 工况(1)下数值解与精确解的对比

2.2 算例2

考虑如下的非线性含源对流扩散反应方程:

精确解取为u(x，t)=tanh( Re(1-2x)/4)et，源项f(x，t)可由原方程确定。

图2给出了Re=1 000，T=0.02时的精确解及CDS、SCD4和本文格式在中心点附近的解(步长取 h=1/320，τ=10-5)。可见，CDS格式和SCD4格式在该网格尺度下均产生明显的数值振荡，而本文格式能够获得很好的数值解，体现了其求解非线性对流扩散反应方程时的良好性能。

图2 算例2不同格式在间断位置附近的解

2.3 算例3

考虑对流扩散方程:

精确解为 u(x，0)=e5x-t(0.25+0.01π2)sin(πx)。

时间步长取τ=0.001，表1给出了T=20时，本文格式(Present)与Crank-Nicolson格式及Ding格式［7］在不同空间网格步长下的L2误差和收敛阶(Rate)的对比情况。收敛阶由下式定义:

其中e1、e2分别表示网格步长为 h1、h2时的 L2误差。

通过对比表明:本文格式与Ding的格式［7］均具有四阶空间精度，高于Crank-Nicolson格式的二阶精度;但Ding格式均存在求解效率不高的问题，即在求解过程中需要进行矩阵求逆和乘积等运算，当网格数很大时，求解效率会明显降低，特别是在求解非线性问题时，由于每个迭代步的迭代矩阵均不相同，这样每个迭代步都必须进行矩阵求逆运算，在高维或多网格节点下，计算效率会明显降低。

表1 算例1中不同网格尺度下L2误差和收敛阶对比

4 结束语

本文给出了一种求解对流扩散反应问题的四阶紧致迎风差分格式，格式构造过程简单易懂。数值算例的验证表明，本文格式能有效控制格式的数值振荡，适用于对流占优问题的模拟。另外，该格式具有构造简单，易于编程等特点，且可直接推广到高维问题。对于定常问题，将∂φ/∂t视为人工时间项，计算程序无需太多修改即可直接求解。

［1］Brandt A，Yavneh I.Inadequacy of first-order upwind difference schemes for some recirculating flows［J］.Journal of Computational Physics，1991，93(1):128-143.

［2］Zhang J.Accelerated multigrid high accuracy solution of the convection-diffusion equation with high Reynolds number［J］.Numerical Methods for Partial Differential E-quations，1997，13(1):77-92.

［3］Noye B J，Tan H H.A third-order semi-implicit finite difference method for solving the one-dimensional convection-diffusion equation［J］.International Journal for Numerical Methods in Engineering，1988，26(7):1615-1629.

［4］Chen G Q，Gao Z，Yang Z F.A perturbational h4exponential finite difference scheme for the convective diffusion equation［J］.Journal of Computational Physics，1993，104(1):129-139.

［5］王彩华.一维对流扩散方程的一类新型高精度紧致差分格式［J］.水动力学研究与进展:A辑，2004，19(5):655-663.

［6］张昆，杨茉.求解对流扩散方程的紧致修正方法［J］.工程热物理学报，2011，32(3):459-461.

［7］Ding H，Zhang Y.A new difference scheme with high accuracy and absolute stability for solving convection-diffusion equations［J］.Journal of Computational and Applied Mathematics，2009，230(2):600-606.

［8］Tian Z F，Yu P X.A high-order exponential scheme for solving 1D unsteady convection-diffusion equations［J］.Journal of Computational and Applied Mathematics，2011，235(8):2477-2491.

［9］Tian Z F，Liang X，Yu P X.A higher order compact finite difference algorithm for solving the incompressible Navier-Stokes equations［J］.International Journal for Numerical Methods in Engineering，2011，88(6):511-532.