本系列将对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 提出的形式。主要的代码如下:1
2
3
4
5
6fvScalarMatrix alphaEqn
(
fvm::ddt(alpha)
+ fvm::div(phic, alpha, scheme)
+ fvm::div(-fvc::flux(-phir, beta, schemer), alpha, schemer)
);
其中phic
与phir
的定义如下:1
2surfaceScalarField phic("phic", phi);
surfaceScalarField phir("phir", phia - phib);
得到分散相体积分率后,连续相体积分率$\alpha_b$则为$1-\alpha_a$,如下1
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
1
2
3
4
5if (kineticTheory.on())
{
kineticTheory.solve(gradUaT);
nuEffa = kineticTheory.mua()/rhoa;
}如果开启
kineticTheory
,则固相粘度nuEffa
是由 KTGF 来计算得到,否则nuEffa
是用一个简单的关联式来计算1
2
3
4else
{
nuEffa = sqr(Ct)*nutb + nua;
}此外,
Rca
项也要加上额外的项1
2
3
4if (kineticTheory.on())
{
Rca -= ((kineticTheory.lambda()/rhoa)*tr(gradUaT))*tensor(I);
}pEqn.H
如果kineticTheory
开启了,则要在压力修正方程中加上额外的固相压力项。1
2
3
4if (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
1
2
3
4if (g0.value() > 0.0)
{
phiDraga -= ppMagf*fvc::snGrad(alpha)*mesh.magSf();
}其中
ppMagf
的定义如下:1
2
3
4
5ppMagf = 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
进行修正:1
2
3
4
5
6
7if (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
进行修正:1
2
3
4
5
6
7
8
9
10
11
12
13
14
15if (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 缺乏理论基础,仅仅是一种不甚高明的数值技巧,其主要的处理是定义了一个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