有限体积法（2）——二维、三维扩散方程的离散推导

稳态扩散方程：
∇⋅(Γ∇ϕ)+Sϕ=0(1)\nabla \cdot ( \Gamma \nabla \phi) + S_\phi =0 \tag{1} ∇⋅(Γ∇ϕ)+Sϕ=0(1)
在有限控制体内积分，并由高斯散度定理有，
∫CV∇⋅(Γ∇ϕ)dV+∫CVSϕdV=∫A~n⋅(Γ∇ϕ)dA+∫CVSϕdV=0(2)\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 \tag{2} \end{aligned}∫CV∇⋅(Γ∇ϕ)dV+∫CVSϕdV=∫A~n⋅(Γ∇ϕ)dA+∫CVSϕdV=0(2)

其中n\bold nn为边界面A~\tilde AA~的法向量，Γ\GammaΓ为扩散系数。

二维扩散方程的离散

根据式(1)(1)(1)，二维模型的扩散方程为，
∂∂x(Γ∂ϕ∂x)+∂∂y(Γ∂ϕ∂y)+Sϕ=0(3)\frac{\partial}{\partial x} \left( \Gamma \frac{\partial \phi}{\partial x} \right) + \frac{\partial }{\partial y} \left( \Gamma \frac{\partial \phi}{\partial y} \right) + S_\phi =0 \tag{3} ∂x∂(Γ∂x∂ϕ)+∂y∂(Γ∂y∂ϕ)+Sϕ=0(3)
与一维扩散方程推导类似，先将二维计算域划分网格，

在二维空间中，梯度矢量有两个分量，∂ϕ∂xi\frac{\partial \phi}{\partial x} \bold i∂x∂ϕi和∂ϕ∂yj\frac{\partial \phi}{\partial y} \bold j∂y∂ϕj，方向分别指向xxx和yyy轴的正方向。

散度的离散

在二维空间，单元P的边界包括w、e、sw、e、sw、e、s和nnn四个边界面，由于网格和坐标轴是平行或垂直的，∂ϕ∂y\frac{\partial \phi}{\partial y}∂y∂ϕ在边界面w、ew、ew、e上的通量为零，同理∂ϕ∂x\frac{\partial \phi}{\partial x}∂x∂ϕ在边界面s、ns、ns、n上的通量为零。所以，

∫CV∇⋅(Γ∇ϕ)dV+∫CVSϕdV=∫A~n⋅(Γ∇ϕ)dA+∫CVSϕdV=∫A~n⋅[(Γ∂ϕ∂x)i+(Γ∂ϕ∂y)j]dA+∫CVSϕdV=[(ΓA∂ϕ∂x)e−(ΓA∂ϕ∂x)w]+[(ΓA∂ϕ∂y)n−(ΓA∂ϕ∂y)s]+SˉϕΔV=0(4)\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 \\ \\ &=\int_{\tilde A} \bold n \cdot \left[ \left( \Gamma \frac{\partial \phi}{\partial x} \right) \mathbf i + \left( \Gamma \frac{\partial \phi}{\partial y} \right) \mathbf j \right] dA + \int_{CV}S_\phi dV \\ \\ &=\left[ \left( \Gamma A \frac{\partial \phi}{\partial x} \right)_e - \left(\Gamma A \frac{\partial \phi}{\partial x} \right)_w \right] + \left[ \left( \Gamma A \frac{\partial \phi}{\partial y} \right)_n - \left(\Gamma A \frac{\partial \phi}{\partial y} \right)_s \right] + \bar S_\phi \Delta V\\ \\ &=0 \tag{4} \end{aligned} ∫CV∇⋅(Γ∇ϕ)dV+∫CVSϕdV=∫A~n⋅(Γ∇ϕ)dA+∫CVSϕdV=∫A~n⋅[(Γ∂x∂ϕ)i+(Γ∂y∂ϕ)j]dA+∫CVSϕdV=[(ΓA∂x∂ϕ)e−(ΓA∂x∂ϕ)w]+[(ΓA∂y∂ϕ)n−(ΓA∂y∂ϕ)s]+SˉϕΔV=0(4)

其中AAA代表边界面的面积，Sˉ\bar{S}Sˉ是控制体单元ΔV\Delta VΔV内的平均源项。边界面www和sss上的通量为负，原因与一维扩散方程情况类似。

梯度的离散

与一维方程相同，梯度项采用中心差分格式离散，
(ΓA∂ϕ∂x)w=ΓwAw(ϕP−ϕW)δxWP(5a)\left ( \Gamma A \frac{\partial \phi}{\partial x} \right)_w = \Gamma_w A_w \frac{ ( \phi_P - \phi_W) }{\delta x_{WP}} \tag{5a} (ΓA∂x∂ϕ)w=ΓwAwδxWP(ϕP−ϕW)(5a)

(ΓA∂ϕ∂x)e=ΓeAe(ϕE−ϕP)δxPE(5b)\left ( \Gamma A \frac{\partial \phi}{\partial x} \right)_e = \Gamma_e A_e \frac{ ( \phi_E - \phi_P) }{\delta x_{PE}} \tag{5b} (ΓA∂x∂ϕ)e=ΓeAeδxPE(ϕE−ϕP)(5b)

(ΓA∂ϕ∂y)s=ΓsAs(ϕP−ϕS)δySP(5c)\left ( \Gamma A \frac{\partial \phi}{\partial y} \right)_s = \Gamma_s A_s \frac{ ( \phi_P - \phi_S) }{\delta y_{SP}} \tag{5c} (ΓA∂y∂ϕ)s=ΓsAsδySP(ϕP−ϕS)(5c)

(ΓA∂ϕ∂y)n=ΓnAn(ϕN−ϕP)δyPN(5d)\left ( \Gamma A \frac{\partial \phi}{\partial y} \right)_n = \Gamma_n A_n \frac{ ( \phi_N - \phi_P) }{\delta y_{PN}} \tag{5d} (ΓA∂y∂ϕ)n=ΓnAnδyPN(ϕN−ϕP)(5d)

将式(5a)(5a)(5a)~(5d)(5d)(5d)带入到式(4)(4)(4)中，有
ΓeAe(ϕE−ϕP)δxPE−ΓwAw(ϕP−ϕW)δxWP+ΓnAn(ϕN−ϕP)δyPN−ΓsAs(ϕP−ϕS)δySP+(Su+SPϕP)=0(6)\begin{aligned} \Gamma_e A_e \frac{(\phi_E - \phi_P)}{\delta x_{PE}} - \Gamma_w A_w \frac{(\phi_P - \phi_W)}{\delta x_{WP}} + \Gamma_n A_n \frac{(\phi_N - \phi_P)}{\delta y_{PN}} \\ \\ - \Gamma_s A_s \frac{(\phi_P - \phi_S)}{\delta y_{SP}} + (S_u + S_P \phi_P) =0 \tag{6} \end{aligned}ΓeAeδxPE(ϕE−ϕP)−ΓwAwδxWP(ϕP−ϕW)+ΓnAnδyPN(ϕN−ϕP)−ΓsAsδySP(ϕP−ϕS)+(Su+SPϕP)=0(6)

式中(Su+SPϕP)=SˉϕΔV(S_u + S_P \phi_P)=\bar S_\phi \Delta V(Su+SPϕP)=SˉϕΔV，整理得，
(ΓwAwδxWP+ΓeAeδxPE+ΓsAsδySP+ΓnAnδyPN−SP)ϕP=(ΓwAwδxWP)ϕW+(ΓeAeδxPE)ϕE+(ΓsAsδySP)ϕS+(ΓnAnδyP)ϕN+Su(7)\begin{aligned} &\left( \frac{\Gamma_w A_w}{\delta x_{WP}} + \frac{\Gamma_e A_e}{\delta x_{PE}} + \frac{\Gamma_s A_s}{\delta y_{SP}} +\frac{\Gamma_n A_n}{\delta y_{PN}} - S_P \right) \phi_P \\ \\ &=\left( \frac{\Gamma_w A_w}{\delta x_{WP}} \right) \phi_W + \left( \frac{\Gamma_e A_e}{\delta x_{PE}} \right) \phi_E + \left( \frac{\Gamma_s A_s}{\delta y_{SP}} \right) \phi_S + \left( \frac{\Gamma_n A_n}{\delta y_{P}} \right) \phi_N + S_u \end{aligned} \tag{7}(δxWPΓwAw+δxPEΓeAe+δySPΓsAs+δyPNΓnAn−SP)ϕP=(δxWPΓwAw)ϕW+(δxPEΓeAe)ϕE+(δySPΓsAs)ϕS+(δyPΓnAn)ϕN+Su(7)

简化之，

aPϕP=aWϕW+aEϕE+aSϕS+aNϕN+Su(8)a_P \phi_P = a_W \phi_W + a_E \phi_E + a_S \phi_S + a_N \phi_N + S_u \tag{8} aPϕP=aWϕW+aEϕE+aSϕS+aNϕN+Su(8)

其中各系数为
aW=ΓwAwδxWP(9a)a_W = \frac{\Gamma_w A_w}{\delta x_{WP}} \tag{9a} aW=δxWPΓwAw(9a)

aE=ΓeAeδxPE(9b)a_E = \frac{\Gamma_e A_e}{\delta x_{PE}} \tag{9b} aE=δxPEΓeAe(9b)

aS=ΓsAsδySP(9c)a_S = \frac{\Gamma_s A_s}{\delta y_{SP}} \tag{9c} aS=δySPΓsAs(9c)

aN=ΓnAnδyPN(9d)a_N = \frac{\Gamma_n A_n}{\delta y_{PN}} \tag{9d} aN=δyPNΓnAn(9d)

aP=aW+aE+aS+aN−SP(9e)a_P = a_W + a_E + a_S + a_N - S_P \tag{9e} aP=aW+aE+aS+aN−SP(9e)
边界处的扩散系数Γ\GammaΓ可以通过线性插值计算，见一维扩散方程推导中公式(12)(12)(12)。

三维离散的扩散方程

三维扩散方程：
∂∂x(Γ∂ϕ∂x)+∂∂y(Γ∂ϕ∂y)+∂∂z(Γ∂ϕ∂z)+Sϕ=0(10)\frac{\partial}{\partial x} \left( \Gamma \frac{\partial \phi}{\partial x} \right) + \frac{\partial}{\partial y} \left( \Gamma \frac{\partial \phi}{\partial y} \right) + \frac{\partial}{\partial z} \left( \Gamma \frac{\partial \phi}{\partial z} \right) + S_\phi = 0 \tag{10}∂x∂(Γ∂x∂ϕ)+∂y∂(Γ∂y∂ϕ)+∂z∂(Γ∂z∂ϕ)+Sϕ=0(10)
三维扩散方程的离散推导和二维的套路是一毛一样的，无非就是多了两个边界面bbb和ttt，散度离散和梯度离散时多了两项，梯度离散方法和方程式(5)(5)(5)一样也用中心差分格式，然后带入，整理，简化之。

散度离散

∫CV∇⋅(Γ∇ϕ)dV+∫CVSϕdV=∫A~n⋅[∂∂x(Γ∂ϕ∂x)i+∂∂y(Γ∂ϕ∂y)j+∂∂z(Γ∂ϕ∂z)k]dA+∫CVSϕdV=[(ΓA∂ϕ∂x)e−(ΓA∂ϕ∂x)w]+[(ΓA∂ϕ∂y)n−(ΓA∂ϕ∂y)s]+[(ΓA∂ϕ∂z)t−(ΓA∂ϕ∂z)b]+SˉϕΔV=0(11)\begin{aligned} &\int_{CV} \nabla \cdot ( \Gamma \nabla \phi ) dV + \int_{CV} S_\phi dV \\ \\ &=\int_{\tilde A} \bold n \cdot \left[ \frac{\partial}{\partial x}\left( \Gamma \frac{\partial \phi}{\partial x} \right) \bold i + \frac{\partial}{\partial y}\left( \Gamma \frac{\partial \phi}{\partial y} \right) \bold j + \frac{\partial}{\partial z}\left( \Gamma \frac{\partial \phi}{\partial z} \right) \bold k \right] dA + \int_{CV} S_\phi dV \\ \\ &=\left[ \left( \Gamma A \frac{\partial \phi}{\partial x} \right)_e - \left(\Gamma A \frac{\partial \phi}{\partial x} \right)_w \right] + \left[ \left( \Gamma A \frac{\partial \phi}{\partial y} \right)_n - \left(\Gamma A \frac{\partial \phi}{\partial y} \right)_s \right] \\ \\ & \qquad + \left[ \left( \Gamma A \frac{\partial \phi}{\partial z} \right)_t - \left(\Gamma A \frac{\partial \phi}{\partial z} \right)_b \right] + \bar S_\phi \Delta V =0 \tag{11} \end{aligned}∫CV∇⋅(Γ∇ϕ)dV+∫CVSϕdV=∫A~n⋅[∂x∂(Γ∂x∂ϕ)i+∂y∂(Γ∂y∂ϕ)j+∂z∂(Γ∂z∂ϕ)k]dA+∫CVSϕdV=[(ΓA∂x∂ϕ)e−(ΓA∂x∂ϕ)w]+[(ΓA∂y∂ϕ)n−(ΓA∂y∂ϕ)s]+[(ΓA∂z∂ϕ)t−(ΓA∂z∂ϕ)b]+SˉϕΔV=0(11)

梯度离散

梯度项采用中心差分格式离散，边界面e、w、s、ne、w、s、ne、w、s、n处的梯度离散和式(5)(5)(5)一样，多出来的两项公式如下，
(ΓA∂ϕ∂z)b=ΓbAb(ϕP−ϕB)δzBP(12a)\left ( \Gamma A \frac{\partial \phi}{\partial z} \right)_b = \Gamma_b A_b \frac{ ( \phi_P - \phi_B) }{\delta z_{BP}} \tag{12a} (ΓA∂z∂ϕ)b=ΓbAbδzBP(ϕP−ϕB)(12a)

(ΓA∂ϕ∂x)t=ΓtAt(ϕT−ϕP)δzPT(12b)\left ( \Gamma A \frac{\partial \phi}{\partial x} \right)_t = \Gamma_t A_t \frac{ ( \phi_T - \phi_P) }{\delta z_{PT}} \tag{12b} (ΓA∂x∂ϕ)t=ΓtAtδzPT(ϕT−ϕP)(12b)
带入式(11)(11)(11)，有
[ΓeAe(ϕE−ϕP)δxPE−ΓwAw(ϕP−ϕW)δxWP]+[ΓnAn(ϕN−ϕP)δyPN−ΓsAs(ϕP−ϕS)δySP]+[ΓtAt(ϕT−ϕP)δzPT−ΓbAb(ϕP−ϕB)δzBP]+(Su+SPϕP)=0(13)\begin{aligned} &\left[ \Gamma_e A_e \frac{(\phi_E - \phi_P)}{\delta x_{PE}} - \Gamma_w A_w \frac{(\phi_P - \phi_W)}{\delta x_{WP}} \right] \\ \\ &+ \left[ \Gamma_n A_n \frac{(\phi_N - \phi_P)}{\delta y_{PN}} - \Gamma_s A_s \frac{(\phi_P - \phi_S)}{\delta y_{SP}} \right ] \\ \\ &+ \left[ \Gamma_t A_t \frac{(\phi_T - \phi_P)}{\delta z_{PT}} - \Gamma_b A_b \frac{(\phi_P - \phi_B)}{\delta z_{BP}} \right ] \\ \\ &+ (S_u + S_P \phi_P) =0 \tag{13} \end{aligned}[ΓeAeδxPE(ϕE−ϕP)−ΓwAwδxWP(ϕP−ϕW)]+[ΓnAnδyPN(ϕN−ϕP)−ΓsAsδySP(ϕP−ϕS)]+[ΓtAtδzPT(ϕT−ϕP)−ΓbAbδzBP(ϕP−ϕB)]+(Su+SPϕP)=0(13)
整理并简化之，
aPϕP=aWϕW+aEϕE+aSϕS+aNϕN+aBϕB+aTϕT+Su(14)a_P \phi_P = a_W \phi_W + a_E \phi_E + a_S \phi_S + a_N \phi_N + a_B \phi_B + a_T \phi_T+ S_u \tag{14} aPϕP=aWϕW+aEϕE+aSϕS+aNϕN+aBϕB+aTϕT+Su(14)
各项系数，
aW=ΓwAwδxWP(15a)a_W = \frac{\Gamma_w A_w}{\delta x_{WP}} \tag{15a} aW=δxWPΓwAw(15a)

aE=ΓeAeδxPE(15b)a_E = \frac{\Gamma_e A_e}{\delta x_{PE}} \tag{15b} aE=δxPEΓeAe(15b)

aS=ΓsAsδySP(15c)a_S = \frac{\Gamma_s A_s}{\delta y_{SP}} \tag{15c} aS=δySPΓsAs(15c)

aN=ΓnAnδyPN(15d)a_N = \frac{\Gamma_n A_n}{\delta y_{PN}} \tag{15d} aN=δyPNΓnAn(15d)

aB=ΓbAbδzBP(15e)a_B = \frac{\Gamma_b A_b}{\delta z_{BP}} \tag{15e} aB=δzBPΓbAb(15e)

aT=ΓtAtδzPT(15f)a_T = \frac{\Gamma_t A_t}{\delta z_{PT}} \tag{15f} aT=δzPTΓtAt(15f)

aP=aW+aE+aS+aN+aB+aT−SP(15h)a_P = a_W + a_E + a_S + a_N + a_B + a_T - S_P \tag{15h} aP=aW+aE+aS+aN+aB+aT−SP(15h)

总结

通过有限体积法离散后的稳态扩散方程可以统一写成如下形式，
aPϕP=Σanbϕnb+Su(16)a_P \phi_P = \Sigma a_{nb}\phi_{nb} + S_u \tag{16} aPϕP=Σanbϕnb+Su(16)
其中，“nbnbnb”代表单元PPP相邻的节点，Σ\SigmaΣ代表所有相邻节点之和。
系数aPa_PaP与anba_{nb}anb之间的关系为：
aP=Σanb−SP(17)a_P = \Sigma a_{nb} - S_P \tag{17} aP=Σanb−SP(17)
相邻节点的系数anba_{nb}anb计算如下表，

	aWa_WaW	aEa_EaE	aSa_SaS	aNa_NaN	aBa_BaB	aTa_TaT
1D	ΓwAwδxWP\frac{\Gamma_w A_w}{\delta x_{WP}}δxWPΓwAw	ΓwAeδxPE\frac{\Gamma_w A_e}{\delta x_{PE}}δxPEΓwAe
2D	ΓwAwδxWP\frac{\Gamma_w A_w}{\delta x_{WP}}δxWPΓwAw	ΓwAeδxPE\frac{\Gamma_w A_e}{\delta x_{PE}}δxPEΓwAe	ΓsAsδySP\frac{\Gamma_s A_s}{\delta y_{SP}}δySPΓsAs	ΓnAnδyPN\frac{\Gamma_n A_n}{\delta y_{PN}}δyPNΓnAn
3D	ΓwAwδxWP\frac{\Gamma_w A_w}{\delta x_{WP}}δxWPΓwAw	ΓwAeδxPE\frac{\Gamma_w A_e}{\delta x_{PE}}δxPEΓwAe	ΓsAsδySP\frac{\Gamma_s A_s}{\delta y_{SP}}δySPΓsAs	ΓnAnδyPN\frac{\Gamma_n A_n}{\delta y_{PN}}δyPNΓnAn	ΓbAbδzBP\frac{\Gamma_b A_b}{\delta z_{BP}}δzBPΓbAb	ΓtAtδzPT\frac{\Gamma_t A_t}{\delta z_{PT}}δzPTΓtAt

参考资料：

Versteeg H K , Malalasekera W . An introduction to computational fluid dynamics : the finite volume method = 计算流体动力学导论[M]. 世界图书出版公司, 2010.