2015-05-17

twoPhaseEulerFoam 全解读之二

本系列将对OpenFOAM-2.1.1 中的 twoPhaseEulerFoam 求解器进行完全解读，共分三部分：方程推导，代码解读，补充说明。本篇对 twoPhaseEulerFoam 中的 UEqn.H 和 pEqn.H 中的代码进行详细地的解读。

2. 代码解读

2.1. UEqn

前一篇导出了分散相的动量守恒方程
$$
\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}
$$
这一篇分析twoPhaseEulerFoam求解器是怎么来对动量方程进行离散的，以及，如果通过构建压力方程来对速度进行修正以保证两相的连续性。

注意上述动量方程中，有两项还需要处理一下：
$$U_a \cdot \nabla U_a=\nabla\cdot(U_aU_a)-U_a(\nabla\cdot U_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)$$
转化前后的形式从数学上来等价的，但是在有限体积离散过程中，转化前的 $U_a \cdot \nabla U_a$ 和 $\frac{\nabla(\alpha_a)}{\alpha_a}\cdot\left[-\nu_{eff} \nabla U_a\right]$ 对于$U_a$来说是非守恒的，转化后的形式是守恒的(参考这篇文献，注意这样转化后，动量方程空间上是守恒的，但时间上仍是不守恒的，这个帖子也有一些有价值的信息)。

这样转化以后，得到的动量方程就跟twoPhaseEulerFoam里的定义是一模一样了:
$$
\begin{aligned}
&(1+\frac{\alpha_b \rho_b}{\rho_a} C_{vm})\left(\frac{\partial U_a}{\partial t} + \nabla\cdot(U_aU_a)-U_a(\nabla\cdot U_a) \right)\\
-&\nabla \cdot \left[ \nu_{eff} \nabla U_a \right ] + \nabla \cdot \left[ R_{c,a}\right] +\nabla \cdot\left[ -\nu_{eff}\frac{(\nabla \alpha_a)}{\alpha_a}(\nabla U_a)\right]-U_a\left(\nabla\cdot(-\nu_{eff}\frac{\nabla \alpha_a}{\alpha_a}) \right) \\
+ & \frac{\nabla(\alpha_a)}{\alpha_a} \cdot \left[ 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}
$$

下面将动量方程的每一项与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)对应$\frac{\partial U_a}{\partial t}$，phia定义为fvc::interpolate(Ua) & mesh.Sf()，于是fvm::div(phia, Ua, "div(phia,Ua)") 和 fvm::Sp(fvc::div(phia), Ua) 便分别对应 $\nabla\cdot(U_aU_a)$ 和 $U_a(\nabla\cdot U_a)$了 [ 注一 ]。

- fvm::laplacian(nuEffa, Ua)
+ fvc::div(Rca)
+ fvm::div(phiRa, Ua, "div(phia,Ua)")
- fvm::Sp(fvc::div(phiRa), Ua)

fvm::laplacian(nuEffa, Ua)对应$\nabla \cdot \left[ \nu_{eff} \nabla U_a \right ]$，fvc::div(Rca)对应$\nabla \cdot \left[ R_{c,a}\right]$。
phiRa的定义是

1 2	-fvc::interpolate(nuEffa)mesh.magSf()fvc::snGrad(alpha) /fvc::interpolate(alpha + scalar(0.001))

相当于$-\nu_{eff}\frac{\nabla \alpha_a}{\alpha_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\left(\nabla\cdot(-\nu_{eff}\frac{\nabla \alpha_a}{\alpha_a}) \right)$。

1	(fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.001)) & Rca)

对应$\frac{\nabla(\alpha_a)}{\alpha_a} \cdot \left[ R_{c,a}\right]$，其中&运算符已重载为计算矢量与张量的点乘积 [ 注三 ]。

 ==
//  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/rhoa*K, Ua)对应 $\frac{\alpha_b}{\rho_a} K U_a$，beta/rhoa*(liftCoeff - Cvm*rhob*DDtUb) 对应
$$
\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\}
$$
其中变量liftCoeff定义为

1	volVectorField liftCoeff(Cl(betarhob + alpharhoa)(Ur ^ fvc::curl(U)));

DDtUb定义为

DDtUb =
        fvc::ddt(Ub)
      + fvc::div(phib, Ub)
      - fvc::div(phib)*Ub;

重力 $g$ 以及曳力的显式项 $\frac{\alpha_b}{\rho_a} K U_b $ 如注释所述，将会在 pEqn 中考虑，压力梯度项 $\frac{\nabla p}{\rho_a}$ 则将在 pEqn 中用来约束两相的连续性。
至此动量方程的每一项都与UEqn.H的代码对应起来了。

2.2. pEqn

压力方程的作用是修正两相速度$U_a$ 和 $U_b$以使速度满足连续性方程。将两相的连续性方程加起来，得到总体的连续性方程如下 [ 注四 ]：
$$\begin{aligned}
& \frac{\partial \alpha_a}{\partial t} + \nabla \cdot (\alpha_a U_a) + \frac{\partial \alpha_b}{\partial t} + \nabla \cdot (\alpha_b U_b) \\
= & \frac{\partial (\alpha_a+\alpha_b)}{\partial t} + \nabla \cdot (\alpha_a U_a+\alpha_b U_b) = 0
\end{aligned}$$
由于$\alpha_a+\alpha_b=1$，于是两相连续性方程等价于
$$
\nabla \cdot (\alpha_a U_a+\alpha_b U_b) = 0
$$

再来看压力方程是如何构建起来的。
完整的动量方程离散后，可以写作如下的统一形式：
$$
a_{p,a}U_{p,a}=H(U_a)-\frac{\nabla p}{\rho_a}+\frac{\alpha_b}{\rho_a} K U_b +g
$$

$$
a_{p,b}U_{p,b}=H(U_b)-\frac{\nabla p}{\rho_b}+\frac{\alpha_a}{\rho_b} K U_a +g
$$
其中$H(U_a)$ 和 $H(U_b)$ 包含了动量方程中除压力梯度项，显式曳力项以及重力项以后所有项的贡献。
由此离散方程可以得到 $U_a$ 和 $U_b$ 的表达式如下：
$$
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
$$

$$
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
$$

如果此 $U_a$ 和 $U_b$ 是方程组的解，那么它们必须满足整体的连续性方程，即
$$\begin{aligned}
& \nabla \cdot \left[ \alpha_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)\right] \\
+ & \nabla \cdot \left[ \alpha_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 ) \right] = 0
\end{aligned}$$
将压力梯度项移到方程的一边，得到
$$\begin{aligned}
& \nabla \cdot \left[ \alpha_a (\frac{1}{a_{p,a}}H(U_a)+\frac{\alpha_b}{ a_{p,a} \rho_a} K U_b +\frac{1}{a_{p,a}} g)\right] + \nabla \cdot \left[ \alpha_b (\frac{1}{a_{p,b}}H(U_b)+\frac{\alpha_a}{ a_{p,b} \rho_b} K U_a +\frac{1}{a_{p,b}} g ) \right] \\
= &\nabla \cdot \left[ (\frac{\alpha_a }{a_{p,a}\rho_a} + \frac{\alpha_b }{a_{p,b}\rho_b}) \nabla p \right ]
\end{aligned}$$
这便是压力修正方程的原型。
在pEqn.H中，压力方程其实修正的是界面通量，压力方程迭代收敛以后能保证界面通量的连续性。所以，散度表达式需要根据高斯定理写成界面通量之和的形式：
$$\begin{aligned}
& \nabla \cdot \left[ \alpha_a (\frac{1}{a_{p,a}}H(U_a)+\frac{\alpha_b}{ a_{p,a} \rho_a} K U_b +\frac{1}{a_{p,a}} g)\right] \\
= & (\alpha_a)_f \left[\sum_f(\frac{1}{a_{p,a}})_f H(U_a)\cdot S_f + \sum_f(\frac{\alpha_b K}{ a_{p,a} \rho_a})_f U_b \cdot S_f +\sum_f(\frac{1}{a_{p,a}})_f g \cdot S_f \right ]
\end{aligned}$$
下标 $_f$ 表示该项将要在代码中用界面上的变量来表示，在OpenFOAM中，即surfaceScalarField，$S_f$ 表示界面的面积矢量，下面的公式里也是一样。
$$\begin{aligned}
& \nabla \cdot \left[ \alpha_b (\frac{1}{a_{p,b}}H(U_b)+\frac{\alpha_a}{ a_{p,b} \rho_b} K U_a +\frac{1}{a_{p,b}} g)\right] \\
= & (\alpha_b)_f \left[\sum_f(\frac{1}{a_{p,b}})_f H(U_b)\cdot S_f + \sum_f(\frac{\alpha_a K}{ a_{p,b} \rho_b})_f U_a \cdot S_f +\sum_f(\frac{1}{a_{p,b}})_f g \cdot S_f \right]
\end{aligned}$$

以及
$$
\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
$$
下面是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 分别对应 $(\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$ 和 $(\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$

由于13-14行的定义，28-29行中的 (fvc::interpolate(Ua) & mesh.Sf()) 和 (fvc::interpolate(Ub) & mesh.Sf()) 便分别对应的是 $(\frac{1}{a_{p,a}})_f H(U_a)\cdot S_f$ 和 $(\frac{1}{a_{p,b}})_f H(U_b)\cdot S_f$ 。

有了上面的定义，可以看出31行定义的phi=alphaf*phia + betaf*phib便表示了压力方程的左边。

再看33-36行定义的Dp，很显然，表示的是压力方程右边的$(\frac{\alpha_a }{a_{p,a}\rho_a} + \frac{\alpha_b }{a_{p,b}\rho_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分别定义成了$\frac{1}{a_{p,a}} H(U_a)$ 和 $\frac{1}{a_{p,b}} H(U_b)$。
回想Ua和Ub的离散方程的统一形式
$$
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
$$

$$
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
$$

会发现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)对应的是 $\nabla \cdot U_a$，根据高斯定理，也就是 $\sum_f (U_a)_f \cdot S_f$，而phia对应着 $(U_a)_f \cdot S_f$，所以，fvc::div(phia)实际进行的运算是将包围每个网格的面上的通量加起来。更详细的说明见我的另一篇博文。
注二：注意这里的$\frac{\nabla \alpha_a}{\alpha_a}$ 在代码中的表示方法，详细说明见我的另一篇博文。
注三：注意这里的$\frac{\nabla \alpha_a}{\alpha_a}$ 在代码中的表示方法，以及与上一个$\frac{\nabla \alpha_a}{\alpha_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

Giskard's CFD Learning Tricks

CFD and Scientific Computing

twoPhaseEulerFoam 全解读之二

2. 代码解读

2.1. UEqn

2.2. pEqn

注释

参考资料