对流占优扩散方程的几种差分解法比较

冯立伟

摘要：对流占优扩散方程是描述对流和扩散过程的模型方程。对于求解对流占优扩散方程的几种差分格式进行分析，并编制了MATLAB程序。使用算例进行了数值试验，差分格式消除了伪数值振荡，实现了对流占优扩散方程的求解，并比较了几种解法的优劣。

关键词： 对流占优扩散方程； 有限差分法； 迎风； 萨马尔斯基；指数

中图分类号：TP18 文献标识码：A 文章编号：1009-3044（2015）20-0155-03

Comparing with Numerical Solution of Several Difference Methods for Convection-dominated Diffusion Equation

FENG Li-wei

（Shenyang University of Chemical Technology ， Shenyang 110142， China）

Abstract： Convection-dominated diffusion equation describes convection and diffusion process. It analyzed finite difference methods solving Convection-dominated diffusion equation， and writes MATLAB programs. It makes numerical experiment using example， these experiments reveals the finite difference methods eliminating false numerical oscillation， And it compares advantages and disadvantages of the several finite difference methods.

Key words： convection-dominated diffusion equation； finite difference method； upwind； samarskii； exponential

对流扩散方程描述了在自然界中大量出现的对流和扩散现象，在流体力学、环境科学以及能源开发等诸多领域有着广泛的应用。因此研究对流扩散问题的数值计算方法就尤为重要。当对流占优即小ε值的情形，传统中心差分格式会产生数值振荡，对流项的处理上以及对小ε值的情形处理，成为求解对流扩散方程的关键问题。在迎风差分格式[1]中对流项采用一阶偏上游的迎风差商离散，修正中心差分格式， 萨马尔斯基格式，指数差分格式[2，3，4]通过添加摄动项来消除数值振荡。

1 对流扩散方程差分格式

一维对流占优扩散方程

[?u?t+a?u?x=ε?2u?x2+fx，t，x，t∈0，L×0，T]

其中[ux，t]为待求量，[a]是对流系数，[ε]是扩散系数，[f（x，t）]为源汇。

1.1 差分格式的建立

首先对[x-t]平面进行网格剖分。分别取[h，τ]为空间步长与时间步长，用两族平行直线[xj=jh ]，[tk=kτ ]将矩形区域[0， T×0， L]分割成矩形网格。采用中心差商对导数进行离散可得到中心差分格式，但在[ε]比较小时扩散效应比较小主要表现为对流效应，该格式会出现伪数值振荡。消除伪数值振荡的一种方法是添加摄动项，不同的摄动项导致了不同的差分格式，另一种消除伪数值振荡，对流项采用偏心差商的离散，从而得到迎风格式。

修正中心差分格式

[Uk+1j-Ukjτ+aUkj+1-Ukj-12h=ε+τa22Ukj+1-2Ukj+Ukj-1h2+fkj]

迎风差分格式

[Uk+1j-Ukjτ+aUkj-Ukj-1h=εUkj+1-2Ukj+Ukj-1h2+fkj]

萨马尔斯基格式

[Uk+1j-Ukjτ+aUkj-Ukj-12h=11+ha/2εεUkj+1-2Ukj+Ukj-1h2+fkj]

指数差分格式

[ Uk+1j-Ukjτ+aUkj+1-Ukj-12h=-ah21+eah/ε1-eah/εUkj+1-2Ukj+Ukj-1h2+fkj] 其中[k=1，2，…，N-1， j=1，2，…，M-1]。

由Taylor公式容易得出：它们都与一维对流扩散方程相容，其截断误差分别为[O（τ+h2）]，[O（τ+h）]，[O（τ+h2）]和[O（τ+h2）]。

1.2 差分格式的稳定性分析

对流扩散方程的中心差分格式的稳定性条件是[ετh2≤12，τ≤2εa2]【1】。对于修正中心差分格式，用[ε+τ2a2]代替中心差分格式稳定性条件中的[ε]得到稳定性条件[ε+τ2a2τh2≤12，]，[τ≤2a2ε+τ2a2]而后式恒成立，所以只需[ε+τ2a2τh2≤12]【1】。迎风差分格式可改写为

[Uk+1j-Ukjτ+aUkj+1-Ukj-12h=ε+ah2Ukj+1-2Ukj+Ukj-1h2+fkj]

用[ε+ah2]代替中心差分格式稳定性条件中的[ε]得到稳定性条件[ε+ah2τh2≤12，][τ≤2a2ε+ah2] 而前式可推出后式，所以稳定性条件为[ε+ah2τh2≤12]。

对于萨马尔斯基格式用[ε1+ha/2ε]代替迎风差分格式稳定性条件中的[ε]得到稳定性为[ah2+2ε22ε+ahτh2≤12]。在指数差分格式中 [1+eah/ε1-eah/ε=σε]，令[σ=ah2εcothah2ε]为拟合因子，指数格式可改写为

[Uk+1j-Ukjτ+aUkj+1-Ukj-12h=σεUkj+1-2Ukj+Ukj-1h2+fkj]，用[σε]代替中心差分格式稳定性条件中的[ε]得到稳定性条件[εστh2≤12，τ≤2εσa2]。

1.3差分格式的求解

令[λ=aτh]，[μ=ετh2]，将上述各个差分格式改写成便于计算的形式修正中心格式

[Uk+1j=Ukj-λ2Ukj+1-Ukj-1+μ+12λ2Ukj+1-2Ukj+Ukj-1+τfkj] 迎风差分格式

[Uk+1j=Ukj-λUkj-Ukj-1+μUkj+1-2Ukj+Ukj-1+τfkj]

萨马尔斯基格式

[Uk+1j=Ukj-λUkj-Ukj-1+11+ha/2εμUkj+1-2Ukj+Ukj-1+τfkj]

指数差分格式

[Uk+1j=Ukj-aτ2hUkj+1-Ukj-1-aτ2h21-eah/ε-1Ukj+1-2Ukj+Ukj-1+τfkj]

其中[k=1，2，…，N-1， j=1，2，…，M-1]。因为各个格式均为显式格式，可以逐层计算。将差分格式与离散化的初边值条件结合，由前一时间层上的结果可以计算下一时间层上的结果，从而得到所有的网格点[xj，tk]上的近似解[Ukj]。

1.4 数值实验

对前面的几种不同差分格式分别编制Matlab程序进行求解初边值问题。定义网格点[xj，tk]上的误差[ekj=u（xj，tk）-Ukj]，使用指标[ek1=1Mj=1Mekj]， [ek2=1Mj=1Mekj2]和[ekc=maxjekj]来度量第k个时间层上的总体误差。

算例1：

[?u?t+?u?x=ε?2u?x2]，相应的初边值条件为

[u（0，t）=exp1-1εx，u（1，t）=exp1ε+1-1εt]，[u（x，0）=expxε]，该问题的解析解为

[u（x，t）=exp1εx+1-1εt]

取[ε=0.1]，空间网格步长为[h=0.1]，时间网格步长为[τ=0.01]，计算结果绘制成图像如下：

图1 t=3时刻的解

图2 2-范数误差

图3 1-范数误差

图4 c-范数误差

取[ε=0.05]，空间网格步长为[h=0.1]，时间网格步长为[τ=0.01]，计算结果绘制成图像如下：

图5 t=3时刻的解

图6 2-范数误差

算例2：

[?u?t+?u?x=ε?2u?x2+πεe-π2tεcos（πxε）]，初边值条件为

[u（x，0）=sin（πx）ε，u（0，t）=0， u（1，t）=0]，该问题的解析解为

[u（x，t）=1εe-π2εtsin（πx）]。

取[ε=0.1]，空间网格步长[h =0.1]，时间网格步长为[τ=0.01]，计算到3秒时结果如图5-图8所示。

图7 t=3时刻的解

图8 2-范数误差

图9 1-范数误差

图10 c-范数误差

取[ε=0.01]，空间网格步长为[h=0.1]，时间网格步长为[τ=0.01]，计算结果绘制成图像如下

图11 t=3时刻的解

图12 2-范数误差

取[ε=0.001]，空间网格步长为[h=0.1]，时间网格步长为[τ=0.01]，计算结果绘制成图像如下。

（下转第215页）

（上接第157页）

图13 t=3时刻的解

图14 2-范数误差

2 结论

通过上述算例实验，可以看到四种差分格式都可抑制数值振荡，实现对流占优扩散方程的求解。在文中定义三种范数度量下，指数差分格式和修正中心差分格式误差最小，其次是萨马尔斯基格式，迎风格式效果最差。

参考文献：

[1] 陆金甫，关治. 偏微分方程数值解法[M].北京： 清华大学出版社，2003：168-170.

[2] 胡健伟，汤怀民. 微分方程数值方法[M].2版.北京：科学出版社，2007：131-138.

[3]Morton K W， Mayers D F， Numerical Solution of Partial Differential Equations[M].2th ed. Cambridge UK： Cambridge University Press，2005：19-26.

[4] 刘国庆， 傅冬生. 非线性奇异摄动问题一致收敛指数型拟合差分格式[J]. 高等学校计算数学学报，1998，10（2）： 173-184.

[5] 陈昊. 对流扩散方程差分格式[D]. 成都： 电子科技大学数学科学学院，2003：12-15.

[6] 王同科. 一维对流扩散方程CRANK-NICOLSON特征差分格式[J]. 应用数学，2001，14（4）： 55-60.

[7] 郑阿奇. MATLAB实用教程[M].北京： 电子工业出版社，2005：175-182.

[8] 彭芳麟. 数学物理方程的MATLAB解法与可视化[M].北京： 清华大学出版社，2004：156-160.