求解变系数对流扩散反应方程的指数型高精度紧致差分方法

求解变系数对流扩散反应方程的指数型高精度紧致差分方法田芳;葛永斌
【摘要】本文给出了一种数值求解变系数对流扩散反应方程的指数型高精度紧致差分方法.我们首先将模型方程变形,借助常系数对流扩散方程的指数型高精度紧致差分格式,采用残量修正法得到变系数对流扩散反应方程的指数型高精度紧致差分格式;并从理论上分析了当Pelect数很大时,本文格式达到四阶计算精度时网格步长的限制条件;离散得到的代数方程组可采用追赶法直接求解.数值实验结果与理论分析完全吻合,表明了本文格式对于边界层问题或大梯度变化的物理量求解问题具有的高精度和鲁棒性的优点.%An exponential high accuracy compact finite difference method is proposed to solve the one-dimension (1D) convection-diffusion-reaction equation with variable coefficients. Fir-stly, the equation is rewritten in the form of convection diffusion equation. Then the exponential high order compact finite difference scheme for the convection diffusion equation with constant coefficients and the remainder term modification approach are utilized to obtain an exponential high accuracy compact finite difference scheme for the 1D convection-diffusion-reaction equation with variable coefficients. Secondly, the necessary condition on grid step length is analyzed theoretically if the scheme in this paper has a fourth-order accuracy when the Peclet number is very high. Lastly, the Thomas approach is applied to deal with the algebraic equations. Numerical examples, mostly with the boundary layer where sharp gradients may appear due to high Peclet number, are
presented to demonstrate the accuracy and robustness of the proposed scheme.
【期刊名称】《工程数学学报》
【年(卷),期】2017(034)003
【总页数】14页(P283-296)
【关键词】对流扩散反应方程;指数型有限差分格式;高精度紧致差分格式;对流占优;边界层
【作者】田芳;葛永斌
【作者单位】宁夏大学数学统计学院,银川 750021;宁夏大学数学统计学院,银川750021
【正文语种】中文
【中图分类】O241.82
1 引言
对流扩散反应方程是自然界和工程应用中的一类非常重要的方程，常常用来描述大气、海洋、河流等污染中污染物的扩散与分布、细菌的浓度分布、核工业中核反应堆的冷却及工业生产中的化学气相沉积等对流扩散反应现象．因此，对于此类方程的数值求解研究具有重要的理论意义和实际应用价值．目前，求解对流扩散反应方程的方法主要有有限元法、有限体积法、边界元法、特征线方法、有限差分法等[1-18]．文献[11]通过网格加密技术与变分多尺度法的耦合，消除了在对流扩散反应方程中由边界层效应和内部层效应引起的数值伪振荡．文献[12]利用微分方程的
通解，提出了数值求解一维定常常系数对流扩散反应方程的一种指数型高精度差分方法；文献[13]分别针对对流扩散方程和对流扩散反应方程采用变量替换消去方程中的对流项，将方程转化为反应扩散方程组进行离散求解；文献[14,15]构造了非
均匀网格上的对流扩散反应方程的多项式型的紧致差分格式，该格式较均匀网格系统下的计算无论在精度上还是分辨率上均有很明显的优势；文献[16]基于泰勒级数展开，结合原方程，发展了求解对流扩散方程的多项式型四阶精度的紧致差分格式；文献[17]发展了求解对流扩散反应方程的多项式型四阶精度的紧致差分格式，并将该格式外推得到多项式型六阶精度的紧致差分格式；文献[18]采用差分修正的思想构造了求解对流扩散方程的指数型差分格式，该格式对于对流占优问题在粗网格上能得到很高的计算精度．本文将针对对流扩散反应模型方程，将文献[18]的差分格式进行推广，发展了求解对流扩散反应方程的指数型高精度紧致差分格式．
2 对流扩散反应方程的指数型高精度紧致差分格式
考虑定常对流扩散反应方程
其中a是扩散系数，一般是常数；c(x)和b(x)分别为对流项系数和反应项系数，可以是常数，亦可以是x的函数；f(x)是x的足够光滑的函数；u(x)是待求未知量；
当b(x)=0时，模型方程(1)为对流扩散方程．
将求解区间[X1,X2]等分为N个子区间
在点xi处由泰勒级数展开得到
其中
2.1 常系数对流扩散反应方程的指数型高精度紧致差分格式
考虑常系数对流扩散反应方程
将方程(4)改写为
对于常系数对流扩散方程
文献[18]中给出了形如
的O(h4)阶精度的指数型紧致差分格式，其中系数为
下面推导常系数对流扩散反应方程(4)的四阶指数型紧致差分格式．对(5)式中的第一个方程，考虑其在点xi处的四阶指数型紧致有限差分格式
由(5)式中的第二个方程直接求导得
将(10)式代入(9)中，整理得
将(2)式代入(11)式右端项，略去高阶项整理得常系数对流扩散反应方程(4)的四阶指数型高精度紧致差分格式
其中
2.2 变系数对流扩散反应方程的指数型高精度紧致差分格式
考虑变系数对流扩散反应方程(1)在点xi处的形如(12)式的差分格式
其中差分系数由(13)给出．由泰勒级数展开得差分格式(14)的修正方程为
由(15)式知，差分格式(14)是方程(1)的O(h2)阶近似差分，要得到其O(h4)的差分格式，需要在原模型方程(1)的左端添加
进行修正，从而得到模型方程(1)的修正微分方程
其中
然后，对修正微分方程(16)考虑在点xi处的形如(12)式的差分格式，有
其中
差分格式(18)和(19)即为对流扩散反应方程(1)的四阶指数型紧致差分格式．(17)和(18)式中的fx,fxx可以直接求导得到，或者采用二阶中心差分算子(3)离散得到，都能保证差分格式(18)的O(h4)精度．将差分算子(3)代入到(18)式中，整理得
(20)式所对应的代数方程组可以采用追赶法直接求解．
3 差分格式的性质
当c̸=0时，记则
由于当时，故当雷诺数很大时，在粗网格差分格式(7)收敛阶也能达到三阶．若要差分格式(7)达到四阶精度，则需要满足
因此，对于差分格式(7)有如下结论：
定理1 1) 如果ac>0，则当网格步长满足
时，差分格式(7)具有四阶精度；
2) 如果ac<0：
(i) 当c>−16a时，则当对任意的网格步长h，差分格式(7)具有四阶精度；(ii) 当c<−16a时，则当网格步长满足时，差分格式(7)具有四阶精度．
证明不妨假设a>0，则分c>0和c<0两种情况证明．由(22)式得
1) 当c>0时，有2ch2+ch−2a<0，由于∆=c2+16ac>0，故
2) 当c<0时，有−2ch2+ch+2a>0，由于∆ =c2+16ac>0，故：当c>−16a时，则对任意的网格步长h，(22)式均成立；当c<−16a时，得
若要差分格式(12)能达到四阶精度，则除(22)式外还需要满足
因此，对于差分格式(12)，我们有如下结论：
定理2 假设a>0，那么：
1) 如果bc>0，则：
(i) 如果b>0,c>0，当网格步长满足
时，差分格式(12)具有四阶精度；
(ii) 如果b>0,c<0，则当网格步长满足时，差分格式(12)具有四阶精度；
2) 如果bc<0，则当网格步长满足时，差分格式(12)具有四阶精度．
证明假设a>0，由(24)式得
1) 如果bc>0，则：
(i) 当b>0,c>0时，有
即4ch2+(b+1)ch−2a(b+1)<0，由于∆=(b+1)2c2+32a(b+1)c>0，故
(ii) 当b<0,c<0时，有
即(b−1)ch+2a(b−1)<0，从而得
2) 如果bc<0，则：
(i) 当b<0,c>0时，有
即(1−b)ch+2a(b−1)<0，从而得
(ii) 当b>0,c<0时，有
即
由于函数是(0,+∞)上关于b的单调递减函数，故有
从而得
综上(i)和(ii)所述，如果bc<0，则当网格步长满足时，差分格式(12)具有四阶精度．
4 数值算例
下面，我们将选取典型算例采用本文格式(12)(简称为EHOC)进行计算，并与精确解和文献[17]中的多项式型格式(简称为FOC)的计算结果进行比较，验证本文格式的精确性和可靠性．收敛阶通过公式计算得到，其中err(N1)和err(N2)分别为网格数为N1和N2时最大绝对误差．
算例1 −εuxx+ux= επ2sin(πx)+ πcos(πx)
精确解为
当ε很小时，该精确解在x=1处有一边界层．此算例中，a=ε,c=1，则当
时，计算收敛阶为四阶．
表1比较了当ε取不同值时，采用EHOC格式和FOC格式计算的最大绝对误差．计算数值结果表明，ε=0.1时，FOC格式的计算误差和EHOC格式的计算误差达到同一个数量级，但当ε=0.01,0.001时，EHOC格式的计算精度明显优于FOC格式，即随着ε取值减小，EHOC格式较FOC格式在计算精度方面的优势明显增加．图1和图2比较了当ε=10−3,10−5时采用EHOC格式和FOC格式的计算结果．从图中明显的看到，在粗网格上(N=16)，EHOC格式的计算解和精确解在计算节点上吻合的很好．
表1: 算例1最大绝对误差比较
图1: 当ε=10−3，网格数=16,128时，采用EHOC、FOC格式的计算解和精确解的比较
图2: 当ε=10−5，网格数=16,128,1024,8192时，采用EHOC、FOC格式的计算解和精确解的比较
算例2 −εuxx+ux+(1+ε)u=0,0<x<1
精确解为
此算例中，a=ε,c=1,b=1+ε，满足bc>0且b>0,c>0，则当
时，计算收敛阶为四阶．
表2比较了当ε取不同值时，采用EHOC格式和FOC格式计算的最大绝对误差．从表格中数值结果我们可以看到，对任何的ε取值，EHOC格式的计算精度
都高于FOC格式的计算精度．
表2: 算例2最大绝对误差比较
算例3 −εuxx−ux+εu=sinx
精确解为
其中
此算例中，a=ε,c=−1,b=ε，满足bc<0，则当时，计算收敛阶为四阶．
表3比较了当ε取不同值时，采用EHOC格式和FOC格式计算的最大绝对误差．从表格中数值结果我们可以看到，对任何的ε取值，EHOC格式的计算精度
都高于FOC格式的计算精度．图3和图4给出了当ε=10−3,10−5时，采用EHOC格式和FOC格式的计算数值结果和精确解的比较．图3中，当ε=10−3时，EHOC格式在粗网格(N=16)上能够很好的逼近精确解(此时最大绝对误差为
2.94×10−4)，而FOC格式若要获得相当量级的计算精度，则至少需要4200个网格点(此时最大绝对误差为
3.18×10−4)．而当ε=10−5时，如图4(d)所示，当网
格数增加到9000时，在靠近边界层的地方，FOC格式的计算解和精确解的误差
(为1.42)很大，而此时EHOC格式的计算误差为5.76×10−10．
表3: 算例3最大绝对误差比较εFOC[17]网格数收敛阶误差误差收敛阶8
9.58(−2)2.65(−3)16 0.11.22(−2)2.983.54(−4)2.90
328.72(−4)3.802.58(−5)3.89 645.22(−5)4.061.56(−6)4.08 128
2.91(−2)8.14(−6)0.012562.43(−3)
3.587.15(−7)3.51 5121.43(−4)EHOC
4.094.25(−8)4.07 10248.93(−6)4.002.68(−9)3.99 512
2.92(−1)7.19(−7)0.00110245.78(−2)2.341.57(−7)2.19
20485.81(−3)3.311.69(−8)3.22 40963.54(−4)4.041.05(−9)4.00
图3: 当ε=10−3，网格数=16,128,4200时，采用EHOC、FOC格式的计算解和
精确解的比较
算例4
精确解为当ε<<1时，该解在x=1处有一边界层．此算例中，满足bc>0且
b>0,c>0，则当
时，计算收敛阶为四阶．
图4: 当ε=10−5，网格数=16,128,9000时，采用EHOC、FOC格式的计算解和
精确解的比较
表4比较了当ε取不同值时，采用EHOC格式和FOC格式计算的最大绝对误差．计算数值结果表明，EHOC格式较FOC格式在计算精度方面具有明显的优势．表4: 算例4最大绝对误差比较
图5和图6给出了当ε=10−3,10−5时，采用EHOC格式和FOC格式的计算数值结果和精确解的比较．从图中明显的看到，在强对流和反应占优的情况下，即使在粗网格上(N=16)，EHOC格式都能很好的逼近精确解．当ε=10−3时，若要获得相当的计算精度，采用FOC格式计算需要在620个网格下的计算，计算量是
EHOC格式计算的38.75倍(EHOC的最大误差为4.33×10−4，FOC的最大误差
为4.48×10−4)．当ε=10−5时，如图6(c)所示，将网格点增加255倍，采用FOC格式计算在光滑区域能很好的和精确解吻合，但在边界处误差很大，而EHOC格式能够很好的逼近精确解：EHOC的最大误差为4.25×10−8，FOC的最大误差为7.45×10−1．
图5: 当ε=10−3，网格数=16,620时，采用EHOC、FOC格式的计算解和精确解的比较
图6: 当ε=10−5，网格数=16,4096时，采用EHOC、FOC格式的计算解和精确
解的比较
5 结论
本文基于常系数对流扩散方程的四阶指数型高精度紧致差分格式，结合残量修正法发展了一种数值求解一维变系数对流扩散反应方程的指数型高精度紧致差分方法，并选取了四个典型算例与文献[18]中的四阶多项式型格式的计算结果进行了比较，所有算例计算结果与表5中的理论分析结论一致，对于对流(反应)占优问题、边界层问题和大梯度变化的物理量问题本文格式的计算精度要明显优于文献[17]的格式．表5: 当ε取不同值时，差分格式达到四阶精度的条件
参考文献：
[1]Radhakrishna Pillai A C.Fourth-order exponential finite difference methods for boundary value problems of convective diffusion
type[J].International Journal for Numerical Methods in
Fluids,2001,37(1):87-106
[2]Chen G Q,Gao Z.A perturbational exponential finite difference scheme
for the convection diffusion equation[J].Journal of Computational Physics,1993,104(1):128-138
[3]Phongthanapanich S,Dechaumphai bined finite volume and finite element method for convectiondiffusion-reaction equation[J].Journal of Mechanical Science and Technology,2009,23(3):790-801
[4]Bause M,Schwegler K.Higher order finite element approximation of systems of convection-diffusion-reaction equations with small
diffusion[J].Journal of Computational and Applied
Mathematics,2013,246:52-64
[5]Ayuso B,Marini L D.Discontinuous Galerkin methods for advection-diffusion-reaction problems[J].SIAM Journal on Numerical
Analysis,2009,47(2):1391-1420
[6]Shih Y,Kellogg R B,Tsai P.A tailored finite point method for convection-diffusion-reaction problems[J].Journal of Scientific
Computing,2010,43(2):239-260
[7]Angelini O,Brenner K,Hihorst D.A finite volume method on general meshes for a degenerate parabolic convection-reaction-diffusion equation[J].Numerische Mathematik,2013,123(2):219-257
[8]张建松，朱江，郭会，等.对流扩散反应方程的特征分裂最小二乘方法[J].高等学校计算数学学报，2012,34(4):289-299 Zhang J S,Zhu J,Guo H,et al.A characteristic splitting least-squares method for convection-diffusionreaction equations[J].Numerical Mathematics:A Journal of Chinese Universities,2012,34(4):289-299
[9]朱国庆，李清善，陈绍春.分层网格上对流扩散反应方程的双二次元逼近[J].河南大学学报(自然科学版)，2007,37(3):221-225 Zhu G Q,Li Q S,Chen S
C.Biquadraic element approximation of convection-diffusion-reaction
equation under graded meshes[J].Journal of Henan University(Natural Science Edition),2007,37(3):221-225
[10]崔翔鹏，贺力平.非线性对流扩散反应方程的预估-校正单调迭代差分方法[J].上海交通大学学报，2007,41(10):1731-1736 Cui X P,He L P.Predictor-corrector monotone iterative difference method for nonlinear convection diffusion reaction equation[J].Journal of Shanghai Jiaotong
University,2007,41(10):1731-1736
[11]朱海涛，欧阳洁.对流-扩散-反应方程的变分多尺度解法[J].工程数学学报，2009,26(6):997-1004 Zhu H T,Ouyang J.Variational multiscale method for the advection-diffusion-reaction equation[J].Chinese Journal of Engineering Mathematics,2009,26(6):997-1004
[12]魏剑英.定常对流扩散反应方程的指数型高阶紧致差分格式[J].宁夏大学学报(自然科学版)，2012,33(2):140-143 Wei J Y.High-order exponential finite difference method for 1D convection-diffusion-reaction equation[J].Journal of Ningxia University(Natural Science Edition),2012,33(2):140-143
[13]杨录峰，李春光.一种求解对流扩散反应方程的高阶紧致差分格式[J].宁夏大学学报(自然科学版)，2013,34(2):101-104 Yang L F,Li C G.A high-order compact finite difference scheme for solving the convection difference reaction equations[J].Journal of Ningxia University(Natural Science Edition),2013,34(2):101-104
[14]田芳，田振夫.定常对流扩散反应方程非均匀网格上高精度紧致差分格式[J].工程数学学报，2009,26(2):219-225 Tian F,Tian Z F.A high accuracy compact difference scheme for convection diffusion reaction equation on non-uniform grid[J].Chinese Journal of Engineering
Mathematics,2009,26(2):219-225
[15]兰斌，薛文强，葛永斌.对流扩散反应方程基于坐标变换的高阶紧致差分格式[J].青岛科技大学学报(自然科学版)，2014,35(1):100-106 Lan B,Xue W Q,Ge Y B.A High-order compact difference scheme based on the coordinate transformation for the convection diffusion reaction equation[J].Journal of Qingdao University of Science and Technology(Natural Science Edition),2014,35(1):100-106
[16]Spotz W F.High-order compact finite difference schemes for computational mechanics[D].Austin:University of Texas at Austin,1995 [17]Sun H W,Zhang J.A high order finite difference discretization strategy based on extrapolation for convection diffusion equations[J].Numerical Methods for Partial Differential Equations,2004,20(1):18-32
[18]Tian Z F,Dai S Q.High-order compact exponential finite difference methods for convection-diffusion type problems[J].Journal of Computational Physics,2007,220(2):952-974。