有限体积法（1）——一维扩散方程的推导

对于变量ϕ\phiϕ的输运方程为：

∂(ρϕ)∂t+∇⋅(ρϕu)=∇⋅(Γ∇ϕ)+Sϕ(1)\frac{\partial ( \rho \phi)}{\partial t} + \nabla \cdot (\rho \phi \bold u) = \nabla \cdot (\Gamma \nabla \phi) + S_\phi \tag{1} ∂t∂(ρϕ)+∇⋅(ρϕu)=∇⋅(Γ∇ϕ)+Sϕ(1)

其中，Γ\GammaΓ为扩散系数。
方程(1)(1)(1)从左到右的各项分别是时间项、对流项、扩散项和源项。将方程(1)(1)(1)中略去时间项和对流项就是稳态扩散方程：
∇⋅(Γ∇ϕ)+Sϕ=0(2)\nabla \cdot (\Gamma \nabla \phi) + S_\phi = 0 \tag{2} ∇⋅(Γ∇ϕ)+Sϕ=0(2)

对方程(2)(2)(2)在有限控制体内积分，并根据高斯散度定理有，
∫CV∇⋅(Γ∇ϕ)dV+∫CVSϕdV=∫A~n⋅(Γ∇ϕ)dA+∫CVSϕdV=0(3)\begin{aligned} &\int_{CV} \nabla \cdot (\Gamma \nabla \phi) dV + \int_{CV} S_\phi dV \\ \\ &= \int_{\tilde A} \bold n \cdot (\Gamma \nabla \phi) dA + \int_{CV} S_\phi dV = 0 \end{aligned} \tag{3}∫CV∇⋅(Γ∇ϕ)dV+∫CVSϕdV=∫A~n⋅(Γ∇ϕ)dA+∫CVSϕdV=0(3)

其中，n\bold nn为边界面A~\tilde AA~的法向量。
考虑一维模型，扩散方程为
ddx(Γdϕdx)+S=0(4)\frac{d}{dx}\left( \Gamma \frac{d\phi}{dx} \right) + S = 0 \tag{4} dxd(Γdxdϕ)+S=0(4)

要离散方程(4)(4)(4)首先要生成网格，如下图，将一维计算域分成5段（即5个网格单元），每个网格单元就是一个有限控制体。两端边界为A和B，每个网格单元有一个节点，如W、P和E，节点一般为网格单元的中心点，也是变量保存的位置。两个相邻节点的中点是网格单元的边界。

网格与边界之间的距离关系如下图，

扩散方程的离散主要包括两个部分：散度的离散和梯度的离散。

散度的离散

根据方程式(3)(3)(3)，一维扩散方程(4)(4)(4)在控制体单元ΔV\Delta VΔV（边界为A~\tilde AA~）上的积分有，
∫ΔVddx(Γdϕdx)dV+∫ΔVSdV=∫A~n⋅(Γdϕdxi)dA+∫ΔVSdV=(ΓAdϕdx)e−(ΓAdϕdx)w+SˉΔV=0(5)\begin{aligned} &\int_{\Delta V} \frac{d}{dx} \left( \Gamma \frac{d \phi}{dx} \right)dV + \int_{\Delta V} SdV \\ \\ &= \int_{\tilde A} \bold n \cdot \left( \Gamma \frac{d \phi}{dx} \bold i \right) dA + \int_{\Delta V}SdV \\ \\ &= \left( \Gamma A \frac{d \phi}{dx} \right)_e - \left( \Gamma A \frac{d \phi}{dx} \right)_w + \bar{S} \Delta V = 0 \end{aligned} \tag{5}∫ΔVdxd(Γdxdϕ)dV+∫ΔVSdV=∫A~n⋅(Γdxdϕi)dA+∫ΔVSdV=(ΓAdxdϕ)e−(ΓAdxdϕ)w+SˉΔV=0(5)

其中，AAA是边界面的面积， Sˉ\bar{S}Sˉ为控制体单元内的平均值。从上图中可见，在边界www处，边界的法向量n\bold nn是指向x轴负向的，所以n⋅i=−1\bold n \cdot \bold i = -1n⋅i=−1。

梯度的离散

方程式(5)(5)(5)仅是半离散方程，需要进一步对梯度项dϕdx\frac{d \phi}{dx}dxdϕ进行离散。从式中看见，梯度项都是位于单元边界面位置的，而边界处是不保存变量值，边界面两边的节点会保存变量值，因此可以用有限差分法求解边界面的导数，比如用二阶中心差分格式，
(Γdϕdx)w=ΓwϕP−ϕWδxWP(6a)\left( \Gamma \frac{d \phi}{dx} \right)_w = \Gamma_w \frac{\phi_P - \phi_W}{\delta x_{WP}} \tag{6a}(Γdxdϕ)w=ΓwδxWPϕP−ϕW(6a)

(Γdϕdx)e=ΓeϕE−ϕPδxPE(6b)\left( \Gamma \frac{d \phi}{dx} \right)_e = \Gamma_e \frac{\phi_E - \phi_P}{\delta x_{PE}} \tag{6b} (Γdxdϕ)e=ΓeδxPEϕE−ϕP(6b)

源项SˉΔV\bar S \Delta VSˉΔV通常是变量ϕP\phi_PϕP的函数，这里假设为线性关系，即
SˉΔV=Su+SPϕP(7)\bar S \Delta V = S_u + S_P \phi_P \tag{7} SˉΔV=Su+SPϕP(7)

将式(6a)(6a)(6a)、(6b)(6b)(6b)和(7)(7)(7)带入到式(5)(5)(5)中，
ΓeAeϕE−ϕPδxPE−ΓwAwϕP−ϕWδxWP+(Su+SPϕP)=0(8)\Gamma_e A_e \frac{\phi_E - \phi_P}{\delta x_{PE}} - \Gamma_w A_w \frac{\phi_P - \phi_W}{\delta x_{WP}} + (S_u + S_P \phi_P) = 0 \tag{8}ΓeAeδxPEϕE−ϕP−ΓwAwδxWPϕP−ϕW+(Su+SPϕP)=0(8)
重新整理方程式(8)(8)(8)得，
(ΓeδxPE+ΓwδxWP−SP)ϕP=(ΓwδxWP)ϕW+(ΓeδxPE)ϕE+Su(9)\left( \frac{\Gamma_e}{\delta x_{PE}} + \frac{\Gamma_w}{\delta x_{WP}} - S_P\right) \phi_P = \left( \frac{\Gamma_w}{\delta x_{WP}} \right) \phi_W + \left( \frac{\Gamma_e}{\delta x_{PE}} \right) \phi_E + S_u \tag{9}(δxPEΓe+δxWPΓw−SP)ϕP=(δxWPΓw)ϕW+(δxPEΓe)ϕE+Su(9)
简化一下，
aPϕP=aWϕW+aEϕE+Su(10)a_P \phi_P = a_W \phi_W + a_E \phi_E + S_u \tag{10}aPϕP=aWϕW+aEϕE+Su(10)
其中各系数为
{aW=ΓwδxWPaE=ΓeδxPEaP=aW+aE−SP(11)\begin{cases} a_W &= \frac{\Gamma_w}{\delta x_{WP}} \\ \\ a_E &= \frac{\Gamma_e}{\delta x_{PE}} \\ \\ a_P &= a_W + a_E - S_P \end{cases}\tag{11}⎩⎨⎧aWaEaP=δxWPΓw=δxPEΓe=aW+aE−SP(11)

方程式(5)(5)(5)中的扩散系数Γ\GammaΓ有两种方法来计算：算术平均法（相当于线性插值）和调和平均法（详见文献[2]）。这里先采用简单的线性插值来表示，即
Γw=ΓW+ΓP2(12a)\Gamma_w = \frac{\Gamma_W + \Gamma_P}{2} \tag{12a} Γw=2ΓW+ΓP(12a)

Γe=ΓP+ΓE2(12b)\Gamma_e = \frac{\Gamma_P + \Gamma_E}{2} \tag{12b} Γe=2ΓP+ΓE(12b)

参考资料：

Versteeg H K , Malalasekera W . An introduction to computational fluid dynamics : the finite volume method = 计算流体动力学导论[M]. 世界图书出版公司, 2010.
陶文铨. 数值传热学-第2版[M]. 西安交通大学出版社, 2001.