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+∇⋅(αaUa)=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+∇⋅(αaU)+∇⋅(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=αaUa+αbUbUr=Ua−Ub
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)
- );
- surfaceScalarField phic("phic", phi);
- surfaceScalarField phir("phir", phia - phib);
- 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
- if (kineticTheory.on())
- {
- kineticTheory.solve(gradUaT);
- nuEffa = kineticTheory.mua()/rhoa;
- }
- else
- {
- nuEffa = sqr(Ct)*nutb + nua;
- }
- if (kineticTheory.on())
- {
- 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 的其他相干论文。
g0
跟 kineticTheory一样，g0 也是只有在模拟气固（液固）两相流时才需要开启，相干的设置在constant/ppProperties里，当g0的值设置为大于零时，则跟g0相干的项会其作用，主要的影响如下：
- pEqn
- if (g0.value() > 0.0)
- {
- phiDraga -= ppMagf*fvc::snGrad(alpha)*mesh.magSf();
- }
- ppMagf = rUaAf*fvc::interpolate
- (
- (1.0/(rhoa*(alpha + scalar(0.0001))))
- *g0*min(exp(preAlphaExp*(alpha - alphaMax)), expMax)
- );
− 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 进行修正：- 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;
- }
- 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)"
- );
- }
在对 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+∇⋅(αaU)+∇⋅(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 缺乏理论底子，仅仅是一种不甚高明的数值本领，其主要的处理是界说了一个
- volScalarField alphaEx(max(alpha - alphaMax, scalar(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企服之家，中国第一个企服评测及商务社交产业平台。 |
