一篇多相流review献给大家，JCP编辑向我邀稿！

乘风破浪_ 2023年9月6日浏览：1475 收藏：1

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

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

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

聊聊多相流的CFD模型

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

如果你需要气泡的形状，波浪的表面，你就需要捕获气泡或者波浪的界面，因此你需要使用界面跟踪类模型。比如下面这个图，就是界面重构出来的气泡界面，明显你可以看出来气泡的形状。
如果你的泡泡特别多，并且你不在乎气泡的形状是什么样的，你只在乎我的水里面有多少空气（气含率），那么这种情况适用于高相分数多相流模型。其不能预测出气泡的形状。下面这个图是用高相分数多相流模型模拟的，很明显，看不出来气泡的形状。这两大模型，前者属于多相流中的直接模拟，后者模化多一些。正是因为前者偏向直接模拟，因此在使用VOF的时候，需要巨量的网格。

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

高相分数多相流模型

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

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

一篇多相流review献给大家，JCP编辑向我邀稿！的图3

连续相模型

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

$\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$ $\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}$

离散相模型

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

把离散相当做连续相，使用NS方程来表述；
把离散相通过拉格朗日粒子进行跟踪，方程通过拉格朗日方程来表述

$\frac{\mathrm{d}\mathbf{X}}{\mathrm{d} t}=\mathbf{U}_ \mathrm{d}$ $m_\mathrm{d}\frac{\mathrm{d}\mathbf{U}}{\mathrm{d} t}=\mathbf{F}$

把离散相通过Generalized Population Balance Equation（GPBE）描述

$\frac{\partial n}{\partial t}+\nabla_\mathbf{X}\cdot(\mathbf{U}_ \mathrm{d}n)+\nabla_{\mathbf{U}_ \mathrm{d}}\cdot(\mathbf{A} n)=S,$

另外简单提一下GPBE，方程是用来描述粒子的大而统一的模型，因此是一个非常底层的模型，在不同的领域，这东西有不同的名字。其实他更广为人知的名字，就是玻尔兹曼方程。

是的，玻尔兹曼方程是比NS方程更底层的模型。NS方程以及拉格朗日方程都可以从GPBE方程推导出来。但是为什么工程上玻尔兹曼方程用得少？因为太底层，难以求解。同时，如果采用矩方法求解GPBE，我们大名鼎鼎的双欧拉模型即对应的双矩GPBE模型。总之，

GPBE这个可以在欧拉下求解，最简单的就是双欧拉模型，
也可以在拉格朗日下求解，演变成为MPPIC。
也可以在欧拉拉格朗日下混合求解，目前用的比较少。
还有更神的，对GPBE的分布函数采用特征线方程求解出来，并求矩来获得宏观方程，就是GKS类方法。

另外一个求解GPBE的方法就是矩方法，矩用的越多，越精确。双欧拉模型就是最简单的二阶矩模型。矩方法这面相当复杂。没准可以奋斗终身...

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

欧拉欧拉模型

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

欧拉欧拉模型相对于VOF模型非常便宜，在文献中我见过模拟气泡的如果采用VOF需要上千万的网格，欧拉欧拉模型几十万网格就搞定；
欧拉欧拉模型界面力模型非常不完善，也就是方程中的。这个包含了曳力、虚拟质量力等，旗下又衍生出来各种人名的模型，这些模型直接影响欧拉欧拉模型的精准度；
相分数这个东西需要保证严格的有界性，也就是，但是在CFD中，这个问题很难保证，但凡存在一点越界行为，直接导致求解发散；
由于NS方程对流项的问题，欧拉欧拉模型难以避免的存在数值耗散的问题；
欧拉欧拉模型在某些情况下不适用，比如颗粒流中的颗粒轨迹交叉（如下图），其中两股颗粒在中心点交叉过去并穿越，欧拉欧拉模型并不能够预测这一类信息。
欧拉欧拉模型中离散相只能给定均一的颗粒属性，比如颗粒都是5mm，传质系数均一；
欧拉欧拉模型无法预测碰撞、聚并、破碎等统计信息，因此大量的文献在欧拉欧拉模型的基础上，和PBE耦合来预测粒子的粒径、传质等信息。但是这种方法依然存在很多的困难点：

欧拉欧拉模型中可以预测离散相的相分数，也即，PBE模型同样也可以预测相分数，二者需要相同才能表明结果是合理的，但是目前这方面的研究很少；
由于上面讨论的第三点，相分数传输方程通常需要特殊的数值格式处理，但是这个PBE中并不一定使用，因此进一步的导致PBE预测的相分数和欧拉欧拉模型预测的相分数存在偏差；
即使结合PBE可以计算不同的颗粒属性（比如有的5mm有的3mm），但是不同颗粒的传输速度在大部分文章中都是给定同样的传输速度。这在有些情况下并不适用；
PBE本身很难求解，目前依然是国际上的研究热点；

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

欧拉拉格朗日模型

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

一篇多相流review献给大家，JCP编辑向我邀稿！的图4

欧拉拉格朗日不仅仅比欧拉欧拉是贵一些，而是贵很多。如果颗粒上万个以上，就明显感觉出来计算要慢。如果颗粒上百万，千万，那欧拉拉格朗日是非常贵了；
但是由于欧拉拉格朗日求解的传输方程的问题，也即方程，拉格朗日方程不存在数值耗散的问题！
虽然方程形式简单，但是复杂的在于颗粒位置跟踪，目前好多的CFD算法在研究如何高效的捕获颗粒位置，否则会产生粒子丢失，计算死循环等问题；
由于源项更新计算问题，欧拉拉格朗日算法要求离散相粒子传输时间步长小于连续相NS方程的时间步长，这进一步的导致计算比较贵；
拉格朗日粒子碰撞模型涉及到了大量的模型，但凡模型的引入，都会导致一定的不确定性，目前软球模型和硬球模型依然是国际上的研究热点；
欧拉拉格朗日可以很好的模拟颗粒轨迹交叉，并且颗粒传输的速度是不相同的；
欧拉拉格朗日的噪音难以消除，需要大量大量的粒子才能减轻；
目前普通的算法需要颗粒直径小于网格直径；

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

欧拉-矩方法

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

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

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

问题的难点在于，如何考虑不同颗粒的传输速度？在很多情况下我们都希望小颗粒有小颗粒的传输速度，大颗粒有大颗粒的上浮速度。很明显，小气泡和大气泡的上浮速度是不同的。很遗憾，PBE不能处理这个问题！。

在这种情况下只能选用GPBE模型。考虑仅仅粒子的速度作为变量，其实就是玻尔兹曼方程，可以采用CM、DSMC、GKS、QBMM等来求解，其中QBMM这面从2008年后得力于Massot、Fox、Marchisio的贡献研究的如火如荼。很遗憾，这方面的算法目前我没见过有一家CFD商软在研究，基本局限于学术界。我目前也在做这个。简单评论一下：

欧拉-矩方法可以非常便宜的处理粒子破碎、聚并等行为，这样就可以很好的预测多相系统中的颗粒粒径大小，不同的传质系数，这对反应流是非常重要的；
欧拉-矩方法可以不依赖于欧拉欧拉模型！因为离散相的传输方程被转化为了混合矩传输方程，这样就不存在欧拉模型预测的相分数和PBM预测的相分数不统一问题；
欧拉-矩方法可以更新不同粒径的颗粒具有不同的传输速度，然而这方面研究就更少了，基本在2015年之后了；但是，确实是可以的！
欧拉-矩方法可以预测颗粒流的颗粒轨迹交叉；
矩传输方程的realizable问题一直是欧拉QBMM的痛点，在对流项的计算中，高阶空间格式的使用会使得计算发散，目前只有2种方法能处理这个realizable的问题（分别在2011和2017年发在了JCP上），Vikas的算法还比较简单易用，非常不错
欧拉-矩方法要比欧拉欧拉模型相对贵一些，但是远远比欧拉拉格朗日便宜。在欧拉欧拉模型中，求解的是相分数方程、一个离散相速度方程、一个连续相速度方程，一个压力方程共4个方程。在欧拉-矩方法中（如果使用3节点QMOM算法），3维问题需要求解一个连续相速度方程，一个压力方程，15个混合矩传输方程，共17个方程。传输方程要多一些，但是依然比上百万个粒子要便宜很多！
欧拉-矩方法模型的精度依赖于调用的节点数量，越多的节点，精度越高，但是稍微要贵一些