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

稳态扩散方程：
? ? ( Γ ? ? ) + S ? = 0 (1) \nabla \cdot ( \Gamma \nabla \phi) + S_\phi =0 \tag{1} ??(Γ??)+S??=0(1)
在有限控制体内积分，并由高斯散度定理有，
∫ C V ? ? ( Γ ? ? ) d V + ∫ C V S ? d V = ∫ A ~ n ? ( Γ ? ? ) d A + ∫ C V S ? d V = 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+∫CV?S??dV=∫A~?n?(Γ??)dA+∫CV?S??dV=0?(2)
其中 n \bold n n为边界面 A ~ \tilde A A~的法向量， Γ \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)
与一维扩散方程推导类似，先将二维计算域划分网格，

文章图片

在二维空间中，梯度矢量有两个分量， ? ? ? x i \frac{\partial \phi}{\partial x} \bold i ?x???i和 ? ? ? y j \frac{\partial \phi}{\partial y} \bold j ?y???j，方向分别指向 x x x和 y y y轴的正方向。
散度的离散
在二维空间，单元P的边界包括 w 、 e 、 s w、e、s w、e、s和 n n n四个边界面，由于网格和坐标轴是平行或垂直的， ? ? ? y \frac{\partial \phi}{\partial y} ?y???在边界面 w 、 e w、e w、e上的通量为零，同理 ? ? ? x \frac{\partial \phi}{\partial x} ?x???在边界面 s 、 n s、n s、n上的通量为零。所以，
∫ C V ? ? ( Γ ? ? ) d V + ∫ C V S ? d V = ∫ A ~ n ? ( Γ ? ? ) d A + ∫ C V S ? d V = ∫ A ~ n ? [ ? ? x ( Γ ? ? ? x ) i + ? ? y ( Γ ? ? ? y ) j ] d A + ∫ C V S ? d V = [ ( Γ 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[ \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 \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+∫CV?S??dV=∫A~?n?(Γ??)dA+∫CV?S??dV=∫A~?n?[?x??(Γ?x???)i+?y??(Γ?y???)j]dA+∫CV?S??dV=[(ΓA?x???)e??(ΓA?x???)w?]+[(ΓA?y???)n??(ΓA?y???)s?]+Sˉ??ΔV=0?(4)
其中 A A A代表边界面的面积， S ˉ \bar{S} Sˉ是控制体单元 Δ V \Delta V ΔV内的平均源项。边界面 w w w和 s s s上的通量为负，原因与一维扩散方程情况类似。
梯度的离散
与一维方程相同，梯度项采用中心差分格式离散，
( Γ A ? ? ? x ) w = Γ w A w ( ? P ? ? W ) δ x W P (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?=Γw?Aw?δxWP?(?P???W?)?(5a)
( Γ A ? ? ? x ) e = Γ e A e ( ? E ? ? P ) δ x P E (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?=Γe?Ae?δxPE?(?E???P?)?(5b)
( Γ A ? ? ? y ) s = Γ s A s ( ? P ? ? S ) δ y S P (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?=Γs?As?δySP?(?P???S?)?(5c)
( Γ A ? ? ? y ) n = Γ n A n ( ? N ? ? P ) δ y P N (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?=Γn?An?δyPN?(?N???P?)?(5d)
将式 ( 5 a ) (5a) (5a)~ ( 5 d ) (5d) (5d)带入到式 ( 4 ) (4) (4)中，有
Γ e A e ( ? E ? ? P ) δ x P E ? Γ w A w ( ? P ? ? W ) δ x W P + Γ n A n ( ? N ? ? P ) δ y P N ? Γ s A s ( ? P ? ? S ) δ y S P + ( S u + S P ? 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} Γe?Ae?δxPE?(?E???P?)??Γw?Aw?δxWP?(?P???W?)?+Γn?An?δyPN?(?N???P?)??Γs?As?δySP?(?P???S?)?+(Su?+SP??P?)=0?(6)
式中 ( S u + S P ? P ) = S ˉ ? Δ V (S_u + S_P \phi_P)=\bar S_\phi \Delta V (Su?+SP??P?)=Sˉ??ΔV，整理得，
( Γ w A w δ x W P + Γ e A e δ x P E + Γ s A s δ y S P + Γ n A n δ y P N ? S P ) ? P = ( Γ w A w δ x W P ) ? W + ( Γ e A e δ x P E ) ? E + ( Γ s A s δ y S P ) ? S + ( Γ n A n δ y P ) ? N + S u (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?Γw?Aw??+δxPE?Γe?Ae??+δySP?Γs?As??+δyPN?Γn?An???SP?)?P?=(δxWP?Γw?Aw??)?W?+(δxPE?Γe?Ae??)?E?+(δySP?Γs?As??)?S?+(δyP?Γn?An??)?N?+Su??(7)
简化之，
a P ? P = a W ? W + a E ? E + a S ? S + a N ? N + S u (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)
其中各系数为
a W = Γ w A w δ x W P (9a) a_W = \frac{\Gamma_w A_w}{\delta x_{WP}} \tag{9a} aW?=δxWP?Γw?Aw??(9a)
a E = Γ e A e δ x P E (9b) a_E = \frac{\Gamma_e A_e}{\delta x_{PE}} \tag{9b} aE?=δxPE?Γe?Ae??(9b)
a S = Γ s A s δ y S P (9c) a_S = \frac{\Gamma_s A_s}{\delta y_{SP}} \tag{9c} aS?=δySP?Γs?As??(9c)
a N = Γ n A n δ y P N (9d) a_N = \frac{\Gamma_n A_n}{\delta y_{PN}} \tag{9d} aN?=δyPN?Γn?An??(9d)
a P = a W + a E + a S + a N = S P (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)
三维扩散方程的离散推导和二维的套路是一毛一样的，无非就是多了两个边界面 b b b和 t t t，散度离散和梯度离散时多了两项，梯度离散方法和方程式 ( 5 ) (5) (5)一样也用中心差分格式，然后带入，整理，简化之。
散度离散
∫ C V ? ? ( Γ ? ? ) d V + ∫ C V S ? d V = ∫ A ~ n ? [ ? ? x ( Γ ? ? ? x ) i + ? ? y ( Γ ? ? ? y ) j + ? ? z ( Γ ? ? ? z ) k ] d A + ∫ C V S ? d V = [ ( Γ 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+∫CV?S??dV=∫A~?n?[?x??(Γ?x???)i+?y??(Γ?y???)j+?z??(Γ?z???)k]dA+∫CV?S??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 、 n e、w、s、n e、w、s、n处的梯度离散和式 ( 5 ) (5) (5)一样，多出来的两项公式如下，
( Γ A ? ? ? z ) b = Γ b A b ( ? P ? ? B ) δ z B P (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?=Γb?Ab?δzBP?(?P???B?)?(12a)
( Γ A ? ? ? x ) t = Γ t A t ( ? T ? ? P ) δ z P T (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?=Γt?At?δzPT?(?T???P?)?(12b)
带入式 ( 11 ) (11) (11)，有
[ Γ e A e ( ? E ? ? P ) δ x P E ? Γ w A w ( ? P ? ? W ) δ x W P ] + [ Γ n A n ( ? N ? ? P ) δ y P N ? Γ s A s ( ? P ? ? S ) δ y S P ] + [ Γ t A t ( ? T ? ? P ) δ z P T ? Γ b A b ( ? P ? ? B ) δ z B P ] + ( S u + S P ? 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} ?[Γe?Ae?δxPE?(?E???P?)??Γw?Aw?δxWP?(?P???W?)?]+[Γn?An?δyPN?(?N???P?)??Γs?As?δySP?(?P???S?)?]+[Γt?At?δzPT?(?T???P?)??Γb?Ab?δzBP?(?P???B?)?]+(Su?+SP??P?)=0?(13)
整理并简化之，
a P ? P = a W ? W + a E ? E + a S ? S + a N ? N + a B ? B + a T ? T + S u (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)
各项系数，
a W = Γ w A w δ x W P (15a) a_W = \frac{\Gamma_w A_w}{\delta x_{WP}} \tag{15a} aW?=δxWP?Γw?Aw??(15a)
a E = Γ e A e δ x P E (15b) a_E = \frac{\Gamma_e A_e}{\delta x_{PE}} \tag{15b} aE?=δxPE?Γe?Ae??(15b)
a S = Γ s A s δ y S P (15c) a_S = \frac{\Gamma_s A_s}{\delta y_{SP}} \tag{15c} aS?=δySP?Γs?As??(15c)
a N = Γ n A n δ y P N (15d) a_N = \frac{\Gamma_n A_n}{\delta y_{PN}} \tag{15d} aN?=δyPN?Γn?An??(15d)
a B = Γ b A b δ z B P (15e) a_B = \frac{\Gamma_b A_b}{\delta z_{BP}} \tag{15e} aB?=δzBP?Γb?Ab??(15e)
a T = Γ t A t δ z P T (15f) a_T = \frac{\Gamma_t A_t}{\delta z_{PT}} \tag{15f} aT?=δzPT?Γt?At??(15f)
a P = a W + a E + a S + a N + a B + a T ? S P (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)
总结 【有限体积法（2）——二维、三维扩散方程的离散推导】通过有限体积法离散后的稳态扩散方程可以统一写成如下形式，
a P ? P = Σ a n b ? n b + S u (16) a_P \phi_P = \Sigma a_{nb}\phi_{nb} + S_u \tag{16} aP??P?=Σanb??nb?+Su?(16)
其中，“ n b nb nb”代表单元 P P P相邻的节点， Σ \Sigma Σ代表所有相邻节点之和。
系数 a P a_P aP?与 a n b a_{nb} anb?之间的关系为：
a P = Σ a n b ? S P (17) a_P = \Sigma a_{nb} - S_P \tag{17} aP?=Σanb??SP?(17)
相邻节点的系数 a n b a_{nb} anb?计算如下表，

	a W a_W aW?	a E a_E aE?	a S a_S aS?	a N a_N aN?	a B a_B aB?	a T a_T aT?
1D	Γ w A w δ x W P \frac{\Gamma_w A_w}{\delta x_{WP}} δxWP?Γw?Aw??	Γ w A e δ x P E \frac{\Gamma_w A_e}{\delta x_{PE}} δxPE?Γw?Ae??
2D	Γ w A w δ x W P \frac{\Gamma_w A_w}{\delta x_{WP}} δxWP?Γw?Aw??	Γ w A e δ x P E \frac{\Gamma_w A_e}{\delta x_{PE}} δxPE?Γw?Ae??	Γ s A s δ y S P \frac{\Gamma_s A_s}{\delta y_{SP}} δySP?Γs?As??	Γ n A n δ y P N \frac{\Gamma_n A_n}{\delta y_{PN}} δyPN?Γn?An??
3D	Γ w A w δ x W P \frac{\Gamma_w A_w}{\delta x_{WP}} δxWP?Γw?Aw??	Γ w A e δ x P E \frac{\Gamma_w A_e}{\delta x_{PE}} δxPE?Γw?Ae??	Γ s A s δ y S P \frac{\Gamma_s A_s}{\delta y_{SP}} δySP?Γs?As??	Γ n A n δ y P N \frac{\Gamma_n A_n}{\delta y_{PN}} δyPN?Γn?An??	Γ b A b δ z B P \frac{\Gamma_b A_b}{\delta z_{BP}} δzBP?Γb?Ab??	Γ t A t δ z P T \frac{\Gamma_t A_t}{\delta z_{PT}} δzPT?Γt?At??