有限体积法求解一维非稳态对流扩散方程

1. 源起

最近想用有限差分法计算二维的顶盖驱动流，在推导过程中遇到了许多问题，例如对流项的离散、“线性化”这一在有限体积法中很常见的操作、1−δ1−δ压差在有限差分法中的实现、压力项的离散等等。以上问题让我感觉自己理论方面实在欠缺，所以还是从最简单的问题练手。

这个系列来自B站Up主“Red-Green鲤鱼”，视频链接

2. 一维非稳态对流扩散方程

∂f∂t+U∂f∂x=D∂2f∂x2

∂f∂t+U∂f∂x=D∂2f∂x2

2.1. 常系数：U和D为定值

方程特解：

f(x,t)=1√4t+1exp(−(x−1−Ut)2D(4t+1))

f(x,t)=14t+1−−−−−√exp(−(x−1−Ut)2D(4t+1))

2.2. 问题描述

用有限差分法求解上述方程，其中U=1.0,D=0.05,x∈[0,9]，并将数值解与t=2.5时的解析解对比

2.2.1. 初始条件

f(x,0)=exp(−(x−1)2D)

2.2.2. 边界条件

直接给出边界的值，为Dirichlet边界条件

f(0,t)=1√4t+1exp(−(1+Ut)2D(4t+1))

f(9,t)=1√4t+1exp(−(8−Ut)2D(4t+1))

3. 方程离散和区域离散

3.1. 区域离散

有限差分法中，将连续的区域离散为一系列结点，每个结点上含有物理量的值。假设结点间等距，距离为δ，那么离散方式可以分为内外两种：

1. 结点和边界重合，为外结点法
2. 结点和边界之间相差δ/2，为内结点法

以下将比较两种结点设置的方式。

1. 外结点法：设置N+1个结点，编号从0~N，结点间距δ=L/N，结点i的坐标为xi=i⋅δ
2. 内结点法：设置N个结点，编号从1~N，结点间距δ=L/N，结点i的坐标为xi=i⋅δ−δ/2

3.2. 方程离散

3.2.1. 非稳态项

对于结点i，采用一阶欧拉

∂f∂t|i=fn+1i−fniΔt

3.2.2. 对流项

一阶迎风，由于给出对流项前置系数U=1.0，因此迎风格式为

∂f∂x|i=fi−fi−1δ

3.2.3. 扩散项

二阶中心差分

∂2f∂x2|i=fi+1−2fi+fi−1δ2

注意：以上离散格式均是针对内结点。

3.3. 差分方程

将以上各项的离散形式组合起来得到差分方程。需要注意的是，如果对流项和扩散项中用到的物理量的值为下一时刻即n+1时间步的值，则为隐式离散；如果为当前时刻即n时间步的值，则为显式离散。

为了表示方便，下文中i用P表示，i−1用W表示，i+1用E表示。

3.3.1. 显式

fn+1P=aPfP+aWfW+aEfE

aP=1−UΔtδ−2DΔtδ2

aW=UΔtδ+DΔtδ2

aE=DΔtδ2

根据Harten定理系数非负原则，该显式格式稳定的条件是aP、aW、aE同时大于0，可见，这三项中只有aP存在小于0的可能，通过调整程序中δ和Δt能够发现这一现象。

3.3.2. 隐式

aPfn+1P=aWfn+1W+aEfn+1E+b

aP=1+UΔtδ+2DΔtδ2

aW=UΔtδ+DΔtδ2

aE=DΔtδ2

b=fnP

3.4. 边界离散

3.4.1. 外结点法

结点编号从0到N，共N+1个结点，对于结点i∈[1,N−1]，均采用上述离散方程求解，对于i=0,N的两个结点，由Dirichlet边界条件直接获得结点上的值。

3.4.2. 内结点法

结点编号从1到N，共N个结点，结点i∈[2,N−1]按照上述离散方程求解。结点i=1,N的对流项以及扩散项分别为

3.4.2.1. 对流项

∂f∂x|i=1=f1−fx=0δ/2,∂f∂x|i=N=fN−fN−1δ

3.4.2.2. 扩散项：

对于左边界：

fx=0=f0=f1−δ2f′+δ28f′′+ε(δ3)f2=f1+δf′+δ22f′′+ε(δ3)

消去f′可得到二阶导数

∂2f∂x2|i=1=2f0+f2−3f1δ2(3/4)+ε(δ)

注意：上式只有一阶精度。同理，右边界按类似的方式处理，令右边界为N+1号结点，则结点N

∂2f∂x2|i=N=2fN+1+fN−1−3fNδ2(3/4)+ε(δ)

使用一阶精度扩散项的显式算法：

fn+11=f1−2UΔtδ(f1−f0)+4DΔt3δ2(2f0+f2−3f1)

fn+1N=fN−UΔtδ(fN−fN−1)+4DΔt3δ2(2fN+1+fN−1−3fN)

使用一阶精度扩散项的隐式算法：

(1+2UΔtδ+4DΔtδ2)fn+11=4DΔt3δ2fn+12+(2UΔtδ+8DΔt3δ2)fn+10+fn1

(1+UΔtδ+4DΔtδ2)fn+1N=(UΔtδ+4DΔt3δ2)fn+1N−1+8DΔt3δ2fn+1N+1+fnN

上式中的fn+10以及fn+1N+1均为已知的边界条件。

3.4.2.3. 扩散项二阶精度的补充结点法

这种方法来自《数值传热学》例4-4

对于内结点法，i=1的结点和左边界距离为δ/2，我们在左边界上设置一个结点，编号为i=0，在左边界往左δ/2处补充结点i=−1，对于结点i=1

∂2f∂x2|i=1=f2+f−1−2f1δ2+ε(δ2)

对于结点i=0

∂2f∂x2|i=0=f1+f−1−2f0(δ/2)2+ε(δ2)

由于结点0处于边界，在满足边界条件的同时应尽量让该结点满足控制方程，所以对结点0的控制方程的显式离散为

fn+10−fn0Δt+Uf0−f−1δ/2=Df1+f−1−2f0(δ/2)2

fn+10以及fn0为Dirichlet边界条件，均已知，因此可求出f−1，将f−1代入结点i=1的扩散项差分方程中，最终可得到满足空间二阶精度的扩散项差分方程。

4. 代数方程组求解方法

显式离散能够直接推导出物理量的更新方程，因此只要给出Δx和Δt便能求解，但因为稳定性条件（如Harten定理）限制了所能采用的Δt

隐式离散的稳定性更好，但是需要联立求解方程组。空间被离散为N个结点，每个结点都对应一个方程，即N个方程N个未知数。方程组求解的快慢直接决定了整个计算流程的效率，因此研究者开发了若干方程组求解方法，可分为两大类，直接求解和迭代求解

4.1. 直接求解

当求解的未知数个数极少时，采用Cramer法则直接求解。本文讨论的一维问题中，每个结点只和周围两个结点相关，即aP只和aW以及aE相关，因此，求解的代数方程组只有三个对角上的元素非零，即三对角矩阵。下面介绍针对该矩阵发展的算法，基于Gauss消元法的Thomas算法：TDMA

4.1.1. TDMA: Tridiagonal matrix method

为了节省篇幅，以下只讨论区域离散方式为外结点法的代数方程组求解。事实上有限体积法的离散方式更像内结点法，因为体心值作为网格平均值有二阶精度（参考链接）。

矩阵形式为AX=B，由于使用了第一类边界条件，x1和xN已知，因此a12=aN,N−1=0

[a11a12a21a22a23⋯ai,i−1aiiai,i+1⋯aN−1,N−2aN−1,N−1aN−1,NaN,N−1aN,N][x1x2⋯xi⋯xN−1xN]=[b1b2⋯bi⋯bN−1bN]

TDMA算法包含两个步骤，消元和回代。消元指的是通过将上一行整体加到下一行从而消去下一行最左边的元素，从矩阵形状上来看就是将三对角矩阵变成了上对角矩阵，此做法的目的是将最后一行的aN,N−1消去从而能直接计算出xN，得到xN后便逐一向上回代解出其他未知数。

为了讨论方便，把通式写成AiTi=BiTi+1+CiTi−1+Di，消元的目的是把该式化为Ti−1=Pi−1Ti+Qi−1，通过整理可以得到系数Pi,Qi和Bi,Ci,Di之间的关系，

Pi=BiAi−CiPi−1,Qi=Di+CiQi−1Ai−CiPi−1

可以看出，Pi和Qi都需要递归求解

P1=B1A1,Q1=D1A1,QN=TN

综上，TDMA计算流程为

1. 计算P1,Q1
2. 更新Pi,Qi，得到PN,QN，其中QN=TN
3. 由Ti−1的表达式计算T1−TN−1

4.2. 迭代求解

待补充

5. 计算结果讨论

5.1. 显式格式

由上述离散过程可见，对流项采用一阶迎风。一阶格式的结果是截断项首项为二阶导数，而二阶导数项有扩散性质，即物理量会向各个方向传播。通常这种由二阶导数引起的扩散现象称为假扩散/人工粘性/数值粘性。具体讨论可参考《数值传热学》第二版5.5节。

下图为显式离散在t=2.5s的计算结果，区域离散方式为外结点和内结点的差别不大。结点间距Δx=0.009，时间步Δt=0.0001s，该步长满足Harten定理。可见，模拟结果与解析解相比，有矮胖特征，意味着存在假扩散现象。

5.2 隐式格式

下图为隐式离散在t=2.5s的计算结果，区域离散采用外结点，结点间距0.009，时间步Δt=0.001s，而在相同时间步下，显式格式的计算结果会发散，比较来看，隐式格式能够适用更大的时间步长。