1.4 多相流
多相流通常包括气-液、液-液、液-固、气-固、气-液-液、气-液-固、气-液-液-固混合物的流动，在机械、能源、动力、核能、石油、化工、冶金、制冷、运输、环境保护及航天技术等许多领域都有其踪迹。多相流的建模计算通常比较难，尤其是追踪流体与流体之间的交界面更难，这里简单介绍COMSOL中多相流的相关内容。
1.4.1 多相流模型的分类
在COMSOL中，多相流模型分为两大类：分离多相流模型与分散多相流模型。
在较小尺度上，我们可以对相边界的形状进行详细建模，并把这种模型称为分离多相流模型。分离多相流模型有清晰的相界面，对此我们一般用相场
描述相界面，
是一相，而
是另一相。分离多相流模型主要用于气泡、液滴和颗粒流的模拟分析，这些相界的尺度量级与流场尺度相当，而且数量较少；也可用于微流体中的多相流、宏观流场中的单相流的自由液面。通常，我们使用表面追踪法来描述此类模型。
在较大尺度上，如果仍详细描述相边界，则模型方程无法求解。这时我们可以使用体积分数场描述不同的相，并把这种模型称为分散多相流模型。在分散多相流模型方程中，相间效应（例如表面张力、浮力和跨越相边界的传递）被视为源和汇。分散多相流模型的相界面不太清晰，为此我们用体积分数场描述相间的关系，其中
。分散多相流模型主要用于气泡数量多而气泡体积小的气泡流，也可用于乳液和气溶胶的模拟，还可用于流场中有大量固体颗粒的情况以及宏观的多相流分析。
如图1-20所示，分离多相流模型详细描述了相边界，分散多相流模型则只考虑分散在连续相中的一个相的体积分数。在分离多相流模型中，不同相之间相互排斥，并存在一个清晰的相边界，在此边界上相场函数发生突变。除了追踪相边界的位置，相场函数没有任何物理意义。在分散多相流模型中，函数描述了气相（分散相）和液相（连续相）的局部平均体积分数。通过平均体积分数可以在该区域的任一点顺利地找到介于0 和1之间的值，这预示着在其他均质域中是存在少量还是大量气泡。也就是说，在分散多相流模型中，可以在同一时间和空间点上定义气相和液相；而在分离多相流模型中，在给定的时间和空间点上，只能定义气相或液相。
（a）分离多相流模型
（b）分散多相流模型
图1-20 多相流模型
1.4.2 分离多相流模型
对于分离多相流模拟，COMSOL软件提供了3种不同的界面追踪方法：相场法、水平集法与移动网格法。
相场法和水平集法都是基于场的方法，其中相之间的界面代表相场或水平集函数的等值面。与上述两种方法完全不同，移动网格法将相界面模拟为分隔两个域的几何表面，每个域对应不同的相。
基于场的问题通常是在固定的网格上解决，而移动网格问题要在移动的网格上解决。在相场法和水平集法中，有限元网格不必与两个相的边界一致，如图1-21所示的搅拌自由液面模拟。
（a）相场法
（b）水平集法
图1-21 相场法与水平集法
对于移动网格，网格与相边界的形状保持一致，并且网格边缘与相边界重合，如图1-22所示的搅拌自由液面模拟，是在单相流中用移动网格模拟自由液面。但是，移动网格模型也有缺点，即目前无法处理拓扑变化（例如界面分离等），而相场法与水平集法不存在这个缺点，可以处理相边界形状的任何变化。
图1-22 移动网格法
1．相场法、水平集法和移动网格法的选择策略
对于给定的网格，移动网格法具有更高的精度。基于这一优势，我们可以直接在相边界上施加力和通量。基于相场的方法需要围绕相边界表面建立密集网格，以解析该表面的等值面。由于很难定义一个精确贴合等值面的自适应网格，因此通常必须在等值面周围建立大量密集网格。在具有相同精度的情况下，与移动网格相比，这样做会降低基于场的方法的效率。
对于不希望发生拓扑变化的微流体系统，通常首选移动网格法；如果需要拓扑变化，则必须使用相场法。如果表面张力的影响较大，则首选相场法；如果可以忽略表面张力，则首选水平集法。
2．分离多相流模型和湍流模型的结合使用
在湍流模型中，由于仅解析平均速度和压力，流体的细节会丢失。从这一点来看，表面张力效应在流体的宏观描述中也变得不那么重要。由于湍流表面的流动比较剧烈，几乎不可能避免拓扑变化，因此对于湍流模型和分离多相流模型的组合，最好使用水平集法。
在COMSOL软件中，所有湍流模型都可以与相场法和水平集法相结合来模拟两相流，例如将水平集法与k-ε湍流模型相结合来模拟反应堆中水和空气的两相流，如图1-23所示。
图1-23 反应堆中水和空气的两相流
1.4.3 分散多相流模型
如果多相流的相边界过于复杂而无法解析，则必须使用分散多相流模型。COMSOL软件的CFD模块提供了4种不同的分散多相流模型（原理上）：①气泡流模型，适合高密度相中包含较小体积分数的低密度相；②混合物模型，适合连续相中包含较小体积分数的分散相（或几个分散相），其密度与一个或多个分散相相近；③欧拉-欧拉模型，适用于任何类型的多相流，可以处理气体中有密集颗粒的多相流，例如流化床；④欧拉-拉格朗日模型，适合包含相对较少（成千上万，而不是数十亿）的气泡、液滴或悬浮颗粒流体，也适合气泡、颗粒、液滴或使用方程模拟的颗粒，该方程假定流体中每个颗粒的力平衡。
1．分散多相流模型的选择策略
（1）气泡流模型显然适用于液体中的气泡。由于忽略了分散相的动量贡献，因此该模型仅在分散相的密度比连续相小几个数量级时才有效。
（2）混合物模型与气泡流模型相似，但考虑了分散相的动量贡献，通常用于模拟分散在液相中的气泡或固体颗粒。混合物模型还可以处理任意数量的分散相。混合物模型和气泡流模型均假设分散相与连续相处于平衡状态，即分散相不能相对于连续相加速。因此，混合物模型无法处理分散在气体中的大固体颗粒。
当多相流混合物被迫通过孔口时，用混合物模型模拟了5种不同大小的气泡，流动中的剪切力导致较大的气泡破裂成较小的气泡，如图1-24所示。
图1-24 混合物模型
（3）欧拉-欧拉模型是最精确的分散多相流模型，也是用途最多的分散多相流模型。该模型可以处理任何类型的分散多相流，允许分散相加速，对不同相的体积分数也没有限制，但是它为每个相定义了一组Navier-Stokes方程。
在实践中，欧拉-欧拉模型仅适用于两相流，并且其计算成本（CPU时间和内存）较高。正因如此，该模型使用起来也相对困难，并且需要良好的初始条件才能在数值解中收敛。使用欧拉-欧拉多相流模型模拟流化床中固体颗粒的体积分数分布如图1-25所示。
图1-25 欧拉-欧拉模型的模拟
（4）如果连续流体中悬浮有一些（成千上万，但不是数十亿）非常小的气泡、液滴或颗粒，则可以使用欧拉-拉格朗日模型模拟多相流系统。该方法的优点是计算成本相对较低。从数值的角度来看，这个方法通常也不错。因此，如果连续流体中分散相的颗粒数量相对较少，那么优选欧拉-拉格朗日模型，其模拟效果如图1-26所示。
图1-26 欧拉-拉格朗日模型的模拟
还有一些方法可以使用欧拉-拉格朗日模型来模拟大量粒子，它们使用的相互作用项和体积分数可以模拟具有数十亿个粒子的系统。这些方法可以在COMSOL软件中实现，但在预定义的物理接口中无法实现，需要自定义或多模块组合来实现，如在COMSOL软件中用附加的CFD模块和粒子追踪模块可实现欧拉-拉格朗日多相流模型。
混合物模型能够处理任何相的组合，并且计算成本较低。在大多数情况下，我们可以使用此模型模拟。对于流化床（具有高密度和高体积分数的大颗粒分散相）之类的系统，只能使用欧拉-欧拉模型模拟。
2．分散多相流模型和湍流模型的结合使用
各种分散多相流模型本质上是近似的，并且也与近似的湍流模型非常吻合。可以在分散相和连续相之间以及在分散相中的气泡、液滴和颗粒之间引入相互作用。这些相互作用的起源可以是用湍流模型模拟的湍流。气泡流、混合物流和欧拉-拉格朗日多相流模型可以与 COMSOL软件中的所有湍流模型结合使用。
1.4.4 水平集法
水平集法是Osher和Sethian于1988年提出的一种用于界面追踪的数值方法。在水平集法中，界面被看作零水平集的光滑函数。由于水平集函数的对流本身很光滑，因此可以代替界面处对流引起的物性陡变梯度。虽然水平集法不像其他某些方法拥有守恒属性，但它的优势在于能够轻易计算界面曲率。水平集法采用连续逼近方法，将表面张力和交界面局部曲率表示为体积力，这简化了在计算中捕捉由表面张力变化引起的拓扑结构变化过程。
在水平集法中，用水平集平滑函数来描述两相交界面。在连续相中水平集函数始终为正，在分散相中始终为负，而相间交界表面是由水平集函数为零的点构成，即
（1-103）
从上面的描述可知，交界面上的单位法线由分散相指向连续相，交界面的曲率可以用水平集函数表示为
（1-104）
（1-105）
交界面的运动可以通过水平集函数的对流来捕获：
（1-106）
流速和压力的控制方程可以用不可压缩N-S方程表示：
（1-107）
（1-108）
其中，F是体积力，包括重力和由于对交界面应力进行水平集处理引入的表面张力项。F的两个分量可以表示为
（1-109）
（1-110）
其中，
为张力系数，
为界面曲率，
函数用来处理交界面处的表面张力项，可以有很多个液-液交界面来标定分散相。描述物性急剧变化的Heaviside函数可以用水平集函数表示为
（1-111）
计算区域流体的密度与黏度可以表示为
（1-112）
（1-113）
其中，
为第一相的密度，
为第二相的密度，
为第一相的黏度，
为第二相的黏度。
在COMSOL软件微流动的不混溶的两相流模拟中，水平集描述的两相流界面输运满足如下方程
（1-114）
其中，第三项为保持数值稳定性所需的项；
为界面厚度控制参数，大多数情况下使用默认值
，
为流域典型网格大小；
为水平集函数重新初始化或稳定性的参数，通常情况下是采用流体流动的最大速度。太小可能会使相界面厚度不再保持恒定，也可能会由于数值不稳定性引起的振荡，而太大会导致相界面的移动不正确。
1.4.5 相场法
相场法是基于Cahn-Hilliard和Ginzburg-Landau方程的改进数值方法，其相场变量由Cahn-Hilliard扩散方程决定，可以用以下两个二阶偏微分方程表示。
（1-115）
（1-116）
其中，为相场变量，取值为（‒1, 1）；λ为混合能密度；ε为界面厚度控制参数，用于评价相界面的厚度，λ与ε两个参数与表面张力系数相关，有
；γ 为迁移率，γ 足够大时，可使界面厚度保持恒定，γ 足够小时，将使得对流项不会被过分抑制，γ=χε2；χ为迁移调节参数；ψ为相场助变量；
为外部自由能（多数情况下为0）。
相场法求解时，主场方程中流体属性的控制方程为
（1-117）
（1-118）
（1-119）
相场法中的参数设置：迁移调节参数
的默认值为1，对于大多数模型来说是一个不错的初始值；界面厚度控制参数
大多数情况下使用默认的，为流域典型网格大小。
1.4.6 混合物模型
混合物模型是一种宏观两相流模型，在许多方面类似于气泡流模型。该模型跟踪平均相浓度或体积分数，并求解混合物速度的单个动量方程，适用于由浸没在液体中的固体颗粒或液滴组成的混合物。
在混合物模型中，颗粒与流体组合被视为具有宏观特性（如密度与黏度）的单个连续流动体，一般是由分散相与连续相组成的两相流。例如，连续相为液体，分散相为固体颗粒、液滴或气泡。然而对于液体中的气泡，气泡流模型更合适些。混合物模型依赖于以下假设：①每相的密度是近似常数；②两相共享相同的压力场；③颗粒的松弛时间少于宏观流动的时间尺度。
混合物模型可以求解连续性方程、动量方程、分散相体积分数输运方程、质量传输方程、湍流方程、滑移速度方程等，这里简单介绍一下前面3个方程，其他的请查阅参考文献[16]。
1．混合物模型的连续性方程
混合物模型的连续性方程为
（1-120）
其中，
是混合密度，计算式为
（1-121）
ս是混合速度，计算式为
（1-122）
其中，
分别是连续相、分散相的体积分数；
分别是连续相、分散相的密度；
分别是连续相、分散相的速度。
2．混合物模型的动量方程
混合物模型的动量方程可以通过对所有相各自的动量方程求和来获得，可表示为
（1-123）
其中，u为速度矢量；
为压力；为减少的密度差，
；
为两相的滑移速度矢量；
为滑移通量，
；
为黏性与湍流应力之和；
为湍流耗散系数；
为分散相到连续相的质量传输率；F为任意外部体积力。
3．分散相体积分数输运方程
分散相体积分数输运方程为
（1-124）
1.4.7 欧拉-欧拉模型
欧拉-欧拉模型界面基于在体积上平均每个当前相的Navier-Stokes方程，该体积与计算区域相比较小，但与分散相（颗粒、液滴或气泡）相比较大。
欧拉-欧拉模型包含连续性方程、动量方程、黏度方程、相间动量输送方程、固体压力方程等，这里简单介绍其中的一些方程，详细内容请查阅参考文献[16]。
1．连续性方程
欧拉-欧拉连续性方程适用于连续相与分散相：
（1-125）
（1-126）
其中，是相体积分数，
，表示密度，u为速度，各项中的下标c和d分别表示连续相与分散相；为分散相到连续相的质量传输率。
2．动量方程
连续相和分散相的动量方程使用非保守形式，即
（1-127）
（1-128）
假设上述方程中的流体相为牛顿流体，黏性应力张量定义为
（1-129）
（1-130）
其中，是混合流体压力，
是每相的黏性应力张量，g为重力加速度，
为相间动量传递项（一相被其他相施加的体积力），F是任何其他体积力，
是相间速度，
为动力黏度，各项中的下标c和d分别表示连续相与分散相。
3．黏度方程
欧拉-欧拉模型中使用混合物黏度的表达式，两个互穿相的动力黏度默认值为
（1-131）
简单的混合物黏度（能覆盖整个颗粒浓度范围）方程可以用Krieger方程表示为
（1-132）
其中，
为最大填充限制，对于固体颗粒，其默认值为0.62。
1.4.8 气泡流模型
双流体欧拉-欧拉模型是两相流体流动的一般宏观模型，它将两相视为互穿介质，跟踪相的平均浓度。其中的速度场与每个相场相关联，动量方程和连续性方程描述每个相的动力学过程。气泡流模型是双流体模型的简化，它依赖于以下假设：①与液体密度相比，气体密度可以忽略不计；②气泡相对于液体的运动由黏性阻力和压力之间的平衡决定；③两个阶段共享相同的压力场。
基于这些假设，列出两相流的动量方程和连续性方程，结合气相输运方程，就可以跟踪气泡的体积分数。
1．动量方程
（1-133）
2．连续方程
（1-134）
上面的方程中，是相体积分数，表示密度，u为速度，p是压力，
是每相的黏性应力张量，
为重力加速度，F是任何额外体积力，
是液体的动力黏度，
为湍流黏度，各项中的下标l和g分别表示液相和气相。
3．气相输运方程
（1-135）
其中，mgl为质气体到液体的质量传输。
在实际应用时，上述模型可能根据实际情况进行变化，限于篇幅，这里不一一说明，详细内容请查阅参考文献[16]。
1.4.9 移动网格模型
对于层流两相流，当关注界面的精确位置时，移动网格模型可用于模拟两种不同的不混溶流体的流动。界面位置由移动网格跟踪，边界条件考虑了表面张力和润湿以及界面上的质量传输。两相流移动网格界面是单相流界面和移动网格界面之间的预定义物理界面耦合。在对应于各个相的区域内，流体流动使用Navier-Stokes方程求解。
我们通过流体流动域内的网格变形，来说明两种流体之间的界面。软件会扰动网格结点，使其与移动界面以及模型中的其他移动或静止边界一致。边界位移在整个区域中传播，以获得平滑的网格变形。这是通过求解网格位移方程（拉普拉斯方程、温斯洛方程或超弹性平滑方程）实现的。通常情况下，在“移动边界平滑”选项中根据式（1-136）平滑法向网格速度。
（1-136）
其中，X为x坐标变化量，n为单位矢量，
为理想的法向网格速度，
是平滑速度，
，
是移动边界平滑调整参数（无量纲），
是网格尺寸（单位：m），
是平均表面曲率（单位：1/m）。
以二维为例，变形网格中的一个位置坐标（x，y）可以与其在原始未变形网格中的坐标（X，Y）相关，写成函数形式为
（1-137）
原始未变形的网格称为材料框架（或参考框架），而变形网格称为空间框架。COMSOL还定义了几何体和网格框架，这些框架与该物理界面的材料框架一致。流体流动方程（以及其他耦合方程，如电场或化学物质传输方程）在网格被扰动的空间框架中求解。因此，在这些界面中考虑了相边界的移动。
用网格位移和流体流动的特定边界条件跟踪两相之间的界面。有两个选项可用：自由表面和流体-流体界面。当外部流体的黏度与内部流体的黏度相比可以忽略时，自由表面边界条件是合适的。在这种情况下，外部流体的压力是建模流体所需的唯一参数，并且在外部流体中不求解流动。对于流体-流体界面，对两个相的流动进行求解。
1．流体-流体界面
两种不混溶流体（流体1和流体2）界面处的边界条件为
（1-138）
（1-139）
（1-140）
其中，
和
分别是流体1和流体2的速度，
是两流体间界面网格速度，
为界面法向（见图1-27），
和
分别为区域1与区域2的总应力张量，
为由界面张力引起的单位面积力，
，
是表面梯度算子，
是穿过界面的质量通量。
图1-27 流体1与流体2界面法向的定义
式（1-139）的切向分量在边界处的流体之间施加无滑移条件。在没有穿过边界的传质的情况下，式（1-138）和式（1-140）确保垂直于边界的流体速度等于界面的速度。当发生传质时，这些方程是质量守恒的结果，在边界静止的框架中很容易导出。
总应力张量的分量
表示垂直于
方向的每单位面积力的第
分量。因此，
（使用求和约定）被解释为作用在边界上的每单位面积的力，通常这不是边界的法线。因此，式（1-139）表示了两种流体之间界面上的力平衡。
2．自由表面
通常情况下，流体1的黏度显著大于流体2的黏度（例如对于气液界面）。在这种情况下，流体2的总应力中的黏度项可以忽略，式（1-139）变为
（1-141）
外部流体（流体2）仅通过压力项进入方程系统，并且系统可以用仅由流体1组成的域表示，该域在流体2的域中具有外部压力
的表达式（或恒定值）。