此文转载自本人在简书平台上发布的文章
1. 源起
最近想用有限差分法计算二维的顶盖驱动流，在推导过程中遇到了许多问题，例如对流项的离散、“线性化”这一在有限体积法中很常见的操作、 1 − δ 1-\delta 1−δ压差在有限差分法中的实现、压力项的离散等等。以上问题让我感觉自己理论方面实在欠缺，所以还是从最简单的问题练手。
这个系列来自B站Up主“Red-Green鲤鱼”，视频链接
2. 一维非稳态对流扩散方程
∂ f ∂ t + U ∂ f ∂ x = D ∂ 2 f ∂ x 2 \frac{\partial f}{\partial t} + U\frac{\partial f}{\partial x} = D \frac{\partial^2 f}{\partial x^2} ∂t∂f+U∂x∂f=D∂x2∂2f
2.1. 常系数：U和D为定值
方程特解：
f ( x , t ) = 1 4 t + 1 e x p ( − ( x − 1 − U t ) 2 D ( 4 t + 1 ) ) f(x,t)=\frac{1}{\sqrt{4t+1}}exp \left( -\frac{(x-1-Ut)^2}{D(4t+1)} \right) f(x,t)=4t+1 1exp(−D(4t+1)(x−1−Ut)2)
2.2. 问题描述
用有限差分法求解上述方程，其中 U = 1.0 , D = 0.05 , x ∈ [ 0 , 9 ] U=1.0, D=0.05, x\in[0,9] U=1.0,D=0.05,x∈[0,9]，并将数值解与 t = 2.5 t=2.5 t=2.5时的解析解对比
2.2.1. 初始条件
f ( x , 0 ) = e x p ( − ( x − 1 ) 2 D ) f(x,0)=exp \left( -\frac{(x-1)^2}{D} \right) f(x,0)=exp(−D(x−1)2)
2.2.2. 边界条件
直接给出边界的值，为Dirichlet边界条件
f ( 0 , t ) = 1 4 t + 1 e x p ( − ( 1 + U t ) 2 D ( 4 t + 1 ) ) f(0,t)=\frac{1}{\sqrt{4t+1}}exp \left( -\frac{(1+Ut)^2}{D(4t+1)} \right) f(0,t)=4t+1 1exp(−D(4t+1)(1+Ut)2)
f ( 9 , t ) = 1 4 t + 1 e x p ( − ( 8 − U t ) 2 D ( 4 t + 1 ) ) f(9,t)=\frac{1}{\sqrt{4t+1}}exp \left( -\frac{(8-Ut)^2}{D(4t+1)} \right) f(9,t)=4t+1 1exp(−D(4t+1)(8−Ut)2)
3. 方程离散和区域离散
3.1. 区域离散
有限差分法中，将连续的区域离散为一系列结点，每个结点上含有物理量的值。假设结点间等距，距离为 δ \delta δ，那么离散方式可以分为内外两种：
结点和边界重合，为外结点法结点和边界之间相差 δ / 2 \delta /2 δ/2，为内结点法
以下将比较两种结点设置的方式。
外结点法：设置 N + 1 N+1 N+1个结点，编号从0~N，结点间距 δ = L / N \delta = L/N δ=L/N，结点 i i i的坐标为 x i = i ⋅ δ x_i=i\cdot \delta xi=i⋅δ内结点法：设置 N N N个结点，编号从1~N，结点间距 δ = L / N \delta = L/N δ=L/N，结点 i i i的坐标为 x i = i ⋅ δ − δ / 2 x_i=i\cdot \delta - \delta /2 xi=i⋅δ−δ/2
3.2. 方程离散
3.2.1. 非稳态项
对于结点 i i i，采用一阶欧拉
∂ f ∂ t ∣ i = f i n + 1 − f i n Δ t \frac{\partial f}{\partial t}\bigg|_{i}=\frac{f^{n+1}_i-f^n_i}{\Delta t} ∂t∂f∣ ∣i=Δtfin+1−fin
3.2.2. 对流项
一阶迎风，由于给出对流项前置系数 U = 1.0 U=1.0 U=1.0，因此迎风格式为
∂ f ∂ x ∣ i = f i − f i − 1 δ \frac{\partial f}{\partial x}\bigg|_{i}=\frac{f_i-f_{i-1}}{\delta} ∂x∂f∣ ∣i=δfi−fi−1
3.2.3. 扩散项
二阶中心差分
∂ 2 f ∂ x 2 ∣ i = f i + 1 − 2 f i + f i − 1 δ 2 \frac{\partial^2 f}{\partial x^2}\bigg|_{i}=\frac{f_{i+1}-2f_i+f_{i-1}}{\delta ^2} ∂x2∂2f∣ ∣i=δ2fi+1−2fi+fi−1
注意：以上离散格式均是针对内结点。
3.3. 差分方程
将以上各项的离散形式组合起来得到差分方程。需要注意的是，如果对流项和扩散项中用到的物理量的值为下一时刻即 n + 1 n+1 n+1时间步的值，则为隐式离散；如果为当前时刻即 n n n时间步的值，则为显式离散。
为了表示方便，下文中 i i i用 P P P表示， i − 1 i-1 i−1用 W W W表示， i + 1 i+1 i+1用 E E E表示。
3.3.1. 显式
f P n + 1 = a P f P + a W f W + a E f E f^{n+1}_P=a_P f_P + a_W f_W + a_E f_E fPn+1=aPfP+aWfW+aEfE
a P = 1 − U Δ t δ − 2 D Δ t δ 2 a_P=1-\frac{U\Delta t}{\delta} -\frac{2D\Delta t}{\delta^2} aP=1−δUΔt−δ22DΔt
a W = U Δ t δ + D Δ t δ 2 a_W=\frac{U\Delta t}{\delta} + \frac{D\Delta t}{\delta^2} aW=δUΔt+δ2DΔt
a E = D Δ t δ 2 a_E=\frac{D\Delta t}{\delta^2} aE=δ2DΔt
根据Harten定理系数非负原则，该显式格式稳定的条件是 a P a_P aP、 a W a_W aW、 a E a_E aE同时大于0，可见，这三项中只有 a P a_P aP存在小于0的可能，通过调整程序中 δ \delta δ和 Δ t \Delta t Δt能够发现这一现象。
3.3.2. 隐式
a P f P n + 1 = a W f W n + 1 + a E f E n + 1 + b a_P f^{n+1}_P=a_W f^{n+1}_W + a_E f^{n+1}_E + b aPfPn+1=aWfWn+1+aEfEn+1+b
a P = 1 + U Δ t δ + 2 D Δ t δ 2 a_P=1+\frac{U\Delta t}{\delta}+\frac{2D\Delta t}{\delta^2} aP=1+δUΔt+δ22DΔt
a W = U Δ t δ + D Δ t δ 2 a_W=\frac{U\Delta t}{\delta} + \frac{D\Delta t}{\delta^2} aW=δUΔt+δ2DΔt
a E = D Δ t δ 2 a_E=\frac{D\Delta t}{\delta^2} aE=δ2DΔt
b = f P n b=f^n_P b=fPn
3.4. 边界离散
3.4.1. 外结点法
结点编号从0到N，共N+1个结点，对于结点 i ∈ [ 1 , N − 1 ] i\in [1,N-1] i∈[1,N−1]，均采用上述离散方程求解，对于 i = 0 , N i=0,N i=0,N的两个结点，由Dirichlet边界条件直接获得结点上的值。
3.4.2. 内结点法
结点编号从1到N，共N个结点，结点 i ∈ [ 2 , N − 1 ] i\in [2,N-1] i∈[2,N−1]按照上述离散方程求解。结点 i = 1 , N i=1,N i=1,N的对流项以及扩散项分别为
3.4.2.1. 对流项
∂ f ∂ x ∣ i = 1 = f 1 − f x = 0 δ / 2 , ∂ f ∂ x ∣ i = N = f N − f N − 1 δ \frac{\partial f}{\partial x}\bigg|_{i=1}=\frac{f_1-f_{x=0}}{\delta/2},\quad \frac{\partial f}{\partial x}\bigg|_{i=N}=\frac{f_N-f_{N-1}}{\delta} ∂x∂f∣ ∣i=1=δ/2f1−fx=0,∂x∂f∣ ∣i=N=δfN−fN−1
3.4.2.2. 扩散项：
对于左边界：
f x = 0 = f 0 = f 1 − δ 2 f ′ + δ 2 8 f ′ ′ + ε ( δ 3 ) f 2 = f 1 + δ f ′ + δ 2 2 f ′ ′ + ε ( δ 3 ) f_{x=0}=f_0=f_1-\frac{\delta}{2}f^{\prime}+\frac{\delta^2}{8}f^{\prime\prime} + \varepsilon(\delta^3) \\ f_2=f_1+\delta f^{\prime} + \frac{\delta^2}{2} f^{\prime\prime} +\varepsilon(\delta^3) \\ fx=0=f0=f1−2δf′+8δ2f′′+ε(δ3)f2=f1+δf′+2δ2f′′+ε(δ3)
消去 f ′ f^{\prime} f′可得到二阶导数
∂ 2 f ∂ x 2 ∣ i = 1 = 2 f 0 + f 2 − 3 f 1 δ 2 ( 3 / 4 ) + ε ( δ ) \frac{\partial^2 f}{\partial x^2}\bigg|_{i=1}=\frac{2f_0+f_{2}-3f_1}{\delta ^2 (3/4)}+\varepsilon(\delta) ∂x2∂2f∣ ∣i=1=δ2(3/4)2f0+f2−3f1+ε(δ)
注意：上式只有一阶精度。同理，右边界按类似的方式处理，令右边界为 N + 1 N+1 N+1号结点，则结点 N N N
∂ 2 f ∂ x 2 ∣ i = N = 2 f N + 1 + f N − 1 − 3 f N δ 2 ( 3 / 4 ) + ε ( δ ) \frac{\partial^2 f}{\partial x^2}\bigg|_{i=N}=\frac{2f_{N+1}+f_{N-1}-3f_N}{\delta ^2 (3/4)}+\varepsilon(\delta) ∂x2∂2f∣ ∣i=N=δ2(3/4)2fN+1+fN−1−3fN+ε(δ)
使用一阶精度扩散项的显式算法：
f 1 n + 1 = f 1 − 2 U Δ t δ ( f 1 − f 0 ) + 4 D Δ t 3 δ 2 ( 2 f 0 + f 2 − 3 f 1 ) f^{n+1}_1=f_1-\frac{2U\Delta t}{\delta}(f_1-f_0)+\frac{4D\Delta t}{3\delta^2}(2f_0+f_2-3f_1) f1n+1=f1−δ2UΔt(f1−f0)+3δ24DΔt(2f0+f2−3f1)
f N n + 1 = f N − U Δ t δ ( f N − f N − 1 ) + 4 D Δ t 3 δ 2 ( 2 f N + 1 + f N − 1 − 3 f N ) f^{n+1}_N=f_N-\frac{U\Delta t}{\delta}(f_N-f_{N-1})+\frac{4D\Delta t}{3\delta^2}(2f_{N+1}+f_{N-1}-3f_N) fNn+1=fN−δUΔt(fN−fN−1)+3δ24DΔt(2fN+1+fN−1−3fN)
使用一阶精度扩散项的隐式算法：
( 1 + 2 U Δ t δ + 4 D Δ t δ 2 ) f 1 n + 1 = 4 D Δ t 3 δ 2 f 2 n + 1 + ( 2 U Δ t δ + 8 D Δ t 3 δ 2 ) f 0 n + 1 + f 1 n \left( 1+\frac{2U\Delta t}{\delta}+\frac{4D\Delta t}{\delta ^2} \right)f_1^{n+1}=\frac{4D\Delta t}{3\delta^2}f_2^{n+1}+ \left( \frac{2U\Delta t}{\delta}+\frac{8D\Delta t}{3\delta^2} \right)f_0^{n+1} + f_1^n (1+δ2UΔt+δ24DΔt)f1n+1=3δ24DΔtf2n+1+(δ2UΔt+3δ28DΔt)f0n+1+f1n
( 1 + U Δ t δ + 4 D Δ t δ 2 ) f N n + 1 = ( U Δ t δ + 4 D Δ t 3 δ 2 ) f N − 1 n + 1 + 8 D Δ t 3 δ 2 f N + 1 n + 1 + f N n \left( 1+\frac{U\Delta t}{\delta}+\frac{4D\Delta t}{\delta^2} \right)f_N^{n+1}=\left( \frac{U\Delta t}{\delta}+\frac{4D\Delta t}{3\delta^2} \right)f_{N-1}^{n+1}+\frac{8D\Delta t}{3\delta^2}f_{N+1}^{n+1} + f_N^n (1+δUΔt+δ24DΔt)fNn+1=(δUΔt+3δ24DΔt)fN−1n+1+3δ28DΔtfN+1n+1+fNn
上式中的 f 0 n + 1 f_0^{n+1} f0n+1以及 f N + 1 n + 1 f_{N+1}^{n+1} fN+1n+1均为已知的边界条件。
3.4.2.3. 扩散项二阶精度的补充结点法
这种方法来自《数值传热学》例4-4
对于内结点法， i = 1 i=1 i=1的结点和左边界距离为 δ / 2 \delta /2 δ/2，我们在左边界上设置一个结点，编号为 i = 0 i=0 i=0，在左边界往左 δ / 2 \delta/2 δ/2处补充结点i = − 1 i=-1 i=−1，对于结点 i = 1 i=1 i=1
∂ 2 f ∂ x 2 ∣ i = 1 = f 2 + f − 1 − 2 f 1 δ 2 + ε ( δ 2 ) \frac{\partial^2 f}{\partial x^2}\bigg|_{i=1}=\frac{f_{2}+f_{-1}-2f_1}{\delta ^2}+\varepsilon(\delta ^2) ∂x2∂2f∣ ∣i=1=δ2f2+f−1−2f1+ε(δ2)
对于结点 i = 0 i=0 i=0
∂ 2 f ∂ x 2 ∣ i = 0 = f 1 + f − 1 − 2 f 0 ( δ / 2 ) 2 + ε ( δ 2 ) \frac{\partial^2 f}{\partial x^2}\bigg|_{i=0}=\frac{f_{1}+f_{-1}-2f_{0}}{(\delta/2) ^2}+\varepsilon(\delta ^2) ∂x2∂2f∣ ∣i=0=(δ/2)2f1+f−1−2f0+ε(δ2)
由于结点0处于边界，在满足边界条件的同时应尽量让该结点满足控制方程，所以对结点0的控制方程的显式离散为
f 0 n + 1 − f 0 n Δ t + U f 0 − f − 1 δ / 2 = D f 1 + f − 1 − 2 f 0 ( δ / 2 ) 2 \frac{f_0^{n+1}-f_0^n}{\Delta t}+U\frac{f_0-f_{-1}}{\delta /2}=D\frac{f_{1}+f_{-1}-2f_{0}}{(\delta/2) ^2} Δtf0n+1−f0n+Uδ/2f0−f−1=D(δ/2)2f1+f−1−2f0
f 0 n + 1 f_0^{n+1} f0n+1以及 f 0 n f_0^n f0n为Dirichlet边界条件，均已知，因此可求出 f − 1 f_{-1} f−1，将 f − 1 f_{-1} f−1代入结点 i = 1 i=1 i=1的扩散项差分方程中，最终可得到满足空间二阶精度的扩散项差分方程。
4. 代数方程组求解方法
显式离散能够直接推导出物理量的更新方程，因此只要给出 Δ x \Delta x Δx和 Δ t \Delta t Δt便能求解，但因为稳定性条件（如Harten定理）限制了所能采用的 Δ t \Delta t Δt
隐式离散的稳定性更好，但是需要联立求解方程组。空间被离散为 N N N个结点，每个结点都对应一个方程，即 N N N个方程 N N N个未知数。方程组求解的快慢直接决定了整个计算流程的效率，因此研究者开发了若干方程组求解方法，可分为两大类，直接求解和迭代求解
4.1. 直接求解
当求解的未知数个数极少时，采用Cramer法则直接求解。本文讨论的一维问题中，每个结点只和周围两个结点相关，即 a P a_P aP只和 a W a_W aW以及 a E a_E aE相关，因此，求解的代数方程组只有三个对角上的元素非零，即三对角矩阵。下面介绍针对该矩阵发展的算法，基于Gauss消元法的Thomas算法：TDMA
4.1.1. TDMA: Tridiagonal matrix method
为了节省篇幅，以下只讨论区域离散方式为外结点法的代数方程组求解。事实上有限体积法的离散方式更像内结点法，因为体心值作为网格平均值有二阶精度（参考链接）。
矩阵形式为 A X = B AX=B AX=B，由于使用了第一类边界条件， x 1 x_1 x1和 x N x_N xN已知，因此 a 12 = a N , N − 1 = 0 a_{12}=a_{N,N-1}=0 a12=aN,N−1=0
[ a 11 a 12 a 21 a 22 a 23 ⋯ a i , i − 1 a i i a i , i + 1 ⋯ a N − 1 , N − 2 a N − 1 , N − 1 a N − 1 , N a N , N − 1 a N , N ] [ x 1 x 2 ⋯ x i ⋯ x N − 1 x N ] = [ b 1 b 2 ⋯ b i ⋯ b N − 1 b N ] \left[ \begin{matrix} a_{11} & a_{12} & & & & \\ a_{21} & a_{22} & a_{23} & & & & \\ & & & \cdots \\ & & a_{i,i-1} & a_{ii} & a_{i,i+1} & & \\ & & & \cdots \\ & & & & a_{N-1,N-2} & a_{N-1,N-1} & a_{N-1,N} \\ & & & & & a_{N,N-1} & a_{N,N} \\ \end{matrix} \right] \left[ \begin{matrix} x_{1} \\ x_{2} \\ \cdots \\ x_{i} \\ \cdots \\ x_{N-1} \\ x_{N} \end{matrix} \right]= \left[ \begin{matrix} b_{1} \\ b_{2} \\ \cdots \\ b_{i} \\ \cdots \\ b_{N-1} \\ b_{N} \end{matrix} \right] ⎣ ⎡a11a21a12a22a23ai,i−1⋯aii⋯ai,i+1aN−1,N−2aN−1,N−1aN,N−1aN−1,NaN,N⎦ ⎤⎣ ⎡x1x2⋯xi⋯xN−1xN⎦ ⎤=⎣ ⎡b1b2⋯bi⋯bN−1bN⎦ ⎤
TDMA算法包含两个步骤，消元和回代。消元指的是通过将上一行整体加到下一行从而消去下一行最左边的元素，从矩阵形状上来看就是将三对角矩阵变成了上对角矩阵，此做法的目的是将最后一行的 a N , N − 1 a_{N,N-1} aN,N−1消去从而能直接计算出 x N x_N xN，得到 x N x_N xN后便逐一向上回代解出其他未知数。
为了讨论方便，把通式写成 A i T i = B i T i + 1 + C i T i − 1 + D i A_iT_i=B_iT_{i+1}+C_iT_{i-1}+D_i AiTi=BiTi+1+CiTi−1+Di，消元的目的是把该式化为 T i − 1 = P i − 1 T i + Q i − 1 T_{i-1}=P_{i-1}T_i+Q_{i-1} Ti−1=Pi−1Ti+Qi−1，通过整理可以得到系数 P i , Q i P_i,Q_i Pi,Qi和 B i , C i , D i B_i,C_i,D_i Bi,Ci,Di之间的关系，
P i = B i A i − C i P i − 1 , Q i = D i + C i Q i − 1 A i − C i P i − 1 P_i=\frac{B_i}{A_i-C_iP_{i-1}},\quad Q_i=\frac{D_i+C_iQ_{i-1}}{A_i-C_iP_{i-1}} Pi=Ai−CiPi−1Bi,Qi=Ai−CiPi−1Di+CiQi−1
可以看出， P i P_i Pi和 Q i Q_i Qi都需要递归求解
P 1 = B 1 A 1 , Q 1 = D 1 A 1 , Q N = T N P_1=\frac{B_1}{A_1},\quad Q_1=\frac{D_1}{A_1},\quad Q_N=T_N P1=A1B1,Q1=A1D1,QN=TN
综上，TDMA计算流程为
计算 P 1 , Q 1 P_1,Q_1 P1,Q1更新 P i , Q i P_i,Q_i Pi,Qi，得到 P N , Q N P_N,Q_N PN,QN，其中 Q N = T N Q_N=T_N QN=TN由 T i − 1 T_{i-1} Ti−1的表达式计算 T 1 − T N − 1 T_1-T_{N-1} T1−TN−1
4.2. 迭代求解
待补充
5. 计算结果讨论
5.1. 显式格式
由上述离散过程可见，对流项采用一阶迎风。一阶格式的结果是截断项首项为二阶导数，而二阶导数项有扩散性质，即物理量会向各个方向传播。通常这种由二阶导数引起的扩散现象称为假扩散/人工粘性/数值粘性。具体讨论可参考《数值传热学》第二版5.5节。
下图为显式离散在 t = 2.5 s t=2.5s t=2.5s的计算结果，区域离散方式为外结点和内结点的差别不大。结点间距 Δ x = 0.009 \Delta x=0.009 Δx=0.009，时间步 Δ t = 0.0001 s \Delta t=0.0001s Δt=0.0001s，该步长满足Harten定理。可见，模拟结果与解析解相比，有矮胖特征，意味着存在假扩散现象。
5.2 隐式格式
下图为隐式离散在 t = 2.5 s t=2.5s t=2.5s的计算结果，区域离散采用外结点，结点间距0.009，时间步 Δ t = 0.001 s \Delta t=0.001s Δt=0.001s，而在相同时间步下，显式格式的计算结果会发散，比较来看，隐式格式能够适用更大的时间步长。
6. 代码
代码链接
