denseParticleFoam解析

引言

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

欧拉拉格朗日模型是可用于描述多相系统的介尺度模型。在欧拉欧拉模型中，连续相由Navier-Stokes方程描述，离散相采用拉格朗日粒子进行跟踪。界面间的动量/能量传递采用模型进行模化。相对于欧拉欧拉模型，欧拉拉格朗日模型可以处理多分散系统，并且具有明细的粒子追踪特征。若不考虑离散相分数对连续相的影响，通常被称之为单向耦合。若考虑离散相和连续相的交互作用，通常被称之为双向耦合。若进一步考虑粒子之间的交互，通常被称之为四向耦合。OpenFOAM中的欧拉拉格朗日模型求解器为denseParticleFoam（老版本中叫做DPMFoam与MPPICFoam），由于植入的模型定位，其主要用于模拟气固流动。本文主要介绍denseParticleFoam/DPMFoam求解器中的数值模型与求解细节。

控制方程与算法

在欧拉拉格朗日框架中，粒子的位置和速度可以通过求解下面的方程获得：

(1)\[\begin{equation} \label{DPM1} \frac{\mathrm{d}\mathbf{X}^p}{\mathrm{d} t}=\mathbf{U}^p , \end{equation}\]

(2)\[ m^p\frac{\mathrm{d}\mathbf{U}^p}{\mathrm{d} t}=\mathbf{F}^p , \]

其中 \(\bfX^p\) 是位置矢量，\(\bfU^p\) 是粒子的速度，\(m^p\) 是粒子质量，\(\bfF^p\) 是作用于离散粒子上的力。这里假设粒子质量是常数，并且没有因为质量传递而产生的净动量交换。对于颗粒流，由于浮力和其他非曳力型力非常小，通常是考虑曳力和重力作用。而对于气泡流，非曳力型力在某些情况下对于求解粒子的横向运动非常关键。

大部分情况下，最重要的力是由连续相作用于粒子表面的应力而产生的曳力和浮力，可以用下面的方程表示：

(3)\[\begin{equation}\label{drag} \bfF_\drag^p=\frac{1}{\tau^p}\left(\bfU_\rc-\bfU^p\right)=\frac{3}{4}\frac{m^p}{\rho^p} \frac{C_\rd}{d^p}\rho_\rc \left|\bfU_\rc-\bfU^p\right| \left(\bfU_\rc-\bfU^p\right) \end{equation}\]

(4)\[\begin{equation} \bfF_\buo^p=m^p\bfg\left(1-\frac{\rho_\rc}{\rho^p}\right), \end{equation}\]

其中 \(\tau^p\) 是曳力特征时间尺度，\(\bfU_\rc\)是定义在网格上的连续相的速度，\(\rho_\rc\) 和 \(\rho^p\) 分别是连续相和离散相的密度；\(C_\rd\) 是曳力系数，可以用下列方程计算：

(5)\[\begin{equation} C_\rd=\left\{\begin{matrix} 24\left(1+0.15\mathrm{Re}^{0.687} \right) /\mathrm{Re}, & \mathrm{Re}\leqslant 1000 \\ 0.44, & \mathrm{Re} > 1000 \end{matrix}\right. \end{equation}\]

其中 \(\mathrm{Re} \) 是粒子雷诺数。除此之外，在这里不考虑其他的力的作用。上述方程都是针对一个粒子的情况，在欧拉拉格朗日耦合框架，一个网格会存在多个粒子，因此定义相分数\(\alpha^p\)：

(6)\[\begin{equation} \alpha^p=\frac{\sum V^p}{V_{\mathrm{cell}}} \end{equation}\]

其中\(V^p\) 是粒子体积，\(V_{\mathrm{cell}}\) 是网格单元体积。将一个网格单元的粒子的力进行加和，即为与欧拉项进行交换的总的力。但在OpenFOAM中并不是采用将每个粒子的力的加和的方式来计算。OpenFOAM通过将粒子的受力进行加和后，通过方程(2)来更新\(\rd\bfU/\rd t\)，因此总网格内部的粒子的总受力为：

(7)\[ \sum \bfF = \sum m^p\frac{\mathrm{d}\mathbf{U}^p}{\mathrm{d} t} \]

在欧拉-拉格朗日方法中，连续相是在欧拉框架下计算。在欧拉框架下，纳维斯托克斯方程是统计平均计算，湍流和相间现象用封闭模型计算。在这种假设下，不可压缩连续相的质量和动量守恒可以用下列方程表示：

(8)\[ \frac{\p\alpha_\rc}{\p t} + \nabla \cdot \left( {{\alpha_\rc }{\bfU_\rc}} \right) = 0, \]

(9)\[ \frac{{\p \left( {{\alpha_\rc}{\bfU_\rc}} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rc}{{\bfU_\rc} {\bfU_\rc}} } \right) - \nabla \cdot \left( {{\alpha_\rc}{\boldsymbol{\tau_\rc}}} \right) = - \alpha_\rc\nabla \tilde{p}_\rc + \alpha_\rc\bfg - \frac{1}{\rho_\rc}\frac{ \bfF}{ V_{cell}} \]

其中\(V_{{cell}}\) 是网格单元体积，\(\tilde{p}_\rc\) 是压力除以密度的值。本研究聚焦于CFD模拟粒子流时欧拉-拉格朗日方法的数值处理，因此，忽略热传递和质量传递。可以看出，欧拉-拉格朗日方法中连续相的控制方程和欧拉-欧拉方法相似。

为了更好的有界性，denseParticleFoam实际求解的方程(2)为非守恒形式，其方程左侧可以写为：

(10)\[ \bfU_\rc\frac{\p\alpha_\rc}{\p t} + \alpha_\rc\frac{\p\bfU_\rc}{\p t} +\alpha_\rc\bfU_\rc\cdot\nabla\bfU_\rc +\bfU_\rc(\nabla\cdot(\alpha\bfU_\rc)) - \nabla \cdot \left( {{\alpha_\rc}{\boldsymbol{\tau_\rc}}} \right) \]

调用方程(3)，其可以化简为

(11)\[ \alpha_\rc\frac{\p\bfU_\rc}{\p t} +\alpha_\rc\bfU_\rc\cdot\nabla\bfU_\rc - \nabla \cdot \left( {{\alpha_\rc}{\boldsymbol{\tau_\rc}}} \right) \]

上述方程可以更好的保证有界。然而非守恒形式的方程由于不存在通量函数，上述方程的对流项无法调用隐性有限体积法离散，因此其继续可以推导为：

(12)\[\begin{split} \frac{{\p \left( {{\alpha_\rc}{\bfU_\rc}} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rc}{{\bfU_\rc} {\bfU_\rc}} } \right) - \bfU_\rc\left(\frac{\p\alpha_\rc}{\p t}+\nabla\cdot(\alpha_\rc\bfU_\rc)\right) - \nabla \cdot \left( {{\alpha_\rc}{\boldsymbol{\tau_\rc}}} \right) \\\\ = - \alpha_\rc\nabla \tilde{p}_\rc + \alpha_\rc\bfg - \frac{1}{\rho_\rc}\frac{ \bfF}{ V_{cell}} \end{split}\]

下面我们对方程(12)进行离散，在离散时，去掉方程右侧的重力项以及压力项，可写为：

(13)\[ \frac{{\p \left( {{\alpha_\rc}{\bfU_\rc}} \right)}}{{\p t}} + \nabla \cdot \left( {{\alpha_\rc}{{\bfU_\rc} {\bfU_\rc}} } \right) - \bfU_\rc\left(\frac{\p\alpha_\rc}{\p t}+\nabla\cdot(\alpha_\rc\bfU_\rc)\right) - \nabla \cdot \left( {{\alpha_\rc}{\boldsymbol{\tau_\rc}}} \right) = - \frac{1}{\rho_\rc}\frac{ \bfF}{ V_{cell}} \]

其离散形式可以写为：

(14)\[ A_{\rc,\rP}\bfU_{\rc,\rP} + \sum{A_{\rc,\mathrm{N}}\bfU_{\rc,\mathrm{N}}} = S_{\rc,\rP}, \]

求解方程(15)有预测速度\(\bfHbyA\)：

(15)\[ \bfHbyA_{\rc,\rP}=\frac{1}{A_{\rc,\rP}}\left( -\sum{A_{\rc,\mathrm{N}}\bfU_{\rc,\mathrm{N}}} + S_{\rc,\rP} \right) \]

在\(\bfHbyA\)的基础上，添加压力项、重力项的贡献：

(16)\[ \bfHbyA_{\rc,\rP}' =\bfHbyA_{\rc,\rP} + \frac{1}{A_{\rc,\rP}} \left(-\alpha_\rc\nabla \tilde{p}_\rc+\alpha_\rc\bfg\right) \]

在收敛的情况下，需满足：

(17)\[ \frac{\p\alpha_\rc}{\p t} + \nabla \cdot \left( \alpha_\rc \bfHbyA_{\rc,\rP}' \right) = 0 \]

也即：

(18)\[ \nabla \cdot \left( \frac{\alpha_\rc^2 }{A_{\rc,\rP}} \nabla \tilde{p}_\rc \right) = \frac{\p\alpha_\rc}{\p t} + \nabla \cdot \left( \frac{\alpha_\rc }{A_{\rc,\rP}}\left( \bfHbyA_{\rc,\rP} +\alpha_\rc\bfg\right) \right) \]

上述方程为最终的压力方程。

关键代码

fvVectorMatrix UcEqn
(
    fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
  - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
  + continuousPhaseTurbulence->divDevRhoReff(Uc)
 ==
    (1.0/rhoc)*cloudSU
);

其中(1.0/rhoc)*cloudSU中的cloudSU为fvScalarMatrix。方程(2)的单位为\(\mathrm{m}\cdot\mathrm{s}^2\)。离散后的fvScalarMatrix的单位发生了统一的变化。上述代码的方程左侧的fvScalarMatrix的单位为\(\mathrm{m}^3/\mathrm{s} \cdot \mathrm{m}/\mathrm{s}\)。cloudSU的单位即为力的单位，是\(\mathrm{kg}\cdot\mathrm{m}/\mathrm{s}^2\)，在数学上对应\(\sum\bfF\)。(1.0/rhoc)*cloudSU对应的数学项为\(- \frac{1}{\rho_\rc}\frac{ \sum\bfF}{ V_{cell}}\)，其看起来缺失了\(V_{cell}\)。其实不然。这是因为对于\(-\frac{ \sum\bfF}{ V_{cell}}\)这一项，其为一个场，在将其离散为fvScalarMatrix的源项的时候，需要乘以网格单元体积。因此在cloudSU中\(-\frac{ \sum\bfF}{ V_{cell}}V_{cell}\)，仅仅剩下\(\sum\bfF\)需要处理。

在代码中，cloudSU表示为：

fvm.source() = -UTrans()/(this->db().time().deltaT());

其中的-UTrans()/(this->db().time().deltaT())的数学公式为\(\sum m^p\frac{\rd\bfU^p}{\dt}\)，即为\(\sum\bfF^p\)。这与上述的分析是一致的。

代码中的cloudVolSUSu定义为：

cloudVolSUSu.primitiveFieldRef() = -cloudSU.source()/mesh.V();

可见，cloudVolSUSu的数学形式即为\(-\frac{ \sum\bfF}{ V_{cell}}\)。

另一方面，从代码来看:

fvScalarMatrix pEqn
(
    fvm::laplacian(alphacf*rAUcf, p)//非平方项
  ==
    fvc::ddt(alphac) + fvc::div(alphacf*phiHbyA)
);

方程(18)的压力梯度相的相分数不是\(\alpha_\rc^2\)而是\(\alpha_\rc\)，这是因为在方程(2)中的压力梯度项并没有乘以相分数。这是两种处理方式。在不同的工作中有不同的处理方式：

在Yuan et al.、 Deen et al.、Li et al.中均采用了乘以相分数的方式；
在这个帖子中OpenFOAM创始人Weller表示OpenFOAM的植入与一些工作存在不一致；
在我的工作Li et al.中，我对乘以相分数以及不乘以相分数的问题作了简单对比，至少从速度与相分数来看，二者差异可以忽略；