return 0;
wmake

多相流的数学模型

如果打不开图像,请右键在新标签页打开图像后刷新几次


引言

多相系统(multiphase system)存在于大量的工业设备中,如鼓泡反应器、流化床、火焰喷射器、萃取塔等。多相系统在日常生活中大量存在且更加直观。瓶装水可以认为是水/气两相系统,桑拿房也可以认为是水/气两相系统。在水中加一点沙子,可以认为是固液多相系统。所有工业设备中的多相系统看起来并没有联系,但是通常我们用连续相(continuous phase)和离散相(disperse phase)对多相系统中的相进行区分。其中的离散相,往往存在粒径的分布。


气泡两相流

段塞流

对多相系统的研究主要分为两个方面:1)离散相的演变,2)多相流体动力学。在离散相演变的研究领域,主要关心均衡状态下的离散相粒子如何聚并破碎、如何进行碰撞,这些特征主要由群体平衡方程(population balance equation,PBE)控制。然而在许多工业设备中,离散相和连续相的耦合至关重要。目前多相流体力学的研究主要通过计算流体力学(computational fluid dynamics,CFD)来进行。多相流体力学由于离散相和连续相的耦合作用,导致其流场演变和单相情况大不相同。相对于单相系统,多相系统的数学模型还远不完善。

多相系统进而可分为单分散系统(monodisperse phase)和多分散系统(polydisperse phase)。单分散系统通常指分散相的物理属性是均一(homogenerous)的,例如水中所有的泡泡是相同的(不论形状还是温度)。多分散系统通常指粒子的某些物理属性不均一,如气泡具有不同的直径大小,颗粒具有不同的温度。另外一种重要的多分散系统为颗粒的其他属性均相同,但传输速度不同,例如大的颗粒移动速度较慢,小的颗粒移动速度较慢。这种速度的非均匀分布,对多分散系统数学模型的描述产生了一定的挑战。

如果打不开图像,请右键在新标签页打开图像后刷新几次
单分散系统以及多分散系统

在多分散系统中,既然存在非均一的粒子属性,必然就存在属性的分布。例如,如果考虑气泡的直径,那么就存在气泡的直径分布。这个属性的分布通常被称之为数量密度函数(number density function,NDF)。其在不同的领域可能具有其他的名字,如颗粒分布函数(particle distribution function)、结晶粒径分布(crystal size distribution)等。需要注意的是,多分散系统中粒子速度的分布,由于具有一定的特殊性,通常不使用NDF来描述,而是使用速度分布函数(velocity distribution function,VDF)来描述。在分子动力学领域,介尺度层级描述的多相模型即为动力学方程(kinetic equation)。

如果打不开图像,请右键在新标签页打开图像后刷新几次
NDF、VDF示意图

现在用两个例子介绍多分散系统。第一个例子为流体中细小颗粒(fine particles)输运的多分散系统。在这个多分散系统中,最重要的假定为粒子移动的速度和流体的速度相同。这个多分散系统中的离散相可以用下面的方程进行描述(Marchisio and Fox, 2013):

\begin{equation}\label{pbe} \frac{\p n}{\p t}+\nabla_\bfx\cdot\left(\bfU_\rc n\right)-\nabla_{\bfx}\cdot\left(\Gamma\nabla_\bfx n\right)=\mathbb{S}, \end{equation}

其中$n$表示NDF,$\bfU_\rc $表示连续相的传输速度,$\nabla_\bfx$表示对空间的散度操作符,$\Gamma$表示扩散系数,$\mathbb{S}$表示聚并破碎增长等源项。熟悉CFD的用户可以看出,方程\eqref{pbe}类似CFD中的标量传输方程,只不过其中的未知量为$n$,其并不是一个标量而是一个分布。如果我们用方程\eqref{pbe}描述粒子直径的变化,那么NDF可以表示为$n\left(t,\bfx,d_\rd\right)$,其中$d_\rd$表示粒子的直径。方程\eqref{pbe}也被称之为群体平衡方程(population balance equation,PBE)(Ramkrishna, 2000)。

第二个例子为考虑粒子的传输速度非均一,在这种情况下,需要考虑粒子的速度分布函数。描述这个多分散系统的数学模型通常称之为动力学方程(参见附录)。

描述分散系统的多尺度模型

微观模型、宏观模型以及介尺度模型

用于描述分散系统模型的模型通常可以分为微观模型(microscale model)、宏观模型(macroscale model),以及介尺度模型(mesoscale model),其对应多相流体力学中的三种尺度(Li, 1994)。在微观模型中,在粒子尺度以下对所有的特征长度以及特征时间进行解析。微观模型是多相流研究领域的直接模拟(DNS)。对于气固多相流动,微观模型甚至解析颗粒表面的流动边界层,以及相应的粒子之间的碰撞。

宏观模型主要调用水动力学特征(hydrodynamic description)来对不同的相进行描述。宏观模型通常被称之为多流体模型(multi-fluid model,MFM)。如果仅有一个离散相则被称之为双流体模型(two fluid model,TFM)。这是因为在宏观模型中,离散的粒子也被当做为一种流体。在宏观模型中,用户无法直接观测到粒子的外部形态。

介尺度模型是介于宏观模型以及微观模型的一种模型。拉格朗日模型,普适性群体平衡模型(generalized population balance equation,GPBE)以及EMMS模型等都是典型的介尺度模型(Marchisio and Fox, 2013; Yang et al., 2003)。方程\eqref{pbe}也是介尺度模型的一种。在拉格朗日模型中,粒子的运动由拉格朗日方程进行控制。实际上,粒子拉格朗日方程的解即为粒子动力学方程的随机解。因此拉格朗日模拟也会产生噪音。

在微观模型上进行体平均(volume-averaging)或者集合平均(emsemble-averaging)可以推导出宏观模型。在介尺度模型的基础上,对动力学方程中的密度函数取矩(moment)也可以推导出宏观模型(Struchtrup, 2005)。各种尺度间的相互推导参见附录。



多相流的多尺度模型(Sundaresan, 2018)

局限性和普适性群体平衡模型

宏观模型必然存在一些无法预测的局限。衡量介尺度模型和宏观模型适用性的方法为评估粒子的克努森数$\mathrm{Kn}$(Knudsen number)。宏观模型只适用于$\mathrm{Kn}<0.05$的情况(Struchtrup, 2005; Marchisio and Fox, 2013; Chalons et al., 2010),在这种情况下,粒子的碰撞较为频繁引致较为均衡的粒子速度分布,满足连续介质假定,也即满足水动力学(hydrodynamic)特征。在很多的多相系统中,离散相单体的运动并不是碰撞主宰的,即$\mathrm{Kn}>0.05$,在极端情况下,可能是无碰撞的,即$\mathrm{Kn}\approx\infty$。一个典型的例子即为较高大气层中的稀薄气体(rarefied gas),其中的分子几乎无碰撞。另外一个例子为无碰撞气固流动,其中的耦合作用只有离散相和连续相之间动量的传递。在这些情况下CFD中的连续介质假定不再适用。因此不建议使用宏观模型。

例如,考虑稀疏的气固两相流,在$\mathrm{Kn}\approx\infty$的情况下,颗粒马赫数趋向于无穷大,在这种情况下使用宏观模型将会出现激波(Chen and Liu, 2003)。这种流动的典型特征为颗粒轨迹交叉(particle trajectory crossing,PTC)(Wells and Stokes, 1983)。这在数学上很容易解释,这是因为双流体模型在一个空间点只能预测一个速度。当无碰撞的颗粒群穿过一个空间点的时候,颗粒的速度被平均为一个平均的速度导致无法穿越奇异点。


通过 MPPIC预测的颗粒轨迹交叉现象(O'Rourke2012 and Snider, 2012)

介尺度模型中的欧拉拉格朗日模型克服了宏观模型的缺陷。但是依然存在一些问题。对于相间耦合不是特别强的情况下,欧拉拉格朗日模拟具有较高的精准性。在某些情况下,用户希望通过增加网格分辨率来获得更精确的流场信息,但是这导致更高的拉格朗日粒子统计误差(statistics error)。在相见耦合作用比较强的时候,统计误差甚至足够大来导致流型的不同(Kaufmann et al., 2008)。另一方面,网格分辨率的降低会减少统计误差,但是却会带来离散误差(discretization error)。同时,相间的耦合需要拉格朗日粒子对于时间非常精准的信息(其并不是时均值),这会导致欧拉拉格朗日的时间步长非常小。一个完美的数学模型,需要在各个情况下都适用。对于多相流领域,本质上应该存在一个最底层的模型,其他相关模型都可由此产生。问题是这样一个模型是否存在?答案是肯定的。

普适性群体平衡模型(generalized population balance equation,GPBE)是描述多分散系统的最基本的方程。其在不同的领域具有不同的名称,例如在颗粒动力学领域又名动力学方程(kinetic equation,KE),在计算化学领域又名玻尔兹曼方程(Boltzmann equation),在气溶胶领域又名普适性气体动力学方程(generalized particle-dynamic equation)、在燃烧喷雾领域又名威廉玻尔兹曼方程(Williams-Boltzmann equation)。GPBE属于一种介尺度模型,其主要依据概率密度函数(probability density function,PDF),在六自由度下来描述粒子的运动状态。GPBE的重要性在于其可通过6D非线性求解器进行直接求解,也可以通过求矩来演化为矩方法、多流体模型甚至MPPIC(multi-phase particle-in-cell)。GPBE目前为认可度最高的用于描述多分散系统最基本的控制方程,其可以表示为(Marchisio and Fox, 2013)

\begin{equation}\label{GPBE} \frac{\p n}{\p t}+\nabla_\bfx\cdot\left(\bfU_\rd n\right)+\nabla_{\bfU_\rd}\cdot\left(\bfA n\right)=\mathbb{C}, \end{equation}

其中$\bfU_\rd$表示粒子传输速度,$\nabla_{\bfU_\rd}$表示在粒子传输速度空间进行的散度操作,$\bfA$表示外力引起的粒子加速,$\mathbb{C}$表示粒子碰撞引起的运动。很显然,方程\eqref{GPBE}即为分子模拟中的玻尔兹曼方程。

数学上,之所以说GPBE为多分散系统的最基本模型,是因为对GPBE使用Chapman-Enskog模型对$\mathrm{Kn}$进行展开,可以分别得到离散相欧拉方程($\mathrm{Kn}=0$,$0$阶展开)、离散相NS方程($\mathrm{Kn}<0.1$,$1$阶展开)以及离散相Burnett方程($\mathrm{Kn}>0.1$,$2$阶展开)。同时,宏观模型的速度、相分数、温度等均为GPBE中概率分布函数的统计学高阶矩(Struchtrup, 2005; Marchisio and Fox, 2013)。例如,在TFM中,可认为是对离散相调用二阶矩模型。

控制方程

离散相模型之拉格朗日模型(介尺度模型)

考虑单一的粒子,其占据的体积为$V_\dpm$,质量为$m_\dpm$,密度为$\rho_\dpm$,瞬时速度为$\bfU_\dpm$。依据牛顿第二定律,粒子运动主要由下面的方程控制

\begin{equation}\label{DPMmass} \mathrm{All.forces.acting.on.particle}+\mathrm{Momentum.change.due.to.mass.transfer}=\frac{\rd }{\rd t}\left(m_\dpm\bfU_\dpm\right). \end{equation}

在这里不考虑传质,那么粒子其受力$\bfF_\dpm$主要有曳力(阻力)$\bfF_{\drag,\dpm}$、重力$\bfF_{\mathrm{g},\dpm}$、以及浮力$\bfF_{\mathrm{buo},\dpm}$。其中曳力、浮力为表面力(surface force),重力属于体积力(body force):

\begin{equation}\label{DPMmass1} \bfF_\dpm=\bfF_{\drag,\dpm}+\bfF_{\mathrm{g},\dpm}+\bfF_{\mathrm{buo},\dpm}. \end{equation}

其中重力可以表示为

\begin{equation}\label{gravity} \bfF_{\mathrm{g},\dpm}\left[\frac{\mathrm{kg\cdot m}}{\mathrm{s}^2}\right]=m_\dpm \bfg\left[\frac{\mathrm{m}}{\mathrm{s}^2}\right]. \end{equation}

其中浮力可以表示为

\begin{equation}\label{buoyant} \bfF_{\mathrm{buo},\dpm}\left[\frac{\mathrm{kg\cdot m}}{\mathrm{s}^2}\right]=-\frac{m_\dpm}{\rho_\dpm}\nabla p. \end{equation}

其中曳力可以表示为

\begin{equation}\label{drag} \bfF_{\drag,\dpm}\left[\frac{\mathrm{kg}\cdot\mathrm{m}}{\mathrm{s}^2}\right]=\frac{m_\dpm}{\rho_\dpm}\frac{3}{4}\frac{C_\rD\rho_\rc}{d_\rd}\left|\bfU_\rc-\bfU_\dpm\right|\left(\bfU_\rc-\bfU_\dpm\right), \end{equation}

$C_\rD$表示曳力系数。方程\eqref{drag}描述的为每单位质量粒子所受的曳力。依据牛顿第二定律,单一粒子的运动遵守下面的方程

\begin{equation}\label{DPM0} \bfF_{\dpm}\left[\frac{\mathrm{kg\cdot m}}{\mathrm{s}^2}\right]=\frac{\rd }{\rd t}\left(m_\dpm\bfU_\dpm\right). \end{equation}

方程\eqref{DPM0}表示粒子在$\bfF_{\drag,\dpm}$的作用下产生的位移。在不考虑传质以及破碎等情况下,$m_\rd$为常量可以提出:

\begin{equation}\label{DPM} m_\dpm\frac{\rd \bfU_\dpm}{\rd t}=-\frac{m_\dpm}{\rho_\dpm}\nabla p+m_\dpm\bfg+\frac{m_\dpm}{\rho_\dpm}\mathrm{Kd}\left(\bfU_\rc-\bfU_\dpm\right)+\mathrm{OtherForces}. \end{equation}

其中$\mathrm{Kd}$为

\begin{equation}\label{Kd} \mathrm{Kd}=\frac{3}{4}\frac{C_\rD\rho_\rc}{d_\rd}\left|\bfU_\rc-\bfU_\dpm\right|. \end{equation}

方程\eqref{DPM}为一个常微分方程,可用常规的方法求解。如果仅仅考虑曳力项,其可以简写为

\begin{align} \frac{\rd \bfU_\dpm}{\rd t}&=\bfg+\frac{1}{\rho_\dpm}\nabla p+\frac{1}{\rho_\dpm}\Kd\left(\bfU_\rc-\bfU_\dpm\right) \label{ODE} \end{align}

对方程\eqref{ODE}进行移项有(略去$_\dpm$下标)

\begin{align}\label{ODE2} \frac{1}{\left(\bfg+\nabla p/\rho+\Kd\left(\bfU_\rc-\bfU\right)/\rho\right)}\rd \bfU &= \rd t \end{align}

对其做积分有:

\begin{align} \mathrm{ln}\left|\bfg+\frac{\nabla p}{\rho}+\frac{\Kd\left(\bfU_\rc-\bfU\right)}{\rho} \right|&=-\frac{\mathrm{Kd}}{\rho}t+C \label{ODE3} \end{align}

当$t=0$的时候,$\bfU=\bfU^t$,有

\begin{equation}\label{C} C=\mathrm{ln}\left|\bfg+\frac{\nabla p}{\rho}+\frac{\Kd\left(\bfU_\rc-\bfU^t\right)}{\rho} \right|. \end{equation}
代入有精确解为:

\begin{equation} \left|\bfg+\frac{\nabla p}{\rho}+\frac{\Kd\left(\bfU_\rc-\bfU^{t+\Delta t}\right)}{\rho} \right|=\exp\left(-\frac{\Kd}{\rho}\Delta t\right)\left|\bfg+\frac{\nabla p}{\rho}+\frac{\Kd\left(\bfU_\rc-\bfU^t\right)}{\rho}\right| \end{equation}

整理后有

\begin{align} \bfU^{t+\Delta t} &=\bfU^t e^{-\Kd\cdot\Delta t/\rho}+\frac{\rho\bfg+\nabla p+\Kd\cdot\bfU_\rc}{\Kd}\left(1-e^{-\Kd\cdot\Delta t/\rho}\right) \notag\\ &=\bfU^t e^{-\Kd\cdot\Delta t/\rho}+\left(\bfU_\rc+\frac{\rho\bfg+\nabla p}{\Kd}\right)\left(1-e^{-\Kd\cdot\Delta t/\rho}\right) \end{align}

此即为求解拉格朗日速度的精确解,对应开源CFD代码OpenFOAM中的analytical时间格式。同时,依据距离等于速度和时间的乘积有

\begin{equation}\label{DPM2} \frac{\rd \bfx}{\rd t}=\bfU_\dpm. \end{equation}

方程\eqref{DPM}和\eqref{DPM2}也即为离散相的拉格朗日模型。拉格朗日模型用于在拉格朗日框架下更新粒子的速度并更新粒子的位置信息。拉格朗日模型在数值上求解比欧拉模型容易,且不存在假扩散现象,但对于耦合插值以及高效率的粒子跟踪算法(并行计算)需要特殊的技巧。

需要提及的是,方程\eqref{DPM}可以进一步推导。依据伯努利定律有

\begin{equation}\label{bnl} \nabla p=\rho_\rc\bfg \end{equation}

将方程\eqref{bnl}代入到\eqref{DPM},有

\begin{equation}\label{DPM3} m_\dpm\frac{\rd \bfU_\dpm}{\rd t}=m_\dpm\bfg\left(1-\frac{\rho_\rc}{\rho_\dpm}\right)+\frac{m_\dpm}{\rho_\dpm}\frac{3}{4}\frac{C_\rD\rho_\rc}{d_\rd}\left|\bfU_\rc-\bfU_\dpm\right|\left(\bfU_\rc-\bfU_\dpm\right)+\mathrm{OtherForces}. \end{equation}

离散相模型之欧拉模型(宏观模型)

离散相如果使用欧拉模型描述,那么多相系统模型可称之为欧拉欧拉模型(eulerian eulerian model)。欧拉模型为一种宏观模型。双流体模型(two fluid model)为一个只有一个离散相一个连续相的宏观模型,主要用来描述高相分数两相混合。另外需要提及的是,对于相分数较低的情况下,可以使用体分数模型(volume of fraction)对相界面直接捕获,其近似为多相流的直接模拟(微观模型)。

离散相欧拉模型的动量方程可以表示为:

\begin{equation}\label{momentum1} \frac{{\p \left( {{\alpha_\rd }{\rho_\rd }{\bfU_\rd }} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rd}{\rho_\rd }\left( {{\bfU_\rd} \otimes {\bfU_\rd}} \right)} \right) - \nabla \cdot \left( {{\alpha_\rd}{ \rho_\rd}{\bfR_\rd}} \right) = - {\alpha_\rd} \nabla p_\rd + {\alpha_\rd}{\rho_\rd} \bfg + {\bfM_\drag}, \end{equation}

其中$\alpha_\rd$表示离散相的相分数,$\rho_\rd$表示密度,$\bfU_\rd$表示速度,$p_\rd$表示压力,$\bfR_\rd$表示雷诺应力,$\bfg$表示重力矢量,$\bfM_\drag$为在欧拉框架下的界面交换力(此处仅仅考虑曳力),其表示每单位体积每单位时间下两相之间的动量传递,单位为$\frac{\mathrm{kg}\cdot\mathrm{m}}{\mathrm{s}}/\mathrm{m}^3/\mathrm{s}=\frac{\mathrm{kg}}{\mathrm{m}^2\mathrm{s}^2}$。

在这里要注意的是,表面张力在双流体模型中并没有考虑,因为其相对于其他力(如曳力)等的作用较小。这也是目前一个研究问题:泡沫模拟。粒子动量变化的产生原因主要是由于体积力以及粒子-液体之间的作用力。在各种体积力中,重力是最为重要的。但是由于重力不仅作用在粒子上,也作用在液体上,因此重力并不会产生动量交换。在界面交换力中,最重要的力为浮力以及曳力。

观察方程\eqref{momentum1},可见其左边和单相的动量方程类似。欧拉模型主要需要处理的即为方程右边的界面交换相$\bfM_\drag$。$\bfM_\drag$以及$\alpha_\rd$负责两相之间的耦合作用。

在欧拉框架下,并不对单个的粒子进行追踪,因此在离散相欧拉模型下无法展现单独的颗粒/气泡/液滴。但方程\eqref{momentum1}中的$- {\alpha_\rd} \nabla p_\rd + {\alpha_\rd}{\rho_\rd} \bfg + {\bfM_\drag}$可以从离散相的拉格朗日模型导出。如果我们考虑每单位体积的粒子所受的力,对方程\eqref{gravity}、\eqref{buoyant}以及\eqref{drag}左右除以$V_\dpm$有:

\begin{equation}\label{dragEuler21} \bfA_{\mathrm{g}}\left[\frac{\mathrm{kg}}{\mathrm{m}^2\mathrm{s}^2}\right]=\rho_\dpm \bfg. \end{equation}
\begin{equation}\label{dragEuler22} \bfA_{\mathrm{buo}}\left[\frac{\mathrm{kg}}{\mathrm{m}^2\mathrm{s}^2}\right]=\nabla p. \end{equation}
\begin{equation}\label{dragEuler2} \bfA_{\drag}\left[\frac{\mathrm{kg}}{\mathrm{m}^2\mathrm{s}^2}\right]=\frac{\bfF_{\drag,\dpm}}{V_\dpm}=\underset{\mathrm{Kd}}{\underbrace{\frac{3}{4}\frac{C_\rD \rho_\rc}{d_\rd}\left|\bfU_\rc-\bfU_\dpm\right|}} \left(\bfU_\rc-\bfU_\dpm\right). \end{equation}

其中$\bfA$的单位为力的单位除以体积的单位,即每单位体积所受的曳力加速项(见附录),其符号标识取决于英文Acceleration;下括号$\mathrm{Kd}\left[\frac{\mathrm{kg}}{\mathrm{m}^3\mathrm{s}}\right]$表示GPBEFoam代码中的Kd。$\bfA$和$\bfM$的单位是一致的。

在这里,从粒子的加速项$\bfA_\drag$可以定义弛豫时间$\tau$:

\begin{equation}\label{time} \bfA_{\drag}=\frac{m_\dpm}{\tau}\left(\bfU_\rc-\bfU_\dpm\right). \end{equation}

在欧拉框架下,$\bfM$定义在每个网格单元上。从物理意义上来解释,其为每个网格单元$V_\cell$内所有粒子粒子的受力加和并除以网格单元的体积,即

\begin{equation}\label{MEuler} \bfM_\drag\left[\frac{\mathrm{kg}}{\mathrm{m}^2\mathrm{s}^2}\right]=\frac{\sum \left(V_\dpm A_\drag\right)}{V_\cell}=A_\drag\frac{\sum V_\dpm }{V_\cell}. \end{equation}

方程\eqref{MEuler}引入了一个重要的假定,即对于不同的粒子,有

\begin{equation}\label{UU} \bfU_\dpm=\bfU_\rd. \end{equation}

这个假定一方面是合理的,因为在宏观模型中,假定颗粒大小为均一的,因此质量为均一的。但这个假定也是宏观模型的缺陷之一。在介尺度模型中,颗粒大小非均一,因此可以具有不同的运动速度。

依据相分数的定义

\begin{equation}\label{alpha} \alpha_\rd=\frac{\sum V_\dpm}{V_\cell}. \end{equation}
\begin{equation}\label{Md} \bfM_\drag=\alpha_\rd \frac{3}{4} \frac{C_\rD\rho_\rc}{d_\rd}|\bfU_\rc-\bfU_\rd|\left(\bfU_\rc-\bfU_\rd\right). \end{equation}

方程\eqref{Md}即为欧拉框架下界面交换相中曳力的定义。同理,欧拉框架下的浮力以及重力项可以从方程\eqref{dragEuler22}和\eqref{dragEuler21}导出,即$\alpha_\rd\nabla p_\rd$和$\alpha_\rd\rho_\rd\bfg$:

在欧拉框架下,$\bfM_\drag$和粒子的弛豫时间$\tau$的关系为

\begin{equation}\label{Frt} \bfM_\drag=\frac{\sum m_\dpm}{V_\mathrm{cell}}\frac{1}{\tau}\left(\bfU_\rc-\bfU_\rd\right)=\rho_\rd\frac{\sum V_\dpm}{V_\mathrm{cell}}\frac{1}{\tau}\left(\bfU_\rc-\bfU_\rd\right)=\frac{\rho_\rd\alpha_\rd}{\tau}\left(\bfU_\rc-\bfU_\rd\right) \end{equation}

结合方程\eqref{Md}和\eqref{Frt},有欧拉框架下粒子的弛豫时间为

\begin{equation}\label{tau} {\tau}=\frac{4}{3}\frac{\rho_\rd d_\rd}{C_\rD \rho_\rc |\bfU_\rc-\bfU_\rd|} \end{equation}

其中的曳力系数$C_\rD$可以通过不同的模型进行模化,比如SchillerNaumann模型定义的$C_\rD$如下

\begin{equation}\label{SN} C_\rd=\left\{\begin{matrix} \frac{24}{\mathrm{Re}}\left(1+0.15\mathrm{Re}^{0.687}\right) & \mathrm{Re}\leq 1000\\ 0.445 & \mathrm{Re}> 1000 \end{matrix}\right.. \end{equation}

其中$\mathrm{Re}$为粒子的雷诺数,定义为

\begin{equation}\label{Re} \mathrm{Re}=\frac{d_\rd\rho_\rc|\bfU_\rc-\bfU_\rd|}{\mu_\rc}. \end{equation}

如果假定粒子的雷诺数非常小(爬流),另一个比较简单的曳力常数$C_\rD$可以这样计算

\begin{equation}\label{CdRe} C_\rD=\frac{24}{\mathrm{Re}}. \end{equation}

将方程\eqref{CdRe},\eqref{Re}代入到方程\eqref{tau}中,有欧拉框架下爬流中的弛豫时间为

\begin{equation}\label{tauCreep} {\tau}=\frac{\rho_\rd d_\rd^2}{18\mu_\rc}. \end{equation}

最后,将方程\eqref{Md}代入到方程\eqref{momentum1},则离散相的动量方程为(区别仅仅在将$\bfM_\drag$展开)

\begin{equation}\label{momentum3} \frac{{\p \left( {{\alpha_\rd\rho_\rd }{\bfU_\rd }} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rd}\rho_\rd\bfU_\rd\otimes\bfU_\rd} \right) - \nabla \cdot \left( {{\alpha_\rd}\rho_\rd{\bfR_\rd}} \right) = - \alpha_\rd \nabla p_\rd + {\alpha_\rd}\rho_\rd \bfg + \alpha_\rd\frac{3}{4} \frac{C_\rD\rho_\rc}{d_\rd}|\bfU_\rc-\bfU_\rd|\left(\bfU_\rc-\bfU_\rd\right), \end{equation}

有关欧拉模型在OpenFOAM中的离散求解,可以参考之前写的一篇文章

连续相模型之欧拉模型(宏观模型)

在微观模型中,连续相的模化不需要模化体积分数。在介尺度模型以及宏观模型中,连续相的模型均使用欧拉模型来描述:

\begin{equation}\label{momentum2} \frac{{\p \left( {{\alpha_\rc}{\rho_\rc}{\bfU_\rc}} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rc}{\rho_\rc}\left( {{\bfU_\rc} \otimes {\bfU_\rc}} \right)} \right) - \nabla \cdot \left( {{\alpha_\rc}{\rho_\rc}{\bfR_\rc}} \right) = - {\alpha_\rc} \nabla p_\rc + {\alpha_\rc}{\rho_\rc} \bfg - {\bfM_\drag}. \end{equation}

方程\eqref{momentum2}即连续相的欧拉模型。进一步的,连续相的方程\eqref{momentum2}结合离散相的方程\eqref{DPM}和\eqref{DPM2}构成多相流的欧拉拉格朗日模型(eulerian-lagrangian method)。连续相的方程\eqref{momentum2}结合离散相的方程\eqref{momentum1}构成多相流的欧拉欧拉模型(eulerian-eulerian method)

速度压力耦合数值求解

在介尺度模型和宏观模型中,连续相和离散相的相互作用可以区分为单向耦合、双向耦合以及四向耦合。在单向耦合中,连续相影响离散相的运动且离散相不影响连续相。在这种情况下,连续相的控制方程演化为单相流的控制方程。在双向耦合中,连续相和离散相相互影响。在四向耦合中,离散相间的影响也需要同时考虑进来。

同时,在上文的推导中可以看出,不管是欧拉拉格朗日模型、还是欧拉欧拉模型,或者其他模型如Eulerian-QBMM等,连续相模型均使用欧拉方法进行描述。在求解中,连续相的求解大同小异,只需要处理耦合项即可。下面讨论欧拉拉格朗日下的单向耦合和双向耦合中的连续相的速度压力更新方法。

单向耦合

在单向耦合中,不考虑相分数以及粒子对连续相的影响,方程\eqref{momentum2}简化为单相流动量方程:

\begin{equation}\label{oneway} \frac{{\p \left(\rho_\rc\bfU_\rc\right) }}{{\p t}} + \nabla \cdot {\left( \rho_\rc{{\bfU_\rc} \otimes {\bfU_\rc}} \right)}- \nabla \cdot {\rho_\rc{\bfR_\rc}} = - \nabla p_\rc. \end{equation}

方程\eqref{oneway}的求解请参考icoFoam解析

双向耦合

考虑不可压缩的方程\eqref{momentum2},并将方程\eqref{Md}代入有(此处采用有相分数的方程,参见讨论欧拉及拉格朗日下的动量方程压力项竟然不一样):

\begin{equation}\label{NS1} \frac{{\p \left( {{\alpha_\rc}\rho_\rc{\bfU_\rc}} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rc\rho_\rc}\left( {{\bfU_\rc} \otimes {\bfU_\rc}} \right)} \right) - \nabla \cdot \left( {{\alpha_\rc\rho_\rc}{\bfR_\rc}} \right) =-\alpha_\rc\nabla p_\rc + {\alpha_\rc}\rho_\rc\bfg - \alpha_\rd\underset{\mathrm{Kd}}{\underbrace{ \frac{3}{4} \frac{C_\rD\rho_\rc}{d_\rd}|\bfU_\rc-\bfU_\rd|}}\left(\bfU_\rc-\bfU_\rd\right). \end{equation}

去掉方程\eqref{NS1}中的阿基米德浮力项、重力项、以及界面交换项中的$\bfU_\rd$项:

\begin{equation}\label{NS2} \frac{{\p \left( {{\alpha_\rc}\rho_\rc{\bfU_\rc}} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rc}\rho_\rc\left( {{\bfU_\rc} \otimes {\bfU_\rc}} \right)} \right) - \nabla \cdot \left( {{\alpha_\rc}\rho_\rc{\bfR_\rc}} \right) =-\alpha_\rd\underset{\mathrm{Kd}}{\underbrace{ \frac{3}{4} \frac{C_\rD\rho_\rc}{d_\rd}|\bfU_\rc-\bfU_\rd|}}\bfU_\rc. \end{equation}

并对方程\eqref{NS2}离散化后(参考icoFoam解析),可组建速度矩阵(在开源CFD代码OpenFOAM中为fvVectorMatrix)。对此速度矩阵求解后,有预测速度$\mathbf{HbyA}$:

\begin{equation} \mathbf{HbyA}=\frac{1}{A_\mathrm{P}\left(\alpha_\rc,\rho_\rc\right)}\left( { - \sum {{A_\mathrm{N}\left(\alpha_\rc,\rho_\rc\right)}\mathbf{U}_\mathrm{N}^r}+...} \right) \label{hbya} \end{equation}

在$\mathbf{HbyA}$的基础上,考虑重力、阿基米德浮力以及界面交换项的影响,同时略去$A_\mathrm{P}\left(\alpha_\rc,\rho_\rc\right)$中的$\alpha_\rc,\rho_\rc$:

\begin{equation} \bfU_\rc=\mathbf{HbyA}+\frac{1}{A_\mathrm{P}}\left(-\alpha_\rc\nabla p_\rc+\alpha_\rc\rho_\rc\bfg+\alpha_\rd\underset{\mathrm{Kd}}{\underbrace{ \frac{3}{4} \frac{C_\rD\rho_\rc}{d_\rd}|\bfU_\rc-\bfU_\rd|}}\bfU_\rd\right) \label{hbyaU} \end{equation}

考虑连续相的连续性方程(Drew, 1983)

\begin{equation}\label{cont1} \frac{\p \alpha_\rc\rho_\rc}{\p t}+\nabla\cdot\left(\alpha_\rc\rho_\rc\bfU_\rc\right)=0. \end{equation}

方程\eqref{hbyaU}中的$\bfU_\rc$应满足离散后的方程\eqref{cont1}。现构建方程\eqref{hbyaU}的通量

\begin{equation} \phi_\rc=\mathbf{HbyA}_f\cdot\bfS_f+\frac{1}{A_{\mathrm{P},f}}\left(-\alpha_{\rc,f}\nabla_f p_\rc+\alpha_{\rc,f}\rho_{\rc,f}\bfg_f+\alpha_{\rd,f}\underset{\mathrm{Kd}_f}{\underbrace{ \frac{3}{4}\frac{C_\rD \rho_{\rc,f}}{d_\rd}|\bfU_{\rc,f}-\bfU_{\rd,f}|}}\bfU_{\rd,f}\right)\cdot\bfS_f, \label{phi} \end{equation}
也即
\begin{equation} \phi_\rc=\mathbf{HbyA}_f\cdot\bfS_f+\mathbf{Pre}\cdot\bfS_f-\frac{\alpha_{\rc,f}}{A_{\mathrm{P},f}}\left(\nabla_f p_\rc\right)\cdot\bfS_f. \label{phi2} \end{equation}
其中
\begin{equation} \mathbf{Pre}=\frac{\alpha_{\rc,f}\rho_{\rc,f}}{A_{\mathrm{P},f}}\bfg_f+\frac{1}{A_{\mathrm{P},f}}\left(\alpha_{\rd,f}\underset{\mathrm{Kd}_f}{\underbrace{ \frac{3}{4}\frac{C_\rD \rho_{\rc,f}}{d_\rd}|\bfU_{\rc,f}-\bfU_{\rd,f}|}}\bfU_{\rd,f}\right). \label{pre} \end{equation}

其中的$\phi_\rc$应满足离散后的方程\eqref{cont1},即对应压力泊松方程:

\begin{equation} \frac{\p \alpha_\rc\rho_\rc}{\p t}+\nabla\cdot\left(\alpha_\rc\rho_\rc\left(\mathbf{Pre}+\mathbf{HbyA}\right)\right)=\nabla\cdot\left(\alpha_\rc\rho_\rc\frac{\alpha_{\rc}}{A_{\mathrm{P}}}\nabla p_\rc\right). \label{pressure} \end{equation}

在这里需要提及一下,对于存在静水压力以及重力的压力方程,直接求解$p$方程容易在计算域底部引起假拟速度的出现(李东岳, 2018)。一种解决办法是定义一个去水压力p_rgh:

\begin{equation} p_\rgh=p_\rc-\rho\bfg\bfh, \label{prgh} \end{equation}
其中$\bfh$表示网格单元的位置矢量。对去水压力进行散度操作有
\begin{equation} \nabla p_{\mathrm{rgh}}=\nabla p_\rc - \rho_\rc \mathbf{g}-\mathbf{g} \mathbf{h} \nabla \rho_\rc, \label{prgh2} \end{equation}
\begin{equation} -\alpha_\rc\nabla p_\rc + \alpha_\rc\rho_\rc \mathbf{g}=-\alpha_\rc\left(\nabla p_{\mathrm{rgh}}+\mathbf{g} \mathbf{h} \nabla \rho_\rc\right). \label{prgh3} \end{equation}

这样,方程\eqref{NS1}变为去水压力方程

\begin{equation}\label{NSprgh} \frac{{\p \left( {{\alpha_\rc}\rho_\rc{\bfU_\rc}} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rc\rho_\rc}\left( {{\bfU_\rc} \otimes {\bfU_\rc}} \right)} \right) - \nabla \cdot \left( {{\alpha_\rc\rho_\rc}{\bfR_\rc}} \right) =-\alpha_\rc\left(\nabla p_{\mathrm{rgh}}+\mathbf{g} \mathbf{h} \nabla \rho_\rc\right) - \alpha_\rd\underset{\mathrm{Kd}}{\underbrace{ \frac{3}{4} \frac{C_\rD\rho_\rc}{d_\rd}|\bfU_\rc-\bfU_\rd|}}\left(\bfU_\rc-\bfU_\rd\right). \end{equation}

对方程\eqref{NSprgh}求解,可得更加符合物理的解。


附录

动力学方程

定义粒子的可能性密度函数为$n(t,\bfx,\bfU)$,那么动力学方程可以表示为

\begin{equation}\label{ke} \frac{\p n}{\p t}+\nabla_\bfx\cdot\left(\bfU_\rd n\right)+\nabla_{\bfU_\rd}\cdot\left(\bfA n\right)=\mathbb{C}, \end{equation}

其中$\bfA$表示加速项,$\mathbb{C}$表示碰撞项。考虑阿基米德浮力、重力、以及曳力的加速项可以表示为

\begin{equation}\label{A} \bfA=-\frac{1}{\rho_\rd}\nabla p_\rc +\bfg +\frac{1}{\tau}\left(\bfU_\rc -\bfU_\rd\right). \end{equation}

碰撞项可以表示为

\begin{equation} \mathbb{C}=\frac{1}{\tau_c}(n_\mathrm{eq}-n), \end{equation}

其中$n_\mathrm{eq}$表示均衡状态的密度函数分布。

将用于描述粒子运动的动力学方程,即方程\eqref{ke},\eqref{A},\eqref{C},结合流体的NS方程,即可用于描述多相系统的流动状态。在这里再次强调的是,动力学方程的求解需要6D求解器。因此,直接求解动力学方程是很困难的。

如果对动力学方程从数学上进行操作,可以演变为多相流求解方法中的矩方法(quadrature based moment methods)、多流体模型。矩方法本文从略,下面简述从介尺度模型之动力学方程推到之宏观模型之NS方程。

在宏观模型中,我们对概率密度函数的信息不感兴趣(并且也无法获知),因此我们只关心粒子概率密度函数的统计学信息。因此我们可以对动力学方程\eqref{ke}左右两边取矩:

\begin{equation}\label{ke1} \int_{-\infty}^\infty\left(\frac{\p n}{\p t}+\nabla_\bfx\cdot\left(\bfU_\rd n\right)+\nabla_{\bfU_\rd}\cdot\left(\bfA n\right)\right)\rd \bfU_\rd=\int_{-\infty}^\infty\mathbb{C}\rd \bfU_\rd, \end{equation}

考虑方程\eqref{ke1}的第一项,其实际上即为粒子体分数的定义

\begin{equation}\label{kealpha} \int_{-\infty}^\infty n \rd\bfU_\rd=\alpha_\rd. \end{equation}

同时依据速度矩的定义有

\begin{equation}\label{vm} \int_{-\infty}^\infty\left(\bfU_\rd n\right) \rd\bfU_\rd=\alpha_\rd\bfU_\rd. \end{equation}

考虑方程\eqref{ke1}的第二项,结合方程\eqref{vm}有

\begin{equation}\label{kealphaU} \int_{-\infty}^\infty \left(\nabla_\bfx\cdot\left(\bfU_\rd n\right)\right)\rd\bfU_\rd=\nabla\cdot\left(\alpha_\rd\bfU_\rd\right). \end{equation}

考虑方程\eqref{ke1}的第三项有

\begin{equation}\label{keA} \int_{-\infty}^\infty\left(\nabla_{\bfU_\rd}\cdot\left(\bfA n\right)\right)\rd \bfU_\rd=An|_{-\infty}^{\infty}=0, \end{equation}

其中调用了隐含的条件即概率密度函数在两端趋向于$0$。考虑方程\eqref{ke1}的第四项,依据质量守恒有

\begin{equation}\label{keC} \int_{-\infty}^\infty\left(\mathbb{C} \right)\rd \bfU_\rd=0. \end{equation}

将方程\eqref{kealpha}、\eqref{kealphaU}、\eqref{keA}、以及\eqref{keC}代入到方程\eqref{ke1}有

\begin{equation}\label{cont} \frac{\p \alpha_\rd}{\p t}+\nabla_\bfx\cdot\left(\alpha_\rd\bfU_\rd\right)=0. \end{equation}

方程\eqref{cont}即为宏观模型中的连续性方程。宏观模型中的动量方程可以参考类似的方法得出。

Eulerian-QBMM

在正文的推导中,可以看出双流体模型只具备一个离散相速度$\bfU_\rd$。这是一个非常严格的假定,在某些情况下并不适用(例如小颗粒可能移动速度要比大颗粒快)。读者可能会问,是否可以存在多个离散相速度?答案是肯定的。这即对应的多流体模型。

另一种更抽象的方法,是直接对动力学方程进行矩封闭,并假定速度基于颗粒大小的多项式拟定,计算不同颗粒大小的不同速度。这种方法也被称之为Eulerian-QBMM。在Eulerian-QBMM中,$\bfM_\drag$项需要进行特殊处理。

考虑方程\eqref{NS1},在Eulerian-QBMM中,不可压缩连续相的动量方程为

\begin{equation}\label{QBMM1} \frac{{\p \left( {{\alpha_\rc}{\bfU_\rc}} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rc}\left( {{\bfU_\rc} \otimes {\bfU_\rc}} \right)} \right) - \nabla \cdot \left( {{\alpha_\rc}{\bfR_\rc}} \right) =-\alpha_\rc\nabla\frac{p_\rc}{\rho_\rc} + {\alpha_\rc}\bfg - \sum_{i=0}^N\left(\frac{w_i \frac{\pi}{6}d_i^3}{\rho_\rc V_\cell} \frac{3}{4} \frac{C_{\rD,i}\rho_\rc}{d_i}|\bfU_\rc-\bfU_{d_i}|\left(\bfU_\rc-\bfU_{d_i}\right)\right). \end{equation}

其中$w_i$表示第$i$组离散相具有的权重(百分比),$d_i$表示第$i$组离散相的粒径,$C_{\rD,i}$表示第$i$组离散相的曳力系数,$\bfU_{d_i}$表示第$i$组离散相的速度。如果各组的粒径是相同的,且速度相同,即:

\begin{equation}\label{QBMM2} d_1=d_2=d_3=...=d_N,\bfU_{d_1}=\bfU_{d_2}=\bfU_{d_3}=...=\bfU_{d_N} \end{equation}

那么方程\eqref{QBMM1}和方程\eqref{NS1}是相同的。

现在对方程\eqref{QBMM1}去掉阿基米德浮力项、重力项,构建速度方程:

\begin{equation}\label{QBMM3} \frac{{\p \left( {{\alpha_\rc}{\bfU_\rc}} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rc}\left( {{\bfU_\rc} \otimes {\bfU_\rc}} \right)} \right) - \nabla \cdot \left( {{\alpha_\rc}{\bfR_\rc}} \right) =- \sum_{i=0}^N\left(\frac{w_i \frac{\pi}{6}d_i^3}{\rho_\rc V_\cell} \frac{3}{4} \frac{C_{\rD,i}\rho_\rc}{d_i}|\bfU_\rc-\bfU_{d_i}|\right)\bfU_\rc. \end{equation}

对其求解后有$\bfHbyA$。然后附加阿基米德浮力项、重力项来构建速度:

\begin{equation}\label{QBMM4} \bfU_\rc=\bfHbyA+\frac{1}{A_\mathrm{P}}\left(-\alpha_\rc\nabla\frac{p_\rc}{\rho_\rc}+\alpha_\rc\bfg+\sum_{i=0}^N\left(\frac{w_i \frac{\pi}{6}d_i^3}{\rho_\rc V_\cell} \frac{3}{4} \frac{C_{\rD,i}\rho_\rc}{d_i}|\bfU_\rc-\bfU_{d_i}|\right)\bfU_{\rd_i}\right) \end{equation}

其对应的通量为

\begin{equation}\label{QBMM5} \phi_\rc=\bfHbyA_f\cdot\bfS_f+\frac{1}{A_{\mathrm{P},f}}\left(-\alpha_{\rc,f}\nabla_f\frac{p_\rc}{\rho_\rc}+\alpha_{\rc,f}\bfg_f+\sum_{i=0}^N\left(\frac{w_i \frac{\pi}{6}d_i^3}{\rho_\rc V_\cell} \frac{3}{4} \frac{C_{\rD,i}\rho_\rc}{d_i}|\bfU_\rc-\bfU_{d_i}|\right)_f\bfU_{\rd_i,f}\right)\cdot\bfS_f \end{equation}

$\phi_\rc$可用来求解压力$p$,随后应通过方程\eqref{QBMM4}进行速度的更新。然后即可进入下一个时间步。

终端速度精确解

考虑一个水中上浮的气泡,其稳定的上浮速度为终端速度(terminal velocity)。下面,以仅考虑曳力和重力的情况下,推导终端速度的精确解。气泡的上浮中,收到的重力和曳力相互平衡,即

\begin{equation}\label{balance} \frac{m_\dpm}{\rho_\dpm}\frac{3}{4}\frac{C_\rD\rho_\rc}{d_\rd}\left|\bfU_\rc-\bfU_\dpm\right|\left(\bfU_\rc-\bfU_\dpm\right)=m_\dpm \bfg. \end{equation}

其中的未知数只有$\bfU_\dpm$,求解后即为终端速度:$\bfU_\dpm=\bfg/\mathrm{Kd}$。

参考文献

Marchisio, D. L., Fox, R. O., 2013. Computational models for polydisperse particulate and multiphase systems. Cambridge University Press.

Ramkrishna, D., 2000. Population balances: theory and applications to particulate systems in engineering. Academic press.

Sundaresan, S., Ozel, A., Kolehmainen, J., 2018. Toward constitutive models for momentum, species, and energy transport in gas–particle flows. To Appear. Annual review of chemical and biomolecular engineering.

Li, J., 1994. Particle-fluid two-phase flow: the energy-minimization multi-scale method. Metallurgical Industry Press.

Yang, N., Wang, W., Ge, W., Li, J., 2003. CFD simulation of concurrent-up gas–solid flow in circulating fluidized beds with structure-dependent drag coefficient. Chemical Engineering Journal, 96, 71-80.

Struchtrup, H., 2005. Macroscopic transport equations for rarefied gas flows. Springer.

Chalons, C., Fox, R. O., Massot, M., 2010. A multi-Gaussian quadrature method of moments for gas-particle flows in a LES framework. In Proceedings of the Summer Program, Center for Turbulence Research, Stanford University.

Chen, G. Q., Liu, H., 2003. Formation of δ-shocks and vacuum states in the vanishing pressure limit of solutions to the Euler equations for isentropic fluids. SIAM journal on mathematical analysis, 34, 925-938.

Wells, M. R., Stock, D. E., 1983. The effects of crossing trajectories on the dispersion of particles in a turbulent flow. Journal of fluid mechanics, 136, 31-62.

O'Rourke, P. J., Snider, D. M., 2012. Inclusion of collisional return-to-isotropy in the MP-PIC method. Chemical engineering science, 80, 39-54.

Kaufmann, A., Moreau, M., Simonin, O., Helie, J., 2008. Comparison between Lagrangian and mesoscopic Eulerian modelling approaches for inertial particles suspended in decaying isotropic turbulence. Journal of Computational Physics, 227, 6448-6472.

Drew, D. A., 1983. Mathematical modeling of two-phase flow. Annual Review of Fluid Mechanics, 15, 261-291.

李东岳, 2018. 如果用DPMFoam求解稀相流会怎么样?

更新历史
2018.04.27将连续相方程变为可压缩形式,增加参考文献 | 2018.04.22增补内容及图片 | 2018.04.12增加曳力、重力、浮力的讨论增加拉格朗日方程精确解 | 2018.04.11更正方程\eqref{NS1} | 2018.04.04新增附录Eulerian-QBMM | 2018.03.21经四川大学李巧列建议修订公式(5|14|18)以及并增加单位量纲 | 2018.03.16创立页面

东岳流体®版权所有
勘误、讨论、补充内容请前往CFD中文网