版本对应：

OpenFOAM-11中的incompressibleDenseParticleFluid模块，对应OpenFOAM-10之前的denseParticleFoam，以及更老版本中的DPMFoam。

CFD: 多相欧拉拉格朗日算法

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

欧拉拉格朗日模型是可用于描述多相系统的介尺度模型。在欧拉欧拉模型中，连续相由Navier-Stokes方程描述，离散相采用拉格朗日粒子进行跟踪。界面间的动量/能量传递采用模型进行模化。相对于欧拉欧拉模型，欧拉拉格朗日模型可以处理多分散系统，并且具有明细的粒子追踪特征。若不考虑离散相分数对连续相的影响，通常被称之为单向耦合。若考虑离散相和连续相的交互作用，通常被称之为双向耦合。若进一步考虑粒子之间的交互，通常被称之为四向耦合。

欧拉拉格朗日算法并没有某一篇sci进行明确的讨论（像PISO、SIMPLE那样）。因为其不需要特殊的算法来处理。只要欧拉相可以求解，拉格朗日相可以求解，耦合策略为一种普适性策略。

控制方程与算法

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

(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=m^p\frac{\bfU_\rc-\bfU^p}{\tau^p}=\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}}\) 是网格单元体积。将一个网格单元的粒子的力进行加和，即为与欧拉项进行交换的总的力。

在欧拉-拉格朗日方法中，连续相是在欧拉框架下计算。在这种假设下，不可压缩连续相的质量和动量守恒可以用下列方程表示：

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

(8)\[ \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{ \sum \bfF^p}{ V_{cell}} \]

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

注意：

源项\(- \frac{1}{\rho_\rc}\frac{ \sum \bfF^p}{ V_{cell}}\)在离散过程中，方程左右两侧需要乘以网格单元体积，因此再fvVectorMatrix层面，源项变为\(- \frac{1}{\rho_\rc} \sum \bfF^p\)

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

(9)\[ \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)，其可以化简为

(10)\[ \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) \]

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

(11)\[\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{\sum \bfF^p}{ V_{cell}} \end{split}\]

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

(12)\[ \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{ \sum\bfF^p}{ V_{cell}} \]

其离散形式可以写为：

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

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

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

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

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

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

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

也即：

(17)\[ \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) \]

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

最终，我们如果考虑其他作用力的情况下，力源项中的\(\sum \bfF^p\)可以写为：

(18)\[ \sum \bfF^p = Sp(\bfU-\bfU^p) +Su \]

反过来，如果仅仅存在曳力的情况下，即为：

(19)\[ Sp(\bfU-\bfU^p) +Su=\sum m^p\frac{\bfU_\rc-\bfU^p}{\tau^p} \]

也即

(20)\[ Sp=\sum \frac{m^p}{\tau^p}, Su=0 \]

OpenFOAM中定义\(\mathrm{Ucoeff}\)为

(21)\[ \mathrm{Ucoeff}=Sp\dt \]

同时定义\(\mathrm{UTrans}\)为

(22)\[ \mathrm{UTrans}=\dt\sum \bfF^p \]

关键代码

欧拉拉格朗日耦合的速度耦合至关重要。在OpenFOAM中，可以调用2种算法：如果耦合项采用完全的显性算法，即在算例层面semiImplicit设置为0，有代码：

fvVectorMatrix cloudSU(kinematicCloud.SU(Uc));

上述代码中的cloudSU表示为：

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

如方程(22)所示，其中的-UTrans()/(this->db().time().deltaT())的表达式为\(\sum\bfF^p\)。同时，cloudSU其仅仅源项存在值（对角线系数为0）。随后，

volVectorField cloudVolSUSu
(
    IOobject
    (
        "cloudVolSUSu",
        runTime.timeName(),
        mesh
    ),
    mesh,
    dimensionedVector
    (
        "0",
        cloudSU.dimensions()/dimVolume,
        Zero
    ),
    zeroGradientFvPatchVectorField::typeName
);

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

cloudVolSUSu.correctBoundaryConditions();

cloudVolSUSu被赋值为\( -\frac{\sum\bfF^p}{V_{cell}}\)。

cloudSU.source() = cloudSU.diag()*Uc();

cloudSU源相置为0，到此为止，cloudSU为一个空的fvVectorMatrix。然后：

surfaceScalarField phicForces
(
   fvc::flux(rAUc*cloudVolSUSu/rhoc) + rAUcf*(g & mesh.Sf())
);

cloudVolSUSu的值通过phicForces进入到了压力方程中。

如果耦合项采用完全的半隐性算法，即在算例层面semiImplicit设置为1，有代码：

const volScalarField::Internal
    Vdt(mesh_.V()*this->db().time().deltaT());

return UTrans()/Vdt - fvm::Sp(UCoeff()/Vdt, U) + UCoeff()/Vdt*U;

首先-fvm::Sp(UCoeff()/Vdt, U)返回矩阵对角线系数：

(23)\[\begin{split} \left(\begin{matrix} -\frac{Sp\dt\dV}{\dV \dt} & .. & ..\\ .. & .. &.. \\ .. & .. & .. \end{matrix}\right) \end{split}\]

UTrans()/Vdt进入矩阵源项，UCoeff()/Vdt*U也进入矩阵源项。其中UCoeff()/Vdt*U与-fvm::Sp(UCoeff()/Vdt, U)相互抵消，只不过一个隐性离散，一个显性离散。这也即为半隐性离散的意义。因此，起作用的即为UTrans()/Vdt。在进入源项的时候，其要乘以网格单位（伪代码为UTrans()/Vdt*V），因此其数学形式为\(\sum \bfF^p\)，也即总受力。可以看出半隐性离散主要是通过在对角线上增加一项，同时在源相上抵消其作用，来增加对角占有特性。

volVectorField cloudVolSUSu
(
    IOobject
    (
        "cloudVolSUSu",
        runTime.timeName(),
        mesh
    ),
    mesh,
    dimensionedVector
    (
        "0",
        cloudSU.dimensions()/dimVolume,
        Zero
    ),
    zeroGradientFvPatchVectorField::typeName
);

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

cloudVolSUSu.correctBoundaryConditions();

随后，cloudVolSUSu同样被赋值为\( -\frac{\sum\bfF^p}{V_{cell}}\)。

cloudSU.source() = cloudSU.diag()*Uc();

cloudSU源相置为\(Sp\bfU_c\)，然后在速度方程中：

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分为对角线贡献与源项贡献，对角线贡献为\(\frac{1}{\rho_\rc}Sp\bfU\)，源项贡献为\(\frac{1}{\rho_\rc}\sum\bfF^p+\frac{1}{\rho_\rc}Sp\bfU\)，其中\(\frac{1}{\rho_\rc}Sp\bfU\)的对角线贡献与源项贡献抵消。

surfaceScalarField phicForces
(
   fvc::flux(rAUc*cloudVolSUSu/rhoc) + rAUcf*(g & mesh.Sf())
);

cloudVolSUSu的值再一次的通过phicForces进入到了压力方程中。

代码层面另外一个关键点是:

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

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

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