twoPhaseEulerFoam 全解读之三（转)

twoPhaseEulerFoam 全解读之三（转)

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

前文提到，经过求解 pEqn，并修正速度Ua和Ub以后，总体的一连性便得到了保证。为了得到分散相的体积分率alpha，还需要使用得到的速度场来求解分散相的一连性方程，即alphaEqn。分散相一连性方程可以表达如下：
                                                                ∂                                               α                                     a                                                                  ∂                                  t                                                 +                            ∇                            ⋅                            (                                       α                               a                                                 U                               a                                      )                            =                            0                                  \frac{\partial \alpha_a}{\partial t}+\nabla \cdot (\alpha_a U_a)=0                     ∂t∂αa​​+∇⋅(αa​Ua​)=0
为了让每一项都写成守恒情势，而且保证                                             α                            a                                       \alpha_a                  αa​的有界性，Weller 将分散相一连性方程写成如下情势
                                                                ∂                                               α                                     a                                                                  ∂                                  t                                                 +                            ∇                            ⋅                            (                                       α                               a                                      U                            )                            +                            ∇                            ⋅                            (                                       U                               r                                                 α                               a                                      (                            1                            −                                       α                               a                                      )                            )                            =                            0                                  \frac{\partial \alpha_a}{\partial t}+\nabla \cdot (\alpha_a U)+\nabla \cdot(U_r\alpha_a(1-\alpha_a))=0                     ∂t∂αa​​+∇⋅(αa​U)+∇⋅(Ur​αa​(1−αa​))=0
其中
                                         U                            =                                       α                               a                                                 U                               a                                      +                                       α                               b                                                 U                               b                                                          U                               r                                      =                                       U                               a                                      −                                       U                               b                                            U=\alpha_a U_a + \alpha_b U_b \\ U_r=U_a-U_b                     U=αa​Ua​+αb​Ub​Ur​=Ua​−Ub​
OpenFOAM-2.1.1 的alphaEqn求解的正是 Weller 提出的情势。主要的代码如下：

1. fvScalarMatrix alphaEqn
2.        (
3.             fvm::ddt(alpha)
4.           + fvm::div(phic, alpha, scheme)
5.           + fvm::div(-fvc::flux(-phir, beta, schemer), alpha, schemer)
6.        );

复制代码

其中phic与phir的界说如下：

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

复制代码

得到分散相体积分率后，一连相体积分率                                             α                            b                                       \alpha_b                  αb​则为                                   1                         −                                   α                            a                                       1-\alpha_a                  1−αa​，如下

1. beta = scalar(1) - alpha;

复制代码

补充说明

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

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

• UEqn.h
  
  1.   if (kineticTheory.on())
  2.   {
  3.       kineticTheory.solve(gradUaT);
  4.       nuEffa = kineticTheory.mua()/rhoa;
  5.   }

  复制代码
  假如开启kineticTheory，则固相粘度nuEffa是由 KTGF 来计算得到，否则nuEffa是用一个简朴的关联式来计算
  
  1. else
  2.   {
  3.       nuEffa = sqr(Ct)*nutb + nua;
  4.   }

  复制代码
  别的，Rca项也要加上额外的项
  
  1. if (kineticTheory.on())
  2.   {
  3.       Rca -= ((kineticTheory.lambda()/rhoa)*tr(gradUaT))*tensor(I);
  4.   }

  复制代码
• pEqn.h
  假如kineticTheory开启了，则要在压力修正方程中加上额外的固相压力项。
  
  1. if (kineticTheory.on())
  2. {
  3.     phiDraga -= rUaAf*fvc::snGrad(kineticTheory.pa()/rhoa)*mesh.magSf();
  4. }

  复制代码

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

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

• pEqn
  1. if (g0.value() > 0.0)
  2.   {
  3.       phiDraga -= ppMagf*fvc::snGrad(alpha)*mesh.magSf();
  4.   }

  复制代码
  其中ppMagf的界说如下：
  1. ppMagf = rUaAf*fvc::interpolate
  2.   (
  3.    (1.0/(rhoa*(alpha + scalar(0.0001))))
  4.    *g0*min(exp(preAlphaExp*(alpha - alphaMax)), expMax)
  5.   );

  复制代码
  这一段的结果，相当于在本系列第一篇最后的分散相动量方程的等式右边额外增加一项：
                                                −                               g                               0                               ∗                               m                               i                               n                               (                                           e                                               p                                     r                                     e                                     A                                     l                                     p                                     h                                     a                                     E                                     x                                     p                                     ∗                                     (                                                   α                                        a                                                  −                                                   α                                                       M                                           a                                           x                                                                )                                                      ,                               e                               x                               p                               M                               a                               x                               )                               ∇                                           α                                  a                                          /                               (                                           α                                  a                                                      ρ                                  a                                          )                                      -g0*min(e^{preAlphaExp*(\alpha_a - \alpha_{Max})},expMax)\nabla \alpha_a / (\alpha_a \rho_a)                        −g0∗min(epreAlphaExp∗(αa​−αMax​),expMax)∇αa​/(αa​ρa​)
  可见，这一项的作用是给固相额外施加了一个力，可以认为是固相压力，这个力与                                                   α                               a                                            \alpha_a                     αa​的梯度有关，且与梯度的方向相反。
  留意这一项的特点：g0, preAlphaExp, alphaMax, expMax 都是需要用户指定的参数，一般将g0, preAlphaExp和expMax设置为一个比力大的正整数（                                        1                                       0                               3                                            10^3                     103的量级），当alpha - alphaMax < 0时，                                                   e                                                       p                                  r                                  e                                  A                                  l                                  p                                  h                                  a                                  E                                  x                                  p                                  ∗                                  (                                               α                                     a                                              −                                               α                                                   M                                        a                                        x                                                           )                                                       e^{\,preAlphaExp*(\alpha_a - \alpha_{Max})}                     epreAlphaExp∗(αa​−α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 进行修正：
  1. if (g0.value() > 0.0)
  2. {
  3.     surfaceScalarField alphaf(fvc::interpolate(alpha));
  4.     surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha)*mesh.magSf());
  5.     phir += phipp;
  6.     phic += fvc::interpolate(alpha)*phipp;
  7. }

  复制代码
  另一处是对 alphaEqn 进行修正：
  1. if (g0.value() > 0.0)
  2. {
  3.     ppMagf = rUaAf*fvc::interpolate
  4.     (
  5.      (1.0/(rhoa*(alpha + scalar(0.0001))))
  6.      *g0*min(exp(preAlphaExp*(alpha - alphaMax)), expMax)
  7.      );
  8.    
  9.     alphaEqn -= fvm::laplacian
  10.     (
  11.      (fvc::interpolate(alpha) + scalar(0.0001))*ppMagf,
  12.      alpha,
  13.      "laplacian(alphaPpMag,alpha)"
  14.      );
  15. }

  复制代码
  为什么当 g0 > 0 时，需要在 alphaEqn 中额外减去一个 laplacian 项呢？这里要联合上面提到的两段代码来进行分析。
  在对 phic 和 phir 进行修正这一段， phic 和 phir 分别加上了一项。将修正过的 phic 和 phir 代入到 alphaEqn 中，会导致 alphaEqn 与上文的分散相一连性方程
                                                                         ∂                                                   α                                        a                                                                        ∂                                     t                                                      +                               ∇                               ⋅                               (                                           α                                  a                                          U                               )                               +                               ∇                               ⋅                               (                                           U                                  r                                                      α                                  a                                          (                               1                               −                                           α                                  a                                          )                               )                               =                               0                                      \frac{\partial \alpha_a}{\partial t}+\nabla \cdot (\alpha_a U)+\nabla \cdot(U_r\alpha_a(1-\alpha_a))=0                        ∂t∂αa​​+∇⋅(αa​U)+∇⋅(Ur​αa​(1−α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] $ 。
多出的这两项加起来，为
                                         ∇                            ⋅                            [                            α                            ∗                                       p                               p                               M                               a                               g                               f                                                                   ∇                            α                            ]                                  \nabla \cdot [\alpha *\mathrm{ppMagf} \ \nabla \alpha]                     ∇⋅[α∗ppMagf ∇α]
对比上面的 alphaEqn 减去的 laplacian 项，会发现这一项跟谁人 laplacian 是一样的。
以是，g0 > 0 时，先将固相压力带来的通量代入到 alphaEqn 的构建中，然后再减去对应的 laplacian 项，这么做的目的，应该是为了计算稳固性以及保证 alpha 的有界性（保证                                    0                         <                         a                         l                         p                         h                         a                         <                         a                         l                         p                         h                         a                         M                         a                         x                              0<alpha<alphaMax                  0<alpha<alphaMax）。但是这背后的数值原理，我也没有完全明白。
最后，有必要说明一下，g0相干的项着实就是在动量方程中增加了一项颗粒压力，想要达到的结果是防止固相过度堆积。kineticTheory开启后将计算颗粒相的应力，其中包括了颗粒相压力。以是 g0可以当作 KTGF 的某种简朴的替换。从原理上讲，kineticTheory和g0不应该同时开启。
packingLimiter

packingLimiter 的作用，从名字可以看出，也是用来处理过度堆积问题的。packingLimiter 缺乏理论底子，仅仅是一种不甚高明的数值本领，其主要的处理是界说了一个

1. 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 方法来求解。这些有待于下一步进行解读。

---

参考资料

• 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
  

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。