多相流与矩方法的数学模型

[李东岳]

首先，多相流模拟相对于单向流模拟一直都不是很成熟。模拟多相流的小伙伴们，这是个好事，这表明，咱们有的玩。要是很成熟了，也基本没什么玩的了。多相流再掺杂湍流，那更是不好搞。不成熟的方面主要有：

• 模型不成熟，比如各种曳力模型，升力模型，破碎模型等；
• 求解算法不成熟，比如在CFD求解多相流模型中如何保证相分数严格的小于1大于0；欧拉拉格朗日那面需要颗粒大小小于网格大小；
• 体系复杂性，气液、液液、气固体系的模型不同，即使气固也分为稠密和稀疏；
• 其他等

由于CFD模型在多相流这面远不成熟，在多相流比较常见的工业（如石油化工、船舶水动力、食品、核工业等），CFD应用的还不如单向流工业（如汽车飞机、建筑、传热等）广泛。这是个坏事，这表示就业机会没有单向流CFD行业多。有好有坏。

多相流模型

多相流中区分为多个相。比如空气是一相、液滴是一相。也可以称之为气相、液相或者固相。参考上图。多相模拟主要分为界面类模型（如VOF、front-tracking）和高相分数多相流模型。二者的区别非常简单：

• 如果你需要气泡的形状，波浪的表面，你就需要捕获气泡或者波浪的界面，因此你需要使用界面跟踪类模型。
  

  比如上面这个图，就是界面重构出来的气泡界面，明显你可以看出来气泡的形状。

• 如果你的泡泡特别多，并且你不在乎气泡的形状是什么样的，你只在乎我的水里面有多少空气（气含率），那么这种情况适用于高相分数多相流模型。其不能预测出气泡的形状。

  

  上面这个图是用高相分数多相流模型模拟的，很明显，看不出来气泡的形状。这两大模型，前者属于多相流中的直接模拟，后者模化多一些。正是因为前者偏向直接模拟，因此在使用VOF的时候，需要巨量的网格。

目前界面捕获类模拟和高相分数模拟基本属于两个研究方向，各自有各自的研究热点和方向，基本上二者不重合。本文主要讨论第二类，高相分数模型。

高相分数多相流模型

在高相分数多相流模型中，一般区分为连续相和离散相（同样参考第一个图）。举例：

1. 水里吹气泡（鼓泡床），气泡是离散相，水是连续相；
2. 沙子里吹气（流化床），空气是连续相，沙子是离散相；
3. 水里面加油（液液分散），水是连续相，油是分散相；

连续相模型

连续相模型一般都使用NS方程来描述

\begin{equation}
\frac{\partial \alpha_\mathrm{c}\rho_\mathrm{c}}{\partial t}+\nabla\cdot\left(\alpha_\mathrm{c}\rho_\mathrm{c} \mathbf{U}_\mathrm{c}\right)=0
\end{equation}

\begin{equation}\label{NS}
\frac{{\partial \left({{\alpha_\mathrm{c}}{\rho_\mathrm{c}}{\mathbf{U}_ \mathrm{c}}} \right)}}{{\partial t}} + \nabla \cdot \left( {{\alpha_\mathrm{c}}{\rho_\mathrm{c} }\left( {{\mathbf{U}_ \mathrm{c}} \otimes {\mathbf{U}_ \mathrm{c}}} \right)} \right) - \ \nabla \cdot \left( {{\alpha_\mathrm{c}}\rho_\mathrm{c}{\tau_\mathrm{c}}} \right)
= - {\alpha_\mathrm{c}} \nabla p_\mathrm{c} + {\alpha_\mathrm{c}}{\rho_\mathrm{c}} \mathbf{g} - {\mathbf{F}}
\end{equation}

离散相模型

离散相模型就比较百花齐放。目前我比较熟悉的主要有：

1. 把离散相当做连续相，使用NS方程来表述；
2. 把离散相通过拉格朗日粒子进行跟踪，方程通过拉格朗日方程来表述
   \begin{equation}
   \frac{\mathrm{d}\mathbf{X}}{\mathrm{d} t}=\mathbf{U}_ \mathrm{d}
   \end{equation}
   \begin{equation}\label{DPM}
   m_\mathrm{d}\frac{\mathrm{d}\mathbf{U}}{\mathrm{d} t}=\mathbf{F}
   \end{equation}
3. 把离散相通过Generalized Population Balance Equation（GPBE）描述
   \begin{equation}
   \frac{\partial n}{\partial t}+\nabla_\mathbf{X}\cdot(\mathbf{U}_ \mathrm{d}n)+\nabla_{\mathbf{U}_ \mathrm{d}}\cdot(\mathbf{A} n)=S,
   \label{GPBE}
   \end{equation}

另外简单提一下GPBE，方程$\eqref{GPBE}$是用来描述粒子的大而统一的模型，因此是一个非常底层的模型，NS方程以及拉格朗日方程都可以从GPBE方程，也即方程$\eqref{GPBE}$，推导出来。但是为什么方程$\eqref{GPBE}$很少有人使用？因为太底层，难以求解。同时，如果采用矩方法求解GPBE，我们大名鼎鼎的双欧拉模型即对应的双矩GPBE模型。MPPIC模型的底层模型也和GPBE模型是相同的。总之，

• GPBE这个可以在欧拉下求解，最简单的就是双欧拉模型，
• 也可以在拉格朗日下求解，演变成为MPPIC。
• 也可以在欧拉拉格朗日下混合求解，目前用的非常少。

另外一个求解GPBE的方法就是矩方法，矩用的越多，越精确。双欧拉模型就是最简单的二阶矩模型。矩方法这面相当复杂。没准可以奋斗终身... 当然还有一些其他方法，如离散速度法啊，气体动理学格式等

下面我们根据这三种不同的离散相模型，和连续相模型耦合。讨论不同的结合情况。

欧拉欧拉模型

欧拉欧拉模型是最为人知的多相流模型，也就是双欧拉模型。其中连续相使用NS方程，离散相还是使用NS方程。欧拉欧拉模型在目前的CFD代码中已经被大量的植入了。在此我们不讨论这个模型的方程，我个人对这个模型有一点评论：

1. 欧拉欧拉模型相对于VOF模型非常便宜，在文献中我见过模拟气泡的如果采用VOF需要上千万的网格，欧拉欧拉模型几十万网格就搞定；

2. 欧拉欧拉模型界面力模型非常不完善，也就是方程$\eqref{NS}$中的$\mathbf{F}$。这个$\mathbf{F}$包含了曳力、虚拟质量力等，旗下又衍生出来各种人名的模型，这些模型直接影响欧拉欧拉模型的精准度；

3. 相分数也即$\alpha$，这个东西需要保证严格的有界性，也就是$0<\alpha<1$，但是在CFD中，这个问题很难保证，但凡存在一点越界行为，直接导致求解发散；

4. 由于NS方程对流项的问题，欧拉欧拉模型难以避免的存在数值耗散的问题；

5. 欧拉欧拉模型在某些情况下不适用，比如颗粒流中的颗粒轨迹交叉（如下图），其中两股颗粒在中心点交叉过去并穿越，欧拉欧拉模型并不能够预测这一类信息。

6. 欧拉欧拉模型中离散相只能给定均一的颗粒属性，比如颗粒都是5mm，传质系数均一；

7. 欧拉欧拉模型无法预测碰撞、聚并、破碎等统计信息，因此大量的文献在欧拉欧拉模型的基础上，和PBE耦合来预测粒子的粒径、传质等信息。但是这种方法依然存在很多的困哪点：

   1. 欧拉欧拉模型中可以预测离散相的相分数，也即$\alpha$，PBE模型同样也可以预测相分数，二者需要相同才能表明结果是合理的，但是目前这方面的研究很少；

   2. 由于上面讨论的第三点，相分数传输方程通常需要特殊的数值格式处理，但是这个PBE中并不一定使用，因此进一步的导致PBE预测的相分数和欧拉欧拉模型预测的相分数存在偏差；

   3. 即使结合PBE可以计算不同的颗粒属性（比如有的5mm有的3mm），但是不同颗粒的传输速度在大部分文章中都是给定同样的传输速度。这在有些情况下并不适用；

   4. PBE本身很难求解，目前依然是国际上的研究热点；

但是，由于欧拉欧拉模型应该是发现起来非常早的模型了（应该起始于上个世纪80年代），目前大部分研究都是基于欧拉欧拉模型，并且被大量的用于工业上的模拟。

欧拉拉格朗日模型

相对于欧拉欧拉模型，欧拉拉格朗日模型要贵很多。因为它需要跟踪每个离散粒子的位置信息。欧拉拉格朗日模型目前在CFD代码中也已经被大量的植入。并且算法个人感觉要比欧拉欧拉模型略简单一些。

1. 欧拉拉格朗日不仅仅比欧拉欧拉是贵一些，而是贵很多。如果颗粒上万个以上，就明显感觉出来计算要慢。如果颗粒上百万，千万，那欧拉拉格朗日是非常贵了；

2. 但是由于欧拉拉格朗日求解的传输方程的问题，也即方程$\eqref{DPM}$，拉格朗日方程不存在数值耗散的问题！

3. 虽然方程$\eqref{DPM}$形式简单，但是复杂的在于颗粒位置跟踪，目前好多的CFD算法在研究如何高效的捕获颗粒位置，否则会产生粒子丢失，计算死循环等问题；

4. 由于源项更新计算问题，欧拉拉格朗日算法要求离散相粒子传输时间步长小于连续相NS方程的时间步长，这进一步的导致计算比较贵；

5. 拉格朗日粒子碰撞模型涉及到了大量的模型，但凡模型的引入，都会导致一定的不确定性，目前软球模型和硬球模型依然是国际上的研究热点；

6. 欧拉拉格朗日可以很好的模拟颗粒轨迹交叉，并且颗粒传输的速度是不相同的；

7. 欧拉拉格朗日的噪音难以消除，需要大量大量的粒子才能减轻；

8. 目前普通的算法需要颗粒直径小于网格直径；

   

总体来说，个人觉得欧拉拉格朗日方法用于模拟气固流动要比欧拉欧拉模型精度高。但是就是太贵了。90年代后期一个稍微便宜一点的欧拉拉格朗日方法MPPIC被研究出来，但是还是略贵。欧拉拉格朗日算法不存在人工粘性、不存在singular方程。但是粒子跟踪的有效性太复杂了。也见过2篇发在JCP上的有效提升粒子跟踪信息的算法。

欧拉-矩方法

在这里强调一下，严格来讲，欧拉-矩方法也算是一种欧拉欧拉方法，因为矩传输是在欧拉框架下，但是为了区别传统的欧拉欧拉，在这里成为欧拉-矩方法。

上面的欧拉欧拉模型，欧拉拉格朗日模型基本都起始于上世纪80,90年代。应该比我们大多在读的CFDer还要老。欧拉QBMM模型的研究起始于2010年左右，这个算法属于新兴的小弟弟了。目前国内，甚至国际上用的也是极为、非常少。从算法来讲，这个方法依然是通过NS方程来模化连续相，但是离散相的模化采用了最为底层的模型，也即方程$\eqref{GPBE}$。用一句话来概括欧拉QBMM，就是对最底层的离散相方程$\eqref{GPBE}$通过矩进行封闭，可以在欧拉、拉格朗日甚至混合框架下自由的求解。

方程$\eqref{GPBE}$除了通过矩进行封闭，也可以通过别的模型来封闭，比如大名鼎鼎的Class Method（CM）。通过CM来进行封闭主要是我同事的工作。我主要用矩进行封闭。但问题是这个方程不管怎么封闭，依然是太复杂，这里讨论的封闭有点像是二阶矩湍流模型，你在封闭湍流的时候，如果使用二阶矩进行封闭，那就会产生三阶矩项。怎么办？模化。如果你把三阶矩处理掉，那就会出现四阶矩，怎么办？模化。不管怎么搞，最后都的模化。

问题的难点在于，如何考虑不同颗粒的传输速度？在很多情况下我们都希望小颗粒有小颗粒的传输速度，大颗粒有大颗粒的上浮速度。很明显，小气泡和大气泡的上浮速度是不同的。很遗憾，PBE不能处理这个问题！。在这种情况下只能选用GPBE模型。考虑仅仅粒子的速度作为变量，就变成了大名鼎鼎的动理学方程。动理学方程可以采用CM或者QBMM来求解，其中QBMM这面从2008年后得力于Massot、Fox、Marchisio的贡献研究的如火如荼。很遗憾，这方面的算法目前我没见过有一家CFD商软在研究，基本局限于学术界。我目前也在做这个。简单评论一下：

1. 欧拉-矩方法可以非常便宜的处理粒子破碎、聚并等行为，这样就可以很好的预测多相系统中的颗粒粒径大小，不同的传质系数，这对反应流是非常重要的；

2. 欧拉-矩方法可以不依赖于欧拉欧拉模型！ 因为离散相的传输方程被转化为了混合矩传输方程，这样就不存在欧拉模型预测的相分数和PBM预测的相分数不统一问题；

3. 欧拉-矩方法可以更新不同粒径的颗粒具有不同的传输速度，然而这方面研究就更少了，基本在2015年之后了；但是，确实是可以的！

4. 欧拉-矩方法可以预测颗粒流的颗粒轨迹交叉；

5. 矩传输方程的realizable问题一直是欧拉QBMM的痛点，在对流项的计算中，高阶空间格式的使用会使得计算发散，目前只有2种方法能处理这个realizable的问题（分别在2011和2017年发在了JCP上），Vikas的算法还比较简单易用，非常不错

6. 欧拉-矩方法要比欧拉欧拉模型相对贵一些，但是远远比欧拉拉格朗日便宜。在欧拉欧拉模型中，求解的是相分数方程、一个离散相速度方程、一个连续相速度方程，一个压力方程共4个方程。在欧拉-矩方法中（如果使用3节点QMOM算法），3维问题需要求解一个连续相速度方程，一个压力方程，15个混合矩传输方程，共17个方程。传输方程要多一些，但是依然比上百万个粒子要便宜很多！

7. 欧拉-矩方法模型的精度依赖于调用的节点数量，越多的节点，精度越高，但是稍微要贵一些

工业应用

那么这一段讨论应该会轻松一些。毕竟不会涉及到各种乱七八糟的模型。下面我只讨论我比较熟悉的工业应用，也就是说，这些模型的应用远不局限于下面的讨论。我没有列举，应为我没有应用过因此不好评价 :cheeky: 同时我在工业应用这面研究的也不是很透彻，这方面的工程师应该比我更有经验。

气液/液液分散体系：那在这方面的应用就很多了。在工业中，气液/液液分散体系一般选用最便宜的欧拉欧拉模型，预测一下气含率、搅拌功率，算一算传质系数，产率多少，费多少电？放大一下试试？

气固体系：气固那面玩的和气液/液液不太一样（比如各种流化床，气力输送等），固体颗粒不同意气泡或者液滴，一般不会破碎。但是交叉，团聚效应某些情况下非常明显。同时，可能还需要加上点传热以及化学反应。dusty-gas那面还需要考虑激波问题，也很复杂。

其他的工业应用也有很多，比如稀薄气体、喷雾、燃烧烫烟、大气污染、石油管道输运、这些文章都见过一些，但是并没有详细去研究过 多相流难搞，我这几年，一直在和收敛作斗争。欢迎大家多多交流、补充自己的见解，互相学习。

[yhdthu]

@李东岳 前辈好总结！可以看到:happy:

[peperpanDUT]

@李东岳 可以

[hurricane007]

我记得之前在群里面讨论的时候有哥们说 MPPIC 模拟工业流化床比E-E快，似乎是 守望？ @散漫守望2016

[anubis]

如果你需要气泡的形状，波浪的表面，你就需要捕获气泡或者波浪的界面，因此你需要使用界面跟踪类模型
这里的图2得多少网格

[李东岳]

@anubis
图二还不是VOF，是discrete bubble model，相对VOF也便宜不少，网格没有VOF多，具体多少不清楚，哈哈。我手头有那个文献，不过人在机场2018年给你发 :cheeky:

[anubis]

@李东岳 好的 谢谢！

[金石为开]

老师 目前在用欧拉欧拉研究气泡减阻 ，想考虑别的模型，尝试过，但失败了，不知道咋办，另外您的euler-dqmm可以用于气液吗？

[李东岳]

@金石为开
可以，气液、液液、气固都可以，但这个算法太复杂了，入行需要谨慎。IATE那个挺好的，算法比较简单，用的也比较少，很适合发文章。

[金石为开]

@李东岳 想用复杂的 opefoam目前好像没有，=哎

[李东岳]

@金石为开
IATE不合适么，我还计划对比一下IATE、QMOM呢。
复杂点就是CFD-PBM耦合，用CM或者矩方法求解。你要复杂的就得自己搞啊。你也可以试试OPENQBMM。

[金石为开]

@李东岳 曾经网上下载过 没研究明白

[李东岳]

@金石为开
Pbm算法很复杂，问题也非常多，入坑需谨慎，已经有前人为止奋斗终生

[李东岳]

我对界面类模型也比较感兴趣。据我所知用的比较多的是

• CICSAM
  * HRIC
  * OpenFOAM植入的，暂时称之为Weller+MULES
  • 网格依附类

* Level-Set类
* SLIC等

但是目前研究的不是很深入，只是了解界面捕获类，感觉写不出来很宏观很整体的介绍。有人对界面那面各种方法都玩的比较深入的么，概括概括让大家学习学习。

[wwzhao]

@李东岳 界面跟踪方法（interface tracking method）要求网格与两相界面完全重合，所以level-set和SLIC应该属于界面捕获方法（interface capturing method）。

[winsway_zero]

总结的非常好，给我刚入门的小白有了一个框架，非常感谢！

[wenjinlv]

@李东岳

[yhdthu]

@wwzhao 前辈这种分类是按照advect的显隐式处理分的类，几何重构属于显式tracking，代数重构属于隐式capturing

[李东岳]

以上主观言论来自于我这几年一行一行拆解开源代码OpenFOAM中VOF求解器interFoam，双流体模型求解器twoPhaseEulerFoam、自己写的Euler-QBMM求解器以及近期研究的欧拉拉格朗日求解器MPPICFoam的经验总结。

我们最终确定将使用Euler-QBMM求解器进行计算的算例以及代码提交给Computer Physics Communication

@wenjinlv 咋了老铁

@anubis 抱歉才想起来，请看这个文章 https://www.sciencedirect.com/science/article/pii/S0009250914004448

[l.j刘侃]

十分赞许李博士在多相流领域所做的开创性工作。我在多相流的工业应用领域也有一定的尝试，主要是针对灰尘堆积部分。自己也发展了一些比较直接的工程化模型。多相流的收敛失败去年达到840多次。个人认为，收敛问题是一个非常关键的问题，如果不能在这个方面上取得突破，难以在真正多相流的大项上推广，可能只能针对某一特殊工艺。如果在碰撞和孔隙结构对工况的影响不能忽略的情况下，采用拉格朗日法是难以避免的。而即便采用目前的MP-PIC方式，在一个全工程尺度上，依然是不实际的。从现在的发展来看，可能近期想可能能够取得的应用还主要是针对特定工况的。PBM对于空间的体分的拟合还能接受，但是对于静态堆积的拟合还是有一定差距。而对于一些非微观模型，如何在工业中直接获取模型中的变量值又是一个头疼的问题。总得来看，针对多相流的工业化应用，发展出一个统一的适用模型，在近期可能还是比较困难的。针对特定工况的，多相流求解器可能是近期的一个热点。