在OpenFOAM中添加新的曳力模型
参考链接：在OpenFOAM中添加新的曳力_姜蜉蝣的博客-CSDN博客
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曳力模型在算得曳力的基础上乘以一个系数HD。ErgunWenYu曳力的求解在程序中是这样实现的：
1 template<class CloudType> 2 Foam::forceSuSp Foam::ErgunWenYuDragForce<CloudType>::calcCoupled 3 ( 4 const typename CloudType::parcelType& p, 5 const scalar dt, 6 const scalar mass, 7 const scalar Re, 8 const scalar muc 9 ) const 10 { 11 scalar alphac(alphac_[p.cell()]); 12 13 if (alphac < 0.8) 14 { 15 return forceSuSp 16 ( 17 vector::zero, 18 (mass/p.rho()) 19 *(150.0*(1.0 - alphac)/alphac + 1.75*Re)*muc/(alphac*sqr(p.d())) 20 ); 21 } 22 else 23 { 24 return forceSuSp 25 ( 26 vector::zero, 27 (mass/p.rho()) 28 *0.75*CdRe(alphac*Re)*muc*pow(alphac, -2.65)/(alphac*sqr(p.d())) 29 ); 30 } 31 }
可以发现算法非常明确，与文献中的曳力公式是完全对应的。在修改此曳力时，只需要在这个成员函数内部计算出EMMS曳力模型中的非均匀因子Hd，并将其乘入ErgunWenYu曳力公式即可。修改方式如下：
1 template<class CloudType> 2 Foam::forceSuSp Foam::EMMSDragForce<CloudType>::calcCoupled 3 ( 4 const typename CloudType::parcelType& p, 5 const scalar dt, 6 const scalar mass, 7 const scalar Re, 8 const scalar muc 9 ) const 10 { 11 scalar alphac(alphac_[p.cell()]); 12 13 // Calculate Hd. Hd is a function of 'Re' and alphac. 14 // scalar Hd = func(Re, alphac); 15 16 if (alphac < 0.8) 17 { 18 return forceSuSp 19 ( 20 vector::zero, 21 (mass/p.rho()) 22 *(150.0*(1.0 - alphac)/alphac + 1.75*Re)*muc/(alphac*sqr(p.d())) 23 *Hd // Add Hd !!! 24 ); 25 } 26 else 27 { 28 return forceSuSp 29 ( 30 vector::zero, 31 (mass/p.rho()) 32 *0.75*CdRe(alphac*Re)*muc*pow(alphac, -2.65)/(alphac*sqr(p.d())) 33 *Hd // Add Hd !!! 34 ); 35 } 36 }
修改构建曳力配置文件
打开文件 src/lagrangian/intermediate/parcels/include/makeParcelForces.H ，
向其头文件中添加新增曳力模型的头文件，即 #include "EMMSDragForce.H" ，
并在其宏部分添加 makeParticleForceModelType(EMMSDragForce, CloudType); \ 这样构建曳力的配置项。
======================================================
=======================注意一下========================
======================================================
在OpenFOAM-7中，发现后面只需要在 src/lagrangian/intermediate 路径下，
先 wclean 之后再 wmake ，即可在DPMFoam求解中直接使用编译的曳力模型，
无需后面的繁琐内容
======================================================
======================================================
重新编译曳力
首先要将新添加的曳力模型头文件及源文件链接到 src/lagrangian/intermediate/lnInclude 文件夹下。进入该文件夹，并执行如下命令：
1 ln -s ../submodels/Kinematic/ParticleForces/Drag/EMMSDrag/* .
之后，在 src/lagrangian/intermediate 重新编译相关的库文件即可:
1 wclean 2 wmake
重新编译求解器
找到MPPICFoam求解器所在的位置，并重新编译该求解器：
1 cd $FOAM_APP/solvers/lagrangian/DPMFoam 2 ./Allwclean 3 ./Allwmake
此时，MPPICFoam（以及DPMFoam）求解器就已经可以使用EMMS曳力模型了。
