OpenFOAM: twoPhaseEulerFoam解读

twoPhaseEulerFoam全解读之一(转载)

本系列将对OpenFOAM-2.1.1 中的 twoPhaseEulerFoam 求解器进行完全解读，共分三部分：方程推导，代码解读，补充说明。本篇进行方程推导，详细介绍如果从双流体模型出发得到 twoPhaseEulerFoam 中的 UEqn.H 对应的模型方程形式。

方程推导

双流体模型方程可以表达成如下形式：

连续性方程：
$$\frac{\partial ({\alpha }_{\varphi }{\rho }_{\varphi })}{\partial t}+\nabla \cdot ({\alpha }_{\varphi }{\rho }_{\varphi }{U}_{\varphi })=0$$
动量守恒方程:
$$\frac{\partial ({\alpha }_{\varphi }{\rho }_{\varphi }{U}_{\varphi })}{\partial t}+\nabla \cdot ({\alpha }_{\varphi }{\rho }_{\varphi }{U}_{\varphi }{U}_{\varphi })+\nabla \cdot ({\alpha }_{\varphi }{\tau }_{\varphi })+\nabla \cdot ({\alpha }_{\varphi }{\rho }_{\varphi }{R}_{\varphi })=-{\alpha }_{\varphi }\nabla p+{\alpha }_{\varphi }{\rho }_{\varphi }g+M_{\varphi }$$
式中，下标$\varphi =a,b$分别代表分散相和连续相，${\tau }_{\varphi }$表示粘性应力项，$R_{\varphi }$表示雷诺应力项，$M_{\varphi }$表示相间作用项。
上述方程是完全守恒形式的，但是注意到上述动量方程的瞬变项是$\frac{\partial ({\alpha }_{\varphi }{\rho }_{\varphi }{U}_{\varphi })}{\partial t}$，等于说解这个方程能得到的是每个时间步的动量，若要转化成速度，则需要用动量除以密度与体积分率的乘积，即$\frac{({\alpha }_{\varphi }{\rho }_{\varphi }{U}_{\varphi })}{{\alpha }_{\varphi }{\rho }_{\varphi }}$。那么当离散相a的体积分率${\alpha }_a\to 0$时，这个除法就要出问题了。于是，Weller [1] 提出通过构造一种”phase-intensive”形式的动量方程来避开这个问题，见下面的详细推导。
Weller提出的方法的核心是将${\alpha }_{\varphi }{\rho }_{\varphi }$从动量方程的瞬变项中剥离出来，以使动量方程直接对速度进行演化，而不是动量。
首先对动量方程的瞬变项和对流项进行如下转化：
$$\frac{\partial ({\alpha }_{\varphi }{\rho }_{\varphi }{U}_{\varphi })}{\partial t}={\alpha }_{\varphi }{\rho }_{\varphi }\frac{\partial (U_{\varphi })}{\partial t}+U_{\varphi }\frac{\partial ({\alpha }_{\varphi }{\rho }_{\varphi })}{\partial t}$$

$$\nabla \cdot ({\alpha }_{\varphi }{\rho }_{\varphi }{U}_{\varphi }{U}_{\varphi })={\alpha }_{\varphi }{\rho }_{\varphi }{U}_{\varphi }\cdot \nabla (U_{\varphi })+U_{\varphi }\nabla \cdot ({\alpha }_{\varphi }{\rho }_{\varphi }{U}_{\varphi })$$
注：这里到了张量运算公式$[\nabla \cdot \mathbf{vw}]=[\mathbf{v}\cdot \nabla \mathbf{w}]+\mathbf{w}(\nabla \cdot \mathbf{v})$，具体可参考 Bird 的 Transport Phenomenon 的 Appendix A。

于是，瞬变项和对流项的加和可以写成如下形式：
$$\frac{\partial ({\alpha }_{\varphi }{\rho }_{\varphi }{U}_{\varphi })}{\partial t}+\nabla \cdot ({\alpha }_{\varphi }{\rho }_{\varphi }{U}_{\varphi }{U}_{\varphi })={\alpha }_{\varphi }{\rho }_{\varphi }\left[\frac{\partial (U_{\varphi })}{\partial t}+U_{\varphi }\cdot \nabla (U_{\varphi })\right]+U_{\varphi }\left[\frac{\partial ({\alpha }_{\varphi }{\rho }_{\varphi })}{\partial t}+\nabla \cdot ({\alpha }_{\varphi }{\rho }_{\varphi }{U}_{\varphi })\right]$$
注意右边第二项的括号里其实就是连续性方程的左边，其值为0，因此得到：
$$\frac{\partial ({\alpha }_{\varphi }{\rho }_{\varphi }{U}_{\varphi })}{\partial t}+\nabla \cdot ({\alpha }_{\varphi }{\rho }_{\varphi }{U}_{\varphi }{U}_{\varphi })={\alpha }_{\varphi }{\rho }_{\varphi }\left[\frac{\partial (U_{\varphi })}{\partial t}+U_{\varphi }\cdot \nabla (U_{\varphi })\right]$$

于是得到第一步转化之后的动量方程：
$${\alpha }_{\varphi }{\rho }_{\varphi }\left[\frac{\partial (U_{\varphi })}{\partial t}+U_{\varphi }\cdot \nabla (U_{\varphi })\right]+\nabla \cdot ({\alpha }_{\varphi }{\tau }_{\varphi })+\nabla \cdot ({\alpha }_{\varphi }{\rho }_{\varphi }{R}_{\varphi })=-{\alpha }_{\varphi }\nabla p+{\alpha }_{\varphi }{\rho }_{\varphi }g+M_{\varphi }$$

下面处理粘性应力项和雷诺应力项。

$$\nabla \cdot ({\alpha }_{\varphi }{\tau }_{\varphi })+\nabla \cdot ({\alpha }_{\varphi }{\rho }_{\varphi }{R}_{\varphi })=\nabla \cdot \left[{\alpha }_{\varphi }{\rho }_{\varphi }(\frac{{\tau }_{\varphi }}{{\rho }_{\varphi }}+R_{\varphi })\right]=\nabla \cdot \left[{\alpha }_{\varphi }{\rho }_{\varphi }{R}_{eff,\varphi }\right]$$

其中 $R_{eff,\varphi }=\frac{{\tau }_{\varphi }}{{\rho }_{\varphi }}+R_{\varphi }$。

根据定义(此处参考BubbleFoam的Wiki页面)：
$${\mathbf{\tau }}_{\varphi }=-{\rho }_{\varphi }{\nu }_{\varphi }\left[\nabla {\mathbf{U}}_{\varphi }+{\nabla }^{\text{T}}{\mathbf{U}}_{\varphi }\right]+\frac{2}{3}{\rho }_{\varphi }{\nu }_{\varphi }\left(\nabla \cdot {\mathbf{U}}_{\varphi }\right)\mathbf{I}$$
以及
$${\mathbf{R}}_{\varphi }=-{\nu }_{\varphi ,\text{t}}\left[\nabla {\mathbf{U}}_{\varphi }+{\nabla }^{\text{T}}{\mathbf{U}}_{\varphi }\right]+\frac{2}{3}{\nu }_{\varphi ,\text{t}}\left(\nabla \cdot {\mathbf{U}}_{\varphi }\right)\mathbf{I}+\frac{2}{3}{k}_{\varphi }\mathbf{I}$$
代入到 $R_{eff,\varphi }$中，得：
$$R_{eff,\varphi }=-({\nu }_{\varphi }+{\nu }_{\varphi ,t})\left[\nabla {\mathbf{U}}_{\varphi }+{\nabla }^{\text{T}}{\mathbf{U}}_{\varphi }\right]+\frac{2}{3}({\nu }_{\varphi }+{\nu }_{\varphi ,t})\left(\nabla \cdot {\mathbf{U}}_{\varphi }\right)\mathbf{I}+\frac{2}{3}{k}_{\varphi }\mathbf{I}$$
令 ${\nu }_{eff}={\nu }_{\varphi }+{\nu }_{\varphi ,t}$ ，则：
$$R_{eff,\varphi }=-{\nu }_{eff}\left[\nabla {\mathbf{U}}_{\varphi }+{\nabla }^{\text{T}}{\mathbf{U}}_{\varphi }\right]+\frac{2}{3}{\nu }_{eff}\left(\nabla \cdot {\mathbf{U}}_{\varphi }\right)\mathbf{I}+\frac{2}{3}{k}_{\varphi }\mathbf{I}=-{\nu }_{eff}\nabla U_{\varphi }+R_{c,\varphi }$$
其中$$R_{c,\varphi }=-{\nu }_{eff}\nabla {\mathbf{U}}_{\varphi }^{\text{T}}+\frac{2}{3}{\nu }_{eff}\left(\nabla \cdot {\mathbf{U}}_{\varphi }\right)\mathbf{I}+\frac{2}{3}{k}_{\varphi }\mathbf{I}$$

于是得到：
$$\begin{array}{rl}\nabla \cdot \left[{\alpha }_{\varphi }{\rho }_{\varphi }{R}_{eff,\varphi }\right]=& \nabla ({\alpha }_{\varphi }{\rho }_{\varphi })\cdot \left[R_{eff,\varphi }\right]+{\alpha }_{\varphi }{\rho }_{\varphi }\nabla \cdot \left[R_{eff,\varphi }\right]\\ =& {\alpha }_{\varphi }{\rho }_{\varphi }\nabla \cdot \left[-{\nu }_{eff}\nabla U_{\varphi }\right]+{\alpha }_{\varphi }{\rho }_{\varphi }\nabla \cdot \left[R_{c,\varphi }\right]+\nabla ({\alpha }_{\varphi }{\rho }_{\varphi })\left[-{\nu }_{eff}\nabla U_{\varphi }+R_{c,\varphi }\right]\end{array}$$

代入到动量方程中，并且方程两边同时除以${\alpha }_{\varphi }{\rho }_{\varphi }$，得到：
$$\frac{\partial U_{\varphi }}{\partial t}+U_{\varphi }\cdot \nabla U_{\varphi }-\nabla \cdot \left[{\nu }_{eff}\nabla U_{\varphi }\right]+\nabla \cdot \left[R_{c,\varphi }\right]+\frac{\nabla ({\alpha }_{\varphi }{\rho }_{\varphi })}{{\alpha }_{\varphi }{\rho }_{\varphi }}\cdot \left[-{\nu }_{eff}\nabla U_{\varphi }+R_{c,\varphi }\right]=-\frac{\nabla p}{{\rho }_{\varphi }}+g+\frac{M_{\varphi }}{{\alpha }_{\varphi }{\rho }_{\varphi }}$$

如果假定两相流体均为不可压缩，密度恒为常数，于是可以得到不可压缩的双流体模型的方程组：

连续性方程
$$\frac{\partial ({\alpha }_{\varphi })}{\partial t}+\nabla \cdot ({\alpha }_{\varphi }{U}_{\varphi })=0$$

动量方程
$$\frac{\partial U_{\varphi }}{\partial t}+U_{\varphi }\cdot \nabla U_{\varphi }-\nabla \cdot \left[{\nu }_{eff}\nabla U_{\varphi }\right]+\nabla \cdot \left[R_{c,\varphi }\right]+\frac{\nabla ({\alpha }_{\varphi })}{{\alpha }_{\varphi }}\cdot \left[-{\nu }_{eff}\nabla U_{\varphi }+R_{c,\varphi }\right]=-\frac{\nabla p}{{\rho }_{\varphi }}+g+\frac{M_{\varphi }}{{\alpha }_{\varphi }{\rho }_{\varphi }}$$

方程中还剩下相间作用项没有处理，对于分散相和连续项形式，相间作用力是大小相等符号想反，这里只考虑分散相的形式，令$\varphi =a$，则得到分散相的动量方程：
$$\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{\nabla p}{{\rho }_a}+g+\frac{M_a}{{\alpha }_a{\rho }_a}$$

相间作用只考虑曳力，升力以及虚拟质量力，即$M,a=M_{drag}+M_{lift}+M_{vm}$，下面分别考虑每一种相间作用力。

曳力
$M_{drag}=-\beta (U_a-U_b)$，其中$\beta$为曳力系数。
升力
$M_{lift}=-{\alpha }_a{\alpha }_b{C}_l({\alpha }_b{\rho }_b+{\alpha }_a{\rho }_a)U_r\times (\nabla \times U)$ ，其中 $U_r=U_a-U_b$，$U={\alpha }_a{U}_a+{\alpha }_b{U}_b$
虚拟质量力
$M_{vm}={\alpha }_a{\alpha }_b{C}_{vm}{\rho }_b\left[\frac{D{U}_b}{Dt}-\frac{D{U}_a}{Dt}\right]$，其中$\frac{D}{Dt}$表示物质导数，$\frac{D{U}_b}{Dt}=\frac{\partial U_b}{\partial t}+U_b\cdot \nabla U_b$，$\frac{D{U}_a}{Dt}=\frac{\partial U_a}{\partial t}+U_a\cdot \nabla U_a$
考虑到形式的统一，令$K=\frac{\beta }{{\alpha }_a{\alpha }_b}$，则曳力可表示为$M_{drag}=-{\alpha }_a{\alpha }_{b}K(U_a-U_b)$

代入到分散相a的动量方程中，得到：
$$\begin{array}{rl}& \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{\nabla p}{{\rho }_a}+g-\frac{{\alpha }_b}{{\rho }_a}K(U_a-U_b)-\frac{{\alpha }_b}{{\rho }_a}{C}_l({\alpha }_b{\rho }_b+{\alpha }_a{\rho }_a)U_r\times (\nabla \times U)\\ +& \frac{{\alpha }_b}{{\rho }_a}{C}_{vm}{\rho }_b\left[\frac{\partial U_b}{\partial t}+U_b\cdot \nabla U_b-(\frac{\partial U_a}{\partial t}+U_a\cdot \nabla U_a)\right]\end{array}$$

将相关的项合并，并调整顺序，便得到与twoPhaseEulerFoam求解器的UEqn.H文件中相同形式的分散相动量方程：
$$\begin{array}{rl}& (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{array}$$

连续相b的动量方程形式相仿，这里就不再重复了。这里有几点注意事项：

此处的双流体模型在推导的过程中，是把a当作分散相，b当作连续相的。分散相的体积分率${\alpha }_a$可以等于0，但是连续项体积分率${\alpha }_b$不能等于0 ，否则会出问题。
曳力系数 $\beta$ 的形式就是文献中常见的形式，比如，WenYu 曳力系数 $\beta =\frac{3}{4}\frac{(1-{\alpha }_b){\alpha }_b}{d_{p,a}}|U_b-U_a|C_{D0}{\alpha }_b^{-2.7}$，Ergun 曳力系数 $\beta =150\frac{(1-{\alpha }_b{)}^2{\mu }_b}{{\alpha }_b{d}_a^2}+1.75\frac{(1-{\alpha }_b){\rho }_b{U}_b-U_a}{d_a}$。而程序中定义的$K=\frac{\beta }{{\alpha }_a{\alpha }_b}$，所以，当${\alpha }_b\to 0$时，如果用WenYu曳力那还不会出错，因为曳力系数中的分子里同时含有${\alpha }_a{\alpha }_b$，运算$K=\frac{\beta }{{\alpha }_a{\alpha }_b}$不会出现除以0的问题；但如果用Ergun曳力，那就要出问题了，因为Ergun曳力系数中两项的分子都没有${\alpha }_b$，所以运算$K=\frac{\beta }{{\alpha }_a{\alpha }_b}$就要出问题了。