对流扩散方程的非振荡守恒特征差分法

Non Oscillation Conservative Characteristic Difference Method for Solving Convection Diffusion Equations

摘要: 当求解对流占优扩散问题时，采用传统的二次拉格朗日插值的特征差分法，会出现数值振荡，且不满足质量守恒。结合算子分裂，本质非震荡和MMOCAA质量校正，提出了求解对流扩散方程的非震荡的守恒特征差分法。首先采用局部一维(LOD)法把一个二维偏微分方程分裂成x方向和y方向的两个一维的偏微分方程组；其次在每个方向上利用二阶本质非振荡和MMOCAA格式进行数值计算。数值实验验证格式满足非震荡和质量守恒，能够有效地解决大型对流占优扩散问题。

Abstract: The characteristic difference method based on quadratic Lagrange interpolation is used to solve the convection dominated diffusion problem, but there will be a larger numerical oscillation and not conservative. Combining the operator splitting, non-oscillatory and mass correction methods, the non-oscillation conservative characteristic difference methods are proposed to solve the convection diffusion equations. Firstly, the partial differential equations in two dimensions are splitting into two one-dimensional partial differential equations along the x-direction and the y-direction, respectively. Secondly, the second-order essentially non-oscillation and MMOCAA schemes are presented to compute the equations. By the numerical results, it shows that the scheme not only meets non-oscillatory and mass conservation, but also effectively solves the convection-dominant diffusion problems.

1. 引言

对流扩散方程 [1] [2] [3] 是一种常见的用来描述海洋、大气等环境和化学中传热、质量输运的等运动现象。众所周知，对流占优导致对流扩散方程呈双曲特性，采用标准有限差分法会导致数值震荡和数值弥散。为有效抑制数值震荡和数值弥散现象，国内外学者提出了求解对流扩散方程的迎风 [4] [5] [6] 和特征格式。Douglas [2] 利用二次Lagrange线性插值和差分法，提出了二阶修正特征线法(MMOC)法，并给出了收敛性分析。为改进二次Lagrange插值带来的数值震荡，陆金甫 [7] 、由同顺 [8] 和张鲁明 [9] 等人分别对特征差分法进行了修正和改进。由于MMOC特征差分不守恒，Douglas [10] 提出了修正对流的质量守恒的MMOCAA差分格式，芮红兴 [11] 提出了质量守恒的特征有限元格式，梁栋 [12] 提出了守恒的高精度质量守恒的特征差分法。

考虑到实际研究问题的规模大，计算时间长和计算存储量大的问题，算子分裂和区域分解算法 [13] [14] 被广泛应用到求解抛物型问题中。算子分裂把一个高维问题分解成多个一维问题，进而可以把一个十分复杂且包含多个变量的方程转化为包含一个或者较少的自变量的多个方程组。柳忠伟 [15] 结合LOD算子分裂、特征差分法和二阶本质非振荡法求解对流方程，但此算法不满足质量守恒。

结合算子分裂，二阶本质非震荡和质量校正法，本文提出了求解对流扩散方程的守恒特征法。在每个时间步，算法分三步，第一步，利用LOD算子分裂将方程分解成沿x方向和y方向的两个一维的偏微分方程组。第二步，在每个方向上，结合二阶本质非振荡和特征差分法，构造特征差分法。第三步，引入保守的MMOCAA格式，提出本质非振荡守恒特征差分格式。最后，通过数值实验验证格式质量守恒。

2. 模型

2.1. 对流扩散方程

对于伴有物质输运、气溶胶运动和分子扩散的物理过程以及粘性流体的流动等流体运动现象，都可用如下的数学模型表述。

$\{\begin{array}{l}\frac{\partial c}{\partial t}+\nabla \left(uc-D\nabla c\right)=f,\left(x,y\right)\in \Omega ,t\in \left[0,T\right],\\ c\left(x,y,t\right)=0,\left(x,y\right)\in \partial \Omega ,t\in \left[0,T\right],\\ c_0\left(x,y\right)=c_0,\left(x,y\right)\in \Omega .\end{array}$ (2.1)

其中有限区域 $c=c\left(x,y,t\right)$ 是流体的浓度， $\partial \Omega$ 为区域的边界，其中D， $u$ 分别表示扩散系数和速度场，不失一般性，假设它们均是常数，即 $u=\left(u_1,u_2\right),D=\left(D_1,D_2\right)$ 。

2.2. LOD算子分裂

局部一维格式

利用LOD算子分裂，将对流扩散方程分别沿x方向和y方向进行分解，得到如下的方程组：

$\{\begin{array}{l}\frac{\partial c}{\partial t}+\frac{\partial }{\partial x}\left(u_{1}c-D_1\frac{\partial c}{\partial x}\right)=f\\ \frac{\partial c}{\partial t}+\frac{\partial }{\partial y}\left(u_{2}c-D_2\frac{\partial c}{\partial y}\right)=0\end{array}$ (2.2)

分裂后的两个一维偏微分方程，若利用时间一阶格式和空间二阶修正迎风格式离散，截断误差可以达到 $O\left(h^2+\Delta t\right)$ ，且每个方向的方程可用追赶法求解。考虑到二阶修正迎风格式尽管可避免数值震荡，但在一定程度上带来了数值弥散，尤其在高对流占优扩散问题时，弥散现象变得更加明显。接下来，将结合特征差分和二阶本质非震荡法，解决对流扩散问题。

3. 算法

为了有效解决对流占优问题，第一步，分别对每个方向的方程采用特征线法处理，构造本质非振荡特征差分格式；第二步，加入质量校正项，引入保守的MMOCAA格式，得到校正后的本质非振荡守恒特征差分格式，并表示为矩阵形式；第三步，对于方程右端的未知项用本质非振荡插值求解，使其转化为已知项；第四步，利用追赶法 [16] 求解对角占优的三对角矩阵。

3.1. 特征线法

Step 1：沿x方向

$\frac{\partial c}{\partial t}+\frac{\partial }{\partial x}\left(u_{1}c-D_1\frac{\partial c}{\partial x}\right)=f,$ (3.1)

采用特征线法

$\frac{\partial c}{\partial {\tau }_1}=\frac{\partial c}{\partial t}\cos \alpha +\frac{\partial c}{\partial x}\cos \beta .$ (3.2)

${\tau }_1$ 为x方向特征线， $\beta$ 和 $\alpha$ 分别为特征线与x和t方向的夹角，则

$\{\begin{array}{l}\cos \alpha =\frac{1}{\sqrt{1+u_1^2}}\\ \cos \beta =\frac{u_1}{\sqrt{1+u_1^2}}\end{array}$

代入(3.2)得

$\frac{\partial c}{\partial t}\frac{1}{\sqrt{1+u_1^2}}+\frac{\partial c}{\partial x}\frac{u_1}{\sqrt{1+u_1^2}}=\frac{\partial c}{\partial {\tau }_1}.$

方程两边同乘 $\sqrt{1+u_1^2}$ ，代入(3.2)式，得

$\frac{\partial c}{\partial {\tau }_1}\sqrt{1+u_1^2}-D_1\frac{{\partial }^{2}c}{\partial x^2}=f,$

离散上式，得

$\frac{C_{ij}^{n-\frac{1}{2}}-{\bar{C}}_{ij}^{n-1}}{\Delta t}-D_1\frac{C_{i\text{+}1，j}^{n-\frac{1}{2}}-2C_{ij}^{n-\frac{1}{2}}+C_{i-1,j}^{n-\frac{1}{2}}}{h^2}=f_{ij}.$ (3.3)

$\left\{C_{ij}\right\}$ 表示 $\left\{c_{ij}\right\}$ 的数值解， $\left\{{\bar{C}}_{ij}\right\}$ 将在后面给出。

Step 2：沿y方向

$\frac{\partial c}{\partial t}+\frac{\partial }{\partial y}\left(u_{2}c-D_2\frac{\partial c}{\partial y}\right)=0,$

令

$\frac{\partial c}{\partial {\tau }_2}=\frac{\partial c}{\partial t}\cos \varphi +\frac{\partial c}{\partial y}\cos \phi ,$

其中 ${\tau }_2$ 是y方向的特征线， $\phi$ ， $\varphi$ 分别为特征线与y，t方向夹角，离散化，得

$\frac{C_{ij}^n-{\bar{C}}_{ij}^{n-\frac{1}{2}}}{\Delta t}-D_2\frac{C_{i,j+1}^n-2C_{ij}^n+C_{i,j-1}^n}{h^2}=0,$

其中， $\left\{C_{ij}\right\}$ 表示 $\left\{c_{ij}\right\}$ 的数值解， $\left\{{\bar{C}}_{ij}^{n-1}\right\}$ 将在后面给出。而此时数值格式质量不守恒。

为解决格式质量不守恒问题，结合MMOCAA格式，现提出守恒的特征差分格式。

(1-a) 沿x方向，令

${\bar{C}}_{ij}^{n-1}=I_1{C}_{ij}^{n-1}\left({\bar{x}}_i\right),{\lambda }_1=-\frac{u_1\Delta t}{h},{\bar{x}}_i=x_i-u_1\Delta t.$

不失一般性，我们假定 $\Delta t=O\left(h^2\right)$ ，即点 ${\bar{x}}_i$ 落在 $x_{i-1}$ 和 $x_{i+1}$ 之间。

(1-b) 引入三个质量通量： $t=t^{n-1}$ 时的初始质量通量

$R^{n-1}={\displaystyle {\sum }_{i=1}^N{\displaystyle {\sum }_{j=1}^N{C}_{ij}^{n-1}{h}^2}}.$

中间质量通量

${\bar{R}}^{n-1}={\displaystyle {\sum }_{i=1}^N{\displaystyle {\sum }_{j=1}^N{\bar{C}}_{ij}^{n-1}{h}^2}}.$

改进的质量通量

${\tilde{R}}^{n-1}={\displaystyle {\sum }_{i=1}^N{\displaystyle {\sum }_{j=1}^N{\tilde{C}}_{ij}^{n-1}{h}^2}}.$

其中，

${\tilde{C}}_{ij}^{n-1}=\{\begin{array}{l}\max \left(I_1{C}_{ij}^{n-1}\left({\tilde{x}}_i^{+}\right),I_1{C}_{ij}^{n-1}\left({\tilde{x}}_i^{-}\right)\right),R^{n-1}>{\bar{R}}^{n-1},\\ \min \left(I_1{C}_{ij}^{n-1}\left({\tilde{x}}_i^{+}\right),I_1{C}_{ij}^{n-1}\left({\tilde{x}}_i^{-}\right)\right),R^{n-1}\le {\bar{R}}^{n-1},\end{array}$

${\tilde{x}}_i^{+}={\bar{x}}_i+{\delta }_1$ ， ${\tilde{x}}_i^{-}={\bar{x}}_i-{\delta }_1$ ， ${\delta }_1=r_1{u}_1{\left(\Delta t\right)}^2$ ， $r_1=\frac{D_1\Delta t}{h^2}$ 且 $0<r_1<2$ 。

(1-c) 质量扰动参数 ${\theta }_1$ 的确定

如果 ${\bar{R}}^{n-1}={\tilde{R}}^{n-1}\ne R^{n-1}$ ，在这个时间步长内，误差是允许接收的。当 ${\bar{R}}^{n-1}\ne {\tilde{R}}^{n-1}$ 时， ${\theta }_1$ 被确定为

$R^{n-1}={\theta }_1{\bar{R}}^{n-1}+\left(1-{\theta }_1\right){\tilde{R}}^{n-1}.$

(1-d) 定义 $\left\{{\hat{C}}_{ij}^{n-1}\right\}$

${\hat{C}}_{ij}^{n-1}={\theta }_1{\bar{C}}_{ij}^{n-1}+\left(1-{\theta }_1\right){\tilde{C}}_{ij}^{n-1}.$

它是容易验证格式沿x方向满足质量守恒

(1-e) 计算 $\left\{C_{ij}^{n-\frac{1}{2}}\right\}$ 。

加入质量校正后，(3.3)式变成

$\frac{C_{ij}^{n-\frac{1}{2}}-{\hat{C}}_{ij}^{n-1}}{\Delta t}-D_1\frac{C_{i+1,j}^{n-\frac{1}{2}}-2C_{ij}^{n-\frac{1}{2}}+C_{i-1,j}^{n-\frac{1}{2}}}{h^2}=f_{ij}.$ (3.4)

两边同乘 $\Delta t$ ，令 $r_1=\frac{D_1\Delta t}{h^2}$ 则

$C_{ij}^{n-\frac{1}{2}}-{\hat{C}}_{ij}^{n-1}-r_1\left(C_{i+1,j}^{n-\frac{1}{2}}-2C_{ij}^{n-\frac{1}{2}}+C_{i-1,j}^{n-\frac{1}{2}}\right)=f_{ij}\Delta t.$

化简得

$-r_1{C}_{i+1,j}^{n-\frac{1}{2}}+\left(2r_1+1\right)C_{ij}^{n-\frac{1}{2}}-r_1{C}_{i-1,j}^{n-\frac{1}{2}}=f_{ij}\Delta t+{\hat{C}}_{ij}^{n-1}.$ (3.5)

$\left(\begin{array}{ccccc}2r_1+1& -r_1& & & \\ -r_1& 2r_1+1& -r_1& & \\ & \ddots & \ddots & \ddots & \\ & & \ddots & \ddots & -r_1\\ & & & -r_1& 2r_1+1\end{array}\right)\left(\begin{array}{c}{c}_1\\ c_2\\ ⋮\\ c_{N-1}\\ c_N\end{array}\right)=\left(\begin{array}{c}{F}_1\\ F_2\\ ⋮\\ F_{N-1}\\ F_N\end{array}\right)$

记为

$AC=F$

其中系数矩

$A=\left(\begin{array}{ccccc}2r_1+1& -r_1& & & \\ -r_1& 2r_1+1& -r_1& & \\ & \ddots & \ddots & \ddots & \\ & & \ddots & \ddots & -r_1\\ & & & -r_1& 2r_1+1\end{array}\right)$

因为 $2r_1+1>\left|r_1\right|+\left|r_1\right|$ ，所以系数矩阵A是一个严格对角占优的三对角矩阵。

(2-a) 沿y方向，

${\bar{C}}_{ij}^{n-\frac{1}{2}}=I_2{C}_{ij}^{n-\frac{1}{2}}\left({\bar{y}}_j\right),{\lambda }_2=-\frac{u_2\Delta t}{h},$

假定 $\Delta t=O\left(h^2\right),{\bar{y}}_j\in \left[y_{j-1},y_{j+1}\right]$ 。

(2-b) 质量通量。 $t=t^{n-\frac{1}{2}}$ 时的质量通量

$R^{n-\frac{1}{2}}={\displaystyle {\sum }_{i=1}^N{\displaystyle {\sum }_{j=1}^N{C}_{ij}^{n-\frac{1}{2}}{h}^2}}.$

中间质量通量

${\bar{R}}^{n-\frac{1}{2}}={\displaystyle {\sum }_{i=1}^N{\displaystyle {\sum }_{j=1}^N{\bar{C}}_{ij}^{n-\frac{1}{2}}{h}^2}}.$

改进的质量通量

${\tilde{R}}^{n-\frac{1}{2}}={\displaystyle {\sum }_{i=1}^N{\displaystyle {\sum }_{j=1}^N{\tilde{C}}_{ij}^{n-\frac{1}{2}}{h}^2}}.$

其中，

${\tilde{C}}_{ij}^{n-\frac{1}{2}}=\{\begin{array}{l}\max \left(I_2{C}_{ij}^{n-\frac{1}{2}}\left({\tilde{y}}_j^{+}\right),I_2{C}_{ij}^{n-\frac{1}{2}}\left({\tilde{y}}_j^{-}\right)\right),R^{n-\frac{1}{2}}>{\bar{R}}^{n-\frac{1}{2}},\\ \min \left(I_2{C}_{ij}^{n-\frac{1}{2}}\left({\tilde{y}}_j^{+}\right),I_2{C}_{ij}^{n-\frac{1}{2}}\left({\tilde{y}}_j^{-}\right)\right),R^{n-\frac{1}{2}}\le {\bar{R}}^{n-\frac{1}{2}},\end{array}$

${\tilde{y}}_j^{+}={\bar{y}}_j+{\eta }_2$ ， ${\tilde{y}}_j^{-}={\bar{y}}_j-{\eta }_2$ 其中， ${\eta }_2=r_2{u}_2{\left(\Delta t\right)}^2$ 。

(2-c) ${\theta }_2$ 的确定

如果 ${\bar{R}}^{n-\frac{1}{2}}={\tilde{R}}^{n-\frac{1}{2}}\ne R^{n-\frac{1}{2}}$ ，误差可以接收。

当 ${\bar{R}}^{n-\frac{1}{2}}\ne {\tilde{R}}^{n-\frac{1}{2}}$ 时， ${\theta }_2$ 被确定为

$R^{n-\frac{1}{2}}={\theta }_2{\bar{R}}^{n-\frac{1}{2}}+\left(1-{\theta }_2\right){\tilde{R}}^{n-\frac{1}{2}}.$

(2-d) 计算 $\left\{C_{ij}^n\right\}$ 。定义 $\left\{{\hat{C}}_{ij}^{n-\frac{1}{2}}\right\}$ 满足

${\hat{C}}_{ij}^{n-\frac{1}{2}}={\theta }_2{\bar{C}}_{ij}^{n-\frac{1}{2}}+\left(1-{\theta }_2\right){\tilde{C}}_{ij}^{n-\frac{1}{2}}.$

(3.4)式加入质量校正后，得

$\frac{C_{ij}^n-{\hat{C}}_{ij}^{n-\frac{1}{2}}}{\Delta t}-D_2\frac{C_{i,j+1}^n-2C_{ij}^n+C_{i,j-1}^n}{h^2}=0,$ (3.6)

令 $r_2=\frac{D_2\Delta t}{h^2}$ ，得

$-r_2{C}_{ij}^n+\left(2r_2+1\right)C_{ij}^n-r_2{C}_{i,j-1}^n={\hat{C}}_{ij}^{n-\frac{1}{2}}.$ (3.7)

表示成矩阵形式，即

$BC=F$

由于A，B均为严格对角占优三对角矩阵，故采用追赶法求解。

3.2. 本质非振荡插值

本部分给出 $\left\{{\bar{C}}_{ij}^{n-1},{\tilde{C}}_{ij}^{n-1},{\bar{C}}_{ij}^{n-\frac{1}{2}},{\tilde{C}}_{ij}^{n-\frac{1}{2}}\right\}$ 的计算方法。为了避免数值震荡，欲通过比较两个二阶差商的绝对值，选择绝对值相对较小的方法给出(文献 [8] )。

1) 当 $u_1$ 为正数时，

(1-1) 当 $\left|C_{i-2,j}^{n-1}-2C_{i-1,j}^{n-1}+C_{ij}^{n-1}\right|>\left|C_{i+1,j}^{n-1}-2C_{ij}^{n-1}+C_{i-1,j}^{n-1}\right|$ 时，利用牛顿插值得

$\begin{array}{c}{\bar{C}}_{ij}^{n-1}=C_{ij}^{n-1}+\frac{C_{ij}^{n-1}-C_{i-1,j}^{n-1}}{h}\left({\bar{x}}_i-x_i\right)+\frac{C_{i+1,j}^{n-1}-2C_{ij}^{n-1}+C_{i-1,j}^{n-1}}{2h^2}\left({\bar{x}}_i-x_{i-1}\right)\left({\bar{x}}_i-x_i\right)\\ =\frac{1}{2h^2}\left({\bar{x}}_i-x_{i-1}\right)\left({\bar{x}}_i-x_i\right)C_{i+1,j}^{n-1}+\left[1+\frac{1}{h}\left({\bar{x}}_i-x_i\right)-\frac{1}{h^2}\left({\bar{x}}_i-x_{i-1}\right)\left({\bar{x}}_i-x_i\right)\right]C_{ij}^{n-1}\\ +\left[-\frac{1}{h}\left({\bar{x}}_i-x_i\right)+\frac{1}{2h^2}\left({\bar{x}}_i-x_{i-1}\right)\left({\bar{x}}_i-x_i\right)\right]C_{i-1,j}^{n-1}\\ =\frac{\left(h-u_1\Delta t\right)\left(-u_1\Delta t\right)}{2h^2}{C}_{i+1,j}^{n-1}+\left(1-\frac{1}{h}{u}_1\Delta t-\frac{1}{h^2}\left(h-u_1\Delta t\right)\left(-u_1\Delta t\right)\right)C_{ij}^{n-1}\\ +\left(\frac{1}{h}{u}_1\Delta t+\frac{\left(h-u_1\Delta t\right)\left(-u_1\Delta t\right)}{2h^2}\right)C_{i-1,j}^{n-1}\end{array}$ (3.8)

同理

${\tilde{C}}_{ij}^{n-1}=C_{ij}^{n-1}+\frac{C_{ij}^{n-1}-C_{i-1,j}^{n-1}}{h}\left({\tilde{x}}_i-x_i\right)+\frac{C_{i+1,j}^{n-1}-2C_{ij}^{n-1}+C_{i-1,j}^{n-1}}{2h^2}\left({\tilde{x}}_i-x_{i-1}\right)\left({\tilde{x}}_i-x_i\right).$

令

$\begin{array}{l}{k}_1=\frac{1}{h^2}\left(h-u_1\Delta t+{\delta }_1\right)\left(-u_1\Delta t+{\delta }_1\right),k_2=\frac{1}{h}\left(-u_1\Delta t+{\delta }_1\right),\\ k_3=\frac{1}{h^2}\left(h-u_1\Delta t-{\delta }_1\right)\left(-u_1\Delta t-{\delta }_1\right),k_4=\frac{1}{h}\left(-u_1\Delta t-{\delta }_1\right).\end{array}$

由上式计算可得，当 ${\tilde{x}}_i^{+}={\bar{x}}_i+{\delta }_1$ 时，

${\tilde{C}}_{ij}^{n-1,+}=\frac{k_1}{2}{C}_{i+1,j}^{n-1}+\left(1+k_2-k_1\right)C_{ij}^{n-1}+\left(-k_2+\frac{k_1}{2}\right)C_{i-1,j}^{n-1}.$

当 ${\tilde{x}}_i^{-}={\bar{x}}_i-{\delta }_1$ 时，

${\tilde{C}}_{ij}^{n-1,-}=\frac{k_3}{2}{C}_{i+1,j}^{n-1}+\left(1+k_4-k_3\right)C_{ij}^{n-1}+\left(-k_4+\frac{k_3}{2}\right)C_{i-1,j}^{n-1}.$

当 $R^{n-1}>{\bar{R}}^{n-1}$ 时， ${\tilde{C}}_{ij}^{n-1}=\max \left\{{\tilde{C}}_{ij}^{n-1,+},{\tilde{C}}_{ij}^{n-1,-}\right\}$ ， $R^{n-1}\le {\bar{R}}^{n-1}$ 时， ${\tilde{C}}_{ij}^{n-1}=\min \left\{{\tilde{C}}_{ij}^{n-1,+},{\tilde{C}}_{ij}^{n-1,-}\right\}$ 。

(1-2) 当 $\left|C_{i-2,j}^{n-1}-2C_{i-1,j}^{n-1}+C_{ij}^{n-1}\right|\le \left|C_{i+1,j}^{n-1}-2C_{ij}^{n-1}+C_{i-1,j}^{n-1}\right|$ 时，得

$\begin{array}{c}{\bar{C}}_{ij}^{n-1}=C_{ij}^{n-1}+\frac{C_{ij}^{n-1}-C_{i-1,j}^{n-1}}{h}\left({\bar{x}}_i-x_i\right)+\frac{C_{ij}^{n-1}-2C_{i-1,j}^{n-1}+C_{i-2,j}^{n-1}}{2h^2}\left({\bar{x}}_i-x_{i-1}\right)\left({\bar{x}}_i-x_i\right)\\ =\left[1+\frac{1}{h}\left({\bar{x}}_i-x_i\right)+\frac{1}{2h^2}\left({\bar{x}}_i-x_{i-1}\right)\left({\bar{x}}_i-x_i\right)\right]C_{ij}^{n-1}+\frac{1}{2h^2}\left({\bar{x}}_i-x_{i-1}\right)\left({\bar{x}}_i-x_i\right)C_{i-2,j}^{n-1}\\ -\left[\frac{1}{h}\left({\bar{x}}_i-x_i\right)+\frac{1}{h^2}\left({\bar{x}}_i-x_{i-1}\right)\left({\bar{x}}_i-x_i\right)\right]C_{i-1,j}^{n-1}\\ =\left(1-\frac{1}{h}{u}_1\Delta t+\frac{1}{2h^2}\left(h-u_1\Delta t\right)\left(-u_1\Delta t\right)\right)C_{ij}^{n-1}+\frac{\left(h-u_1\Delta t\right)\left(-u_1\Delta t\right)}{2h^2}{C}_{i-2,j}^{n-1}\\ -\left(-\frac{1}{h}{u}_1\Delta t+\frac{\left(h-u_1\Delta t\right)\left(-u_1\Delta t\right)}{h^2}\right)C_{i-1,j}^{n-1}\end{array}$ (3.9)

${\tilde{C}}_{ij}^{n-1}=C_{ij}^{n-1}+\frac{C_{ij}^{n-1}-C_{i-1,j}^{n-1}}{h}\left({\tilde{x}}_i-x_i\right)+\frac{C_{ij}^{n-1}-2C_{i-1,j}^{n-1}+C_{i-2,j}^{n-1}}{2h^2}\left({\tilde{x}}_i-x_{i-1}\right)\left({\tilde{x}}_i-x_i\right)$

由上式计算可得，当 ${\tilde{x}}_i^{+}={\bar{x}}_i+{\delta }_1$ 时，

${\tilde{C}}_{ij}^{n-1,+}=\left(1+k_2+\frac{k_1}{2}\right)C_{ij}^{n-1}+\frac{k_1}{2}{C}_{i-2,j}^{n-1}-\left(k_1+k_2\right)C_{i-1,j}^{n-1}.$

当 ${\tilde{x}}_i^{-}={\bar{x}}_i-{\delta }_1$ 时，

${\tilde{C}}_{ij}^{n-1,-}=\left(1+k_4+\frac{k_3}{2}\right)C_{ij}^{n-1}+\frac{k_3}{2}{C}_{i-2,j}^{n-1}-\left(k_3+k_4\right)C_{i-1,j}^{n-1}.$

当 $R^{n-1}>{\bar{R}}^{n-1}$ 时， ${\tilde{C}}_{ij}^{n-1}=\max \left({\tilde{C}}_{ij}^{n-1,+},{\tilde{C}}_{ij}^{n-1,-}\right)$

当 $R^{n-1}\le {\bar{R}}^{n-1}$ 时， ${\tilde{C}}_{ij}^{n-1}=\min \left({\tilde{C}}_{ij}^{n-1,+},{\tilde{C}}_{ij}^{n-1,-}\right)$

2) 当 $u_1$ 为负数时，利用牛顿插值， $p_1=\frac{1}{h^2}\left(-h-u_1\Delta t+{\delta }_1\right)\left(-u_1\Delta t+{\delta }_1\right)$ ，

$p_2=\frac{1}{h}\left(-u_1\Delta t+{\delta }_1\right)$ ， $p_3=\frac{1}{h^2}\left(-h-u_1\Delta t-{\delta }_1\right)\left(-u_1\Delta t-{\delta }_1\right)$ ， $p_4=\frac{1}{h}\left(-u_1\Delta t-{\delta }_1\right)$

(2-1) $\left|C_{i-1,j}^{n-1}-2C_{ij}^{n-1}+C_{i+1,j}^{n-1}\right|>\left|C_{ij}^{n-1}-2C_{i+1,j}^{n-1}+C_{i+2,j}^{n-1}\right|$ 时，得

$\begin{array}{c}{\bar{C}}_{ij}^{n-1}=C_{ij}^{n-1}+\frac{C_{ij}^{n-1}-C_{i+1,j}^{n-1}}{-h}\left({\bar{x}}_i-x_i\right)+\frac{C_{ij}^{n-1}-2C_{i+1,j}^{n-1}+C_{i+2,j}^{n-1}}{2h^2}\left({\bar{x}}_i-x_i\right)\left({\bar{x}}_i-x_{i+1}\right)\\ =\frac{\left(-h-u_1\Delta t\right)\left(-u_1\Delta t\right)}{2h^2}{C}_{i+2,j}^{n-1}+\left(1+\frac{1}{h}{u}_1\Delta t+\frac{1}{2h^2}\left(-h-u_1\Delta t\right)\left(-u_1\Delta t\right)\right)C_{ij}^{n-1}\\ +\left(-\frac{1}{h}{u}_1\Delta t-\frac{\left(-h-u_1\Delta t\right)\left(-u_1\Delta t\right)}{h^2}\right)C_{i+1,j}^{n-1}\end{array}$

$\begin{array}{c}{\bar{C}}_{ij}^{n-1}=C_{ij}^{n-1}+\frac{C_{ij}^{n-1}-C_{i+1,j}^{n-1}}{-h}\left({\tilde{x}}_i-x_i\right)+\frac{C_{ij}^{n-1}-2C_{i+1,j}^{n-1}+C_{i+2,j}^{n-1}}{2h^2}\left({\tilde{x}}_i-x_i\right)\left({\tilde{x}}_i-x_{i+1}\right)\\ =\frac{1}{2h^2}\left({\tilde{x}}_i-x_i\right)\left({\tilde{x}}_i-x_{i+1}\right)C_{i+2,j}^{n-1}+\left[1-\frac{1}{h}\left({\tilde{x}}_i-x_i\right)+\frac{1}{2h^2}\left({\tilde{x}}_i-x_i\right)\left({\tilde{x}}_i-x_{i+1}\right)\right]C_{ij}^{n-1}\\ +\left[\frac{1}{h}\left({\tilde{x}}_i-x_i\right)-\frac{1}{h^2}\left({\tilde{x}}_i-x_i\right)\left({\tilde{x}}_i-x_{i+1}\right)\right]C_{i+1,j}^{n-1}\end{array}$

由上式计算可得，当 ${\tilde{x}}_i^{+}={\bar{x}}_i+{\delta }_1$ 时，

${\tilde{C}}_{ij}^{n-1,+}=\left(1+p_2+\frac{p_1}{2}\right)C_{ij}^{n-1}+\frac{p_1}{2}{C}_{i+2,j}^{n-1}+\left(p_2-p_1\right)C_{i+1,j}^{n-1}.$

当 ${\tilde{x}}_i^{-}={\bar{x}}_i-{\delta }_1$ 时，

${\tilde{C}}_{ij}^{n-1,-}=\left(1+p_4+\frac{p_3}{2}\right)C_{ij}^{n-1}+\frac{p_3}{2}{C}_{i+2,j}^{n-1}+\left(p_4-p_3\right)C_{i+1,j}^{n-1}.$

当 $R^{n-1}>{\bar{R}}^{n-1}$ 时， ${\tilde{C}}_{ij}^{n-1}=\max \left\{{\tilde{C}}_{ij}^{n-1,+},{\tilde{C}}_{ij}^{n-1,-}\right\}$ ， $R^{n-1}\le {\bar{R}}^{n-1}$ 时， ${\tilde{C}}_{ij}^{n-1}=\min \left\{{\tilde{C}}_{ij}^{n-1,+},{\tilde{C}}_{ij}^{n-1,-}\right\}$

(2-2) $\left|C_{i-1,j}^{n-1}-2C_{ij}^{n-1}+C_{i+1,j}^{n-1}\right|\le \left|C_{ij}^{n-1}-2C_{i+1,j}^{n-1}+C_{i+2,j}^{n-1}\right|$ 时，得

$\begin{array}{c}{\bar{C}}_{ij}^{n-1}=C_{ij}^{n-1}+\frac{C_{ij}^{n-1}-C_{i+1,j}^{n-1}}{-h}\left({\bar{x}}_i-x_i\right)+\frac{C_{i-1j}^{n-1}-2C_{i,j}^{n-1}+C_{i+1,j}^{n-1}}{2h^2}\left({\bar{x}}_i-x_i\right)\left({\bar{x}}_i-x_{i+1}\right)\\ =\left(1+\frac{1}{h}{u}_1\Delta t+\frac{1}{h^2}\left(-h-u_1\Delta t\right)\left(-u_1\Delta t\right)\right)C_{ij}^{n-1}+\frac{\left(-h-u_1\Delta t\right)\left(-u_1\Delta t\right)}{2h^2}{C}_{i-1,j}^{n-1}\\ +\left(-\frac{1}{h}{u}_1\Delta t+\frac{\left(-h-u_1\Delta t\right)\left(-u_1\Delta t\right)}{2h^2}\right)C_{i+1,j}^{n-1}\end{array}$

$\begin{array}{c}{\tilde{C}}_{ij}^{n-1}=C_{ij}^{n-1}+\frac{C_{ij}^{n-1}-C_{i+1,j}^{n-1}}{-h}\left({\tilde{x}}_i-x_i\right)+\frac{C_{i-1,j}^{n-1}-2C_{ij}^{n-1}+C_{i+1,j}^{n-1}}{2h^2}\left({\tilde{x}}_i-x_i\right)\left({\tilde{x}}_i-x_{i+1}\right)\\ =\left[1-\frac{1}{h}\left({\tilde{x}}_i-x_i\right)+\frac{1}{h^2}\left({\tilde{x}}_i-x_i\right)\left({\tilde{x}}_i-x_{i+1}\right)\right]C_{ij}^{n-1}+\frac{1}{2h^2}\left({\tilde{x}}_i-x_i\right)\left({\tilde{x}}_i-x_{i+1}\right)C_{i-1,j}^{n-1}\\ +\left[\frac{1}{h}\left({\tilde{x}}_i-x_i\right)+\frac{1}{2h^2}\left({\tilde{x}}_i-x_i\right)\left({\tilde{x}}_i-x_{i+1}\right)\right]C_{i+1,j}^{n-1}\end{array}$

计算可得，当 ${\tilde{x}}_i^{+}={\bar{x}}_i+{\delta }_1$ 时，

${\tilde{C}}_{ij}^{n-1,+}=\left(1-p_2+p_1\right)C_{ij}^{n-1}+\frac{p_1}{2}{C}_{i-1,j}^{n-1}+\left(p_2+\frac{p_1}{2}\right)C_{i+1,j}^{n-1}.$

当 ${\tilde{x}}_i^{-}={\bar{x}}_i-{\delta }_1$ 时，

${\tilde{C}}_{ij}^{n-1,-}=\left(1-p_4+p_3\right)C_{ij}^{n-1}+\frac{p_3}{2}{C}_{i-1,j}^{n-1}+\left(p_4+\frac{p_3}{2}\right)C_{i+1,j}^{n-1}.$

当 $R^{n-1}>{\bar{R}}^{n-1}$ 时， ${\tilde{C}}_{ij}^{n-1}=\max \left\{{\tilde{C}}_{ij}^{n-1,+},{\tilde{C}}_{ij}^{n-1,-}\right\}$ ， $R^{n-1}\le {\bar{R}}^{n-1}$ 时， ${\tilde{C}}_{ij}^{n-1}=\min \left\{{\tilde{C}}_{ij}^{n-1,+},{\tilde{C}}_{ij}^{n-1,-}\right\}$ 。

同理，可得沿y方向。

4. 数值实验

方程(1)中的速度场系数 $u_1=u_2=0.5$ ，扩散系数 $d_1=d_2=0.02$ ，有限区域 $\Omega =\left[0,5\right]\times \left[0,5\right]$ ，时间剖分 $\Delta t=\frac{1}{100}$ 和空间剖分 $h=\frac{1}{10}$ 。假定有初始时刻浓度只在 $\left(\frac{5}{3},\frac{5}{3}\right)$ ， $\left(\frac{5}{3},\frac{10}{3}\right)$ ， $\left(\frac{10}{3},\frac{5}{3}\right)$ 和 $\left(\frac{10}{3},\frac{10}{3}\right)$ 位置为1，其他位置浓度值为0。假定右端项 $f=0$ 。

图1~图6展示了浓度在不同时刻的表面图和等值线图，由于对流项的存在，浓度会随着时间，沿着对角线方向行进，同时，由于扩散项的存在，浓度会不断的减小，这与物理现象相吻合。因此，我们的格式可很好地描述对流扩散问题，而且我们空间步长比较大时( $h=0.1$ )，没有出现数值震荡问题。

验证格式的质量误差，令时间剖分 $\Delta t=0.01$ 和空间步长分别取 $h=0.1,0.05,0.02,0.01$ 。并与文献 [15] 进行比较。

Figure 1. T = 0.0 s time concentration surface map and contour map

图1. T = 0.0 s时刻浓度表面图和等值线图

Figure 2. T = 0.1 s time concentration surface map and contour map

图2. T = 0.1 s时刻浓度表面图和等值线图

Figure 3. T = 0.2 s time concentration surface map and contour map

图3. T = 0.2 s时刻浓度表面图和等值线图

Figure 4. T = 0.5 s time concentration surface map and contour map

图4. T = 0.5 s时刻浓度表面图和等值线图

Figure 5. T = 0.8 s time concentration surface map and contour map

图5. T = 0.8 s时刻浓度表面图和等值线图

Figure 6. T = 1.0 s time concentration surface map and contour map

图6. T = 1.0 s时刻浓度表面图和等值线图

Table 1. Mass errors at the different time

表1. 不同时刻的质量误差

由表1中数据，可知在加入质量校正前，格式质量之差数值比较大，当加入质量校正后，格式质量之差可以达到 ${10}^{-11}$ ，可以认为是质量守恒的。

基金项目

国家自然科学基金资助项目(61703250, 61503227)，山东自然科学基金资助项目(ZR2017BA029, ZR2017BF002)。

NOTES

^*通讯作者。