-
@hurricane007
我对DEM算法不是很熟悉，不过我在单相流和多相流中添加过单向耦合的粒子。如果你需要耦合，无非就是在DPMFoam类求解器下添加反应和能量，要么就是在反应和能量求解器中添加粒子。我看你目前的做法是基于DPMFoam添加反应和能量，是否可以换一种思路？比如在reactingFoam里面添加粒子。 -
这个帖子发出来2年多，2015年的这个时候我主要做欧拉欧拉模型，2017年年初开始做一种欧拉欧拉和欧拉拉格朗日的混合方法，最近转了纯欧拉拉格朗日，加深了一点理解。
多相系统
欧拉欧拉和欧拉拉格朗日主要用于描述多相系统，欧拉欧拉是两相均用欧拉NS方程描述，欧拉拉格朗日是一相用欧拉，一相用拉格朗日方程描述。对于两相系统，欧拉欧拉方法模型为双流体模型（Two fluid model，TFM），欧拉拉格朗日方法为Discrete Particle Model（DPM）。
从数学形式上，DPM的控制方程要比TFM中的离散相方程简单，比如这个是DPM中描述粒子的速度方程（牛顿第二定律）：
\begin{equation}\label{DPM1}
m\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}=\mathbf{F}
\end{equation}
\begin{equation}\label{DPM2}
\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\mathbf{U}
\end{equation}
这个是描述粒子的欧拉方程（2018.01.04更新：经石油大学某同学建议改正符号错误）：
\begin{equation}\label{TFM1}
\frac{\partial \alpha\mathbf{U}}{\partial t}+\nabla\cdot\left(\alpha\mathbf{U}\mathbf{U} \right )-\nabla\cdot\tau=-\alpha\nabla p+\mathbf{g}+\mathbf{M}
\end{equation}
\begin{equation}\label{TFM2}
\frac{\partial \alpha}{\partial t}+\nabla\cdot\left(\alpha\mathbf{U} \right )=0
\end{equation}
很明显，拉格朗日方程是一个单纯的ODE，在将方程右边某些力简化的情况下，这个ODE是存在解析解的。欧拉方程需要数值求解，并且耗散性比较严重。如果只是单纯的求解方程\eqref{DPM1}和\eqref{DPM2}，你都不需要网格，直接求解就可以，这只是单纯的ODE。
个人感觉TFM要比DPM的求解数值上问题比较多，比如比较出名的相分数有界、高梯度区域力的模化、相分离导致的奇异、高阶格式震荡等等等等。DPM都不存在这些问题。
但是如果看DPM的算法，如何跟踪粒子的位置信息占据了大头。目前的算法，即使将网格拆分成为四面体，依然会发生粒子丢失的现象，在某些网格较差的情况下，会陷入死循环。
求解器不发散，就是一直loop :big_mouth:
DPM，DEM，MPPIC
在考虑多相流的情况下，一般颗粒相用方程\eqref{DPM1}和\eqref{DPM2}来描述，连续相用NS方程来描述，也就是同时求解\eqref{DPM1}，\eqref{DPM2}，\eqref{TFM1}和\eqref{TFM2}。
同时求解这四个方程最重要的问题就是效率问题。
-
求解不需要网格的ODE方程需要高效的处理粒子跟踪位置问题，这个问题处理不好很容易导致粒子算着算着丢失了。一个一个粒子去找当然可以，但是如果有10e7个粒子肯定没办法处理。在OpenFOAM中，通过将网格四面体化来处理，由于商软代码不开源，不知道商软如何处理粒子位置跟踪的效率问题。
-
粒子间作用力的处理，在这里，如果考虑OpenFOAM中的DPM模型（商软中的DPM不清楚是否也是这样），方程\eqref{DPM1}中右边的力主要包含了曳力、虚拟质量力等非常主要的力。但是在DEM中，除了上述的力，还有旋转、摩擦、凝聚等，比较重要的颗粒碰撞也通过软球模型来处理（下图是google随便找的）。
所以说，DEM 处理的力更加全面，并且在DEM中对于颗粒的形状也会进行考虑，形状会进一步影响力的作用。 -
MPPIC是为了更加的追求高效的一套算法。在MPPIC中，多个particle被打包到一起并被称之为parcel，在计算parcel运动的同时，插值出欧拉场。欧拉场的作用在于在MPPIC的控制方程中添加一个固相应力。在TFM中，也存在这样一个固相应力，因为TFM是欧拉框架下的。也就是说，MPPIC中，碰撞模型通过欧拉场的固相应力进行表征。为什么这么做？大量的节省用于计算颗粒碰撞导致特别小的时间步长。这在OpenFOAM的代码中很直接的表现出来：
#ifdef MPPIC #include "basicKinematicMPPICCloud.H" #define basicKinematicTypeCloud basicKinematicMPPICCloud #else #include "basicKinematicCollidingCloud.H" #define basicKinematicTypeCloud basicKinematicCollidingCloud #endif如果是MPPIC方法，调用MPPIC的碰撞模型（
basicKinematicMPPICCloud.H），如果是DPM，调用碰撞模型（basicKinematicCollidingCloud.H）。
Misleading
另外好像现在各大厂里面的各种算法虽然是一个名字，但是使用的算法不同，比如按照 @散漫守望2016 所说：
Fluent中的DPM模型只能处理稀相的气固（气液）两相流，一般我们称它为“双向耦合”，颗粒颗粒间不发生碰撞以及在流体方程中不考虑空隙率。
但是在OpenFOAM中的
DPMFoam也可以处理碰撞以及体积分数耦合问题。同时DPMFoam也可以使用DEM中的软球模型，虽然他被叫成了DPMFoam，并且DPMFoam也可以调用MPPIC中parcel的概念。据我所知，parcel这个概念最早出现在Andrews and O'Rourke 1995年的文章中。同时，MPPIC在ANSYS Fluent中被称之为DDPM，在称之为CPFD。
最近我在用欧拉拉格朗日和精确解对比然后模拟3D鼓泡床 :cheeky:
-
-
东岳兄这个回复很详细，我有一些补充。
- 我感觉DPM的控制方程虽然是个ODE，但是在运算的过程中其实 F 是在变的，理论上说解析解存在，但是在一个实际系统里面由于F在变，所以其实并没有意义。
- 跟踪粒子的位置信息占了大头，其实应该说是碰撞探测占了大头，在MPPIC里面也跟踪粒子，但是速度快得多；
- MPPIC 高效在用一个 solid stress 来模拟碰撞，因此不需要确定颗粒碰撞了，这一点省了很多时间，对于parcel 用一个点来跟踪多个颗粒，在DPM也有做，叫做coarsed grain DPM； 但是做parcel都需要一个子模型来模拟这个parcel的行为，我测试下来，对于一个 N 颗粒系统，如果一个parcel包含nParcel 个颗粒，那么 计算时间的比值大概是 T(single particle) / T(parcel) = (nParcel)^1/3 。
-
-
@马乔 MPPIC 精髓还是在于那个solid stress 来模拟particle-particle force，而DPM/DEM 要直接计算颗粒间作用力。我觉得MP-PIC 在某种程度上是 "TFM DPM 混合"的方法，TFM来算 particle particle force，DPM 来追踪颗粒。而 Coarsed grain DPM 最初的概念是 particle-particle force 用大颗粒算，gas-particle force 用小颗粒算，似乎和PIC也没有关系啊。但是到后面发现 1. 一堆小颗粒在一起应该会“变软”；2. 转动惯量（也就是切向力）有影响，还有一些不记得了，所以后面就加了一些修正了，当然可能后面也有用PIC的思想来修正的
-
关于在 DPMFoam 添加化学反应，CFD Online 上有一个关于 coalCollidingChemistryFoam 的讨论（https://www.cfd-online.com/Forums/openfoam-programming-development/158458-colliding-coal-cloud-coalcollidingchemistryfoam.html ）， 即在煤燃烧的算例中添加 DPMFoam 中碰撞项的处理。在该帖子中提供的求解器可以在 OpenFOAM 2.3.x 版本中编译并运行。
但是对该模型中添加 submodel 的部分不太理解，直接将该方法用于向 MPPICFoam 中添加化学反应则一直失败。不知道有没有人有类似需求的，是否可以讨论下？
-
我是直接根据上面提到的帖子修改的。在向 OpenFOAM 自带的 coalCombustion 库添加碰撞时，上述帖子中直接在
coalCloud.H的定义中将typedef ReactingMultiphaseCloud<ReactingCloud<ThermoCloud<KinemaitcCloud<Cloud<coalParcel>>>>> coalCloud;修改为
typedef ReactingMultiphaseCloud<ReactingCloud<ThermoCloud<CollidingCloud<KinemaitcCloud<Cloud<coalParcel>>>>>> coalCloud;对颗粒的定义也采用了类似的方法。
我将该例子中的CollidingCloud均改为MPPICCloud，但是编译无法通过。错误出在 MPPICCloud.C 中td.updateAverages(*this);这一行，编译器提示这里的updateAverages找不到定义，应该是按上面方法修改coalCloud和coalParcel定义后，*this所代表的具体定义跟函数里参数形式不一样导致的。我现在在想上面的思路是不是可行，如果说需要在
thermoparcel中修改代码，加判断的话，是不是意味着我不应该采取上述的方法，而应自己将MPPICParcel和thermoparcel以及反应相关的部分重新写成一个新的颗粒类？非常感谢大家的回复，有必要的话我会想办法把编译错误贴出来。
