控制方程与算法
在双流体模型中，离散相和连续相的动量方程可以表示为：
(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\)不会对数值产生任何问题，其贡献仅仅在速度方程中，因此在下文中将其略去。\(\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)进行整理有：
(5)\[
\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
\]
(6)\[
\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
\]
其中
(7)\[
\Kd=\frac{3}{4}\alpha_\rd\rho_\rc C_\rD\frac{1}{d_\rd} \left|\bfU_\rc-\bfU_\rd\right|.
\]
定义\(\mathrm{Re}\)：
(8)\[
\mathrm{Re}=\frac{\rho_c \left|\bfU_\rc-\bfU_\rd\right| d }{\mu_c}
\]
有：
(9)\[
\Kd=\frac{3}{4}C_\rD \mathrm{Re} \frac{ \rho_\rc \nu_c }{d^2_\rd} \alpha_\rd
\]
可以看出曳力项中仅仅保留\(\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)
\]
对方程(5)和(6)进行有限体积法离散有：
(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
\]
相比较于方程(5)与(6)，很明显方程(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)可用于组建压力泊松方程。
