twoPhaseEulerFoam全解读之一(转载)
本系列将对OpenFOAM-2.1.1 中的 twoPhaseEulerFoam 求解器进行完全解读,共分三部分:方程推导,代码解读,补充说明。本篇进行方程推导,详细介绍如果从双流体模型出发得到 twoPhaseEulerFoam 中的 UEqn.H 对应的模型方程形式。
方程推导
双流体模型方程可以表达成如下形式:
连续性方程:
∂
(
α
ϕ
ρ
ϕ
)
∂
t
+
∇
⋅
(
α
ϕ
ρ
ϕ
U
ϕ
)
=
0
\frac{\partial(\alpha_\phi\rho_\phi)}{\partial t}+\nabla\cdot(\alpha_\phi\rho_\phi U_\phi)=0
∂t∂(αϕρϕ)+∇⋅(αϕρϕUϕ)=0
动量守恒方程:
∂
(
α
ϕ
ρ
ϕ
U
ϕ
)
∂
t
+
∇
⋅
(
α
ϕ
ρ
ϕ
U
ϕ
U
ϕ
)
+
∇
⋅
(
α
ϕ
τ
ϕ
)
+
∇
⋅
(
α
ϕ
ρ
ϕ
R
ϕ
)
=
−
α
ϕ
∇
p
+
α
ϕ
ρ
ϕ
g
+
M
ϕ
\frac{\partial(\alpha_\phi\rho_\phi U_\phi)}{\partial t}+\nabla\cdot(\alpha_\phi\rho_\phi U_\phi U_\phi)+\nabla\cdot(\alpha_\phi\tau_\phi)+\nabla\cdot(\alpha_\phi\rho_\phi R_\phi )=-\alpha_\phi\nabla p+\alpha_\phi\rho_\phi g+M_\phi
∂t∂(αϕρϕUϕ)+∇⋅(αϕρϕUϕUϕ)+∇⋅(αϕτϕ)+∇⋅(αϕρϕRϕ)=−αϕ∇p+αϕρϕg+Mϕ
式中,下标
ϕ
=
a
,
b
\phi=a,b
ϕ=a,b分别代表分散相和连续相,
τ
ϕ
\tau_\phi
τϕ表示粘性应力项,
R
ϕ
R_\phi
Rϕ表示雷诺应力项,
M
ϕ
M_\phi
Mϕ表示相间作用项。
上述方程是完全守恒形式的,但是注意到上述动量方程的瞬变项是
∂
(
α
ϕ
ρ
ϕ
U
ϕ
)
∂
t
\frac{\partial(\alpha_\phi\rho_\phi U_\phi)}{\partial t}
∂t∂(αϕρϕUϕ),等于说解这个方程能得到的是每个时间步的动量,若要转化成速度,则需要用动量除以密度与体积分率的乘积,即
(
α
ϕ
ρ
ϕ
U
ϕ
)
α
ϕ
ρ
ϕ
\frac{(\alpha_\phi\rho_\phi U_\phi)}{\alpha_\phi\rho_\phi}
αϕρϕ(αϕρϕUϕ)。那么当离散相a的体积分率
α
a
→
0
\alpha_a\to0
αa→0时,这个除法就要出问题了。于是,Weller [1] 提出通过构造一种”phase-intensive”形式的动量方程来避开这个问题,见下面的详细推导。
Weller提出的方法的核心是将
α
ϕ
ρ
ϕ
\alpha_\phi\rho_\phi
αϕρϕ从动量方程的瞬变项中剥离出来,以使动量方程直接对速度进行演化,而不是动量。
首先对动量方程的瞬变项和对流项进行如下转化:
∂
(
α
ϕ
ρ
ϕ
U
ϕ
)
∂
t
=
α
ϕ
ρ
ϕ
∂
(
U
ϕ
)
∂
t
+
U
ϕ
∂
(
α
ϕ
ρ
ϕ
)
∂
t
\frac{\partial(\alpha_\phi\rho_\phi U_\phi)}{\partial t}=\alpha_\phi\rho_\phi\frac{\partial( U_\phi)}{\partial t}+U_\phi\frac{\partial(\alpha_\phi\rho_\phi )}{\partial t}
∂t∂(αϕρϕUϕ)=αϕρϕ∂t∂(Uϕ)+Uϕ∂t∂(αϕρϕ)
∇
⋅
(
α
ϕ
ρ
ϕ
U
ϕ
U
ϕ
)
=
α
ϕ
ρ
ϕ
U
ϕ
⋅
∇
(
U
ϕ
)
+
U
ϕ
∇
⋅
(
α
ϕ
ρ
ϕ
U
ϕ
)
\nabla\cdot(\alpha_\phi\rho_\phi U_\phi U_\phi)= \alpha_\phi\rho_\phi U_\phi\cdot \nabla( U_\phi) + U_\phi\nabla\cdot(\alpha_\phi\rho_\phi U_\phi)
∇⋅(αϕρϕUϕUϕ)=αϕρϕUϕ⋅∇(Uϕ)+Uϕ∇⋅(αϕρϕUϕ)
注:这里到了张量运算公式
[
∇
⋅
v
w
]
=
[
v
⋅
∇
w
]
+
w
(
∇
⋅
v
)
[\nabla\cdot \mathbf{vw}]=[\mathbf{v}\cdot\nabla\mathbf{w}]+\mathbf{w}(\nabla\cdot\mathbf{v})
[∇⋅vw]=[v⋅∇w]+w(∇⋅v),具体可参考 Bird 的 Transport Phenomenon 的 Appendix A。
于是,瞬变项和对流项的加和可以写成如下形式:
∂
(
α
ϕ
ρ
ϕ
U
ϕ
)
∂
t
+
∇
⋅
(
α
ϕ
ρ
ϕ
U
ϕ
U
ϕ
)
=
α
ϕ
ρ
ϕ
[
∂
(
U
ϕ
)
∂
t
+
U
ϕ
⋅
∇
(
U
ϕ
)
]
+
U
ϕ
[
∂
(
α
ϕ
ρ
ϕ
)
∂
t
+
∇
⋅
(
α
ϕ
ρ
ϕ
U
ϕ
)
]
\frac{\partial(\alpha_\phi\rho_\phi U_\phi)}{\partial t}+\nabla\cdot(\alpha_\phi\rho_\phi U_\phi U_\phi)=\alpha_\phi\rho_\phi\left[\frac{\partial( U_\phi)}{\partial t}+U_\phi \cdot \nabla( U_\phi)\right]+U_\phi \left[ \frac{\partial(\alpha_\phi\rho_\phi )}{\partial t}+\nabla\cdot(\alpha_\phi\rho_\phi U_\phi)\right]
∂t∂(αϕρϕUϕ)+∇⋅(αϕρϕUϕUϕ)=αϕρϕ[∂t∂(Uϕ)+Uϕ⋅∇(Uϕ)]+Uϕ[∂t∂(αϕρϕ)+∇⋅(αϕρϕUϕ)]
注意右边第二项的括号里其实就是连续性方程的左边,其值为0,因此得到:
∂
(
α
ϕ
ρ
ϕ
U
ϕ
)
∂
t
+
∇
⋅
(
α
ϕ
ρ
ϕ
U
ϕ
U
ϕ
)
=
α
ϕ
ρ
ϕ
[
∂
(
U
ϕ
)
∂
t
+
U
ϕ
⋅
∇
(
U
ϕ
)
]
\frac{\partial(\alpha_\phi\rho_\phi U_\phi)}{\partial t}+\nabla\cdot(\alpha_\phi\rho_\phi U_\phi U_\phi) = \alpha_\phi\rho_\phi\left[\frac{\partial( U_\phi)}{\partial t}+U_\phi \cdot \nabla( U_\phi)\right]
∂t∂(αϕρϕUϕ)+∇⋅(αϕρϕUϕUϕ)=αϕρϕ[∂t∂(Uϕ)+Uϕ⋅∇(Uϕ)]
于是得到第一步转化之后的动量方程:
α
ϕ
ρ
ϕ
[
∂
(
U
ϕ
)
∂
t
+
U
ϕ
⋅
∇
(
U
ϕ
)
]
+
∇
⋅
(
α
ϕ
τ
ϕ
)
+
∇
⋅
(
α
ϕ
ρ
ϕ
R
ϕ
)
=
−
α
ϕ
∇
p
+
α
ϕ
ρ
ϕ
g
+
M
ϕ
\alpha_\phi\rho_\phi\left[\frac{\partial( U_\phi)}{\partial t}+U_\phi\cdot\nabla( U_\phi)\right] + \nabla\cdot(\alpha_\phi\tau_\phi) + \nabla\cdot(\alpha_\phi\rho_\phi R_\phi ) = -\alpha_\phi\nabla p + \alpha_\phi\rho_\phi g + M_\phi
αϕρϕ[∂t∂(Uϕ)+Uϕ⋅∇(Uϕ)]+∇⋅(αϕτϕ)+∇⋅(αϕρϕRϕ)=−αϕ∇p+αϕρϕg+Mϕ
下面处理粘性应力项和雷诺应力项。
∇ ⋅ ( α ϕ τ ϕ ) + ∇ ⋅ ( α ϕ ρ ϕ R ϕ ) = ∇ ⋅ [ α ϕ ρ ϕ ( τ ϕ ρ ϕ + R ϕ ) ] = ∇ ⋅ [ α ϕ ρ ϕ R e f f , ϕ ] \nabla\cdot(\alpha_\phi\tau_\phi) + \nabla\cdot(\alpha_\phi\rho_\phi R_\phi )=\nabla\cdot\left[\alpha_\phi\rho_\phi(\frac{\tau_\phi}{\rho_\phi}+R_\phi)\right] = \nabla\cdot\left[\alpha_\phi\rho_\phi R_{eff,\phi}\right ] ∇⋅(αϕτϕ)+∇⋅(αϕρϕRϕ)=∇⋅[αϕρϕ(ρϕτϕ+Rϕ)]=∇⋅[αϕρϕReff,ϕ]
其中 R e f f , ϕ = τ ϕ ρ ϕ + R ϕ R_{eff,\phi}=\frac{\tau_\phi}{\rho_\phi}+R_\phi Reff,ϕ=ρϕτϕ+Rϕ。
根据定义(此处参考BubbleFoam的Wiki页面):
τ
ϕ
=
−
ρ
ϕ
ν
ϕ
[
∇
U
ϕ
+
∇
T
U
ϕ
]
+
2
3
ρ
ϕ
ν
ϕ
(
∇
⋅
U
ϕ
)
I
\boldsymbol{\tau}_{\phi} = - \rho_{\phi} \nu_{\phi} \left[\nabla \mathbf{U}_{\phi} + \nabla^{\textrm{T}} \mathbf{U}_{\phi} \right] + \frac{2}{3}\rho_{\phi}\nu_{\phi} \left( \nabla \cdot \mathbf{U}_{\phi} \right) \mathbf{I}
τϕ=−ρϕνϕ[∇Uϕ+∇TUϕ]+32ρϕνϕ(∇⋅Uϕ)I
以及
R
ϕ
=
−
ν
ϕ
,
t
[
∇
U
ϕ
+
∇
T
U
ϕ
]
+
2
3
ν
ϕ
,
t
(
∇
⋅
U
ϕ
)
I
+
2
3
k
ϕ
I
\mathbf{R}_{\phi} = - \nu_{\phi,\textrm{t}} \left[ \nabla \mathbf{U}_{\phi} +\nabla^{\textrm{T}} \mathbf{U}_{\phi} \right] + \frac{2}{3} \nu_{\phi,\textrm{t}} \left( \nabla \cdot \mathbf{U}_{\phi} \right) \mathbf{I} + \frac{2}{3} k_{\phi} \mathbf{I}
Rϕ=−νϕ,t[∇Uϕ+∇TUϕ]+32νϕ,t(∇⋅Uϕ)I+32kϕI
代入到
R
e
f
f
,
ϕ
R_{eff,\phi}
Reff,ϕ中,得:
R
e
f
f
,
ϕ
=
−
(
ν
ϕ
+
ν
ϕ
,
t
)
[
∇
U
ϕ
+
∇
T
U
ϕ
]
+
2
3
(
ν
ϕ
+
ν
ϕ
,
t
)
(
∇
⋅
U
ϕ
)
I
+
2
3
k
ϕ
I
R_{eff,\phi}=-(\nu_\phi+\nu_{\phi , t})\left[ \nabla \mathbf{U}_{\phi} +\nabla^{\textrm{T}} \mathbf{U}_{\phi} \right]+\frac{2}{3}(\nu_\phi+\nu_{\phi , t}) \left (\nabla \cdot \mathbf{U}_{\phi}\right ) \mathbf{I} + \frac{2}{3} k_{\phi} \mathbf{I}
Reff,ϕ=−(νϕ+νϕ,t)[∇Uϕ+∇TUϕ]+32(νϕ+νϕ,t)(∇⋅Uϕ)I+32kϕI
令
ν
e
f
f
=
ν
ϕ
+
ν
ϕ
,
t
\nu_{eff}=\nu_\phi+\nu_{\phi , t}
νeff=νϕ+νϕ,t ,则:
R
e
f
f
,
ϕ
=
−
ν
e
f
f
[
∇
U
ϕ
+
∇
T
U
ϕ
]
+
2
3
ν
e
f
f
(
∇
⋅
U
ϕ
)
I
+
2
3
k
ϕ
I
=
−
ν
e
f
f
∇
U
ϕ
+
R
c
,
ϕ
R_{eff,\phi}=-\nu_{eff}\left[ \nabla \mathbf{U}_{\phi} +\nabla^{\textrm{T}} \mathbf{U}_{\phi} \right]+\frac{2}{3}\nu_{eff} \left (\nabla \cdot \mathbf{U}_{\phi}\right ) \mathbf{I} + \frac{2}{3} k_{\phi} \mathbf{I} = -\nu_{eff}\nabla U_\phi + R_{c,\phi}
Reff,ϕ=−νeff[∇Uϕ+∇TUϕ]+32νeff(∇⋅Uϕ)I+32kϕI=−νeff∇Uϕ+Rc,ϕ
其中
R
c
,
ϕ
=
−
ν
e
f
f
∇
U
ϕ
T
+
2
3
ν
e
f
f
(
∇
⋅
U
ϕ
)
I
+
2
3
k
ϕ
I
R_{c,\phi}=-\nu_{eff} \nabla \mathbf{U}^\textrm{T}_{\phi}+\frac{2}{3}\nu_{eff} \left (\nabla \cdot \mathbf{U}_{\phi}\right ) \mathbf{I} + \frac{2}{3} k_{\phi} \mathbf{I}
Rc,ϕ=−νeff∇UϕT+32νeff(∇⋅Uϕ)I+32kϕI
于是得到:
∇
⋅
[
α
ϕ
ρ
ϕ
R
e
f
f
,
ϕ
]
=
∇
(
α
ϕ
ρ
ϕ
)
⋅
[
R
e
f
f
,
ϕ
]
+
α
ϕ
ρ
ϕ
∇
⋅
[
R
e
f
f
,
ϕ
]
=
α
ϕ
ρ
ϕ
∇
⋅
[
−
ν
e
f
f
∇
U
ϕ
]
+
α
ϕ
ρ
ϕ
∇
⋅
[
R
c
,
ϕ
]
+
∇
(
α
ϕ
ρ
ϕ
)
[
−
ν
e
f
f
∇
U
ϕ
+
R
c
,
ϕ
]
\begin{aligned} \nabla\cdot\left[\alpha_\phi\rho_\phi R_{eff,\phi}\right ] = & \nabla(\alpha_\phi\rho_\phi)\cdot\left[ R_{eff,\phi}\right] + \alpha_\phi\rho_\phi\nabla\cdot \left [ R_{eff,\phi}\right ]\\ =& \alpha_\phi\rho_\phi\nabla\cdot\left[ -\nu_{eff}\nabla U_\phi\right] + \alpha_\phi\rho_\phi\nabla\cdot\left[ R_{c,\phi}\right] + \nabla(\alpha_\phi\rho_\phi)\left[ -\nu_{eff}\nabla U_\phi + R_{c,\phi}\right] \end{aligned}
∇⋅[αϕρϕReff,ϕ]==∇(αϕρϕ)⋅[Reff,ϕ]+αϕρϕ∇⋅[Reff,ϕ]αϕρϕ∇⋅[−νeff∇Uϕ]+αϕρϕ∇⋅[Rc,ϕ]+∇(αϕρϕ)[−νeff∇Uϕ+Rc,ϕ]
代入到动量方程中,并且方程两边同时除以
α
ϕ
ρ
ϕ
\alpha_\phi\rho_\phi
αϕρϕ,得到:
∂
U
ϕ
∂
t
+
U
ϕ
⋅
∇
U
ϕ
−
∇
⋅
[
ν
e
f
f
∇
U
ϕ
]
+
∇
⋅
[
R
c
,
ϕ
]
+
∇
(
α
ϕ
ρ
ϕ
)
α
ϕ
ρ
ϕ
⋅
[
−
ν
e
f
f
∇
U
ϕ
+
R
c
,
ϕ
]
=
−
∇
p
ρ
ϕ
+
g
+
M
ϕ
α
ϕ
ρ
ϕ
\frac{\partial U_\phi}{\partial t} + U_\phi\cdot\nabla U_\phi -\nabla \cdot \left[ \nu_{eff} \nabla U_\phi \right ] + \nabla \cdot \left[ R_{c,\phi}\right] + \frac{\nabla(\alpha_\phi\rho_\phi)}{\alpha_\phi\rho_\phi}\cdot \left[ -\nu_{eff}\nabla U_\phi + R_{c,\phi}\right] = -\frac{\nabla p}{\rho_\phi} + g + \frac{M_\phi}{\alpha_\phi\rho_\phi}
∂t∂Uϕ+Uϕ⋅∇Uϕ−∇⋅[νeff∇Uϕ]+∇⋅[Rc,ϕ]+αϕρϕ∇(αϕρϕ)⋅[−νeff∇Uϕ+Rc,ϕ]=−ρϕ∇p+g+αϕρϕMϕ
如果假定两相流体均为不可压缩,密度恒为常数,于是可以得到不可压缩的双流体模型的方程组:
连续性方程
∂
(
α
ϕ
)
∂
t
+
∇
⋅
(
α
ϕ
U
ϕ
)
=
0
\frac{\partial(\alpha_\phi)}{\partial t}+\nabla\cdot(\alpha_\phi U_\phi)=0
∂t∂(αϕ)+∇⋅(αϕUϕ)=0
动量方程
∂
U
ϕ
∂
t
+
U
ϕ
⋅
∇
U
ϕ
−
∇
⋅
[
ν
e
f
f
∇
U
ϕ
]
+
∇
⋅
[
R
c
,
ϕ
]
+
∇
(
α
ϕ
)
α
ϕ
⋅
[
−
ν
e
f
f
∇
U
ϕ
+
R
c
,
ϕ
]
=
−
∇
p
ρ
ϕ
+
g
+
M
ϕ
α
ϕ
ρ
ϕ
\frac{\partial U_\phi}{\partial t} + U_\phi\cdot\nabla U_\phi -\nabla \cdot \left[ \nu_{eff} \nabla U_\phi \right ] + \nabla \cdot \left[ R_{c,\phi}\right] + \frac{\nabla(\alpha_\phi)}{\alpha_\phi} \cdot \left[ -\nu_{eff}\nabla U_\phi + R_{c,\phi}\right] = -\frac{\nabla p}{\rho_\phi} + g + \frac{M_\phi}{\alpha_\phi\rho_\phi}
∂t∂Uϕ+Uϕ⋅∇Uϕ−∇⋅[νeff∇Uϕ]+∇⋅[Rc,ϕ]+αϕ∇(αϕ)⋅[−νeff∇Uϕ+Rc,ϕ]=−ρϕ∇p+g+αϕρϕMϕ
方程中还剩下相间作用项没有处理,对于分散相和连续项形式,相间作用力是大小相等符号想反,这里只考虑分散相的形式,令
ϕ
=
a
\phi=a
ϕ=a,则得到分散相的动量方程:
∂
U
a
∂
t
+
U
a
⋅
∇
U
a
−
∇
⋅
[
ν
e
f
f
∇
U
a
]
+
∇
⋅
[
R
c
,
a
]
+
∇
(
α
a
)
α
a
⋅
[
−
ν
e
f
f
∇
U
a
+
R
c
,
a
]
=
−
∇
p
ρ
a
+
g
+
M
a
α
a
ρ
a
\frac{\partial U_a}{\partial t} + U_a\cdot\nabla U_a -\nabla \cdot \left[ \nu_{eff} \nabla U_a \right ] + \nabla \cdot \left[ R_{c,a}\right] + \frac{\nabla(\alpha_a)}{\alpha_a} \cdot \left[ -\nu_{eff}\nabla U_a + R_{c,a}\right] = -\frac{\nabla p}{\rho_a} + g + \frac{M_a}{\alpha_a\rho_a}
∂t∂Ua+Ua⋅∇Ua−∇⋅[νeff∇Ua]+∇⋅[Rc,a]+αa∇(αa)⋅[−νeff∇Ua+Rc,a]=−ρa∇p+g+αaρaMa
相间作用只考虑曳力,升力以及虚拟质量力,即 M , a = M d r a g + M l i f t + M v m M,a=M_{drag}+M_{lift}+M_{vm} M,a=Mdrag+Mlift+Mvm,下面分别考虑每一种相间作用力。
曳力
M
d
r
a
g
=
−
β
(
U
a
−
U
b
)
M_{drag}=-\beta(U_a-U_b)
Mdrag=−β(Ua−Ub),其中
β
\beta
β为曳力系数。
升力
M
l
i
f
t
=
−
α
a
α
b
C
l
(
α
b
ρ
b
+
α
a
ρ
a
)
U
r
×
(
∇
×
U
)
M_{lift}=-\alpha_a\alpha_b C_l (\alpha_b \rho_b + \alpha_a \rho_a)U_r \times (\nabla \times U)
Mlift=−αaαbCl(αbρb+αaρa)Ur×(∇×U) ,其中
U
r
=
U
a
−
U
b
U_r=U_a-U_b
Ur=Ua−Ub,
U
=
α
a
U
a
+
α
b
U
b
U=\alpha_a U_a + \alpha_b U_b
U=αaUa+αbUb
虚拟质量力
M
v
m
=
α
a
α
b
C
v
m
ρ
b
[
D
U
b
D
t
−
D
U
a
D
t
]
M_{vm}=\alpha_a\alpha_b C_{vm}\rho_b\left[ \frac{DU_b}{Dt}-\frac{DU_a}{Dt}\right]
Mvm=αaαbCvmρb[DtDUb−DtDUa],其中
D
D
t
\frac{D}{Dt}
DtD表示物质导数,
D
U
b
D
t
=
∂
U
b
∂
t
+
U
b
⋅
∇
U
b
\frac{DU_b}{Dt}=\frac{\partial U_b}{\partial t} + U_b \cdot \nabla U_b
DtDUb=∂t∂Ub+Ub⋅∇Ub,
D
U
a
D
t
=
∂
U
a
∂
t
+
U
a
⋅
∇
U
a
\frac{DU_a}{Dt}=\frac{\partial U_a}{\partial t}+U_a \cdot \nabla U_a
DtDUa=∂t∂Ua+Ua⋅∇Ua
考虑到形式的统一,令
K
=
β
α
a
α
b
K=\frac{\beta}{\alpha_a\alpha_b}
K=αaαbβ,则曳力可表示为
M
d
r
a
g
=
−
α
a
α
b
K
(
U
a
−
U
b
)
M_{drag}=-\alpha_a\alpha_b K(U_a-U_b)
Mdrag=−αaαbK(Ua−Ub)
代入到分散相a的动量方程中,得到:
∂
U
a
∂
t
+
U
a
⋅
∇
U
a
−
∇
⋅
[
ν
e
f
f
∇
U
a
]
+
∇
⋅
[
R
c
,
a
]
+
∇
(
α
a
)
α
a
⋅
[
−
ν
e
f
f
∇
U
a
+
R
c
,
a
]
=
−
∇
p
ρ
a
+
g
−
α
b
ρ
a
K
(
U
a
−
U
b
)
−
α
b
ρ
a
C
l
(
α
b
ρ
b
+
α
a
ρ
a
)
U
r
×
(
∇
×
U
)
+
α
b
ρ
a
C
v
m
ρ
b
[
∂
U
b
∂
t
+
U
b
⋅
∇
U
b
−
(
∂
U
a
∂
t
+
U
a
⋅
∇
U
a
)
]
\begin{aligned} &\frac{\partial U_a}{\partial t} + U_a\cdot \nabla U_a -\nabla \cdot \left[ \nu_{eff} \nabla U_a \right ] + \nabla \cdot \left[ R_{c,a}\right] + \frac{\nabla(\alpha_a)}{\alpha_a} \cdot \left[ -\nu_{eff}\nabla U_a + R_{c,a}\right] \\ = & -\frac{\nabla p}{\rho_a} + g - \frac{\alpha_b}{\rho_a} K (U_a-U_b) - \frac{\alpha_b}{\rho_a} C_l (\alpha_b \rho_b + \alpha_a \rho_a) U_r \times (\nabla \times U) \\ +& \frac{\alpha_b}{\rho_a} C_{vm}\rho_b\left[ \frac{\partial U_b}{\partial t} + U_b \cdot \nabla U_b - (\frac{\partial U_a}{\partial t}+U_a \cdot \nabla U_a)\right] \end{aligned}
=+∂t∂Ua+Ua⋅∇Ua−∇⋅[νeff∇Ua]+∇⋅[Rc,a]+αa∇(αa)⋅[−νeff∇Ua+Rc,a]−ρa∇p+g−ρaαbK(Ua−Ub)−ρaαbCl(αbρb+αaρa)Ur×(∇×U)ρaαbCvmρb[∂t∂Ub+Ub⋅∇Ub−(∂t∂Ua+Ua⋅∇Ua)]
将相关的项合并,并调整顺序,便得到与twoPhaseEulerFoam求解器的UEqn.H文件中相同形式的分散相动量方程:
(
1
+
α
b
ρ
b
ρ
a
C
v
m
)
(
∂
U
a
∂
t
+
U
a
⋅
∇
U
a
)
−
∇
⋅
[
ν
e
f
f
∇
U
a
]
+
∇
⋅
[
R
c
,
a
]
+
∇
(
α
a
)
α
a
⋅
[
−
ν
e
f
f
∇
U
a
+
R
c
,
a
]
=
−
α
b
ρ
a
K
U
a
−
α
b
ρ
a
{
C
l
(
α
b
ρ
b
+
α
a
ρ
a
)
U
r
×
(
∇
×
U
)
−
C
v
m
ρ
b
[
∂
U
b
∂
t
+
U
b
⋅
∇
U
b
]
}
−
∇
p
ρ
a
+
g
+
α
b
ρ
a
K
U
b
\begin{aligned} &(1+\frac{\alpha_b \rho_b}{\rho_a} C_{vm})(\frac{\partial U_a}{\partial t} + U_a\cdot \nabla U_a ) -\nabla \cdot \left[ \nu_{eff} \nabla U_a \right ] + \nabla \cdot \left[ R_{c,a}\right] + \frac{\nabla(\alpha_a)}{\alpha_a} \cdot \left[ -\nu_{eff}\nabla U_a + R_{c,a}\right] \\ = & -\frac{\alpha_b}{\rho_a} K U_a - \frac{\alpha_b}{\rho_a} \left\{ {C_l (\alpha_b \rho_b + \alpha_a \rho_a) U_r \times (\nabla \times U) - C_{vm}\rho_b\left[ {\frac{\partial U_b}{\partial t} + U_b \cdot \nabla U_b }\right] } \right\} \\ &- \frac{\nabla p}{\rho_a} + g + \frac{\alpha_b}{\rho_a} K U_b \end{aligned}
=(1+ρaαbρbCvm)(∂t∂Ua+Ua⋅∇Ua)−∇⋅[νeff∇Ua]+∇⋅[Rc,a]+αa∇(αa)⋅[−νeff∇Ua+Rc,a]−ρaαbKUa−ρaαb{Cl(αbρb+αaρa)Ur×(∇×U)−Cvmρb[∂t∂Ub+Ub⋅∇Ub]}−ρa∇p+g+ρaαbKUb
连续相b的动量方程形式相仿,这里就不再重复了。这里有几点注意事项:
此处的双流体模型在推导的过程中,是把a当作分散相,b当作连续相的。分散相的体积分率
α
a
\alpha_a
αa可以等于0,但是连续项体积分率
α
b
\alpha_b
αb不能等于0 ,否则会出问题。
曳力系数
β
\beta
β 的形式就是文献中常见的形式,比如,WenYu 曳力系数
β
=
3
4
(
1
−
α
b
)
α
b
d
p
,
a
∣
U
b
−
U
a
∣
C
D
0
α
b
−
2.7
\beta=\frac{3}{4}\frac{(1-\alpha_b)\alpha_b}{d_{p,a}}|U_b-U_a|C_{D0}\alpha_b^{-2.7}
β=43dp,a(1−αb)αb∣Ub−Ua∣CD0αb−2.7,Ergun 曳力系数
β
=
150
(
1
−
α
b
)
2
μ
b
α
b
d
a
2
+
1.75
(
1
−
α
b
)
ρ
b
U
b
−
U
a
d
a
\beta=150\frac{(1-\alpha_b)^2\mu_b}{\alpha_b d_a^2}+1.75\frac{(1-\alpha_b)\rho_b{U_b-U_a}}{d_a}
β=150αbda2(1−αb)2μb+1.75da(1−αb)ρbUb−Ua。而程序中定义的
K
=
β
α
a
α
b
K=\frac{\beta}{\alpha_a\alpha_b}
K=αaαbβ,所以,当
α
b
→
0
\alpha_b\to 0
αb→0时,如果用WenYu曳力那还不会出错,因为曳力系数中的分子里同时含有
α
a
α
b
\alpha_a\alpha_b
αaαb,运算
K
=
β
α
a
α
b
K=\frac{\beta}{\alpha_a\alpha_b}
K=αaαbβ不会出现除以0的问题;但如果用Ergun曳力,那就要出问题了,因为Ergun曳力系数中两项的分子都没有
α
b
\alpha_b
αb,所以运算
K
=
β
α
a
α
b
K=\frac{\beta}{\alpha_a\alpha_b}
K=αaαbβ就要出问题了。