twoPhaseEulerFoam 全解读之三

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

2.3. alphaEqn

前文提到，经过求解 pEqn，并修正速度Ua和Ub以后，总体的连续性便得到了保证。为了得到分散相的体积分率alpha，还需要利用得到的速度场来求解分散相的连续性方程，即alphaEqn。分散相连续性方程可以表达如下：
$$
\frac{\partial \alpha_a}{\partial t}+\nabla \cdot (\alpha_a U_a)=0
$$

为了让每一项都写成守恒形式，并且保证$\alpha_a$的有界性，Weller 将分散相连续性方程写成如下形式
$$
\frac{\partial \alpha_a}{\partial t}+\nabla \cdot (\alpha_a U)+\nabla \cdot(U_r\alpha_a(1-\alpha_a))=0
$$
其中
$$
U=\alpha_a U_a + \alpha_b U_b \\
U_r=U_a-U_b
$$
OpenFOAM-2.1.1 的alphaEqn求解的正是 Weller 提出的形式。主要的代码如下：

fvScalarMatrix alphaEqn
       (            fvm::ddt(alpha)          + fvm::div(phic, alpha, scheme)          + fvm::div(-fvc::flux(-phir, beta, schemer), alpha, schemer)       );
            fvm::ddt(alpha)
          + fvm::div(phic, alpha, scheme)
          + fvm::div(-fvc::flux(-phir, beta, schemer), alpha, schemer)
       )

其中phic与phir的定义如下：

surfaceScalarField phic("phic", phi);
surfaceScalarField phir("phir", phia - phib);

得到分散相体积分率后，连续相体积分率$\alpha_b$则为$1-\alpha_a$，如下

beta = scalar(1) - alpha;

3. 补充说明

想必读者肯定发现了，我前面的代码解读相比于twoPhaseEulerFoam的源码其实省略了很多，总体上来讲，省略了三大块，一是跟kineticTheory相关的，二是跟g0相关的，三是packingLimiter，下面对这三部分进行一些补充说明。

3.1 kineticTheory

kineticTheory是 KTGF(Kinetic Theory of Granular Flow) 方法在OpenFOAM-2.1.1里的实现。KTGF 的主要作用是计算固相压力和固相粘度，以封闭前述的分散相动量方程，所以kineticTheory只有在模拟气固（液固）两相流时才需要开启。 kineticTheory 是否开启以及 KTGF 模型的参数需要在算例的constant/kineticTheoryProperties文件里进行设置。如果开启了kineticTheory，主要的影响如下：

• UEqn.H

  if (kineticTheory.on())
  {
      kineticTheory.solve(gradUaT);
      nuEffa = kineticTheory.mua()/rhoa;
  }

  如果开启kineticTheory，则固相粘度nuEffa是由 KTGF 来计算得到，否则nuEffa是用一个简单的关联式来计算

  else
  {
      nuEffa = sqr(Ct)*nutb + nua;
  }

  此外，Rca项也要加上额外的项

  if (kineticTheory.on())
  {    Rca -= ((kineticTheory.lambda()/rhoa)*tr(gradUaT))*tensor(I);}
      Rca -= ((kineticTheory.lambda()/rhoa)*tr(gradUaT))*tensor(I);
  }

• pEqn.H
  如果kineticTheory开启了，则要在压力修正方程中加上额外的固相压力项。

  if (kineticTheory.on())
  {
      phiDraga -= rUaAf*fvc::snGrad(kineticTheory.pa()/rhoa)*mesh.magSf();
  }

有关 OpenFOAM 中使用的 KTGF 模型的理论可以参考Derivation, Implementation, and Validation of Computer Simulation Models for Gas-Solid Fluidized Beds，以及B.G.M. van Wachem 的其他相关论文。

3.2 g0

跟 kineticTheory一样，g0 也是只有在模拟气固（液固）两相流时才需要开启，相关的设置在constant/ppProperties里，当g0的值设置为大于零时，则跟g0相关的项会其作用，主要的影响如下：

• pEqn

  if (g0.value() > 0.0)
  {
      phiDraga -= ppMagf*fvc::snGrad(alpha)*mesh.magSf();
  }

  其中ppMagf的定义如下：

  ppMagf = rUaAf*fvc::interpolate
  (
   (1.0/(rhoa*(alpha + scalar(0.0001))))
   *g0*min(exp(preAlphaExp*(alpha - alphaMax)), expMax)
  );

  这一段的效果，相当于在本系列第一篇最后的分散相动量方程的等式右边额外增加一项：
  $$
  -g0*min(e^{preAlphaExp*(\alpha_a - \alpha_{Max})},expMax)\nabla \alpha_a / (\alpha_a \rho_a)
  $$
  可见，这一项的作用是给固相额外施加了一个力，可以认为是固相压力，这个力与$\alpha_a$的梯度有关，且与梯度的方向相反。
  注意这一项的特点：g0, preAlphaExp, alphaMax, expMax 都是需要用户指定的参数，一般将g0, preAlphaExp和expMax设置为一个比较大的正整数（$10^3$的量级），当alpha - alphaMax < 0时，$e^{\,preAlphaExp*(\alpha_a - \alpha_{Max})}$ 将是一个很小的数，此时整个增加的一项也将是一个较小的数，所以，alpha - alphaMax < 0时额外增加的那一项对动量方程的贡献很小。但是，当alpha - alphaMax > 0时，随着 alpha 偏离 alphaMax 越来越远，exp(preAlphaExp*(alpha - alphaMax)将迅速增大，min(exp(preAlphaExp*(alpha - alphaMax)), expMax)的值也将迅速增大，直到设定值expMax。可见，g0项的作用可以理解为为了防止固相过度堆积。气固系统中，固相的体积分率是有上限的，可以将alphaMax设置为这个上限值，当某个网格里的alpha超过设定的alphaMax时，就让固相迅速弹开，以防止固相体积分率过大。

• alphaEqn.H
  alphaEqn.H 里有两处跟 g0 有关的，
  一处是对 phic 和 phir 进行修正：

  if (g0.value() > 0.0)
  {
      surfaceScalarField alphaf(fvc::interpolate(alpha));
      surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha)*mesh.magSf());
      phir += phipp;
      phic += fvc::interpolate(alpha)*phipp;
  }

另一处是对 alphaEqn 进行修正：

if (g0.value() > 0.0)
{
    ppMagf = rUaAf*fvc::interpolate
    (
     (1.0/(rhoa*(alpha + scalar(0.0001))))
     *g0*min(exp(preAlphaExp*(alpha - alphaMax)), expMax)
     );
    
    alphaEqn -= fvm::laplacian
    (
     (fvc::interpolate(alpha) + scalar(0.0001))*ppMagf,
     alpha,
     "laplacian(alphaPpMag,alpha)"
     );
}

为什么当 g0 > 0 时，需要在 alphaEqn 中额外减去一个 laplacian 项呢？这里要结合上面提到的两段代码来进行分析。
在对 phic 和 phir 进行修正这一段， phic 和 phir 分别加上了一项。将修正过的 phic 和 phir 代入到 alphaEqn 中，会导致 alphaEqn 与上文的分散相连续性方程
$$
\frac{\partial \alpha_a}{\partial t}+\nabla \cdot (\alpha_a U)+\nabla \cdot(U_r\alpha_a(1-\alpha_a))=0
$$
不一致，其中 alphaEqn 多出了几项：
其中
fvm::div(phic, alpha, scheme) 比 $ \nabla \cdot (\alpha_a U) $ 多了 fvm::div(alphaf*phipp, alpha, scheme) 一项，写成公式，就是 $ \nabla \cdot [\alpha (\alpha *\mathrm{ppMagf}) \nabla \alpha] $ 。

而 fvm::div(-fvc::flux(-phir, beta, schemer), alpha, schemer) 则比 $ \nabla \cdot (\alpha_a U)+\nabla \cdot(U_r\alpha_a(1-\alpha_a) $ 多出了 fvm::div(-fvc::flux(-phipp, beta, schemer), alpha, schemer) 一项，写成公式就是 $ \nabla \cdot [\alpha (\beta *\mathrm{ppMagf}) \nabla \alpha] $ 。
多出的这两项加起来，为
$$
\nabla \cdot [\alpha *\mathrm{ppMagf} \ \nabla \alpha]
$$

对比上面的 alphaEqn 减去的 laplacian 项，会发现这一项跟那个 laplacian 是一样的。

所以，g0 > 0 时，先将固相压力带来的通量代入到 alphaEqn 的构建中，然后再减去对应的 laplacian 项，这么做的目的，应该是为了计算稳定性以及保证 alpha 的有界性（保证 $0<alpha<alphaMax$）。但是这背后的数值原理，我也没有完全理解。

最后，有必要说明一下，g0相关的项其实就是在动量方程中增加了一项颗粒压力，想要达到的效果是防止固相过度堆积。kineticTheory开启后将计算颗粒相的应力，其中包括了颗粒相压力。所以 g0可以当作 KTGF 的某种简单的替代。从原理上讲，kineticTheory和g0不应该同时开启。

3.3 packingLimiter

packingLimiter 的作用，从名字可以看出，也是用来处理过度堆积问题的。packingLimiter 缺乏理论基础，仅仅是一种不甚高明的数值技巧，其主要的处理是定义了一个

volScalarField alphaEx(max(alpha - alphaMax, scalar(0)));

当alpha - alphaMax > 0时，alphaEx = alpha - alphaMax > 0，否则alphaEx = 0。
packingLimiter 本质操作是，当某个网格的固相体积分率超过设定的最大值时，就让该网格的固相往它的邻居网格匀一匀，就是这么简(ren)单(xing)。所以，一般情况下，不建议开启packingLimiter。

至此，OpenFOAM-2.1.1 版本的twoPhaseEulerFoam便解读完了。最近的 OpenFOAM-2.3 系列中的twoPhaseEulerFoam 变化很大，求解的是全守恒形式的动量方程了（而不是 phase-intensive形式的），耦合了传热模型，考虑了可压缩效应，而且alphaEqn的求解不再是利用隐式迭代的方式，而是改成了用 MULES 方法来求解。这些有待于下一步进行解读。

参考资料

1. 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
2. https://openfoamwiki.net/index.php/BubbleFoam