控制方程与算法
在双流体模型中，离散相和连续相的动量方程可以表示为：
(1)\[
\frac{{\p \left( {{\alpha_\rd }{\rho_\rd }{\bfU_\rd }} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rd}{\rho_\rd } {{\bfU_\rd} {\bfU_\rd}} } \right) + \nabla \cdot \left( {{\alpha_\rd}{ \rho_\rd}{\tau_\rd}} \right)
=
- {\alpha_\rd} \nabla p_\rd + {\alpha_\rd}{\rho_\rd} \bfg + {\bfM},
\]
(2)\[
\frac{{\p \left( {{\alpha_\rc}{\rho_\rc}{\bfU_\rc}} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rc}{\rho_\rc} {{\bfU_\rc} {\bfU_\rc}} } \right) + \nabla \cdot \left( {{\alpha_\rc}{\rho_\rc}{\tau_\rc}} \right)
=
- {\alpha_\rc} \nabla p_\rc + {\alpha_\rc}{\rho_\rc} \bfg - {\bfM}.
\]
其中下标\(_\rd\)和\(_\rc\)表示离散相变量和连续相变量，\(\alpha\)表示相分数，\(\rho\)表示密度，\(\bfU\)表示速度，\(p\)表示压力，\(\tau\)表示有效应力，\(\bfg\)表示重力矢量，\(\bfM\)为在欧拉框架下的界面交换力，其表示每单位体积的两相之间的动量传递。连续性方程可以表示为
(3)\[
\frac{{\p \left( {{\alpha _\rc }{\rho_\rc }} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rc }{\rho_\rc}{\bfU_\rc}} \right) = 0,
\]
(4)\[
\frac{{\p \left( {{\alpha _\rd}{\rho_\rd}} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rd}{\rho_\rd}{\bfU_\rd}} \right) = 0,
\]
考虑牛顿流体，有效应力\(\tau\)可以表示为
(5)\[
\tau_\rd=-\nu_\mathrm{d,eff}\left(\nabla \bfU_\rd+\nabla^\rT \bfU_\rd\right)+\frac{2}{3}\nu_\mathrm{d,eff}\nabla \cdot \left(\bfU_\rd \cdot\bfI\right),
\]
(6)\[
\tau_\rc=-\nu_\mathrm{c,eff}\left(\nabla \bfU_\rc+\nabla^\rT \bfU_\rc\right)+\frac{2}{3}\nu_\mathrm{c,eff}\nabla \cdot \left(\bfU_\rc \cdot\bfI\right),
\]
其中\(\nu_\mathrm{eff}\)表示有效粘度。其需要不同的湍流模型进行封闭，如kEpsilon模型或kOmega模型。\(\bfM\)通常被分解为不同的贡献。其主要可以分为轴向作用力曳力\(\bfM_\drag\)，径向作用力\(\bfM_\lift\)，湍流分散力\(\bfM_\turb\)。具体的，曳力可以表示为
\[
\bfM_{\mathrm{drag}}=\frac{3}{4}\alpha_\rd\rho_\rc C_\rD\frac{1}{d_\rd} \left|\bfU_\rc-\bfU_\rd\right| \left(\bfU_\rc-\bfU_\rd\right),
\]
其中\(C_\rD\)表示曳力系数。升力模型可以表示为
\[
\bfM_\lift=\alpha_\rd C_\rL\rho_\rc\bfU_\rr\times\left(\nabla\times\bfU_\rc\right),
\]
其中\(C_\rL\)表示升力系数，\(\bfU_\rr\)表示相对速度\(\bfU_\rc-\bfU_\rd\)。湍流分散力可以表示为
\[
\bfM_\turb=C_\rT\rho_\rc k_\rc\nabla\alpha_\rd,
\]
其中\(C_\rT\)表示湍流分散力系数。壁面润滑力可以表示为
\[
\bfM_\wall=C_\wall\rho_\rc\alpha_\rd|\bfU_\rc-\bfU_\rd|^2\cdot\bfn
\]
其中\(C_\rT\)表示壁面润滑力系数，\(\bfn\)表示壁面法向矢量。在这些作用力中，曳力被认为是最重要的界面交换力，曳力系数主要取决于流动形态以及连续性属性。升力来源于离散相在连续相运动的过程中的剪切作用，其作用于离散相运动方向的垂直方向。升力系数可正可负，其取决于不同的数学模型。不同符号的升力系数会引起离散相反向的运动。湍流分散力来源于湍流涡旋对离散相的随机扰动，类似于一种扩散行为。因此，湍流分散力的大小取决于相分数的梯度以及湍流分散力系数的大小。壁面润滑力仅仅作用于壁面附件的部分区域，作用于靠近壁面的气泡使其不会过于依附壁面。离散相和连续相之间的交互除了通过界面交换力传递之外，研究表明离散相（例如气泡）的存在会对连续相产生附加的湍流作用，其通常被称之为气泡引致湍流。气泡引致湍流可以通过1）在湍流变量上引入附加的源项进行处理，2）直接对连续相的有效粘度进行修正。
下面我们对方程进行离散，首先，略去浮力项、压力项，对方程(1)与(2)进行整理有：
(7)\[
\frac{{\p \left( {{\alpha_\rd }{\rho_\rd }{\bfU_\rd }} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rd}{\rho_\rd } {{\bfU_\rd} {\bfU_\rd}} } \right) - \nabla \cdot \left( {{\alpha_\rd}{ \rho_\rd}{\tau_\rd}} \right)
= -\Kd\bfU_\rd+\bfM_\lift+\bfM_\turb+\bfM_\wall,
\]
(8)\[
\frac{{\p \left( {{\alpha_\rc}{\rho_\rc}{\bfU_\rc}} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rc}{\rho_\rc} {{\bfU_\rc} {\bfU_\rc}} } \right) - \nabla \cdot \left( {{\alpha_\rc}{\rho_\rc}{\tau_\rc}} \right)
= -\Kd\bfU_\rc-\bfM_\lift-\bfM_\turb-\bfM_\wall
\]
其中
(9)\[
\Kd=\frac{3}{4}\alpha_\rd\rho_\rc C_\rD\frac{1}{d_\rd} \left|\bfU_\rc-\bfU_\rd\right|.
\]
读者可以看出曳力项中仅仅保留\(-\Kd\bfU_\rd\)和\(-\Kd\bfU_\rc\)。这是因为曳力中的负系数可以提高动量方程的对角占优特性，因此将其保留。对方程(7)和(8)进行有限体积法离散有：
(10)\[
{A_{\rd,\mathrm{P}}}\mathbf{U}_{\rd,\mathrm{P}}{\rm{ + }}\sum {A_{\rd,\mathrm{N}}\mathbf{U}_{\rd,\mathrm{N}}} = S_{\rd,\mathrm{P}},
\]
(11)\[
{A_{\rd,\mathrm{P}}}\mathbf{U}_{\rc,\mathrm{P}}{\rm{ + }}\sum {A_{\rc,\mathrm{N}}\mathbf{U}_{\rc,\mathrm{N}}} = S_{\rc,\mathrm{P}},
\]
其中下标\(_\rP\)表示当前网格点，下标\(_\rN\)表示相邻网格点，\(A_\mathrm{P}\)表示矩阵对角系数，\(A_\mathrm{N}\)表示临点系数，\(S_\mathrm{P}\)表示离散后产生的源项。求解方程(10)和(11)有预测速度：
(12)\[
\mathbf{HbyA}_{\rd,\mathrm{P}} = \frac{1}{{{A_{\rd,\mathrm{P}}}}}\left( { - \sum {{A_{\rd,\mathrm{N}}}\mathbf{U}_{\rd,\mathrm{N}}} + S_{\rd,\mathrm{P}}} \right),
\]
(13)\[
\mathbf{HbyA}_{\rc,\mathrm{P}} = \frac{1}{{{A_{\rc,\mathrm{P}}}}}\left( { - \sum {{A_{\rc,\mathrm{N}}}\mathbf{U}_{\rc,\mathrm{N}}} + S_{\rc,\mathrm{P}}} \right),
\]
下面组建压力泊松方程。在组建压力泊松方程之前需要进行一些必要的数值处理。首先，双流体模型中通常求解的为共享压力，即假定\(p=p_\rc=p_\rd\)。另外为简化压力边界条件的指定，同时减少非正交网格下静水压力引起的虚假速度，压力项与浮力项通常进行整合。首先定义\(p_\rgh\)为
(14)\[
p_\mathrm{rgh}=p-\rho \bfg \cdot\bfh
\]
(15)\[
\nabla p_\mathrm{rgh}=\nabla p -\rho\bfg-\bfg\cdot\bfh\nabla\rho
\]
其中\(\bfh\)表示网格单元体心的位置矢量。同时有
(16)\[\begin{split}
-\alpha_\rd\nabla p + \alpha_\rd\rho_\rd\bfg=-\alpha_\rd\nabla p_\mathrm{rgh}-\alpha_\rd\rho\bfg-\alpha_\rd\bfg\cdot\bfh\nabla\rho+\alpha_\rd\rho_\rd\bfg
\\
=\alpha_\rd\nabla p_\mathrm{rgh}-\alpha_\rd\alpha_\rc\left(\rho_\rc-\rho_\rd\right)\bfg-\alpha_\rd\bfg\cdot\bfh\nabla\rho ,
\end{split}\]
(17)\[\begin{split}
-\alpha_\rc\nabla p + \alpha_\rc\rho_\rc\bfg=-\alpha_\rc\nabla p_\mathrm{rgh}-\alpha_\rc\rho\bfg-\alpha_\rc\bfg\cdot\bfh\nabla\rho+\alpha_\rc\rho_\rc\bfg
\\
=\alpha_\rc\nabla p_\mathrm{rgh}-\alpha_\rc\alpha_\rd\left(\rho_\rd-\rho_\rc\right)\bfg-\alpha_\rc\bfg\cdot\bfh\nabla\rho
\end{split}\]
其中关于\(p_\mathrm{rgh}\)的项可用于组建压力泊松方程，其它浮力项将用于连续性方程来修正速度。现将连续性方程(3)与(4)中的密度提出有：
(18)\[
\frac{\p \alpha_\rc}{\p t}+\nabla\cdot\left(\alpha_\rc\bfU_\rc\right)=-\frac{\alpha_\rc}{\rho_\rc}\frac{\rD\rho_\rc}{\rD t},
\]
(19)\[
\frac{\p \alpha_\rd}{\p t}+\nabla\cdot\left(\alpha_\rd\bfU_\rd\right)=-\frac{\alpha_\rd}{\rho_\rd}\frac{\rD\rho_\rd}{\rD t}.
\]
其中\(\frac{\rD}{\rD t}\)表示物质导数。对方程(18)与(19)加和有：
(20)\[
\nabla\cdot\left(\alpha_\rd\bfU_\rd+\alpha_\rc\bfU_\rc\right)=-\frac{\alpha_\rd}{\rho_\rd}\frac{\rD\rho_\rd}{\rD t}-\frac{\alpha_\rc}{\rho_\rc}\frac{\rD\rho_\rc}{\rD t} .
\]
方程(20)右侧表示考虑两相流中较大密度差异引致的对\(\nabla\cdot\bfU=0\)的可压缩项修正（并不同于传统的可压缩算法）。若不考虑方程(20)右侧的可压缩修正，有离散形式的散度限制性条件：
(21)\[
\sum\left(\alpha_{\rc,f}\bfU_{\rc,f}+\alpha_{\rd,f}\bfU_{\rd,f}\right)\cdot\bfS_f=0.
\]
在收敛的情况下，流动变量应满足方程(21)。因此，其可用于组建压力泊松方程。在方程(12)与(13)中，求解的预测速度并没有考虑压力梯度以及浮力项。若考虑浮力、压力、显性曳力项有：
(22)\[
\bfU_{\rd,\rP} = \bfHbyA_{\rd,\rP}+\frac{\alpha_{\rd,\rP}}{A_{\rd,\rP}}\left(\nabla p_{\mathrm{rgh},\rP}-\alpha_{\rc,\rP}\left(\rho_\rc-\rho_\rd\right)\bfg-\bfg\cdot\bfh_\rP\nabla\rho_\rP\right)+\frac{\Kd_\rP}{A_{\rd,\rP}}\bfU_{\rc,\rP},
\]
(23)\[
\bfU_{\rc,\rP} = \bfHbyA_{\rc,\rP}+\frac{\alpha_{\rc,\rP}}{A_{\rc,\rP}}\left(\nabla p_{\mathrm{rgh},\rP}-\alpha_{\rd,\rP}\left(\rho_\rd-\rho_\rc\right)\bfg-\bfg\cdot\bfh_\rP\nabla\rho_\rP\right)+\frac{\Kd_\rP}{A_{\rc,\rP}}\bfU_{\rd,\rP},
\]
将速度对面进行插值，并将其代入到方程(21)有
(24)\[
\alpha_{\rd,f}\phi_{\rd}+\alpha_{\rc,f}\phi_{\rc}=\nabla\cdot\left(\left(\alpha_{\rd,\rP}\frac{\alpha_{\rd,\rP}}{A_{\rd,\rP}}+\alpha_{\rd,\rP}\frac{\alpha_{\rd,\rP}}{A_{\rd,\rP}}
\right)\nabla p_{\mathrm{rgh},\rP}\right),
\]
其中
\[
\phi_{\rd}=\left(\bfHbyA_{\rd,f}+\frac{\alpha_{\rd,f}}{A_{\rd,f}}\left(-\alpha_{\rc,f}\left(\rho_\rc-\rho_\rd\right)\bfg-\bfg\cdot\bfh_f\nabla\rho_f\right)+\frac{\Kd_f}{A_{\rd,f}}\bfU_{\rc,f}\right)\cdot\bfS_f
\]
\[
\phi_{\rc}=\left(\bfHbyA_{\rc,f}+\frac{\alpha_{\rc,f}}{A_{\rc,f}}\left(-\alpha_{\rd,f}\left(\rho_\rd-\rho_\rc\right)\bfg-\bfg\cdot\bfh_f\nabla\rho_f\right)+\frac{\Kd_f}{A_{\rc,f}}\bfU_{\rd,f}\right)\cdot\bfS_f
\]
方程(24)即不可压缩双流体模型中压力泊松方程的最终形式。求解压力之后，压力、浮力等带入到方程(22)和(23)有收敛的速度。
上文中在处理曳力的时候采用一种半隐式的算法，在曳力系数\(\Kd\)远大于对角系数\(A\)的时候会导致相对速度过大，进而相对速度库郎数过大，收敛较慢。Weller在OpenFOAM中采用另外一种更稳健的方法来处理压力速度耦合。在求解压力泊松方程后，定义
(25)\[
\bfU_{\rd,\rP}^s = \bfHbyA_{\rd,\rP}+\frac{\alpha_{\rd,\rP}}{A_{\rd,\rP}}\left(\nabla p_{\mathrm{rgh},\rP}-\alpha_{\rc,\rP}\left(\rho_\rc-\rho_\rd\right)\bfg-\bfg\cdot\bfh_\rP\nabla\rho_\rP\right),
\]
(26)\[
\bfU_{\rc,\rP}^s = \bfHbyA_{\rc,\rP}+\frac{\alpha_{\rc,\rP}}{A_{\rc,\rP}}\left(\nabla p_{\mathrm{rgh},\rP}-\alpha_{\rd,\rP}\left(\rho_\rd-\rho_\rc\right)\bfg-\bfg\cdot\bfh_\rP\nabla\rho_\rP\right),
\]
有：
\[
\bfU_{\rd,\rP} = \bfU_{\rd,\rP}^s+\rD_{\rd,\rP}\bfU_{\rc,\rP}
\]
\[
\bfU_{\rc,\rP} = \bfU_{\rc,\rP}^s+\rD_{\rc,\rP}\bfU_{\rd,\rP}
\]
\[
\bfU_{\rP} = \alpha_{\rd,\rP}\left(\bfU_{\rd,\rP}^s+\rD_{\rd,\rP}\bfU_{\rc,\rP}\right)+\alpha_{\rc,\rP}\left(\bfU_{\rc,\rP}^s+\rD_{\rc,\rP}\bfU_{\rd,\rP}\right)
\]
其中\(\rD_{\rd,\rP}=\frac{\Kd_\rP}{A_{\rd,\rP}}\)，\(\rD_{\rc,\rP}=\frac{\Kd_\rP}{A_{\rc,\rP}}\)。定义
(27)\[
\bfU_{\rP}^r = \bfU_{\rd,\rP}-\bfU_{\rc,\rP}=\bfU_{\rd,\rP}^s+\rD_{\rd,\rP}\bfU_{\rc,\rP}-\bfU_{\rc,\rP}^s-\rD_{\rc,\rP}\bfU_{\rd,\rP}
\]
方程左右两边减去\(\rD_{\rd,\rP}\rD_{\rc,\rP}\left(\bfU_{\rd,\rP}-\bfU_{\rc,\rP}\right)\)有：
\[\begin{split}
\begin{split}
\left(1-\rD_{\rd,\rP}\rD_{\rc,\rP}\right)\left(\bfU_{\rd,\rP}-\bfU_{\rc,\rP}\right)=\qquad\qquad\qquad\qquad\qquad\qquad\qquad
\\
\bfU_{\rd,\rP}^s+\rD_{\rd,\rP}\bfU_{\rc,\rP}-\bfU_{\rc,\rP}^s-\rD_{\rc,\rP}\bfU_{\rd,\rP}-\rD_{\rd,\rP}\rD_{\rc,\rP}\left(\bfU_{\rd,\rP}-\bfU_{\rc,\rP}\right)
\end{split}
\end{split}\]
整理后有：
(28)\[
\begin{split}
\bfU_{\rP}^r=\frac{\left(1-\rD_{\rc,\rP}\right)\bfU_{\rd,\rP}^s-\left(1-\rD_{\rd,\rP}\right)\bfU_{\rc,\rP}^s}{1-\rD_{\rc,\rP}\rD_{\rd,\rP}},
\end{split}
\]
在获得相对速度\(\bfU_{\rP}^r\)后，离散相和连续相速度可通过下式计算：
(29)\[
\bfU_{\rd,\rP}=\bfU_\rP+\alpha_{\rc,\rP}\bfU_\rP^r
\]
(30)\[
\bfU_{\rc,\rP}=\bfU_\rP-\alpha_{\rd,\rP}\bfU_\rP^r
\]
可以看出Weller算法（方程(28)）与传统算法（方程(27)）最重要的数值操作即针对相对速度的不同计算方法。考虑曳力系数相对于对角系数较大的情况（\(\rD \gg 1\)），传统方法计算的相对速度为
(31)\[
\bfU_{\rd,\rP}=\rD_{\rd,\rP}\bfU_{\rc,\rP}-\rD_{\rc,\rP}\bfU_{\rd,\rP}
\]
Weller算法计算的相对速度为
(32)\[
\bfU_{\rd,\rP}=\frac{1}{\rD_{\rd,\rP}}\bfU_{\rd,\rP}^s-\frac{1}{\rD_{\rc,\rP}}\bfU_{\rd,\rP}^s
\]
和明显，Weller算法的相对速度更加趋向于相对速度的定义，并且不会导致过大的相对速度库郎数。
除了相对速度会导致库郎数较大，较大的体积力源项也可能会导致求解结果产生数值震荡。例如在某些情况下（如垂直管道流），升力、壁面润滑力的作用相对重要同时数值较大。这些较大的体积力可能会导致相分数产生震荡。双流体模型中需要进行特殊的数值处理处理防止振荡与不稳定性。这种数值方法的思想即将体积力源项添加到压力方程中。依据这种思想，动量方程可以重新整理为：
(33)\[
\frac{{\p \left( {{\alpha_\rd }{\rho_\rd }{\bfU_\rd }} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rd}{\rho_\rd } {{\bfU_\rd} {\bfU_\rd}} } \right) - \nabla \cdot \left( {{\alpha_\rd}{ \rho_\rd}{\tau_\rd}} \right)
= -\Kd\bfU_\rd,
\]
(34)\[
\frac{{\p \left( {{\alpha_\rc}{\rho_\rc}{\bfU_\rc}} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rc}{\rho_\rc} {{\bfU_\rc} {\bfU_\rc}} } \right) - \nabla \cdot \left( {{\alpha_\rc}{\rho_\rc}{\tau_\rc}} \right)
= -\Kd\bfU_\rc,
\]
相比较于方程(7)与(8)，很明显方程(33)与(34)的右侧并没有考虑除了曳力外的体积力源项。因此，在计算预测速度\(\bfHbyA\)之后，应考虑此处略去的体积力贡献：
(35)\[\begin{split}
\begin{split}
\bfU_{\rd,\rP} = \bfHbyA_{\rd,\rP}+\frac{\alpha_{\rd,\rP}}{A_{\rd,\rP}}\left(\nabla p_{\mathrm{rgh},\rP}-\alpha_{\rc,\rP}\left(\rho_\rc-\rho_\rd\right)\bfg-\bfg\cdot\bfh_\rP\nabla\rho_\rP\right)+\qquad\qquad\qquad
\\
\frac{\Kd_\rP}{A_{\rd,\rP}}\bfU_{\rc,\rP}+\bfM_{\lift,\rP}+\bfM_{\turb,\rP}+\bfM_{\wall,\rP},
\end{split}
\end{split}\]
(36)\[\begin{split}
\begin{split}
\bfU_{\rc,\rP} = \bfHbyA_{\rc,\rP}+\frac{\alpha_{\rc,\rP}}{A_{\rc,\rP}}\left(\nabla p_{\mathrm{rgh},\rP}-\alpha_{\rd,\rP}\left(\rho_\rd-\rho_\rc\right)\bfg-\bfg\cdot\bfh_\rP\nabla\rho_\rP\right)+
\qquad\qquad\qquad
\\
\frac{\Kd_\rP}{A_{\rc,\rP}}\bfU_{\rd,\rP}-\bfM_{\lift,\rP}-\bfM_{\turb,\rP}-\bfM_{\wall,\rP},
\end{split}
\end{split}\]
相对应的通量为：
(37)\[\begin{split}
\begin{split}
\phi_{\rd}=\left(\bfHbyA_{\rd,f}+\frac{\alpha_{\rd,f}}{A_{\rd,f}}\left(-\alpha_{\rc,f}\left(\rho_\rc-\rho_\rd\right)\bfg-\bfg\cdot\bfh_f\nabla\rho_f\right)+\right.
\qquad\qquad\qquad
\\
\left.\frac{\Kd_f}{A_{\rd,f}}\bfU_{\rc,f}+\bfM_{\lift,f}+\bfM_{\turb,f}+\bfM_{\wall,f}\right)\cdot\bfS_f
\end{split}
\end{split}\]
(38)\[\begin{split}
\begin{split}
\phi_{\rc}=\left(\bfHbyA_{\rc,f}+\frac{\alpha_{\rc,f}}{A_{\rc,f}}\left(-\alpha_{\rd,f}\left(\rho_\rd-\rho_\rc\right)\bfg-\bfg\cdot\bfh_f\nabla\rho_f\right)+\right.
\qquad\qquad\qquad
\\
\left.\frac{\Kd_f}{A_{\rc,f}}\bfU_{\rd,f}-\bfM_{\lift,f}-\bfM_{\turb,f}-\bfM_{\wall,f}\right)\cdot\bfS_f
\end{split}
\end{split}\]
方程(37)和(38)可用于组建压力泊松方程。
