版本对应:

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

CFD:多相+不可压+EL算法

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

欧拉拉格朗日模型是可用于描述多相系统的介尺度模型。在欧拉欧拉模型中,连续相由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.中,我对乘以相分数以及不乘以相分数的问题作了简单对比,至少从速度与相分数来看,二者差异可以忽略;