在OpenFOAM中添加新的曳力

OpenFOAM中使用欧拉－拉格朗日方法对气固两相流进行模拟时，需要指定颗粒的曳力模型。有时候我们通过参考文献或者自行推导得到的曳力模型并不在OpenFOAM的默认曳力模型列表中，这时候就需要修改OpenFOAM的源代码，将新的曳力模型添加进去。

这里以EMMS曳力模型为例，简要介绍一下在OpenFOAM中添加新曳力模型的过程。

建立新曳力的类

对于OpenFOAM-2.3.x版本而言，颗粒曳力模型的存放位置为：

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag

可以看到，此文件夹下存在五个曳力模型，包括：

• ErgunWenYuDrag
• NonSphereDrag
• PlessisMasliyahDrag
• SphereDrag
• WenYuDrag

单独进入某个曳力模型的文件夹，例如ErgunWenYuDrag，发现其中包含了两个文件：ErgunWenYuDragForce.H与ErgunWenYuDragForce.C。这两个文件即曳力模型的声明与定义。

在添加EMMS曳力模型时，为了简单起见，可以将ErgunWenYuDrag文件夹复制出一份，并重命名为EMMSDrag，并将此文件夹中包含的ErgunWenYuDragForce.H与ErgunWenYuDragForce.C分别重命名为EMMSDragForce.H与EMMSDragForce.C。打开这两个文件，将其中类名及构造函数等中包含的“ErgunWenYu”字符串替换为“EMMS”字符串即可。

修改曳力模型

根据新的曳力公式，修改曳力模型中的曳力表达。

EMMS曳力模型实际上相当于ErgunWenYu曳力模型在算得曳力的基础上乘以一个系数$H_D$。ErgunWenYu曳力的求解在程序中是这样实现的：

template<class CloudType>
Foam::forceSuSp Foam::ErgunWenYuDragForce<CloudType>::calcCoupled
(
    const typename CloudType::parcelType& p,
    const scalar dt,
    const scalar mass,
    const scalar Re,
    const scalar muc
) const
{
    scalar alphac(alphac_[p.cell()]);

    if (alphac < 0.8)
    {
        return forceSuSp
        (
            vector::zero,
            (mass/p.rho())
           *(150.0*(1.0 - alphac)/alphac + 1.75*Re)*muc/(alphac*sqr(p.d()))
        );
    }
    else
    {
        return forceSuSp
        (
            vector::zero,
            (mass/p.rho())
           *0.75*CdRe(alphac*Re)*muc*pow(alphac, -2.65)/(alphac*sqr(p.d()))
        );
    }
}

可以发现算法非常明确，与文献中的曳力公式是完全对应的。在修改此曳力时，只需要在这个成员函数内部计算出EMMS曳力模型中的非均匀因子Hd，并将其乘入ErgunWenYu曳力公式即可。修改方式如下：

template<class CloudType>
Foam::forceSuSp Foam::EMMSDragForce<CloudType>::calcCoupled
(
    const typename CloudType::parcelType& p,
    const scalar dt,
    const scalar mass,
    const scalar Re,
    const scalar muc
) const
{
    scalar alphac(alphac_[p.cell()]);

    // Calculate Hd. Hd is a function of 'Re' and alphac.
    // scalar Hd = func(Re, alphac);

    if (alphac < 0.8)
    {
        return forceSuSp
        (
            vector::zero,
            (mass/p.rho())
           *(150.0*(1.0 - alphac)/alphac + 1.75*Re)*muc/(alphac*sqr(p.d()))
           *Hd  // Add Hd !!!
        );
    }
    else
    {
        return forceSuSp
        (
            vector::zero,
            (mass/p.rho())
           *0.75*CdRe(alphac*Re)*muc*pow(alphac, -2.65)/(alphac*sqr(p.d()))
           *Hd  // Add Hd !!!
        );
    }
}

修改构建曳力配置文件

打开文件src/lagrangian/intermediate/parcels/include/makeParcelForces.H，向其头文件中添加新增曳力模型的头文件，即#include "EMMSDragForce.H"，并在其宏部分添加makeParticleForceModelType(EMMSDragForce, CloudType); \这样构建曳力的配置项。

重新编译曳力

首先要将新添加的曳力模型头文件及源文件链接到src/intermediate/lnInclude文件夹下。进入该文件夹，并执行如下命令：

ln -s ../submodels/Kinematic/ParticleForces/Drag/EMMSDrag/* .

之后，在src/intermediate重新编译相关的库文件即可:

wclean
wmake

重新编译求解器

找到MPPICFoam求解器所在的位置，并重新编译该求解器：

cd $FOAM_APP/solvers/lagrangian/DPMFoam
./Allwclean
./Allwmake

此时，MPPICFoam（以及DPMFoam）求解器就已经可以使用EMMS曳力模型了。