基于coalChemistryFoam的稠密气固两相求解器开发
前文链接：这里请插入上一篇的链接（基于coalChemistryFoam的稠密气固两相求解器开发)
今天和大家继续学习出自Chalmers University of Technology的OpenFOAM课程学生报告
在完成连续相体积分数的定义之后，我们需要修改原先coalChemistryFoam求解器中的湍流模型，在createFields.H文件中定义了求解器的湍流模型：
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
compressible::turbulenceModel::New
(
rho,
U,
phi,
thermo
)
);
coalChemistryFoam求解器使用的湍流模型是compressible::turbulenceModel::New，其定义在 $FOAM_SRC/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModel.H（因为在求解器主文件coalChemistryFoam.C中该头文件被包括）
namespace compressible
{
typedef ThermalDiffusivity<CompressibleTurbulenceModel<fluidThermo>> turbulenceModel;
该湍流模型是继承CompressibleTurbulenceModel，其中连续相密度仍用描述，并没有包含体积分数
的信息，所以我们需要建立一个新的，继承PhaseCompressibleTurbulenceModel的湍流模型，并使用fluidThermo作为transportModel，仿照turbulentFluidThermoModels文件夹中的文件，我们可以定义一个新的coalChemistryTurbulenceModel模型，这里作者将文件结构做了一些修改：
. // 原turbulentFluidThermoModels 目录
├── derivedFvPatchFields // 边界条件不变
├── makeTurbulenceModel.H // 不变
├── turbulentFluidThermoModel.C
├── turbulentFluidThermoModel.H
├── turbulentFluidThermoModels.C
└── turbulentFluidThermoModels.H
. // coalChemistryFoam目录
├── coalChemistryTurbulenceModelFwd.H // 可以不用
├── coalChemistryTurbulenceModel.H
├── coalChemistryTurbulenceModels.C
作者将turbulentFluidThermoModels.C 与 turbulentFluidThermoModels.H 合并为coalChemistryTurbulenceModels.C，并删掉了coalChemistryTurbulenceModel.C，至于coalChemistryTurbulenceModel.C的作用目前还是不很清楚，因为我尝试了在OF源代码中删去这一部分，编译下来也没有问题，有了解的小伙伴可以后台留言讨论。
coalChemistryTurbulenceModel.H 内容如下：
#ifndef coalChemistryTurbulenceModel_H
#define coalChemistryTurbulenceModel_H
#include "coalChemistryTurbulenceModelFwd.H" // 这个文件可以不用
#include "PhaseCompressibleTurbulenceModel.H"
#include "ThermalDiffusivity.H"
#include "fluidThermo.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
typedef ThermalDiffusivity<PhaseCompressibleTurbulenceModel<fluidThermo>>
coalChemistryTurbulenceModel;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
可以看出现在的coalChemistryTurbulenceModel是PhaseCompressibleTurbulenceModel的子类，然后作者在coalChemistryTurbulenceModels.C中定义了湍流子模型，包括各种RANS模型和LES模型。
通过对比魔改前后的***Models.C文件可以发现：
makeBaseTurbulenceModel // 魔改前
(
geometricOneField, // 返回单位 1
volScalarField,
compressibleTurbulenceModel,
CompressibleTurbulenceModel,
ThermalDiffusivity,
fluidThermo
);
makeBaseTurbulenceModel // 魔改后
(
volScalarField,
volScalarField,
compressibleTurbulenceModel,
PhaseCompressibleTurbulenceModel,
ThermalDiffusivity,
fluidThermo
);
修改后的湍流模型中的增加了一个volScalarField参数，就是连续相的体积分数场。
修改后湍流模型之后，剩下的任务就是魔改控制方程的求解了，在上一期文章中已给出了控制方程的具体表达式，下面以较为复杂的动量方程为例，分析代码的实现方法
// Uc rhoc分别为定义在creatFields.H中的连续相密度和速度
UcEqn =
(
fvm::ddt(alphac, rhoc, Uc) // 时间项
+ fvm::div(alphaRhoPhic, Uc) // 对流项
+ MRF.DDt(alphac*rhoc, Uc) // 旋转体的附加项
+ turbulence->divDevRhoReff(Uc) // 湍流模型
==
alphac*rhoc()*g // 重力
+ coalParcels.SU(Uc) // 颗粒coalParcels与气相的动量交换
+ limestoneParcels.SU(Uc) // 颗粒limestoneParcels与气相的动量交换
+ fvOptions(rhoc, Uc) // fvOptions表示附加源项，可用作点火源
);
UcEqn.relax();
fvOptions.constrain(UcEqn);
if (pimple.momentumPredictor())
{
solve(UcEqn == -alphac*fvc::grad(p)); // 压力梯度项
fvOptions.correct(Uc);
Kc = 0.5*magSqr(Uc);
}
值得注意的一点是，除了颗粒的源项外，动量方程所有的项都经过了体积分数的修正，因为作者没有提供lagrangian类的具体实现方法，所以具体颗粒源项有没有修改不得而知，但根据lagrangian类的源项计算方法
离散相作用于连续相的源项等于网格内所有颗粒产生源项的和对时间的导数，再除以网格的体积，但如果颗粒的体积分数不能忽略，那么应该除以网格当地的连续相体积分数与网格体积之积，即为
最后小结一下，这篇学生报告在coalChemistryFoam求解器的基础上，通过增加离散相体积分数使修改后的求解器能适用于稠密气固两相反应流，我们可以借鉴这种思路去二次开发诸如针对稠密多相流的求解器或者针对某些颗粒并不稠密但本身体积相比网格来说不能忽略的求解器，但这里有一些问题需要注意
- 如果颗粒体积过大，可能拉格朗日的基本假设就不再成立
- 在某一网格内如果连续相的体积分数过小，则可能会造成由于颗粒源项过大而发散的问题
- 颗粒尺寸和网格尺寸的相关性问题
关于以上问题，感兴趣的小伙伴可以参考这篇报告的作者后来发表在顶刊Chemical Engineering Journal上的文章。
