稳态扩散方程:
∇ ⋅ ( Γ ∇ ϕ ) + 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+∫CVSϕdV=∫A~n⋅(Γ∇ϕ)dA+∫CVSϕ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 ) i + ( Γ ∂ ϕ ∂ 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[ \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)
其中 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=ΓwAwδ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=ΓeAeδxPE(ϕE−ϕP)(5b)