JCP：在几何流体体积框架下的多组分液滴蒸发
摘要：这项工作提出了一个创新的多组分相变模型在界面解析模拟。该两相体系采用几何流体体积(VOF)方法描述，并考虑了非等温环境下的多组分，打破了文献中通常研究的纯液滴假设。该模型包括Stefan流，并针对研究液体混合物时出现的复杂性实现了以下解决方案:1)采用耦合方法求解界面跳跃条件;Ii)体积分数场平流的液体速度获取策略合适，同样适用于密度比较大的静态液滴;Iii)和离散标量场输运方程的几何方法。该模型在开源代码Basilisk中实现，并在固定通量蒸发、Stefan问题、Epstein Plesset和Scriven问题等多个基准相变问题上进行了测试。这些算例证明了数值方法对解析解的收敛性。采用假设球对称的数值基准解，对多组分等温液滴和非等温液滴等更为复杂的构型进行了比较。代码和模拟设置在Basilisk网站上发布，使其成为VOF框架中多组件相变的第一个模型和开源实现。
多组分液滴在自然界中广泛存在，并在许多工程应用中进行了多层次的详细研究。大多数文献模型依赖于球对称的近似[1,2]，忽略了表面张力和液滴变形的影响。这种简化对于微重力条件下的液滴是准确的，在微重力条件下，液滴变形和浮力对传热传质现象的影响可以忽略不计[3]。因此，它们可用于分析复杂的化学现象，如燃烧过程[4]，但由于其局限性，需要描述气液体系相变的多维模型。
几何流体体积（VOF）方法[5]是文献中用于描述多相流系统的最流行的模型之一。它通过求解连续性方程和使用保持界面锐度的详细离散化方法，使界面得以保守传输。VOF模型已被用于各种涉及气液系统的模拟中，例如自由表面流、雾化、液体射流以及液滴聚结和破碎等[6-8]。
最近，VOF方法也被考虑用于模拟带有相变的两相系统。相变带来了几个数值挑战，这些挑战主要源于气液界面处存在的膨胀项。这种膨胀导致了一个额外的对流通量，被称为斯蒂芬流，它在速度场中引入了一个强烈且局部的间断。因此，为了避免由于膨胀项导致的界面模糊，必须仔细解决界面和标量场的传输问题。
首次尝试在流体体积模型中加入相变的是Welch和Wilson[9]，他们修改了Juric和Tryggvason[10]以及Son和Dhir[11]提出的用于模拟沸腾流的前缘追踪模型，使其适应于VOF框架，以研究膜态沸腾构型[12,13]。随后，Hardt和Wondra[14]以及Kunkelmann[15]采用了类似的方法，他们使用了一种特定的技术来处理由相变引起的源项模糊，从而消除了界面单元的体积膨胀效应。因此，相变膨胀项不会干扰界面对流过程[15]。他们的沸腾模型被用于基准沸腾蒸发问题，如Stefan问题和膜态沸腾[14]，以及用于模拟成核气泡[15]。Schlottke和Weigand[16]提出了一种迭代方法，以在控制方程中施加速度跳跃条件。该方法用于模拟在高雷诺数和韦伯数下具有大变形的汽化液滴。Scapin等人[17]和Malan等人[18]的近期工作提出了类似的方法，通过求解额外的泊松方程，从单场速度中获得保守的对流速度。这些模型侧重于热诱导[18]或化学物质诱导的相变[17]，这简化了汽化速率的计算，并已在诸如Stefan问题和Scriven问题以及d²定律等基准蒸发问题上进行了测试。其他文献模型，如Palmore和Desjardins[19]以及Zhao等人[20]提出的模型，解决了耦合的热质传递问题，侧重于纯液体液滴，并忽略了气相和液相中普遍存在的多种化学物质。将这些相变模型扩展到多物质系统并不简单，并面临多个挑战：
1）多种化学物质增加了界面跳跃条件的非线性。气液两相的界面温度和质量分数是未知的，必须在一个非线性方程组中耦合能量和质量平衡来获得。
2）必须考虑速度场中的不连续性，并在气液界面处强制执行边界条件，以在气相和液相中传输标量场。
3）不能使用单一场方法来进行物质传输，因为轻质物质的蒸发必然导致重质物质的积累，必须遵守两个不同相中的质量守恒（图1）。
图1. 球形液滴中的物质质量分数分布图。绿线表示液相中轻质组分的质量分数分布，红线表示气相中相同轻质组分的质量分数，而蓝线表示液相中重质化学物质的质量分数分布。（有关图中颜色的解释，请参阅本文的网页版。）
在这项工作中提出的数值模型使用多组分相变公式解决了这些挑战，该公式适用于气相和液相中的任意数量的化学物质。本文提出的方法可以获得体积分数和标量场输运的扩展液体速度，也适用于静态液滴，这是最具挑战性的相变构型。对于静态液滴，速度场由Stefan流主导，必须仔细处理，以避免非物理界面传输。数学模型在第2节中描述，而第3节描述了输运方程的数值离散化。所得到的模型在Basilisk框架中实现[21]，并使用基准沸腾和液滴蒸发问题对其进行验证(第4节)。
图2. 通用域表示。改编自[22]。
图2描述了传输方程所描述的控制体积。该体积包括两个不相混溶的相（𝑉𝑙 和 𝑉𝑔，其中𝑉𝑔重复可能是个错误，应只保留一个以区分两相），它们被一个零厚度的界面Γ分隔开。在同一个控制体积𝑉内的不同相是通过在整个域（𝑉 = 𝑉𝑙 ∪ 𝑉𝑔）中引入一个标量标记函数来识别的，该函数被定义为海维赛德函数：
利用标记函数的定义，可以在整个控制体积上定义一个通用的物理属性𝜙：
其中，𝜙𝑙 和 𝜙𝑔 分别是液相和气相的相属性（即密度、粘度、热导率、比热以及化学物质的扩散率）。本工作假设在整个模拟过程中，气相和液相的属性是恒定的。因此，通用属性𝜙仅根据标记函数𝐻(𝐱,𝑡)变化。
2.1 相变界面平流
在计算上，标记函数使用颜色函数𝑐[23]来近似，该函数定义为控制体积𝑉中参考相（液相）的体积分数：
该函数在液相中取值为1，在气相中取值为0。根据VOF（Volume of Fluid）方法，界面的对流可以通过求解体积分数的守恒方程来实现。该方程可以从液相控制体积𝑉𝑙的质量平衡中推导出来：
其中，右侧项是由于相变而穿过界面的液相质量通量，它可以使用Rankine-Hugoniot关系[23]作为汽化速率的函数来表示：
这是一个移动的气液界面上的质量平衡。等式（4）右侧的最后一项𝑚̇̇是单位表面积上的总汽化速率。将方程（5）代入（4），并假设液相密度仅取决于𝑐，我们得到液相体积分数的传输方程：
等式右侧的项是根据相变现象得到的液相体积的源（或汇）。在给定液相速度场𝐮𝑙和单位表面积总蒸发率𝑚̇的情况下，可以求解该方程。项𝛿Γ仅在界面处应用源，并且定义为控制体积中的界面表面积除以体积本身，结果是面积除以体积的单位。
2.2 具有相变的不可压缩Navier-Stokes方程
通过求解相变不可压缩流的Navier-Stokes方程，得到了流体的速度场和压力场。
流场的不可压缩性在控制体积的两个相中是合理的，因为本工作感兴趣的系统的马赫数较低(Ma << 0.3)。因此，动量方程采用所谓的单场公式编写
这意味着在整个控制体积𝑉上，对在气液界面上急剧变化的物理性质进行积分（有关详细推导，请参阅[23]）。应力张量定义为𝝉 = 𝜇(∇𝐮 + ∇𝐮𝑇) − 𝑝𝐈，其中𝐈是单位张量，而材料属性𝜌和𝜇分别是密度和动力粘度，如方程（2）中定义。
将气相和液相的连续性方程(式(4))分别写成，将两式相加即可得到相变的单场连续性方程[9]。考虑不可压缩流动，利用界面质量平衡(式(5))，得到单场连续性方程:
该方程右侧是气液界面上的局部源项，从物理上讲，它源于相变过程中液相到气相（或反之）的转换伴随着体积的显著变化。这一过程导致了一个膨胀项，使得界面上的一场速度𝐮不是无散度的。
2.3 组分方程
化学物质通过求解守恒方程来传递，守恒方程用双场公式表示质量分数。使用这种方法，两个不同的控制方程求解相同的物质，一个为液相，一个为气相。选择这种方法是因为，当求解多种化学物质时，单一物质的质量分数在两相中呈反比变化。这种效应不能用带源项的单场公式来描述，源项只能增加或减少两相的质量分数。此外，场的分离限制了数值扩散在界面上的分布。通过对化学物质的质量对控制体积Vk的积分(𝑘=𝑙，𝑔)，推导出相𝑘中单个物质的组分方程:
其中，𝜔𝑖,𝑘 表示第 𝑘 相中第 𝑖 种物质的质量分数，而 D𝑖,𝑘 表示扩散系数。类似于体积分数输运方程，界面表面积的积分被界面质量平衡所替代，该平衡针对单一化学物质进行书写：
将方程（10）引入（9）中，相变贡献被收集到项𝑚̇𝑖中，这是对于单一化学物质𝑖的单位表面积的汽化率：
其中，体积分数𝑐𝑘用于在整个控制体积上对方程（9）进行积分。此方程分别对气相和液相进行求解，并最终使用基于颜色函数的算术平均值重建单场质量分数，以便进行可视化。
2.4 能量方程
使用双场公式，根据化学物质质量分数传输所采用的相同方法，对温度进行能量传输方程求解。温度方程可以从控制体积𝑉𝑘上的能量平衡导出，其中𝑘=𝑙,𝑔，并假设流体为牛顿流体，粘性效应和压力功项可忽略不计[24]：
图3. 气液界面处的质量和能量通量。𝑞̇𝐿和𝑞̇𝐺为导热通量，而𝑗̇𝐿和𝑗̇𝐺为物质扩散通量。
其中，𝑇𝑘是第k相的温度，𝑐𝑝,𝑘是定压热容，而𝑘𝑘是热导率。按照[20]中所述的相同方法，将方程(5)代入(13)，并对整个控制体积进行积分，双场能量平衡方程可以改写为：
式中𝑞Γ，𝑘为界面上的热传导，当界面温量值由跳变条件得到时，直接由两侧界面梯度计算得到。
2.5 多组分界面跳跃条件
方程(6)、(8)和(11)中出现的汽化率，以及方程(14)中的界面温度都必须根据界面跳跃条件进行计算。为此，在图3所示的气液界面Γ的一部分对应的控制体积上应用质量和能量平衡。界面处的质量平衡是针对每种化学物质编写的，它要求对流和扩散通量守恒。相反，界面处的能量平衡将传导通量之差等同于蒸发过程所需的能量[3]。因此，考虑到在零厚度界面内不允许质量和能量的积累，并假设界面处发生热力学平衡条件，得到以下方程组：
符号上方的帽子用于表示界面量：x̂i,g 和 x̂i,l 分别是化学物种 i 的气相和液相摩尔分数，ω̂i,g 和 ω̂i,l 分别是气相和液相的质量分数，而
是界面温度。相平衡通过拉乌尔定律计算得出，该定律将热力学平衡常数定义为气相和液相中界面摩尔分数的比值：
其中，P 是系统的热力学压力，
是通用物种 i 的蒸汽压，可以使用不同的气液平衡关系来计算。此处使用了安托万定律：
其中，A、B 和 C 是取决于所研究化学物种的经验参数。在这项工作中，热力学平衡常数是在理想液体混合物、理想气体和可忽略的Poynting修正的假设下获得的[25]。请注意，此公式可以轻松扩展到非理想混合物和更复杂的状态方程。
通过求解界面跳跃条件产生的非线性方程组，可以确定每种化学物种的蒸发速率 𝑚̇𝑖、界面质量分数（𝜔̂𝑖,𝑙 和 𝜔̂𝑖,𝑔）以及界面温度
。该方程组的解提供了表征界面和计算控制方程中源项所需的所有信息。
第2节中描述的控制方程在自适应笛卡尔网格上进行离散化。域被划分为多个方形单元（在三维中为立方体），并使用有限体积法和同位网格排列对控制方程进行数值近似。本文讨论了二维配置，其中Δ是通用单元的长度，下标f表示单元的通用边（或面），而上标n表示模拟时间步。此处使用的时间离散化遵循[26–28]中采用的相同形式，该形式假定体积分数和相关示踪剂（质量分数和温度）在n-1/2时间级别上定义，即比速度和压力场滞后半个时间步。将提出的数值方案扩展到三维是直接的。
3.1 流体体积的几何平流
使用VOF方法，在区域的每个细胞中计算第2节中描述的颜色函数，得到纯液体细胞(𝑐= 1)，纯气体细胞(𝑐= 0)和具有中间体积分数值的界面细胞。使用液相体积分数的平流方程(式(6))保守求解界面的平流。在这项工作中，采用了方向分裂几何VOF方案[29,30]。该方法适用于直角网格，在无散度速度场条件下是保守的。该方法包括两个步骤:重建步骤，即根据分段线性界面构造(Piecewise Linear interface Construction, PLIC)[29]在每个界面单元中对界面进行几何重建;通量计算步骤，即沿每个方向几何计算跨单元表面的液体体积通量，并在下一个时间级别更新每个单元中的体积分数:
方程（18）、（19）中的F和G是跨面的几何通量，而右侧积分是膨胀项，用于纠正由维度分裂引起的散度误差。具体实施细节可参见[22,23,29]。
为了避免界面模糊，用于界面平流的速度𝐮Γ必须在界面上保持连续。因此，我们不能使用单场速度𝐮来近似𝐮Γ。相反，界面速度是通过方程（5）获得的：
其中𝐮𝐸,𝑓是一个无散度的液体速度，其计算方式如第3.3节所述，而𝐮𝑝𝑐是由于相变引起的速度贡献，它会使液相收缩或膨胀。引入界面速度的这种定义允许将方程（6）中的相变源项直接包含在平流步骤中。这一操作并不直接，因为界面梯度是在单元中心计算的，蒸发率𝑚̇也是如此。因此，需要通过插值将蒸发率转换为在单元面上定义的界面回归速度𝐮𝑝𝑐。所提出的插值方案如图4所示：
图4. 界面回归速度从单元中心到单元面的插值（来源：http://basilisk.fr/sandbox/ecipriano/test/interfaceregression.c）。
1. 如果该面将一个界面单元与一个纯单元连接起来，则该面上的界面回归速度取自界面单元，并由相应的界面法向分量加权。
2. 如果表面连接两个界面单元，则通过两个连续蒸发速率之间的线性插值计算界面回归速度，并由界面法向加权。
使用这种方法，得到的界面平流速度𝐮Γ以非保守的方式传输体积分数，但添加或移除的液体体积量与相变现象的效果相匹配。与其他文献方法[18,31]相比，这种方法的优势在于平流和相变是在一个步骤中应用的，从而限制了相变作为显式源项时容易出现的欠射和过射的可能性。
3.2 压力速度耦合
Navier-Stokes方程采用近似投影法进行离散化，其中速度和压力都位于单元中心。Popinet在[26]和[27]中详细描述了这种方法及其在四叉树/八叉树网格中的数值实现。在预测步骤中，动量方程的平流-扩散部分的求解是通过将平流步与扩散步分开来进行的。平流部分采用Bell-Colella-Glaz方法[28]进行估算，该方法是一种二阶未分裂迎风格式，在CFL条件小于1时稳定[27]。而扩散部分则被重新排列为Poisson-Helmholtz方程，并使用多重网格求解器以隐式方式求解[32]。该方法与Crank-Nicholson离散化相结合，具有二阶精度，并且在CFL数大于1时也能保持稳定[27]。时间交错方案简化为以下离散化：
其中下标𝑐→𝑓表示变量从单元中心到单元面的线性插值，反之亦然（𝑓→𝑐）。
一旦计算出临时速度ρ *，投影步骤允许得到一个压力场，这保证了𝑛+ 1时刻的散度符合连续性方程。本文对投影步骤进行了修改，加入了体积展开项，使得速度场散度非零(式(8))。投影步骤简化为泊松方程，使用多网格求解器求解[32]:
其中𝛿Γ近似为离散界面面积（通过PLIC重构获得）与单元体积之比。在投影步骤之后，对同位速度进行重构，加入加速度项和压力梯度：
表面张力f𝜎=𝜎𝜅∇𝑐是通过将高度函数法用于曲率𝜅，以及对表面张力和压力梯度的良好平衡离散化来计算的，这可以恢复静态液滴的平衡解（拉普拉斯方程）[27,33]。在本研究中，我们忽略了相变（反冲压力）引起的额外压力跳跃。根据Mialhe等人的说法[34]，该项通过抵消相变引起的动量损失来强制动量守恒。虽然该项在保持动量的方案中很重要，但其贡献并不影响本文讨论的结果。
3.3 平流速度
在使用VOF方法模拟相变时，主要问题之一是获得扩展速度，该速度代表液相的速度。当使用单场公式对动量和连续性方程进行离散化时，界面单元会受到Stefan流贡献的影响。这意味着气相速度会干扰界面附近液相体积分数的传输。为了克服这个问题，文献中提出了不同的解决方案[14,15,17–19,31]，但尚未找到适用于每个具有相变系统的唯一解决方案。最近的方法基于求解速度势𝜙的附加Poisson方程：
该势用于分离Stefan速度𝐮𝑆的贡献：
因此，可以通过从场速度中减去Stefan速度来获得扩展速度场，从而得到一个从构造上来说是无散度的扩展速度场[18]：
该方法的不同变体被采用。Scapin等人[17]在整个域上求解方程（26），而Malan等人[18]在靠近界面的狭窄单元格带中求解相同的方程。Palmore等人[19]将基于PDE的Aslam外推技术与方程（26）的解相结合，以强制液体速度的连续性，并在外推过程后校正速度散度的误差。Gennari等人[31]提出了一种不同且更新的方法，他们将投影步骤中的扩展项移向靠近界面的纯液体或纯气体单元。这种方法不需要求解额外的Poisson方程，并且在沸腾模拟中表现出了良好的性能[31]。
所有这些方法的一个共同问题是，对于具有强密度比的静态液滴蒸发问题，它们会导致非物理的界面变形。为了克服这个问题，我们提出了一种新方法，该方法基于Gennari等人[31]提出的移动方法和一种新的双压力-速度耦合技术的组合。根据这种方法，将方程（8）中的扩展项移向最近的纯单元，使得得到的速度𝐮包含Stefan流，该流在界面单元上是连续的。在这一步之后，求解第二组Navier-Stokes方程，但不包括扩展项：
求解这组Navier-Stokes方程以找到扩展速度𝐮𝐸，该速度用于界面对流过程。这组方程在整个域上求解，并使用与速度和压力𝐮和𝑝相同的边界条件。采用相同的时间交错离散化方法：
在模拟开始时，𝐮和𝐮𝐸被初始化为相同的值。在时间积分循环中，这两个速度独立求解，但使用相同的离散化方法。唯一的区别在于投影步骤：虽然速度𝐮的投影包含体积膨胀项，并移向最近的纯单元（方程（23）），但速度𝐮𝐸的投影步骤强制执行无散度约束（方程（33））。积分顺序在算法1中进行了总结。通过这种方法，我们获得了两个不同的速度场，𝐮和𝐮𝐸，它们都在界面上连续且无散度，分别用于界面对流和标量场的传输。该方法将Stefan流的影响与界面对流解耦，并在Stefan速度远大于背景速度的情况下（即具有强密度比的静态液滴情况）证明是有效的。
3.4 标量场的输运
标量场的输运方程采用算子分裂法求解。首先，执行对流过程的求解，然后考虑相变贡献来求解扩散。根据本工作中实现的双场方法，在液相和气相中定义标量场（𝑠）的示踪剂形式（𝑡𝑟）是很方便的：在𝑠𝑙中为𝑡𝑟𝑙，在𝑠𝑔中为𝑡𝑟𝑔
除了界面单元外，示踪剂场的值与标量场值完全相同。使用两个单独的示踪剂，可以限制数值扩散，并且很容易在标量场和示踪剂场形式之间切换。在转换过程中，体积分数值较小（𝑐<𝜖）的单元被视为体积分数为零，并且这些单元中相关联的示踪剂场也被设置为零；这里使用了Wenzel和Arienti[36]使用的相同容差（𝜖 < 1 × 10−10），因为这个值足够小，可以避免质量守恒问题。为了后处理的目的，可以简单地通过两相中示踪剂场的和来恢复单场公式：𝑡𝑟 = 𝑡𝑟𝑙 + 𝑡𝑟𝑔。
1）对流运输
方程(11)和(14)的对流部分通过在域的单个单元上积分方程来求解，假设物理性质恒定，并使用体积分数𝑐对积分体积进行校正，以便只考虑正在求解的标量场相𝑘的体积:
图5. 具有相变的扩散项的离散化。
该方程使用与体积分数对流方程（第3.1节）相同的方向分裂过程进行离散化。按照这种方法，使用以下方程考虑对流传输来更新标量场的值：
其中𝐹𝑓是待求解相在通用面𝑓上的体积通量，𝑠𝑛𝑘,𝑓是单元通用面上的标量场值。该值使用二阶迎风Bell-Colella-Glaz方案[28]获得。体积通量使用与VOF场（第3.3节）相同的几何通量计算方法进行计算。对于液滴蒸发问题，气相通量使用包含Stefan流的速度𝐮进行计算，而液相通量则利用扩展速度𝐮𝐸进行计算。在对流步骤之后，每个单元中标量场的值被重构为：
2）扩散输运和相变
求解扩散过程时包括表示相变的源项。在有限体积离散化的背景下，将物质方程（方程（11））的扩散部分重写为积分形式，忽略对流项并假设材料属性是恒定的：
该方程的时间导数使用反向欧拉方案进行离散化，而面梯度则使用共享同一面的两个连续单元之间的有限差分来近似：
该方程被重新排列成泊松-亥姆霍兹方程，并使用多网格求解器[32]以隐式方式求解，如下所示：
其中，相变贡献包含在等式右侧的最后两项中。源项描述了蒸发的物质i的量𝑚̇ 𝑖（显式源），以及由于相变引起的体积变化（隐式源）。由于显式源项的存在，扩散方程的解并非完全隐式，因此可能需要限制时间步长以确保稳定性。准确求解方程（43）需要表面分数值𝑐𝑘|𝑓，该值通过PLIC界面重构[37]进行几何计算得出。使用这种方法，扩散过程被限制在待求解的特定相中，并且简单地忽略穿过空面的扩散通量，而不引入人工扩散（图5）。
对温度方程（方程（14））的扩散部分进行离散化时，也采用了相同的方法，从而得到以下形式：
这可以根据界面温度和沿界面的温度梯度来求解，计算方法将在下一节中解释。
3.5 界面跳变条件下的汽化速率
通过求解每个界面单元中的非线性方程组（方程（15））的界面跳跃条件，可以量化界面属性（温度和质量分数）以及每种化学物质的蒸发速率。为此，开发了一种用于非线性系统解耦求解的策略，以找到可靠的第一猜测值，然后将其用于求根算法以求解非线性系统。通过假设两个相中之一的演变通常较慢，可以以解耦的方式解决该系统。在液滴蒸发过程中，由于扩散率较低，液相中物质质量分数的演变比气相慢。因此，假设液相界面质量分数和温度等于界面单元处场的值。下面报告了系统解耦求解的过程。在每个时间步开始时计算每种化学物质的汽化速率和界面属性，因此使用第𝑛 − 1∕2时间级的体积分数和标量场。
1）在每个界面单元中，界面温度和液相质量分数计算如下：
其中𝜔̂i,l是界面单元周围3×3模板内液相质量分数的体积分数加权平均值。
2）使用质量分数计算界面摩尔分数𝑥̂i,l，并使用前一点中的界面温度计算热力学平衡常数，从而获得气相中的界面摩尔分数：
最后，重新构建了气相中的界面质量分数𝜔̂𝑖,𝑔。
3）一旦知道界面质量分数，就可以根据气液界面的质量平衡计算总蒸发速率，并对液相中的总物种数（𝑁𝐿𝑆）进行求和:
其中，根据定义[24]，液相中所有液相物质的扩散通量之和必须等于零，而同一相中的质量分数之和必须接近1。因此，得到了总蒸发速率的显式方程：
这使得可以显式计算每种化学物质的汽化速率：
4）最后，根据蒸发速率的知识，调整惰性（仅气体）物质的质量分数，并使用求根算法细化界面温度值：
使用此过程，即可获得界面跳跃条件值。然后将此初步猜测解提供给非线性方程组，该方程组通过耦合解来完善这些值。这种方法允许将界面上的传递现象与域中的其余部分解耦。然后将从跳跃条件获得的信息以源项的形式包含在传递方程的解中。
图6. 用于计算界面梯度的嵌入式边界方法。改编自[38]。
为了正确求解界面跳跃条件，必须准确地进行界面梯度计算。在这项工作中，我们利用在笛卡尔网格中为复杂固体边界开发的方法（如Johansen和Colella提出的嵌入式边界方法[38]）来计算界面梯度。采用这种方法，利用PLIC界面重构将气液界面近似为固体边界，并利用Fleckenstein和Bothe提出的VOF平均二阶方案计算界面梯度[39]：
其中，𝑓̂是界面处一般标量场的值，而𝑓1和𝑓2是该场沿界面法线方向两点（图6）的值，使用双二次插值从相邻单元中插值得到[31]。
3.6 时间稳定性和解决方案概述
数值求解算法总结在算法1中。首先，在时间步的初始时刻求解界面跳跃条件来计算汽化速率，以获得𝑚̇𝑖𝑛−1∕2和所有其他界面属性。利用这些信息，以及在时间𝑛时已知的速度场，执行VOF界面平流以将界面推进到下一个时间级别𝑐𝑛+1∕2，并通过算术平均值更新材料属性。使用双场方法求解组分和温度方程，在求解标量传递方程后，获得速度和压力场。数值解的时间步长由表面张力力时间显式离散化的稳定性所需的时间步长控制[33]。但是，由于组分和温度方程中存在显式源项，因此在强蒸发率的情况下，可能需要进一步减小时间步长。
上一节描述的数值模型在Basilisk框架[21]中实现。选择该代码是因为它实现了有效的自适应网格细化系统，以及对表面张力的良好平衡离散化，这对于小液滴和气泡来说非常重要。本工作中开发的模型通过简单配置中的基准测试用例进行了验证，这些测试允许将数值结果与解析解进行比较。更复杂的配置与基准数值模拟进行了比较，讨论了系统的定性行为和质量守恒。不同的测试用例将在接下来的部分中报告，而本工作中开发的所有代码和模拟设置都可在Basilisk沙盒[40]中免费获得。
4.1 固定流量蒸发
将一个二维液滴放置在方形域的中心。液滴以固定的汽化速率蒸发，模拟一直进行到液体完全消耗为止。通过对液相进行简单的质量平衡，可以计算出液滴半径随时间变化的分析解：
在这个测试用例中，施加𝑚(𝑚)使液滴在4秒内完全消耗。本测试用例采用Malan等人[18]报告的物性:𝜌𝑙∕𝜌𝑔= 2，𝑚∕𝑚∕𝜌𝑙= 0.05 m s-1，𝑅0 = 0.23 m。与Malan等[18]不同的是，本次模拟完全去除了表面张力，重点关注相变，避免了对液滴形状的调整。方形区域具有单位长度，施加了流出边界条件，忽略了重力。模拟在5个不同的细化级别上进行，从4到8，这意味着，对于级别4的情况，沿着每个域方向的细胞数量等于24。图7 (a)为液滴体积随时间的变化趋势。模拟结果提供了一个很好的近似解析解，除了最粗糙的细化水平。恢复了收敛趋势，在液滴体积上计算的相对误差显示为二阶收敛速率。图7 (c)显示了三种不同模拟时间(1秒、2秒、3秒)下PLIC界面与液滴形状解析解的对比。尽管没有表面张力，液滴在整个生命周期内都保持球形。从这个测试案例中，我们可以得出结论，双压力-速度耦合的使用为VOF输运方程提供了可靠的速度。这个速度不包含不连续，因此，它不会导致界面涂抹。此外，界面回归速度(式(20))的引入使得液滴收缩过程可以在单个vof平流步骤中以二阶精度进行模拟。Malan等人[18]没有报告该测试用例的收敛率。然而，这两种模型都能很好地模拟，得到的结果与理论解很接近。
图7. 固定流量液滴蒸发结果。液滴体积随时间的变化：不同细化级别的数值模拟与解析解的比较（a）；基于液体体积相对误差计算的收敛速度（b）；模拟时间1秒、2秒和3秒时的液滴球形度（c）（http://basilisk.fr/sandbox/ecipriano/run/fixedflux.c）。
4.2 Stefan问题
Stefan问题是经典的基准测试用例，不同的作者已经使用它们来验证他们的相变模型。这里考虑的Stefan问题包括由气相和液相之间的温度梯度引起的液体平面的蒸发，液相保持在饱和温度。在模拟开始时，一个过热壁加热蒸汽层，使其温度高于液相，导致相变。该系统的分析解由[18]报道，它描述了蒸汽层厚度𝛿(𝑡)随时间的变化以及温度分布。对于这次模拟，我们采用了Malan等人[18]使用的相同物理属性：𝜌𝑙 = 958 kg m-3，𝜌𝑔 = 0.6 kg m-3，𝜎 = 0.059 N m-1，𝜇𝑙 = 2.82 × 10−4 Pa s，𝜇𝑔 = 1.23 × 10−5 Pa s，𝑘𝑙 = 0.68 W m-1 K-1，𝑘𝑔 = 0.025 W m-1 K-1，𝑐𝑝𝑙 = 4216 J kg-1 K-1，𝑐𝑝𝑔 = 2080 J kg-1 K-1，Δℎ𝑒𝑣 = 2.256× 106 J kg-1。模拟的特征是Jacob数为29.84，并使用边长为0.01米的正方形域进行（图8）。蒸汽层厚度的初始值为322.5微米，这样选择是为了在模拟开始时，即使在最粗糙的细化级别下，也至少存在一层蒸汽单元。在域的顶部和底部施加对称边界条件；在域的右侧施加入口零速度条件，温度固定为383 K；而在左侧壁上施加流出边界条件。模拟在三个不同的细化级别上进行：5、6和7。
图8. 10秒时温度场图及界面处速度跳跃的详细情况。白线是体积分数场的0.5等值线。
图9.汽相界面厚度(a);界面厚度相对误差(b);不同精炼水平下气体层温度分布(c) (http://basilisk.fr/sandbox/ecipriano/run/stefanproblem.c)。
图8显示了总模拟时间为10秒时域内的温度场。在气层中，温度从过热壁上的边界条件变化到气液界面处的饱和温度。界面附近区域的放大图显示了由于局部体积膨胀引起的强烈速度跳跃。这种速度跳跃通过双重压力-速度耦合技术得到了正确处理，该技术允许标量场和体积分数通过界面处无散度的速度进行传输。图9中的曲线图显示了在不同细化级别下模拟结果与解析解的比较。即使在最粗糙的网格上，蒸汽层的厚度也能很好地逼近解析解。对不同细化级别下蒸汽层厚度的相对误差进行分析，发现即使在最粗糙的细化级别下，相对误差也约为0.05%，并且随着网格分辨率的增加，相对误差趋于以二阶精度收敛。与Malan等人[18]的比较表明，使用本工作描述的模型获得的相对误差较小，并且比[18]中报告的大约一阶收敛速度更快地收敛到解析解，[18]中的相对误差分别为细化级别6、7和8下的0.64%、0.39%和0.23%。更快的收敛速度可能是由于对流速度的计算和界面梯度的更精确评估相结合的结果。
4.3 Epstein-Plesset问题
这个测试用例描述了在不饱和液体环境中球形气泡的溶解过程。相变过程是由气泡界面与液体主体之间的化学物种浓度梯度驱动的。这里考虑的是纯扩散条件，即不考虑由于Stefan流引起的气相中的对流传输。Epstein和Plesset[41]在准静态条件下获得了气泡半径随时间变化的分析解，即评估汽化速率时忽略了界面运动的影响。在这一假设下，浓度场在任何时间瞬间都满足稳态扩散方程。如果扩散过程的特征时间远小于界面回归的特征时间，则这种近似是合理的[41]。从质量平衡的角度来看，可以得到气泡的半径为：
其中，
和 𝐶𝑏𝑢𝑙𝑘 分别是化学物种在界面处的浓度和液体主体相的浓度。浓度分布可以通过径向坐标下物种扩散方程的解析解获得：
模拟设置借鉴了Farsoiya等人的工作[42]，将两相的密度设置为相同的值，等于1 kg m-3，以便从连续性方程中移除体积膨胀项，从而也移除了Stefan对流。化学物质的扩散率设置为1 m2 s-1，而表面张力𝜎等于0.1 N m-1。在界面气相一侧，浓度设置为固定值，等于0.8 mol m-3，斯坦顿数等于8 × 10-4[42]。模拟在二维轴对称域中进行，气泡在域的左下角初始化，使用最大细化级别为10的自适应网格。在液相边界处，浓度设置为体相值（𝐶𝑏𝑢𝑙𝑘 = 0）。图10（a）中的结果显示了蒸汽气泡的消耗动态，以气泡的归一化半径随时间的变化来表示。在整个气泡寿命期间，这一趋势很好地近似了Epstein-Plesset解。图10（b）中的时间浓度分布图是通过在域中的一个单点（对应于𝑅0 + 0.2 m）采样化学物质浓度的值来绘制的。从这一分布图中可以观察到准静态近似的有效性，因为液相中的浓度从体相值快速过渡到由相变现象确定的浓度值，这是一个非常短暂的瞬态过程。将本工作与[42]中报告的模拟进行比较表明，两个模型都能够准确地预测气泡半径和化学物质浓度的理论分布，除了由于方形域和准静态假设引起的小偏差外。Farsoiya等[42]提出的模型在计算时间方面具有更好的性能，而本工作中描述的模型需要更小的时间步长来保证扩散步长的稳定性(第3.4节)。尽管有这个缺点，这里开发的模型是一个更一般的公式:它既可以用于斯特凡对流起主要作用的蒸发问题，也可以用于纯扩散条件，例如这个测试用例。
图10. Epstein Plesset问题(http://basilisk.fr/ sandbox/ ecpriano /run/epsteinplesset.c)的归一化直径剖面(a)和坐标𝑅0 + 0.2 m处的化学物质浓度(b)。
4.4 Scriven 问题
Scriven测试用例分析了在过热环境中膨胀的气泡。相变过程的驱动力是气泡与体相界面之间的温度梯度。EpsteinPlesset案例的主要区别在于，Scriven[43]提出的解决方案包含了Stefan对流。解析解在[43,44]中有报道，它描述了气泡半径随时间的演变，以及温度场。
图11. Scriven测试案例的温度场和速度场图。该图展示了在模拟的最高细化级别和最后一个时间点上的结果（来源：http://basilisk.fr/sandbox/ecipriano/run/scrivenproblem.c）。
图12. 气泡半径随时间变化的曲线图（a）；气泡半径的相对误差（b）；在四个不同细化级别下，沿域半径的温度分布图（c）（来源：http://basilisk.fr/sandbox/ecipriano/run/scrivenproblem.c）。
这个测试用例在二维轴对称配置上进行模拟，其中初始半径为1毫米的气泡被放置在一个长度为12毫米的正方形域中（图11）。模拟在四个不同的细化级别上进行，从6到9。使用自适应网格细化来最小化模拟的计算成本，模拟运行了0.5秒，对应于最终气泡半径等于初始半径的两倍。与解析解不同，初始半径不能设置为零，因为这将需要模拟气泡成核过程，而这超出了本工作的目的。模拟设置借鉴了[44]，其特点是大密度比，𝜌𝑙/𝜌𝑔 = 1623，以及Jacob数等于3。以下是使用的物理性质：𝜌𝑙 = 958 kg m-3，𝜌𝑔 = 0.59 kg m-3，𝜎 = 0.001 N m-1，𝜇𝑙 = 2.81× 10-4 Pa s，𝜇𝑔 = 1.26× 10-6 Pa s，𝑘𝑙 = 0.6 W m-1 K-1，𝑘𝑔 = 0.026 W m-1 K-1，𝑐𝑝𝑙 = 4216 J kg-1 K-1，𝑐𝑝𝑔 = 2034 J kg-1 K-1，Δℎ𝑒𝑣 = 2.256 × 106 J kg-1。
图12（a）显示，尽管在参考解和最粗糙细化级别的数值结果之间存在位移，但当增加网格分辨率时，气泡半径的动态变化趋于收敛到解析解。数值结果与Scriven解之间的气泡半径相对误差以大约二阶精度收敛。图12（c）绘制了沿径向坐标的温度分布图，它表明从数值解中获得的轮廓也收敛到分析轮廓。正如Tanguy等人[44]所报道的，对于这个测试用例，获得单调收敛解取决于用于获得用于传输界面的扩展速度的方法以及温度场的正确解。Tanguy等人[44]使用Level Set方法对其最有效的方法获得的误差分别为9.5%、1.4%和1.0%，分别对应于7、8和9级。将这些结果与使用本文中描述的模型获得的结果进行比较表明，气泡半径的相对误差是可比较的。在不使用计算成本高昂的外推程序的情况下，在这个测试用例上获得收敛是一个显著的成就：尽管表面张力值较低，但气泡始终保持球形，并且获得了球对称温度分布以及平滑的速度场。
4.5 等温液滴蒸发
根据Pathak等[45]的数值结果，对多组分液滴的蒸发进行了测试。这些结果考虑了一个纯液滴，由于界面和气相之间化学物质浓度的梯度而蒸发。假设界面温度恒定，则界面质量分数和气液界面蒸汽压也恒定。[45]中提出的基准模型采用球面对称和非定常条件。气相和液相用两个贴体网格来描述，它们由界面点连接。由于Stefan流的额外传输包括在内，并且没有使用准静态近似。
图13. 在三个不同模拟时间点的质量分数分布图。白色线条代表界面，是体积分数场的0.5等值线（来源：http://basilisk.fr/sandbox/ecipriano/run/pureisothermal.c）。
图14. 归一化液滴直径平方的时间演变（a）；在三个不同模拟时间点上气相质量分数场的半径分布图（b）。图中的线条是本文描述的模型产生的数值结果，而标记则是Pathak等人[45]给出的基准结果（来源：http://basilisk.fr/sandbox/ecipriano/run/pureisothermal.c）。
一个初始直径为0.4毫米的液滴被放置在正方形域的左下角，该正方形域的尺寸等于液滴直径的两倍（图13）。本模拟使用了以下物理属性：𝜌𝑙 = 10 kg m-3，𝜌𝑔 = 1 kg m-3，𝜎 = 0.01 N m-1，𝜇𝑙 = 1 × 10-4 Pa s，𝜇𝑔 = 1 × 10-5 Pa s，𝐷𝑙 = 2 × 10-3 m2 s-1。尽管这个测试用例是为纯液体液滴设计的，但在这项工作中，它使用了多组分模型进行了运行，其中包含了四种不同的伪物种。按照这种方法，每种化学物质在液相中的质量分数初始化为0.25，在气相中为0。假设所有物种具有相同的物理属性，则液体液滴的行为应等同于纯液滴的行为。
在与液滴接触的边界处设置了对称边界条件。在上边界和右边界采用零压力流出条件，而化学物质的质量分数采用齐次Dirichlet边界条件求解。界面气相侧的质量分数恒定，为0.667，模拟在5 ~ 8四个不同的细化级别上运行。图14为基准测试结果与本文所述数值模型的对比。结果表明，平方直径衰减的动力学与d2定律描述的稳态行为有很大不同，这表明平方直径衰减的变化是恒定的。可以确定三个不同的区域:在第一次瞬态(𝑡< 0.15 × 10−4 s)期间，由于模拟开始时界面和气相之间的强梯度，液滴经历了快速消耗。在𝑡> 0.15 × 10−4 s和𝑡< 0.7× 10−4 s之间，气相质量分数的浓度随着消耗速率的增加而增加，直到𝑡> 0.7× 10−4 s达到稳态汽化速率。随着液滴的收缩，该模型得到的结果与基准溶液之间的位移增大。考虑到液滴每直径的细胞数量随着液滴的消耗而减少，这并不奇怪。实际上，细化级别越高，平方直径衰减越收敛于基准解决方案。这项工作的结果与Pathak等[45]得到的剖面吻合得很好，并且很好地捕捉到了液滴消耗的瞬态。图14 (b)报告了在三个不同的模拟时间瞬间的最高精细化水平的物种质量分数分布图。数值模拟得到的趋势与Pathak等[45]得到的剖面较为接近。由于采用了带有伪组分的多组分模型，因此图14 (b)中绘制的质量分数曲线为每个蒸发组分的气相质量分数之和。用这种方法，总质量分数与纯物质的质量分数一致。
4.6 非等温液滴蒸发
图15.正庚烷质量分数和温度场图 (http://basilisk.fr/sandbox/ecipriano/run/c7pathak.c)。
在非等温环境中模拟了一个液体液滴，这意味着界面质量分数不是像第4.5节那样为常数，而是界面温度的函数。该液滴是纯的，系统压力为28.6bar，且忽略重力。液滴的初始直径为5微米，而域长度（二维轴对称）是液滴直径的两倍，遵循[45]中报告的相同设置。由于问题的复杂性和所涉及的众多物理现象，因此无法获得解析解。因此，使用本工作中开发的模型得到的模拟结果与Pathak等人[45]提供的数值基准结果进行了比较。模拟在三个不同的细化级别上进行，从6到8。为了模拟氮气中高压下正庚烷的性质，选择了该模拟的物理性质：𝜌𝑙 = 626.7 kg m^-3，𝜌𝑔 = 17.51 kg m-3，𝜎 = 0.01 N m-1，𝜇𝑙 = 1 × 10-4 Pa s，𝜇𝑔 = 1 × 10-5 Pa s，𝑘𝑙 = 0.1121 W m-1 K-1，𝐷𝑔 = 6.77 × 10-7 m2 s-1，𝑘𝑔 = 0.4428 W m-1 K-1，𝑐𝑝𝑙 = 2505 J kg-1 K-1，𝑐𝑝𝑔 = 1053 J kg-1 K-1，Δℎ𝑒𝑣 = 3.23 × 105 J kg-1。正庚烷和氮气的分子量分别为100和29 kg kmol-1。液滴的初始温度为363 K，而环境温度为563 K。在液滴内部，正庚烷的质量分数为1，在气相中为零，因为模拟开始时气相中只存在纯氮气。在域的左侧和底部施加对称边界条件，而在顶部和右侧对速度和压力施加流出边界条件。顶部和右侧的正庚烷质量分数条件为Dirichlet边界条件，质量分数设置为0，而气相温度则设置为一个恒定的值，等于初始气体温度。为了与[45]中用于评估汽液平衡的热力学保持一致，使用来自NIST数据库[46,47]的系数，通过Antoine方程计算气液界面处的热力学平衡。图15显示了正庚烷质量分数和温度场的演变。在模拟开始时，由于液滴较冷，蒸发速率较小。经过一段时间后，环境与液滴之间交换的热量加热了液相，提高了温度，从而增加了液滴的蒸发速率，液滴开始被消耗。图16展示了整个模拟过程中正方形直径的衰减和界面温度分布。与等温蒸发 情况（图14）不同，模拟开始时正方形直径分布不会突然衰减，因为在这种情况下，界面质量分数不是固定的，并且随着液滴从热环境中吸热而增加。随着模拟的进行，界面温度增加，直至接近液滴的湿球温度，这是由环境加热和蒸发过程（使界面冷却）之间的相互作用造成的。使用本工作中开发的模型获得的模拟结果与Pathak等人[45]提供的数值结果（包括正方形直径衰减和界面温度）一致。图16还报告了三个不同模拟时间步长的温度和正庚烷质量分数的径向分布。除了模拟结果与基准结果之间的良好一致性外，还观察到液滴的温度增加并在液滴半径范围内保持不变。同时，正庚烷的界面质量分数随时间增加，并从液滴内部的质量分数（等于1）开始，质量分数降至界面质量分数，并遵循气相中的基准分布。
4.7 静态多组分液滴蒸发
研究了二元液滴的等温蒸发，假设液滴由两种具有相同物理性质但挥发性不同的化学物质组成，以便关注这两种物质的不同蒸发速率。最初，两种化学物质的质量分数均等于0.5，而气相中存在一种惰性化合物。初始液滴直径为400微米，域的长度是液滴直径的4倍，忽略重力，并使用了以下物理性质：𝜌𝑔 = 1 kg m-3，𝜎 = 0.03 N m-1，𝜇𝑙 = 1 × 10-4 Pa s，𝜇𝑔 = 1 × 10-5 Pa s，𝐷𝑙 = 4 × 10-6 m2 s-1，𝐷𝑔 = 8 × 10-5 m2 s-1。液相中两种化学物质的相对挥发性等于2。特别是，对于轻组分和重组分，其蒸气压与热力学压力之比分别为0.8和0.4。密度比等于10，以最小化模拟所需的计算时间，从而能够在三个不同的细化级别（7、8和9）上研究液滴消耗趋势和质量守恒。液相中的扩散率经过调整，使得𝜌𝑙和𝐷𝑙的乘积能够代表类似正庚烷的烃类物质的值。在域的左侧和底部边缘施加对称边界条件，而在顶部和右侧对速度、压力和质量分数施加流出边界条件。
图17. 静态二元液滴蒸发的质量分数和网格细化的演变：轻组分（左上角）、重组分（右上角）、惰性组分（右下角）（http://basilisk.fr/sandbox/ecipriano/run/staticbi.c）。
图18. 静态二元液滴的方形直径衰减（a）；气相中每种组分的质量与液相中相同物种的初始质量之比（b）；液相中每种组分的质量与液相中相同物种的初始质量之比（c）。未明确指出的地方，数值结果是在最高细化级别下获得的（http://basilisk.fr/sandbox/ecipriano/run/staticbi.c）。
图17为轻、重、惰性化学物质质量分数随时间的演化情况。轻组分首先开始蒸发，其在液滴内的质量分数在靠近气液界面处优先减小。重物质积累，增加了其在液相中的质量分数，因此，其汽化速率也随时间而增加。与纯等温液滴(第4.5节)不同，由于液相的组成变化，二元液滴的界面质量分数不是恒定的。
图18显示了液滴的消耗动态。平方直径衰减的演变与等温纯液滴蒸发试验类似。在三个不同的细化层次上，数值解收敛于同一个解。图18 (a)中的标记不代表解析解，本测试用例无法获得解析解;相反，他们绘制了从这项工作中开发的模型得到的结果，但没有非线性方程组的解。在这些条件下，可以观察到，求非线性系统的第一猜测值的过程与求非线性系统在界面处的解几乎相同。这证明了该算法能够找到第一猜测的界面值，从而提供可靠的结果，从而导致液滴的消耗动力学相同。图18也报道了质量守恒试验。如果蒸发的每种化学物质的量在气相中被回收并从液相中除去，则蒸发模型被称为保守模型。实际上，蒸发物种的总量:
与液相中同种物质的质量比较:
并与气相中物质的质量比较：
其包括跨越边界的对流质量通量和扩散质量通量的积分。图18 (b)和(c)分别给出了蒸发组分组分的质量与气相和液相中相同组分的初始质量之比。图中表明，相对于重质组分，轻质组分在气相中的质量迅速增加，液相的消耗也遵循相同的动力学规律。轻分量的蒸发质量𝑚E𝑣𝑎𝑝和𝑚𝐿i𝑞之间的差是用相对误差来估计的，定义如下:
对于三个细化级别，质量守恒误差分别为0.186、0.047和0.017。随着网格分辨率的增加，误差得到了改善。
4.8 对流多组分液滴蒸发
研究了在强制对流条件下的二元液滴。液滴的初始直径为400微米，并放置在长度为液滴直径36倍的正方形区域（二维）上。最初，液滴放置在距离左壁6倍初始直径和距离底部18倍初始直径的位置。液滴由两种具有相同物理性质但挥发性不同的化学物质组成，并采用与第4.7节中相同的相对挥发性。使用四叉树离散化进行了从9到12四个不同细化级别的模拟，这对于优化此测试案例的计算时间至关重要。虽然忽略了重力，但从域的左侧注入了气体流量，以在模拟的初始属性下达到Re = 160。因此，在左侧施加了速度和压力的入口边界条件，而对右侧壁使用了流出边界条件。使用以下物理属性进行了此模拟：𝜌𝑙 = 800 kg m-3，𝜌𝑔 = 5 kg m-3，𝜎 = 0.073 N m-1，𝜇𝑙 = 1.138 × 10-3 Pa s，𝜇𝑔 = 1.78 × 10-5 Pa s，𝐷𝑙 = 1.4 × 10-7 m2 s-1，𝐷𝑔 = 1.25 × 10-5 m2 s-1。在静态参考系中运行模拟，直到液滴到达域的右侧（对应于初始液滴体积的近20%被消耗）。
图19. 轻组分质量分数的演变，以及三个不同时间点的液滴内部轻组分和重组分质量分数的详细情况（http://basilisk.fr/sandbox/ecipriano/run/forcedbi.c）。
图20. 在四个不同细化级别下，方形直径衰减的演变（a）；在最大细化级别下，域内每种物种的液体质量变化（实线）与归一化到初始液体质量的蒸发物种总量（标记）之间的比较（b）（http://basilisk.fr/sandbox/ecipriano/run/forcedbi.c）。
图19（a）显示了三个不同时间点的轻组分质量分数。模拟开始时，轻组分蒸发得更快，蒸发的物质在强制对流流量的作用下被输送到域的右侧。所选的雷诺数范围最终导致形成类似卡门涡街的结构，这可以从输运的质量分数场中观察到。图19显示了液滴内两种物质的质量分数场的演变。通过数值模型可以预测化学物质的非球对称分布，观察到涡旋结构在液相内部引起再循环。这一现象直接影响热质传递和液滴的消耗时间。此配置没有分析或数值基准解决方案，但方形直径随时间衰减的趋势（图20（a））表明，随着网格分辨率的增加，消耗曲线趋于相同结果。在图20（b）中报告了质量守恒测试，其中将域中液体质量的变化（按上一节所述计算）与蒸发的总质量进行比较。此分析针对液相中的两种化学物质和总液体质量进行，三种情况下均获得了良好的质量守恒。
本文提出了一种新的VOF多分量相变数值模型。该模型描述了在非等温环境下具有多种化学物质的系统，并考虑了Stefan流。我们提出了一种求解界面跳跃条件的新方法，以及标量输运方程的几何离散化。由相变引起的体积膨胀通过提出膨胀项移位和双压力速度耦合的组合来解决。采用经典相变问题:固定通量蒸发问题、Stefan问题和Scriven问题对所实现的数值方法进行了验证，得到了二阶收敛速率。该模型是通用的，它允许模拟没有Stefan流的纯扩散系统，就像Epstein-Plesset测试用例一样。更复杂的配置，如多组分液滴的蒸发，可以很好地再现一维数值模型的基准结果。质量守恒误差较小，且随着网格的细化而收敛。
该模型的优点是能够模拟包括热效应在内的混合物蒸发，并且能够正确地管理具有强密度比的静态液滴的界面速度跳变。在这项工作中开发的代码和模拟设置被记录并发布在Basilisk沙盒上[40]，允许其他对相变模拟感兴趣的研究人员从现有框架开始开发和测试新方法。
翻译转载自：《Journal of Computational Physics》：“Multicomponent droplet evaporation in a geometric volume-of-fluid framework”
