多相系统存在于大量的工业设备中，如鼓泡反应器、流化床、火焰喷射器、萃取塔等。多相系统在日常生活中大量存在且更加直观。瓶装水可以认为是水/气两相系统，桑拿房也可以认为是水/气两相系统。在水中加一点沙子，可以认为是固液多相系统。所有工业设备中的多相系统看起来并没有联系，但是通常我们用连续相和离散相对多相系统中的相进行区分。多相系统进而可分为单分散系统和多分散系统。单分散系统通常指分散相的物理属性是均一的，例如水中所有的泡泡是相同的（不论形状还是温度）。多分散系统通常指粒子的某些物理属性不均一，如气泡具有不同的直径大小，颗粒具有不同的温度。另外一种重要的多分散系统为颗粒的其他属性均相同，但传输速度不同，例如大的颗粒移动速度较慢，小的颗粒移动速度较慢。这种速度的非均匀分布，对多分散系统数学模型的描述产生了一定的挑战。
双流体模型也被称之为欧拉欧拉模型，是可用于描述多相系统的宏观模型。在双流体模型中，离散相和连续相均由Navier-Stokes方程描述，界面间的动量/能量传递采用模型进行模化。相对于欧拉拉格朗日模型，双流体模型得益于其计算资源耗费低的优点被大量的用于高相分数的工业过程模拟。OpenFOAM中的双流体模型求解器最早被称之为bubbleFoam，仅仅能模拟气液多相流动。随后bubbleFoam演变为twoPhaseEuerFoam，其可模拟不可压缩气液/气固流动。compressibleTwoPhaseEulerFoam则为twoPhaseEuerFoam可压缩版本。在最新版本的OpenFOAM中，multiphaseEuler模块演变的非常普适性，一个附加能量方程的可用于模拟传热传质相变的多流体模型求解器。本文主要介绍双流体模型求解器的数值离散过程，读者也可参考文献[LC17] 了解基本的离散过程和数值问题。
控制方程
在双流体模型中，离散相和连续相的动量方程可以表示为：
(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}{\boldsymbol{\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}{\boldsymbol{\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,
\]
由于有效应力\(\boldsymbol{\tau}\)不会对数值产生任何问题，其贡献仅仅在速度方程中，因此在下文中将其略去。
\(\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）直接对连续相的有效粘度进行修正。
继续考虑\(\bfM_{\mathrm{drag}}\)，定义\(\Kd\)：
(5)\[
\Kd=\frac{3}{4}\alpha_\rd\rho_\rc C_\rD\frac{1}{d_\rd} \left|\bfU_\rc-\bfU_\rd\right|.
\]
同时\(\mathrm{Re}\)：
(6)\[
\mathrm{Re}=\frac{\rho_c \left|\bfU_\rc-\bfU_\rd\right| d }{\mu_c}
\]
有：
(7)\[
\Kd=\frac{3}{4}C_\rD \mathrm{Re} \frac{ \rho_\rc \nu_c }{d^2_\rd} \alpha_\rd
\]
在这种情况下，\(\bfM_{\mathrm{drag}}\)可以简化为：
\[
\bfM_{\mathrm{drag}}=\Kd \left(\bfU_\rc-\bfU_\rd\right)
\]
离散
下面我们对方程进行离散，由于升力项、湍流分散力项等对数值不产生影响，因此将其略去。再次，略去浮力项、压力项等，对方程(1)与(2)进行整理有：
\[
\frac{{\p \left( {{\alpha_\rd }{\rho_\rd }{\bfU_\rd }} \right)}}{{\p t}}
+ \nabla \cdot \left( {{\alpha_\rd}{\rho_\rd } {{\bfU_\rd} {\bfU_\rd}} } \right)
=
\Kd \left(\bfU_\rc-\bfU_\rd\right)
\]
\[
\frac{{\p \left( {{\alpha_\rc}{\rho_\rc}{\bfU_\rc}} \right)}}{{\p t}}
+ \nabla \cdot \left( {{\alpha_\rc}{\rho_\rc} {{\bfU_\rc} {\bfU_\rc}} } \right)
=
-\Kd \left(\bfU_\rc-\bfU_\rd\right)
\]
将其重写为（曳力项处理一部分）
(8)\[
\frac{{\p \left( {{\alpha_\rd }{\rho_\rd }{\bfU_\rd }} \right)}}{{\p t}}
+ \nabla \cdot \left( {{\alpha_\rd}{\rho_\rd } {{\bfU_\rd} {\bfU_\rd}} } \right)
+\Kd\bfU_\rd
=
0
\]
(9)\[
\frac{{\p \left( {{\alpha_\rc}{\rho_\rc}{\bfU_\rc}} \right)}}{{\p t}}
+ \nabla \cdot \left( {{\alpha_\rc}{\rho_\rc} {{\bfU_\rc} {\bfU_\rc}} } \right)
+\Kd\bfU_\rc
=
0
\]
可以看出曳力项中仅仅保留\(\Kd\bfU_\rd\)和\(\Kd\bfU_\rc\)。这是因为曳力中的负系数可以提高动量方程的对角占优特性，因此将其保留。并且将其放在了方程的左边。同时，在这里为了防止出现奇异问题，OpenFOAM将\(\Kd\)进行数值处理：
(10)\[
\Kd=\frac{3}{4}C_\rD \mathrm{Re} \frac{ \rho_\rc \nu_c }{d^2_\rd} \alpha^{small}_d
\]
(11)\[
\alpha^{small}_d=\max\left({\alpha_\rd,\alpha^{small}}\right)
\]
对方程(8)和(9)进行有限体积法离散有：
(12)\[
(A_{\rd,\mathrm{P}}+f(\Kd))\mathbf{U}_{\rd,\mathrm{P}}{\rm{ + }}\sum {A_{\rd,\mathrm{N}}\mathbf{U}_{\rd,\mathrm{N}}} = S_{\rd,\mathrm{P}},
\]
(13)\[
(A_{\rd,\mathrm{P}}+f(\Kd))\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}\)表示矩阵对角系数，\(f(\Kd)\)表示曳力项离散产生的对角阵系数，\(A_\mathrm{N}\)表示临点系数，\(S_\mathrm{P}\)表示离散后产生的源项。求解方程(12)和(13)有预测速度：
(14)\[
\mathbf{HbyA}_{\rd,\mathrm{P}} = \frac{1}{{{A_{\rd,\mathrm{P}}+f(\Kd)}}}\left( { - \sum {{A_{\rd,\mathrm{N}}}\mathbf{U}_{\rd,\mathrm{N}}} + S_{\rd,\mathrm{P}}} \right),
\]
(15)\[
\mathbf{HbyA}_{\rc,\mathrm{P}} = \frac{1}{{{A_{\rc,\mathrm{P}}+f(\Kd)}}}\left( { - \sum {{A_{\rc,\mathrm{N}}}\mathbf{U}_{\rc,\mathrm{N}}} + S_{\rc,\mathrm{P}}} \right),
\]
下面组建压力泊松方程。在组建压力泊松方程之前需要进行一些必要的数值处理。首先，双流体模型中通常求解的为共享压力，即假定\(p=p_\rc=p_\rd\)。另外为简化压力边界条件的指定，同时减少非正交网格下静水压力引起的虚假速度，压力项与浮力项通常进行整合。首先定义\(p_\rgh\)为
(16)\[
p_\mathrm{rgh}=p-\rho \bfg \cdot\bfh
\]
(17)\[
\nabla p_\mathrm{rgh}=\nabla p -\rho\bfg-\bfg\cdot\bfh\nabla\rho
\]
(18)\[
\rho =\alpha_\rc\rho_\rc+\alpha_\rd\rho_\rd
\]
其中\(\bfh\)表示网格单元体心的位置矢量。同时有
(19)\[\begin{split}
-\alpha_\rd^{small} (\nabla p + \rho_\rd\bfg)
=-\alpha_\rd^{small} (\nabla p_\mathrm{rgh}- \rho\bfg- \bfg\cdot\bfh\nabla\rho+ \rho_\rd\bfg)
\\
=-\alpha_\rd^{small} (\nabla p_\mathrm{rgh}- \alpha_\rc\left(\rho_\rc-\rho_\rd\right)\bfg- \bfg\cdot\bfh\nabla\rho)
\end{split}\]
(20)\[\begin{split}
-\alpha_\rc^{small} (\nabla p + \rho_\rc\bfg)
=-\alpha_\rc^{small} (\nabla p_\mathrm{rgh}- \rho\bfg- \bfg\cdot\bfh\nabla\rho+ \rho_\rc\bfg )
\\
=-\alpha_\rc^{small} (\nabla p_\mathrm{rgh}- \alpha_\rd\left(\rho_\rd-\rho_\rc\right)\bfg- \bfg\cdot\bfh\nabla\rho )
\end{split}\]
其中关于\(p_\mathrm{rgh}\)的项可用于组建压力泊松方程，其它浮力项将用于连续性方程来修正速度。现将连续性方程(3)与(4)中的密度提出有：
(21)\[
\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},
\]
(22)\[
\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}\)表示物质导数。对方程(21)与(22)加和有：
(23)\[
\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} .
\]
方程(23)右侧表示考虑两相流中较大密度差异引致的对\(\nabla\cdot\bfU=0\)的可压缩项修正（并不同于传统的可压缩算法）。若不考虑方程(23)右侧的可压缩修正，有离散形式的散度限制性条件：
(24)\[
\sum\left(\alpha_{\rc,f}\bfU_{\rc,f}+\alpha_{\rd,f}\bfU_{\rd,f}\right)\cdot\bfS_f=0.
\]
在收敛的情况下，流动变量应满足方程(24)。因此，其可用于组建压力泊松方程。在方程(14)与(15)中，求解的预测速度并没有考虑压力梯度以及浮力项。若考虑浮力、压力、显性曳力项有：
(25)\[\begin{split}
\bfU_{\rd,\rP} = \bfHbyA_{\rd,\rP}+\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\\\\
\frac{\alpha^{small}_{\rd,\rP}}{A_{\rd,\rP}+f(\Kd)}\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}+f(\Kd)}\bfU_{\rc,\rP},
\end{split}\]
(26)\[\begin{split}
\bfU_{\rc,\rP} = \bfHbyA_{\rc,\rP}+\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\\\\
\frac{\alpha^{small}_{\rc,\rP}}{A_{\rc,\rP}+f(\Kd)}\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}+f(\Kd)}\bfU_{\rd,\rP},
\end{split}\]
将速度对面进行插值，并将其代入到方程(24)有
(27)\[
\alpha_{\rd,f}\phi_{\rd}+\alpha_{\rc,f}\phi_{\rc}
=
\nabla\cdot
\left(
\left(\alpha_{\rd,\rP}\frac{\alpha^{small}_{\rd,\rP}}{A_{\rd,\rP}}+\alpha_{\rc,\rP}\frac{\alpha^{small}_{\rc,\rP}}{A_{\rc,\rP}}
\right)\nabla p_{\mathrm{rgh},\rP}
\right),
\]
其中
(28)\[\begin{split}
\phi_{\rd}=\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\\\\
\left(\bfHbyA_{\rd,f}+\frac{\alpha^{small}_{\rd,f}\left(-\alpha_{\rc,f}\left(\rho_\rc-\rho_\rd\right)\bfg-\bfg\cdot\bfh_f\nabla\rho_f\right)}{A_{\rd,f}+f(\Kd)}+\frac{\Kd_f}{A_{\rd,f}+f(\Kd)}\bfU_{\rc,f}\right)\cdot\bfS_f
\end{split}\]
(29)\[\begin{split}
\phi_{\rc}=\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\\\\
\left(\bfHbyA_{\rc,f}+\frac{\alpha^{small}_{\rc,f}\left(-\alpha_{\rd,f}\left(\rho_\rd-\rho_\rc\right)\bfg-\bfg\cdot\bfh_f\nabla\rho_f\right)}{A_{\rc,f}+f(\Kd)}+\frac{\Kd_f}{A_{\rc,f}+f(\Kd)}\bfU_{\rd,f}\right)\cdot\bfS_f
\end{split}\]
方程(27)即不可压缩双流体模型中压力泊松方程的最终形式。求解压力之后，压力、浮力等带入到方程(25)和(26)有收敛的速度。
强耦合算法
上文中在处理曳力的时候采用一种半隐式的算法，在曳力系数\(\Kd\)远大于对角系数\(A\)的时候会出现问题。很明显，曳力系数\(\Kd\)特别大的时候，方程(28)与(29)会趋向于：
(30)\[
\phi_{\rd}=\bfU_{\rc,f}\cdot\bfS_f
\]
(31)\[
\phi_{\rc}=
\bfU_{\rd,f}\cdot\bfS_f
\]
求解方程(30)与(31)会产生停止，收敛异常缓慢。为了处理这个问题，Weller在OpenFOAM中采用另外一种更稳健的方法。在求解压力泊松方程后，定义
(32)\[
\bfU_{\rd,\rP}^s = \bfHbyA_{\rd,\rP}+\frac{\alpha_{\rd,\rP}}{A_{\rd,\rP}+f(\Kd)}\left(\nabla p_{\mathrm{rgh},\rP}-\alpha_{\rc,\rP}\left(\rho_\rc-\rho_\rd\right)\bfg-\bfg\cdot\bfh_\rP\nabla\rho_\rP\right),
\]
(33)\[
\bfU_{\rc,\rP}^s = \bfHbyA_{\rc,\rP}+\frac{\alpha_{\rc,\rP}}{A_{\rc,\rP}+f(\Kd)}\left(\nabla p_{\mathrm{rgh},\rP}-\alpha_{\rd,\rP}\left(\rho_\rd-\rho_\rc\right)\bfg-\bfg\cdot\bfh_\rP\nabla\rho_\rP\right),
\]
有：
(34)\[
\bfU_{\rd,\rP} = \bfU_{\rd,\rP}^s+\rD_{\rd,\rP}\bfU_{\rc,\rP}
\]
(35)\[
\bfU_{\rc,\rP} = \bfU_{\rc,\rP}^s+\rD_{\rc,\rP}\bfU_{\rd,\rP}
\]
(36)\[
\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}+f(\Kd)}\)，\(\rD_{\rc,\rP}=\frac{\Kd_\rP}{A_{\rc,\rP}+f(\Kd)}\)。定义
(37)\[
\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}\]
整理后有：
(38)\[
\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\)后，离散相和连续相速度可通过下式计算：
(39)\[
\bfU_{\rd,\rP}=\bfU_\rP+\alpha_{\rc,\rP}\bfU_\rP^r
\]
(40)\[
\bfU_{\rc,\rP}=\bfU_\rP-\alpha_{\rd,\rP}\bfU_\rP^r
\]
可以看出Weller算法（方程(38)）与传统算法（方程(37)）最重要的数值操作即速度的更新算法不同。考虑曳力系数\(f(\Kd)\)相对于对角系数较大的情况，上文讨论的传统方法会失效，但是在OpenFOAM中预测的\(\bfU^r\)趋向于0进而离散相与连续相会倾向为获得一个相同的速度。这也是符合物理的（因为曳力系数很大，耦合作用非常明显）。
体积力插值算法
传统算法除了相对速度会导致库郎数较大，较大的体积力源项也可能会导致求解结果产生数值震荡。例如在某些情况下（如垂直管道流），升力、壁面润滑力的作用相对重要同时数值较大。这些较大的体积力可能会导致相分数产生震荡。双流体模型中需要进行特殊的数值处理处理防止振荡与不稳定性。这种数值方法的思想即将体积力源项添加到压力方程中。依据这种思想，动量方程可以重新整理为：
(41)\[
\frac{{\p \left( {{\alpha_\rd }{\rho_\rd }{\bfU_\rd }} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rd}{\rho_\rd } {{\bfU_\rd} {\bfU_\rd}} } \right)
+\Kd\bfU_\rd = 0
\]
(42)\[
\frac{{\p \left( {{\alpha_\rc}{\rho_\rc}{\bfU_\rc}} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rc}{\rho_\rc} {{\bfU_\rc} {\bfU_\rc}} } \right)
+\Kd\bfU_\rc = 0
\]
相比较于方程(8)与(9)，很明显方程(41)与(42)的右侧并没有考虑除了曳力外的体积力源项。因此，在计算预测速度\(\bfHbyA\)之后，应考虑此处略去的体积力贡献：
(43)\[\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}\]
(44)\[\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}\]
相对应的通量为：
(45)\[\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}\]
(46)\[\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}\]
方程(45)和(46)可用于组建压力泊松方程。
相堆积问题与动理学理论
将双流体模型用于模拟固相的情况下，尤其是高相分数的固相的时候，因为固体颗粒会频繁的进行碰撞，因此模型需要进行特殊的处理：
在这里介绍相压力方法。在算例中，相压力方法需要调用phasePressure湍流模型，其本质是在固相动量方程右侧添加压力梯度项（注意固相压力梯度相之前没有乘以相分数）：
(54)\[
-\nabla p_{s}= - g_0 \min (e^{500(\alpha_{s} - \alpha_{s,\max})},\; 1000) \nabla \alpha_{s}
\]
同时固相粘度的计算方法与流体相相同。在这里注意，若使用动理学理论，固相粘度则通过动理学理论来更新。
另一方面，在模拟高相分数颗粒流动的情况下，在局部区域可能会达到最大填充密度，例如针对某种材料，\(\alpha_{s,\max}\)可能为0.63。因此求解器不能预测高于0.63的相分数。相压力方法通过在动量方程中添加源项来进行处理，然而由于固相压力梯度\(-\nabla p_{s}\)实际为相分数的梯度函数，但是其却出现在动量方程中，其存在一个解耦问题以及潜在的数值震荡。OpenFOAM在相分数方程中通过数值操作，抹平潜在的数值震荡并且使得相分数更加有界。
方程(54)在接近\(\alpha_{s,\max}\)的情况下，会出现突然增大到非常大的情况，即使在相梯度为\(1\times 10^6\)的情况下，其值也会非常大。在其他情况则是一个非常非常小的值。因此在数值上可以理解为，在某些区域达到相堆积的情况下，存在一个非常大的扩散效应将相分散开来。这与湍流分散力源项是类似的，但是湍流分散力源项与\(\alpha_{s,\max}\)无关，仅仅与相分数梯度有关。
注意：
在算例层面可以通过fvSolution中的implicitPressure开关来进行控制，在开启的情况下，固相压力梯度在相方程中也有所体现。
类似动量方程中其他体积力的处理方式，固相压力梯度项可以添加至压力方程中。在这一个步骤中参考其他体积力将\(-\nabla p_{s}\)添加至压力方程即可。同时注意到\(-\nabla p_{s}\)不同与其他体积力，其为\(\nabla\alpha_s\)的函数，因此其也可以添加至相分数方程中。参考速度与体积力的更新方法，定义一个新的速度为：
(55)\[
\bfU_{s}^{new}=\bfU_{s}+\frac{|p'|}{A}\nabla\alpha_{s}
\]
其中\(A\)表示速度方程的对角线系数，\(|p'|/A\)表示扩散系数，将在后面进行描述。使用\(\bfU_s^{new}\)来替换\(\bfU_s\)，有：
(56)\[
\bfU_r^{new}=\bfU_r +\frac{|p'|}{A}\nabla\alpha_{s}
\]
(57)\[
\bfU^{new}
= \bfU_{s}\alpha_{s}+\bfU_\rc\alpha_\rc +\alpha_{s}\frac{|p'|}{A}\nabla\alpha_{s}
=\bfU+\alpha_{s}\frac{|p'|}{A}\nabla\alpha_{s}
\]
将方程(56)跟(57)带入到(53)左边有：
(58)\[
\frac{\p\alpha_s}{\p t}
+\nabla\cdot\left( \bfU^{new} \alpha_s\right)
+\nabla\cdot\left( \bfU_r^{new} \alpha_s \alpha_\rc\right)
\]
因为其中使用了\(^{new}\)速度进行了替换，因此其并不等于0。展开有：
(59)\[\begin{split}
\begin{multline}
\frac{\p\alpha_s}{\p t}
+\nabla\cdot\left( \left( \bfU+\alpha_{s}\frac{|p'|}{A}\nabla\alpha_{s} \right) \alpha_s\right)
+\nabla\cdot\left( \left(\bfU_r +\frac{|p'|}{A}\nabla\alpha_{s}\right) \alpha_s \alpha_\rc\right)
\\\\
=\frac{\p\alpha_s}{\p t}+\nabla\cdot(\bfU\alpha_s)+\nabla\cdot(\bfU_r\alpha_s \alpha_\rc)
+\nabla\cdot \left( \alpha_{s}^2\frac{|p'|}{A}\nabla\alpha_{s} \right)
+\nabla\cdot \left( \alpha_{s}\alpha_\rc\frac{|p'|}{A}\nabla\alpha_{s} \right)
\end{multline}
\end{split}\]
再次强调，由于使用了\(^{new}\)速度进行了替换，方程(59)并不等于0。同时注意到：
(60)\[\begin{equation}
\nabla\cdot \left( \alpha_{s}^2\frac{|p'|}{A}\nabla\alpha_{s} \right) +
\nabla\cdot \left( \alpha_{s}\alpha_\rc\frac{|p'|}{A}\nabla\alpha_{s} \right)
=\nabla\cdot \left( \alpha_s\frac{|p'|}{A}\nabla\alpha_{s} \right)
\end{equation}\]
因此有：
(61)\[\begin{multline}
\frac{\p\alpha_s}{\p t}+\nabla\cdot(\bfU\alpha_s)+\nabla\cdot(\bfU_r\alpha_s \alpha_\rc)+
\nabla\cdot \left( \alpha_{s}^2\frac{|p'|}{A}\nabla\alpha_{s} \right) +
\\\\
\nabla\cdot \left( \alpha_{s}\alpha_\rc\frac{|p'|}{A}\nabla\alpha_{s} \right)
-\nabla\cdot \left( \alpha_s\frac{|p'|}{A}\nabla\alpha_{s} \right) =0
\end{multline}\]
(62)\[
\frac{\p\alpha_s}{\p t}
+\nabla\cdot\left( \bfU^{new} \alpha_s\right)
+\nabla\cdot\left( \bfU_r^{new} \alpha_s \alpha_\rc\right)
-\nabla\cdot \left( \alpha_s\frac{|p'|}{A}\nabla\alpha_{s} \right)=0
\]
方程(62)为考虑固相压力梯度后的相传输方程。
回头看\(\frac{|p'|}{A}\)系数，在方程(62)可以更明显的看出，其为一个扩散系数。在使用相压力方法的情况下，其分为2个贡献，第一个贡献来自于固体压力：
(63)\[
|p'|^1=\left|
g_0 \min (e^{500(\alpha_{s} - \alpha_{s,\max})},\; 1000)
\right|
\]
在使用动理学理论的情况下，其通过颗粒温度来进行计算。第二个贡献来自于湍流分散力（若开启这个模型）：
(64)\[
|p'|^2=D
\]
其中\(D\)表示湍流分散力系数。
[LC17]
D. Li and H. Christian. Simulation of bubbly flows with special numerical treatments of the semi-conservative and fully conservative two-fluid model. Chemical Engineering Science, 174:25–39, 2017.
