多相流模拟的数学模型


  • 网格教授 OpenFOAM教授 管理员

    多相流远比单向流复杂,这也导致多相流目前在工业上应用的远比单向流少。上图我列举了针对高相分数两相体系的数学模型。每一种方法下都有大量的人在研究验证,很遗憾,每一种模型都有大量的问题存在。


  • 网格教授 OpenFOAM教授 管理员

    首先,这篇文章是我随便写的,11月份随便写了一点,见一楼,现整体移动到本楼。发在论坛上,不会那么严谨,请各位同行多多包涵!但当然还是尽可能的通俗的总结,学术讨论请参考期刊。我也不会引用各类文献,排版太费事。某些中文译名可能不准确,但是大家都懂的。

    多相流模拟

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

    • 模型不成熟,比如各种曳力模型,升力模型,破碎模型等;
    • 求解算法不成熟,比如在CFD求解多相流模型中如何保证相分数严格的小于1大于0;
    • 其他;这方面我还没怎么研究过

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

    多相流模型

    0_1513856589293_图片1.jpg

    多相流中区分为多个相。比如空气是一相、液滴是一相。也可以称之为气相、液相或者固相。

    参考上图。多相模拟主要分为界面捕获类模型(如VOF)和高相分数多相流模型。二者的区别非常简单:

    • 如果你需要气泡的形状,波浪的表面,你就需要捕获气泡或者波浪的界面,因此你需要使用界面跟踪类模型。
      0_1513850530854_微信图片_20171221105717.jpg

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

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

      0_1513850660103_捕获.JPG

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

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

    高相分数多相流模型

    在高相分数多相流模型中,一般区分为连续相离散相(同样参考第一个图)。举例:

    1. 水里吹气泡(鼓泡床),气泡是离散相,水是连续相;
    2. 沙子里吹气(流化床),空气是连续相,沙子是离散相;
    3. 水里面加油(液液分散),水是连续相,油是分散相;

    连续相模型

    连续相模型一般都使用NS方程来描述(2018.01.04更新:石油大学某同学看的非常仔细,方程2有个正负号错误)。

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

    简单评论一下,上面三种模型都是模化离散相的数学模型。如果连续相的流场给定的话,可以进行单相流模拟

    1. 在使用NS方程模拟离散相,并且给定连续相流场,这是单向耦合
    2. 在使用拉格朗日方程模拟离散相,并且给定连续相流场,这是DPM
    3. 在GPBE模拟离散相,并且给定连续相流场;

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

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

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

    欧拉欧拉模型,评论A

    欧拉欧拉模型是最为人知的多相流模型,也就是双欧拉模型。其中连续相使用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年代),目前大部分研究都是基于欧拉欧拉模型,并且被大量的用于工业上的模拟。

    欧拉拉格朗日模型,评论B

    替代文字

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

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

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

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

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

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

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

      0_1513854644748_捕获.JPG

    总体来说,个人觉得欧拉拉格朗日方法用于模拟气固流动要比欧拉欧拉模型精度高。但是就是太贵了。90年代后期一个稍微便宜一点的欧拉拉格朗日方法MPPIC被研究出来,但是还是略贵,有关MPPIC的看这个贴,目前我也在研究中。不能发表更深的见解。

    对了,最近也在研究OpenFOAM里面植入的欧拉拉格朗日算法,深深地感觉到拉格朗日框架的求解要比NS方程简单的多。起码欧拉拉格朗日算法不存在人工粘性、不存在singular方程。但是粒子跟踪的有效性太复杂了。也见过2篇发在JCP上的有效提升粒子跟踪信息的算法。

    欧拉QBMM,评论C

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

    方程$\eqref{GPBE}$除了通过矩进行封闭,也可以通过别的模型来封闭,比如大名鼎鼎的Class Method(CM)。通过CM来进行封闭主要是我同事的工作。我主要用矩进行封闭。但问题是这个方程不管怎么封闭,依然是太复杂,因此在求解的时候还可以进行进一步的简化。在最简单最简单的情形,GPBE方程化为了单变量的PBE方程。通过把PBE方程和欧拉欧拉模型耦合,就是我们评论A里面提到的一种算法。

    这里讨论的封闭有点像是二阶矩湍流模型,你在封闭湍流的时候,如果使用二阶矩进行封闭,那就会产生三阶矩项,怎么办?模化。如果你把三阶矩处理掉,那就会出现四阶矩,怎么办?模化。不管怎么搞,最后都的模化。

    问题的难点在于,如何考虑不同颗粒的传输速度?在很多情况下我们都希望小颗粒有小颗粒的传输速度,大颗粒有大颗粒的上浮速度。很明显,小气泡和大气泡的上浮速度是不同的。很遗憾,PBE不能处理这个问题!。在这种情况下只能选用GPBE模型。GPBE理论上也可以采用CM或者QBMM来求解,其中QBMM这面从2011年后得力于ISU的Fox的贡献研究的如火如荼。很遗憾,这方面的算法目前我没见过有一家CFD大厂在研究,基本局限于学术界。我目前也在做这个。简单评论一下:

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

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

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

    4. 欧拉QBMM可以预测颗粒流的颗粒轨迹交叉

    5. 矩传输方程的realizable问题一直是欧拉QBMM的痛点,在对流项的计算中,高阶空间格式的使用会使得计算发散,目前只有2种方法能处理这个realizable的问题(分别在2011和2017年发在了JCP上),我也在测试中;个人感觉realizable这个问题,基本上属于QBMM的核心问题,这个问题不解决掉,高阶格式没得用。

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

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

    工业应用

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

    气液/液液分散体系

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

    据我了解,巴斯夫等很多化工500强大厂每年都投入大量的研究经费来进行CFD模拟他们的反应器效率,国内也有些化工大厂想要通过CFD来优化产品设计,减少设计成本。据说一次中试放大有时候都要上亿 :wocao: 那CFD就便宜太多了。但是玩真格的,欧拉欧拉就略心有余力不足。比如我们需要预测气泡的粒径和传质系数,再加点化学反应,传个热,欧拉欧拉模型就不能预测了。需要调用欧拉欧拉-多变量PBE耦合,或者欧拉QBMM。再来点相变、沸腾、结晶?很复杂。

    气固体系

    气固那面玩的和气液/液液不太一样(比如各种流化床,气力输送等),固体颗粒不同意气泡或者液滴,一般不会破碎。但是固体颗粒会交叉!也就是颗粒轨迹交叉。导致欧拉欧拉完全预测不了PTC。如果PTC很重要的话,就需要用欧拉拉格朗日,虽然很贵。也可以用欧拉QBMM。目前使用欧拉QBMM进行气固模拟,基本上工作也只局限于最近2 3年。气固这面在锅炉那面可能还要加上点传热以及化学反应,也很复杂。

    其他工业应用

    其他的工业应用也有很多,比如稀薄气体、喷雾、燃烧烫烟、大气污染、石油管道输运、这些文章都见过一些,但是并没有详细去研究过 :quwan:

    多相流难搞,我这几年,一直在和收敛作斗争。欢迎大家多多交流、补充自己的见解,互相学习。


  • 网格教授 OpenFOAM教授 管理员

    国内GIF能看到么?


  • OpenFOAM讲师

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



  • @李东岳 可以



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



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


  • 网格教授 OpenFOAM教授 管理员

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



  • @李东岳 好的 谢谢!



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


  • 网格教授 OpenFOAM教授 管理员

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



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


  • 网格教授 OpenFOAM教授 管理员

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



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


  • 网格教授 OpenFOAM教授 管理员

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


  • 网格教授 OpenFOAM教授 管理员

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

    • CICSAM
    • HRIC
    • OpenFOAM植入的,暂时称之为Weller+MULES
    • 网格依附类
    • Level-Set类
    • SLIC等

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


  • OpenFOAM教授

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



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



  • @李东岳


  • OpenFOAM讲师

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


  • 网格教授 OpenFOAM教授 管理员

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

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

    @wenjinlv 咋了老铁

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