准备工作

拷贝求解器源文件：cp -r $FOAM_APP/solvers/multiphase/twoPhaseEulerFoam $WM_PROJECT_USER_DIR/twoPhaseEulerFoam-1
拷贝算例：cp -r $FOAM_TUTORIALS/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed $FOAM_RUN/fluidisedBed
生成算例网格：blockMesh -case $FOAM_RUN/fluidisedBed
修改求解器make文件$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/Make/file:

EXE = ($FOAM_RUN)/twoPhaseEulerFoam

运行$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/Allwmake，编译求解器。

新建曳力模型

进入目录cd $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/interfacialModels/dragModels，复制其中一个曳力模型文件夹并命名为EMMScp -r WenYu/ EMMS。然后进入到EMMS文件夹下，把文件名和文件内部涉及到模型名称的地方全部替换成EMMS。
打开vim $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/interfacialModels/Make/files，向其中添加一行dragModels/EMMS/EMMS.C，并将生成目标进行修改：

//- 修改前
LIB = $(FOAM_LIBBIN)/libcompressibleEulerianInterfacialModels
//- 修改后
LIB = $(FOAM_USER_LIBBIN)/libcompressibleEulerianInterfacialModels_hemm

在目录$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/interfacialModels下运行编译生成链接库，运行wmake。生成的链接库可以在目录$FOAM_USER_LIBBIN下：

[mmhe@K226 interfacialModels]$ ll $FOAM_USER_LIBBIN
总用量 1.9M
drwxr-xr-x 2 mmhe ipe   61 10月  9 2019 .
drwxr-xr-x 3 mmhe ipe   24 9月  27 15:30 ..
-rwxr-xr-x 1 mmhe ipe 1.9M 10月  9 2019 libcompressibleEulerianInterfacialModels_hemm.so

修改vim $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/Make/options，在其中将默认的链接库位置换成刚才生成的链接库：

//- 修改前
-lcompressibleEulerianInterfacialModels
//- 修改后
-L$(FOAM_USER_LIBBIN) \
-lcompressibleEulerianInterfacialModels_hemm

运行$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/Allwmake，编译求解器。修改了求解器的Make/options文件后，必须要重新编译，否则无法更新求解器。

曳力模型算例测试

修改vim $FOAM_RUN/fluidisedBed/constant/phaseProperties，将曳力模型设置成EMMS：

drag
(
    (particles in air)
    {
        type            EMMS;
        residualAlpha   1e-6;
        residualRe      1e-3;
        swarmCorrection
        {
            type        none;
        }
    }
);

修改计算时间vim $FOAM_RUN/fluidisedBed/system/controlDict，使算例向前运行一步即可：

startTime       0;
stopAt          endTime;
endTime         2e-4;
deltaT          2e-4;

运行算例：$FOAM_RUN/twoPhaseEulerFoam -case $FOAM_RUN/fluidisedBed，正常运行：

Calculating face flux field phi.air
Selecting diameterModel for phase air: constant
Selecting turbulence model type laminar
Selecting default blending method: none
Selecting dragModel for (particles in air): EMMS
Selecting swarmCorrection for (particles in air): none
Selecting virtualMassModel for (particles in air): constantCoefficient
Selecting heatTransferModel for (particles in air): RanzMarshall
Calculating field DDtU1 and DDtU2

新建固相应力模型

拷贝一个固相应力模型：cp -r $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/EMMSStressModel。
将所有模型关键字phasePressure替换为EMMSStress。
修改vim $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phaseCompressibleTurbulenceModels.C，在最后加上EMMSStressModel模型的相关信息：

#include "EMMSStressModel.H"
#include "EMMSStressModel.H"
makeTurbulenceModel
(phaseModelPhaseCompressibleTurbulenceModel, RAS, EMMSStressModel);

在$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/Make/files文件中添加一行EMMSStressModel/EMMSStressModel.C，并将目标链接库生成配置改为：

//- 修改前
LIB = $(FOAM_LIBBIN)/libphaseCompressibleTurbulenceModels
//- 修改后
LIB = $(FOAM_USER_LIBBIN)/libphaseCompressibleTurbulenceModels_hemm

在目录$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels下运行wmake编译链接库。
修改vim $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/Make/options，在其中将默认的链接库位置换成刚才生成的链接库：

-lphaseCompressibleTurbulenceModels_hemm

运行$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/Allwmake，编译求解器，成功！

应力模型算例测试

修改算例配置文件vim $FOAM_RUN/fluidisedBed/constant/turbulenceProperties.particles，将RAS模型设置成EMMSStress，然后将phasePressure模型的输入参数拷贝成EMMSStress模型的输入参数：

RAS
{
    RASModel EMMSStress;
    EMMSStressCoeffs
    {
        preAlphaExp     500;
        expMax          1000;
        alphaMax        0.62;
        g0              1000;
    }
}

运行算例：$FOAM_RUN/twoPhaseEulerFoam -case $FOAM_RUN/fluidisedBed，正常运行：

Calculating face flux field phi.particles
Selecting diameterModel for phase particles: constant
Selecting turbulence model type RAS
Selecting RAS turbulence model EMMSStress
EMMSStressCoeffs
{
    preAlphaExp     500;
    expMax          1000;
    alphaMax        0.62;
    g0              1000;
}
Calculating face flux field phi.air
Selecting diameterModel for phase air: constant
Selecting turbulence model type laminar
Selecting default blending method: none
Selecting dragModel for (particles in air): EMMS
Selecting swarmCorrection for (particles in air): none
Selecting virtualMassModel for (particles in air): constantCoefficient
Selecting heatTransferModel for (particles in air): RanzMarshall
Calculating field DDtU1 and DDtU2

修改曳力模型

打开文件vim $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/interfacialModels/dragModels/EMMS/EMMS.C，只需要修改其中的函数Foam::tmp<Foam::volScalarField> Foam::dragModels::EMMS::CdRe() const即可。
CdRe()函数返回有效曳力系数与雷诺数的乘积 $C_d\times Re$ ，其中的变量CdsRe表示标准曳力系数与雷诺数的乘积 $C_{d0}\times Re$ 。

Foam::tmp<Foam::volScalarField> Foam::dragModels::EMMS::CdRe() const
{
    Info <<"EMMS曳力模型的CdRe函数！！！"<<endl;
    volScalarField alpha2(max(scalar(1) - pair_.dispersed(), residualAlpha_));
    volScalarField Re(pair_.Re());
    volScalarField CdsRe
    (
        neg(Re - 1000)*24.0*(1.0 + 0.15*pow(Re, 0.687))
      + pos(Re - 1000)*0.44*max(Re, residualRe_)
    );
    //- 获取EMMS模型输入参数
    const phaseModel& dispersed = pair_.dispersed();//离散相
    const phaseModel& continuous = pair_.continuous();//连续相
    const volVectorField Up = dispersed.U();//颗粒速度
    const volVectorField Ug = continuous.U();//流体速度
    const volScalarField Ur = Up.component(vector::Z)-Ug.component(vector::Z);//相对速度
    const volScalarField epsl_g = pair_.continuous();//平均空隙率

    //- 非均匀结构因子
    volScalarField Hd = Ur*epsl_g;//拟合关联式
    Hd.dimensions().reset(dimensionSet(0,0,0,0,0));//设置单位为无量纲数
    

    return
        Hd
       *CdsRe
       *pow(alpha2, -2.65)
       *max(pair_.continuous(), residualAlpha_);
}

在目录$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/interfacialModels下运行编译生成链接库，运行wmake。
运行算例，从输出信息中可以确定已经能够执行该函数，但是因为不收敛而提前终止。

修改应力模型

devRhoReff()函数：
$\begin{equation} -\mathbf{\tau}=-\mu(\nabla\mathbf{U}+\nabla\mathbf{U}^T)-(\rho\lambda-\frac23\mu)(\nabla\cdot\mathbf{U})\mathbf{I} \end{equation}$
divDevRhoReff()函数：
$\begin{equation} -\nabla\cdot\mathbf{\tau}=-\nabla\cdot(\mu\nabla\mathbf{U})-\nabla\cdot(\mu\nabla\mathbf{U}^T)-\nabla(\rho\lambda-\frac23\mu)(\nabla\cdot\mathbf{U}) \end{equation}$
在phasePressure模型中，是不考虑固相应力的，因此这两个函数返回的都是0值。

第一步要改的是四个函数。首先是pPrime函数：

修改vim $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/EMMSStressModel/EMMSStressModel.H：

//- 删除私有变量-Kenetic Theory Model coefficients
/*
        // Kinetic Theory Model coefficients

            //- Maximum packing phase-fraction
            scalar alphaMax_;

            //- Pre-exponential factor
            scalar preAlphaExp_;

            //- Maximum limit of the exponential
            scalar expMax_;

            //- g0
            dimensionedScalar g0_;
*/

修改vim $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/EMMSStressModel/EMMSStressModel.C：

//- 构造函数-删除没有定义的变量初始化
//- read()函数-删除文件数据读取
//- pPrime()和pPrimef()函数-替换掉返回表达式为拟合关联式
//- devRhoReff()函数-复制kineticTheoryModels中对应的部分到此处
//- divDevRhoReff()函数-复制kineticTheoryModels中对应的部分到此处

编译链接库，成功！运行算例，报错！

PIMPLE: iteration 1
new cannot satisfy memory request.
This does not necessarily mean you have run out of virtual memory.
It could be due to a stack violation caused by e.g. bad use of pointers or an out of date shared library

经过检查，是我定义Ur的时候触动的错误。问题解决，应该是我使用了返回的引用，因此导致了运行时发生段错误。

OpenFOAM-twoPhaseEulerFoam修改