reactingTwoPhaseEulerFoam解析

引言

多相系统存在于大量的工业设备中,如鼓泡反应器、流化床、火焰喷射器、萃取塔等。多相系统在日常生活中大量存在且更加直观。瓶装水可以认为是水/气两相系统,桑拿房也可以认为是水/气两相系统。在水中加一点沙子,可以认为是固液多相系统。所有工业设备中的多相系统看起来并没有联系,但是通常我们用连续相和离散相对多相系统中的相进行区分。多相系统进而可分为单分散系统和多分散系统。单分散系统通常指分散相的物理属性是均一的,例如水中所有的泡泡是相同的(不论形状还是温度)。多分散系统通常指粒子的某些物理属性不均一,如气泡具有不同的直径大小,颗粒具有不同的温度。另外一种重要的多分散系统为颗粒的其他属性均相同,但传输速度不同,例如大的颗粒移动速度较慢,小的颗粒移动速度较慢。这种速度的非均匀分布,对多分散系统数学模型的描述产生了一定的挑战。

双流体模型也被称之为欧拉欧拉模型,是可用于描述多相系统的宏观模型。在双流体模型中,离散相和连续相均由Navier-Stokes方程描述,界面间的动量/能量传递采用模型进行模化。相对于欧拉拉格朗日模型,双流体模型得益于其计算资源耗费低的优点被大量的用于高相分数的工业过程模拟。OpenFOAM中的双流体模型求解器最早被称之为bubbleFoam,仅仅能模拟气液多相流动。随后bubbleFoam演变为twoPhaseEuerFoam,其可模拟不可压缩气液/气固流动。compressibleTwoPhaseEulerFoam则为twoPhaseEuerFoam可压缩版本。在最新版本的OpenFOAM中,reactingTwoPhaseEulerFoam演变的非常普适性,一个附加能量方程的可用于模拟传热传质相变的双流体模型求解器。本文主要介绍双流体模型求解器的数值离散过程,读者也可参考这篇文章了解基本的离散过程和数值问题。

_images/bubblyFlow.png

气泡流

_images/slugFlow.png

段塞流

控制方程与算法

在双流体模型中,离散相和连续相的动量方程可以表示为:

(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)可用于组建压力泊松方程。

相分离问题

在多相流中可能会存在一种现象,即某些网格区域内,某一相彻底的消失,\(\alpha=0\)。这会对方程带来数值上的问题。考虑双流体模型中的动量方程,如果\(\alpha_\rd=0\),很明显,方程(1)左右两边均为0,在数学上这被称之为奇异问题。在老版本的OpenFOAM中,并没有对方程(1)和方程(2)原本的形式进行求解,而是将方程左右两边的相分数去掉来处理。这是一种非守恒形式,在OpenFOAM研究领域通常被称之为Phase Intensive形式。在历次的OpenFOAM版本迭代中,处理的方法更加的符合物理。在这里进行简述。假定在某些网格区域,\(\alpha_\rc=0\),因此在这一部分区域,可以认为存在非常非常少量的相,OpenFOAM认为在这种情况下,\(\alpha_\rc\)变为离散相,\(\alpha_\rd\)变为连续相。同时,若连续相为静止状态,则\(\alpha_\rc\)近似以终端速度来进行移动。最终\(\bfU_\rc\)的值为终端速度。

OpenFOAM通过将动量方程简化为终端速度方程来处理的。考虑一个静水中上升的气泡或者液滴,终端速度方程即这个气泡在无限大领域,上升/下降的稳定速度。其仅仅与浮力/重力有关:

(47)\[ 0 = - {\alpha^{small}_\rc} \nabla p+ {\alpha^{small}_\rc}{\rho_\rc} \bfg - \frac{3}{4}\alpha^{small}_\rd\rho_\rc C_\rD\frac{1}{d_\rd} \left|\bfU_\rc-\bfU_\rd\right| \left(\bfU_\rc-\bfU_\rd\right), \]

进一步考虑简单的一维情况,且\(\alpha^{small}_\rc\)为用户给定的一个很小的定值,且另一相的速度为0,方程(47)可以简化为:

(48)\[ \frac{3}{4} \rho_\rc C_\rD\frac{1}{d_\rd} \left|U\right| U = 9.81\rho_\rc -\Delta p \]
(49)\[ U = \frac{ 9.81\rho_\rc -\Delta p } {\frac{3}{4} \rho_\rc C_\rD\frac{1}{d_\rd} \left|U\right| } \]

方程(49)中的\(C_\rD\)取决于速度\(U\),因此其预测的速度\(U\)会随着时间进行推进并最终稳定至终端速度。