twoPhaseEulerFoam 全解读之二（转）
本系列将对OpenFOAM-2.1.1 中的 twoPhaseEulerFoam 求解器进行完全解读，共分三部分：方程推导，代码解读，补充说明。本篇对 twoPhaseEulerFoam 中的 UEqn.H 和 pEqn.H 中的代码进行详细地的解读。
代码解读
UEqn
前一篇导出了分散相的动量守恒方程
(
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
这一篇分析twoPhaseEulerFoam求解器是怎么来对动量方程进行离散的，以及，如果通过构建压力方程来对速度进行修正以保证两相的连续性。
注意上述动量方程中，有两项还需要处理一下：
U
a
⋅
∇
U
a
=
∇
⋅
(
U
a
U
a
)
−
U
a
(
∇
⋅
U
a
)
U_a \cdot \nabla U_a=\nabla\cdot(U_aU_a)-U_a(\nabla\cdot U_a)
Ua⋅∇Ua=∇⋅(UaUa)−Ua(∇⋅Ua)
∇
(
α
a
)
α
a
⋅
[
−
ν
e
f
f
∇
U
a
]
=
∇
⋅
[
−
ν
e
f
f
(
∇
α
a
)
α
a
U
a
]
−
U
a
(
∇
⋅
(
−
ν
e
f
f
∇
α
a
α
a
)
)
\frac{\nabla(\alpha_a)}{\alpha_a}\cdot\left[-\nu_{eff} \nabla U_a\right] = \nabla \cdot\left[ -\nu_{eff}\frac{(\nabla \alpha_a)}{\alpha_a} U_a \right]-U_a\left(\nabla\cdot(-\nu_{eff}\frac{\nabla \alpha_a}{\alpha_a}) \right)
αa∇(αa)⋅[−νeff∇Ua]=∇⋅[−νeffαa(∇αa)Ua]−Ua(∇⋅(−νeffαa∇αa))
转化前后的形式从数学上来等价的，但是在有限体积离散过程中，转化前的
U
a
⋅
∇
U
a
U_a \cdot \nabla U_a
Ua⋅∇Ua 和
∇
(
α
a
)
α
a
⋅
[
−
ν
e
f
f
∇
U
a
]
\frac{\nabla(\alpha_a)}{\alpha_a}\cdot\left[-\nu_{eff} \nabla U_a\right]
αa∇(αa)⋅[−νeff∇Ua] 对于
U
a
U_a
Ua来说是非守恒的，转化后的形式是守恒的(参考这篇文献，注意这样转化后，动量方程空间上是守恒的，但时间上仍是不守恒的，这个帖子也有一些有价值的信息)。
这样转化以后，得到的动量方程就跟twoPhaseEulerFoam里的定义是一模一样了:
下面将动量方程的每一项与twoPhaseEulerFoam的UEqn.H的代码一一对应。
(scalar(1) + Cvm*rhob*beta/rhoa)*
(
fvm::ddt(Ua)
+ fvm::div(phia, Ua, "div(phia,Ua)")
- fvm::Sp(fvc::div(phia), Ua)
)
fvm::ddt(Ua)对应 ∂ U a ∂ t \frac{\partial U_a}{\partial t} ∂t∂Ua，phia定义为fvc::interpolate(Ua) & mesh.Sf()，于是fvm::div(phia, Ua, “div(phia,Ua)”) 和 fvm::Sp(fvc::div(phia), Ua) 便分别对应 ∇ ⋅ ( U a U a ) \nabla\cdot(U_aU_a) ∇⋅(UaUa) 和 U a ( ∇ ⋅ U a ) U_a(\nabla\cdot U_a) Ua(∇⋅Ua)了 [ 注一 ]。
- fvm::laplacian(nuEffa, Ua)
+ fvc::div(Rca)
+ fvm::div(phiRa, Ua, "div(phia,Ua)")
- fvm::Sp(fvc::div(phiRa), Ua)
fvm::laplacian(nuEffa, Ua)对应
∇
⋅
[
ν
e
f
f
∇
U
a
]
\nabla \cdot \left[ \nu_{eff} \nabla U_a \right ]
∇⋅[νeff∇Ua]，fvc::div(Rca)对应
∇
⋅
[
R
c
,
a
]
\nabla \cdot \left[ R_{c,a}\right]
∇⋅[Rc,a]。
phiRa的定义是
-fvc::interpolate(nuEffa)*mesh.magSf()*fvc::snGrad(alpha)
/fvc::interpolate(alpha + scalar(0.001))
相当于
−
ν
e
f
f
∇
α
a
α
a
-\nu_{eff}\frac{\nabla \alpha_a}{\alpha_a}
−νeffαa∇αa [ 注二 ]。
于是 fvm::div(phiRa, Ua, “div(phia,Ua)”) 和 fvm::Sp(fvc::div(phiRa), Ua) 便分别对应 $\nabla \cdot\left[ -\nu_{eff} \frac{(\nabla \alpha_a)}{\alpha_a}(\nabla U_a)\right] $ 和
U
a
(
∇
⋅
(
−
ν
e
f
f
∇
α
a
α
a
)
)
U_a\left(\nabla\cdot(-\nu_{eff}\frac{\nabla \alpha_a}{\alpha_a}) \right)
Ua(∇⋅(−νeffαa∇αa))。
(fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.001)) & Rca)
对应 ∇ ( α a ) α a ⋅ [ R c , a ] \frac{\nabla(\alpha_a)}{\alpha_a} \cdot \left[ R_{c,a}\right] αa∇(αa)⋅[Rc,a]，其中&运算符已重载为计算矢量与张量的点乘积 [ 注三 ]。
==
// g // Buoyancy term transfered to p-equation
- fvm::Sp(beta/rhoa*K, Ua)
//+ beta/rhoa*K*Ub // Explicit drag transfered to p-equation
- beta/rhoa*(liftCoeff - Cvm*rhob*DDtUb)
);
fvm::Sp(beta/rhoaK, Ua)对应
α
b
ρ
a
K
U
a
\frac{\alpha_b}{\rho_a} K U_a
ρaαbKUa，beta/rhoa(liftCoeff - CvmrhobDDtUb) 对应
α
b
ρ
a
{
C
l
(
α
b
ρ
b
+
α
a
ρ
a
)
U
r
×
(
∇
×
U
)
−
C
v
m
ρ
b
[
∂
U
b
∂
t
+
U
b
⋅
∇
U
b
]
}
\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\}
ρaαb{Cl(αbρb+αaρa)Ur×(∇×U)−Cvmρb[∂t∂Ub+Ub⋅∇Ub]}
其中变量liftCoeff定义为
volVectorField liftCoeff(Cl*(beta*rhob + alpha*rhoa)*(Ur ^ fvc::curl(U)));
DDtUb定义为
DDtUb =
fvc::ddt(Ub)
+ fvc::div(phib, Ub)
- fvc::div(phib)*Ub;
重力
g
g
g 以及曳力的显式项 $\frac{\alpha_b}{\rho_a} K U_b $ 如注释所述，将会在 pEqn 中考虑，压力梯度项
∇
p
ρ
a
\frac{\nabla p}{\rho_a}
ρa∇p 则将在 pEqn 中用来约束两相的连续性。
至此动量方程的每一项都与UEqn.H的代码对应起来了。
pEqn
压力方程的作用是修正两相速度
U
a
U_a
Ua 和
U
b
U_b
Ub以使速度满足连续性方程。将两相的连续性方程加起来，得到总体的连续性方程如下 [ 注四 ]：
由于
α
a
+
α
b
=
1
\alpha_a+\alpha_b=1
αa+αb=1，于是两相连续性方程等价于
∇
⋅
(
α
a
U
a
+
α
b
U
b
)
=
0
\nabla \cdot (\alpha_a U_a+\alpha_b U_b) = 0
∇⋅(αaUa+αbUb)=0
再来看压力方程是如何构建起来的。
完整的动量方程离散后，可以写作如下的统一形式：
a
p
,
a
U
p
,
a
=
H
(
U
a
)
−
∇
p
ρ
a
+
α
b
ρ
a
K
U
b
+
g
a_{p,a}U_{p,a}=H(U_a)-\frac{\nabla p}{\rho_a}+\frac{\alpha_b}{\rho_a} K U_b +g
ap,aUp,a=H(Ua)−ρa∇p+ρaαbKUb+g
a
p
,
b
U
p
,
b
=
H
(
U
b
)
−
∇
p
ρ
b
+
α
a
ρ
b
K
U
a
+
g
a_{p,b}U_{p,b}=H(U_b)-\frac{\nabla p}{\rho_b}+\frac{\alpha_a}{\rho_b} K U_a +g
ap,bUp,b=H(Ub)−ρb∇p+ρbαaKUa+g
其中
H
(
U
a
)
H(U_a)
H(Ua) 和
H
(
U
b
)
H(U_b)
H(Ub) 包含了动量方程中除 压力梯度项，显式曳力项以及重力项以后所有项的贡献。
由此离散方程可以得到
U
a
U_a
Ua 和
U
b
U_b
Ub 的表达式如下：
U
a
=
1
a
p
,
a
H
(
U
a
)
−
∇
p
a
p
,
a
ρ
a
+
α
b
a
p
,
a
ρ
a
K
U
b
+
1
a
p
,
a
g
U_{a}=\frac{1}{a_{p,a}}H(U_a)-\frac{\nabla p}{a_{p,a}\rho_a}+\frac{\alpha_b}{ a_{p,a} \rho_a} K U_b +\frac{1}{a_{p,a}} g
Ua=ap,a1H(Ua)−ap,aρa∇p+ap,aρaαbKUb+ap,a1g
U b = 1 a p , b H ( U b ) − ∇ p a p , b ρ b + α a a p , b ρ b K U a + 1 a p , b g U_{b}=\frac{1}{a_{p,b}}H(U_b)-\frac{\nabla p}{a_{p,b}\rho_b}+\frac{\alpha_a}{ a_{p,b} \rho_b} K U_a +\frac{1}{a_{p,b}} g Ub=ap,b1H(Ub)−ap,bρb∇p+ap,bρbαaKUa+ap,b1g
如果此 U a U_a Ua 和 U b U_b Ub 是方程组的解，那么它们必须满足整体的连续性方程，即
将压力梯度项移到方程的一边，得到
这便是压力修正方程的原型。
在pEqn.H中，压力方程其实修正的是界面通量，压力方程迭代收敛以后能保证界面通量的连续性。所以，散度表达式需要根据高斯定理写成界面通量之和的形式：
下标
f
_f
f 表示该项将要在代码中用界面上的变量来表示，在OpenFOAM中，即surfaceScalarField，
S
f
S_f
Sf 表示界面的面积矢量，下面的公式里也是一样。
以及
∇
⋅
[
(
α
a
a
p
,
a
ρ
a
+
α
b
a
p
,
b
ρ
b
)
∇
p
]
=
∑
f
(
α
a
a
p
,
a
ρ
a
+
α
b
a
p
,
b
ρ
b
)
f
(
∇
p
)
⋅
S
f
\nabla \cdot \left[ (\frac{\alpha_a }{a_{p,a}\rho_a} + \frac{\alpha_b }{a_{p,b}\rho_b}) \nabla p \right ] = \sum_f (\frac{\alpha_a }{a_{p,a}\rho_a} + \frac{\alpha_b }{a_{p,b}\rho_b})_f (\nabla p) \cdot S_f
∇⋅[(ap,aρaαa+ap,bρbαb)∇p]=f∑(ap,aρaαa+ap,bρbαb)f(∇p)⋅Sf
下面是pEqn定义了几个跟界面通量有关的变量：
surfaceScalarField alphaf(fvc::interpolate(alpha));
surfaceScalarField betaf(scalar(1) - alphaf);
volScalarField rUaA(1.0/UaEqn.A());
volScalarField rUbA(1.0/UbEqn.A());
phia == (fvc::interpolate(Ua) & mesh.Sf());
phib == (fvc::interpolate(Ub) & mesh.Sf());
rUaAf = fvc::interpolate(rUaA);
surfaceScalarField rUbAf(fvc::interpolate(rUbA));
Ua = rUaA*UaEqn.H();
Ub = rUbA*UbEqn.H();
surfaceScalarField phiDraga
(
fvc::interpolate(beta/rhoa*K*rUaA)*phib + rUaAf*(g & mesh.Sf())
);
surfaceScalarField phiDragb
(
fvc::interpolate(alpha/rhob*K*rUbA)*phia + rUbAf*(g & mesh.Sf())
);
phia = (fvc::interpolate(Ua) & mesh.Sf()) + fvc::ddtPhiCorr(rUaA, Ua, phia) + phiDraga;
phib = (fvc::interpolate(Ub) & mesh.Sf()) + fvc::ddtPhiCorr(rUbA, Ub, phib) + phiDragb;
phi = alphaf*phia + betaf*phib;
surfaceScalarField Dp
(
"(rho*(1|A(U)))",
alphaf*rUaAf/rhoa + betaf*rUbAf/rhob
);
phiDraga 和 phiDragb 分别对应 ( α b K a p , a ρ a ) f U b ⋅ S f + ( 1 a p , a ) f g ⋅ S f (\frac{\alpha_b K}{ a_{p,a} \rho_a})_f U_b \cdot S_f +(\frac{1}{a_{p,a}})_f g \cdot S_f (ap,aρaαbK)fUb⋅Sf+(ap,a1)fg⋅Sf 和 ( α a K a p , b ρ b ) f U a ⋅ S f + ( 1 a p , b ) f g ⋅ S f (\frac{\alpha_a K}{ a_{p,b} \rho_b})_f U_a \cdot S_f +(\frac{1}{a_{p,b}})_f g \cdot S_f (ap,bρbαaK)fUa⋅Sf+(ap,b1)fg⋅Sf
由于13-14行的定义，28-29行中的 (fvc::interpolate(Ua) & mesh.Sf()) 和 (fvc::interpolate(Ub) & mesh.Sf()) 便分别对应的是 ( 1 a p , a ) f H ( U a ) ⋅ S f (\frac{1}{a_{p,a}})_f H(U_a)\cdot S_f (ap,a1)fH(Ua)⋅Sf 和 ( 1 a p , b ) f H ( U b ) ⋅ S f (\frac{1}{a_{p,b}})_f H(U_b)\cdot S_f (ap,b1)fH(Ub)⋅Sf 。
有了上面的定义，可以看出31行定义的phi=alphafphia + betafphib便表示了压力方程的左边。
再看33-36行定义的Dp，很显然，表示的是压力方程右边的 ( α a a p , a ρ a + α b a p , b ρ b ) f (\frac{\alpha_a }{a_{p,a}\rho_a} + \frac{\alpha_b }{a_{p,b}\rho_b})_f (ap,aρaαa+ap,bρbαb)f。
有了以上的定义，便可以构建用于修正界面通量的压力方程了：
fvScalarMatrix pEqn
(
fvm::laplacian(Dp, p) == fvc::div(phi)
);
如上所述，pEqn收敛以后，得到的就是满足连续性的界面通量了，然后再利用求得的界面通量来修正两相的速度，便得到了满足两相连续性的速度：
Ua += fvc::reconstruct(phiDraga - rUaAf*SfGradp/rhoa);
Ua.correctBoundaryConditions();
Ub += fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob);
Ub.correctBoundaryConditions();
U = alpha*Ua + beta*Ub;
注意，Ua和Ub为什么是这样来修正呢？回想上面变量定义那个代码段的13-14行，这两行将Ua和Ub分别定义成了
1
a
p
,
a
H
(
U
a
)
\frac{1}{a_{p,a}} H(U_a)
ap,a1H(Ua) 和
1
a
p
,
b
H
(
U
b
)
\frac{1}{a_{p,b}} H(U_b)
ap,b1H(Ub)。
回想Ua和Ub的离散方程的统一形式
U
a
=
1
a
p
,
a
H
(
U
a
)
−
∇
p
a
p
,
a
ρ
a
+
α
b
a
p
,
a
ρ
a
K
U
b
+
1
a
p
,
a
g
U_{a}=\frac{1}{a_{p,a}}H(U_a)-\frac{\nabla p}{a_{p,a}\rho_a}+\frac{\alpha_b}{ a_{p,a} \rho_a} K U_b +\frac{1}{a_{p,a}} g
Ua=ap,a1H(Ua)−ap,aρa∇p+ap,aρaαbKUb+ap,a1g
U b = 1 a p , b H ( U b ) − ∇ p a p , b ρ b + α a a p , b ρ b K U a + 1 a p , b g U_{b}=\frac{1}{a_{p,b}}H(U_b)-\frac{\nabla p}{a_{p,b}\rho_b}+\frac{\alpha_a}{ a_{p,b} \rho_b} K U_a +\frac{1}{a_{p,b}} g Ub=ap,b1H(Ub)−ap,bρb∇p+ap,bρbαaKUa+ap,b1g
会发现13-14行定义的Ua和Ub都少了几项，所以缺了的这几项的贡献需要在速度修正步骤加回来，而Ua+=后面的fvc::reconstruct(phiDraga - rUaAf*SfGradp/rhoa)刚好就对应着Ua缺少的那几项。因为经过压力方程修正以后，界面通量是连续的，所以，将缺失的几项对应的界面通量通过reconstruct函数从界面通量重构从对体中心的速度的贡献，便得到了满足连续性的体中心速度了。对Ub也是同样的。
经过以上步骤，便能得到满足整体连续性的两相速度Ua 和 Ub了。
注释
【注一】：OpenFOAM 里的div函数，字面意义上看起来好像是散度的意思，实际上，div函数执行的是加和运算。举例说，对于fvc::div(phia)，phia是surfaceScalarField，其值为(fvc::interpolate(Ua) & mesh.Sf())，即将存储在体中心的Ua插值到每个网格对应的面的面心，然后用面心的速度与该面的面积矢量点乘。从代码中看，fvc::div(phia)对应的是
∇
⋅
U
a
\nabla \cdot U_a
∇⋅Ua，根据高斯定理，也就是
∑
f
(
U
a
)
f
⋅
S
f
\sum_f (U_a)_f \cdot S_f
∑f(Ua)f⋅Sf，而phia对应着
(
U
a
)
f
⋅
S
f
(U_a)_f \cdot S_f
(Ua)f⋅Sf，所以，fvc::div(phia)实际进行的运算是将包围每个网格的面上的通量加起来。更详细的说明见我的另一篇博文。
【注二】：注意这里的
∇
α
a
α
a
\frac{\nabla \alpha_a}{\alpha_a}
αa∇αa 在代码中的表示方法，详细说明见我的另一篇博文。
【注三】：注意这里的
∇
α
a
α
a
\frac{\nabla \alpha_a}{\alpha_a}
αa∇αa 在代码中的表示方法，以及与上一个
∇
α
a
α
a
\frac{\nabla \alpha_a}{\alpha_a}
αa∇αa 的区别，详细说明见我的另一篇博文。
【注四】：这里说的总体的连续性方程指的是总体体积的守恒，而不是总体质量的守恒，这二者的差异见Henrik Rusche 的 PHD 论文 P112 的说明。
参考资料
- Henrik Rusche， PHD Thesis， Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions， Imperial College of Science, Technology & Medicine, Department of Mechanical Engineering, 2002
- https://openfoamwiki.net/index.php/BubbleFoam
- http://www.cfd-online.com/Forums/openfoam-solving/71141-rewriting-twophaseeulerfoam-conservative-form.html
![[Modbus] Modbus协议开发-基本概念（一）](https://i-blog.csdnimg.cn/direct/62093e7cf3ea4e75b920bb47a2937d0e.png)
![NSSCTF练习记录：[SWPUCTF 2021 新生赛]jicao](https://i-blog.csdnimg.cn/direct/9728c9e91c384e7dbc1e453b7617a2cc.png)