JFM: 欧拉/拉格朗日框架中液体燃料射流与雾化蒸发

14天前浏览1782

本文摘要：（由ai生成）

本研究通过高精度数值模拟，探讨了韦伯数对液体燃料射流雾化-蒸发过程的影响。研究发现，韦伯数增加会减小燃料液滴的平均直径，且燃料蒸气质量分数在上游区域较高，下游降低。液柱后温度不因蒸发降低。该研究验证了欧拉/拉格朗日框架在模拟气液界面和分散颗粒蒸发中的准确性，并强调了在预测燃料与氧化剂混合和点火中考虑蒸发过程的重要性。

通过高精度数值模拟研究了空气动力学韦伯数𝑊𝑒_𝑔对液体燃料射流雾化-蒸发过程的影响，并考虑了欧拉/拉格朗日框架中气液界面和分散颗粒的蒸发。通过对单个液滴的蒸发进行二维数值模拟，并将计算结果与理论解和湿度数据进行比较，验证了欧拉方式处理的气液界面蒸发模型的准确性。结果表明，本高精度数值模拟合理地再现了不同破碎模式下错流中液体射流的一次和二次雾化特性。此外，考虑欧拉气液界面和拉格朗日分散颗粒的蒸发对于精确预测燃料与氧化剂的混合和点火至关重要。随着𝑊𝑒_𝑔的增加，整个区域产生的燃料液滴的Sauter平均直径（SMD）减小。此外，随着𝑊𝑒_𝑔的增加，上游地区燃料蒸气质量分数较高而下游地区降低的趋势更加明显。然而，尽管由于蒸发而增加热量损失，但液柱后面的温度不会降低。这些行为可以很好地解释为液体射流和热空气横流对𝑊𝑒_𝑔之间的雾化、蒸发和传热差异。

SIMPOP

绪论

近年来，减少燃气涡轮发动机和工业应用中的二氧化碳和氮氧化物等排放已成为实现温室气体排放量减少到净零的碳中和社会的紧迫问题（联合国环境规划署，2022）。因此，迫切需要为燃气涡轮发动机开发和选择最佳运行条件，以进一步提高其效率并减少排放。为此，需要阐明燃烧器中发生的燃烧机理，并高度准确地预测燃烧行为。然而，燃气轮机中的燃烧器通常是使用液体燃料的喷雾燃烧系统。喷雾燃烧是一种高度复杂的现象，涉及多个基本过程，例如液体燃料的雾化、蒸发、与空气混合和化学反应。因此，许多实验（Faeth, 1996; Linne, 2013; Thomas et al., 2019; Lowe et al., 2019）和数值模拟（Jiang et al., 2010; Jenny et al., 2012; Tachibana et al., 2015;Pillai et al., 2020）阐明了控制液体喷雾和喷雾燃烧的机制。然而，实验存在各种条件下重复测量成本高的局限性。此外，在发动机的高温高压运行环境中只能测量有限数量的物理量（Shinjo，2018）。此外，由于计算性能的快速进步，喷雾燃烧的高精度数值模拟一直是人们关注的焦点。

由于过高的计算成本和建模问题，在之前大多数关于喷雾燃烧的数值研究中都没有直接模拟雾化过程（Wen et al., 2022）。取而代之的是，雾化后的液滴直径分布通常是通过根据流入计算域的燃料液滴的实验数据或理论解调整参数来给出的（Tachibana et al., 2015; Moriai et al., 2013; Kitano et al., 2016;Kahila et al., 2018）。然而，实验液滴直径分布并不总是可用的，理论解决方案通常包含简化和假设，这使得它们在各种进样条件下的适用性值得怀疑。因此，非常需要对整个喷雾燃烧进行数值模拟，特别是液体燃料的蒸发、燃料雾化、雾化液滴的蒸发和汽化燃料的燃烧。

通常，欧拉/欧拉模型用于雾化两相流的数值模拟。在这个框架中，使用诸如流体体积（VOF）方法（Hirt and Nichols, 1981），水平集方法（Yue et al.,2003）以及两者的组合（称为耦合的水平集和流体体积（CLSVOF）方法（Sussman and Puckett, 2000）等技术捕获气液界面。然而，这些方法在气液界面处需要较高的计算网格分辨率，这在应用于雾化过程时导致计算成本非常高，因为雾化过程的许多液滴分布在很宽的区域内。因此，欧拉/拉格朗日框架（Zuzio et al., 2018; Herrmann, 2010）已被提出以降低计算成本。在这个框架中，液体燃料柱和大液滴在计算网格上解析，以欧拉方式捕获界面。相反，生成的分散液滴以拉格朗日方式进行跟踪。此前，Li等（2011）应用欧拉/拉格朗日框架研究了液体燃料射流在错流中的雾化-蒸发过程。Wen等（2020）还使用欧拉/拉格朗日框架研究了空气动力学韦伯数W的变化对错流中液体燃料射流雾化蒸发过程的影响。然而，这些研究只考虑了拉格朗日粒子的蒸发，忽略了连续液相（如液体燃料柱和大液滴）的气液界面的蒸发，这些液相以欧拉方式处理。这是因为对气液界面的蒸发进行建模具有挑战性，这可能导致雾化过程的数值模拟中计算不稳定，涉及气液界面的大变形。近年来，只有少数研究考虑了这些蒸发效应。Li and Soteriou（2015）和 Huang and Zhao（2019）分别研究了环境温度对错流和直流液体燃料射流雾化蒸发过程的影响。然而，蒸发模型的精度尚未得到详细研究，讨论仅限于单个𝑊𝑒𝑔条件。

因此，本研究旨在通过考虑欧拉/拉格朗日模型中气液界面和分散颗粒蒸发的详细数值模拟，研究𝑊𝑒𝑔等对液体燃料射流雾化-蒸发过程的影响。使用CLSVOF方法捕获气液界面。此外，通过对单个液滴蒸发现象进行二维数值模拟，并将计算结果与理论解和湿度数据进行比较，详细验证了欧拉方式处理的气液界面蒸发模型的准确性。

数值模拟

1.欧拉/欧拉框架

●1.1界面捕获方法

对于自由表面流动的模拟，例如错流中的雾化过程，需要以良好的精度和质量守恒来捕获具有大变形的气液界面。有几种类型的界面捕获方法，例如流体体积（VOF）方法（Hirt和Nichols，1981），水平集方法（Yue等，2003），以及两者的组合，称为耦合水平集和流体体积（CLSVOF）方法（Sussman和Puckett，2000）。使用VOF方法，界面可以平流，质量守恒良好，但不能被锐利捕获。相比之下，使用水平集方法，可以更清晰地捕获界面，但可能会违反质量守恒。在这项研究中，采用简单耦合水平集和流体体积（S-CLSVOF）方法（Albadawi等，2013）来利用VOF和水平集方法。使用S-CLSVOF方法，首先平流VOF函数。接下来，从VOF函数获取Level-set函数的初始值。最后，将水平集函数重新拉开距离，成为有符号距离函数，用于计算界面曲率。

VOF函数由以下公式控制。

这里，𝐶是VOF函数，𝜌𝑙是液体密度。𝒖𝒍和𝑚(时间导数)分别是扩展液速和蒸发速率，它们都在第2.1.3节中描述。式（1）的右侧（RHS）显示了由于相变引起的界面偏移的源项，并且仅在使用三维（3D）狄拉克delta函数𝛿𝛤的界面处考虑。对于VOF函数的平流，必须特别注意这样一个事实，即由于相变，单流体速度不会保持界面上的无发散条件，尽管已经为无发散速度场开发了传统的几何平流方法。因此，式（1）分为以下两个步骤。

首先，在式（2）中，VOF函数与扩展的液速𝒖𝒍平流，从而在界面上保持无发散条件。方程（2）可以使用定向分裂方法（Ii等，2012）在时间上离散化，其中数值通量以半解析方式确定，如下所示。

这里，𝛥𝑡是时间步长，𝛥𝑥，𝛥𝑦和𝛥𝑧是每个方向的网格间距。𝑢𝑙、𝑣𝑙和𝑤𝑙分别是扩展液速矢量的𝑥-、𝑦-和𝑧-分量。𝑓𝑛，𝑔∗和ℎ∗∗是数值VOF通量，用分段线性界面计算（PLIC）方案的欧拉平流法计算（Gueyffier et al., 1999; Rider and Kothe, 1998）。通过Youngs方法（Youngs，1982）计算界面正态向量，并使用寻根算法Dekker方法通过界面段的迭代放置（Rider和Kothe，1998）重建界面。使用方程（2）更新VOF函数后，求解方程（3）。如果方程（3）中的源项为负值𝐶，则将其设置为零。这种操作似乎违反了质量守恒，但整个界面捕获方法在第 3.1 节中得到了验证方程（3）中𝛿𝛤的值存储在单元中心，计算公式为

这里，𝛥𝑉是计算单元的体积。𝑆𝑖𝑛𝑡是界面段的面积，表示气液界面与电池的交点，界面段的示意图如图1所示。𝑆𝑖𝑛𝑡的计算方法如下。首先，计算使用PLIC方案重构的单元边交点和界面段的位置，即界面段的顶点;对于3D网格，界面段可以是三角形、四边形、五边形或六边形。其次，为了一致地计算这种多边形的面积，将界面段划分为三角形。最后，𝑆𝑖𝑛𝑡可以计算为三角形的总面积。相反，对于二维网格，𝑆𝑖𝑛𝑡的计算归结为线段长度的计算。

图1 界面段的草图，表示气液界面与计算单元的交点；界面段被划分为三角形，用于计算其面积𝑆𝑖𝑛𝑡

VOF函数平流后，它与Levelset函数耦合。水平集函数𝜙0的初始值计算为

其中𝛤是1.5𝛥的表面厚度，𝛥是网格间距。然后，将初始Level-set函数重新拉开距离，以保持|∇𝜙|=1在界面附近的区域中，通过迭代求解以下方程，这称为重新初始化过程。

这里，𝜏是人工时间步长。方程（7）左侧（LHS）的最后一个项是在高阶约束重新初始化（HCR）方案中引入的强制项（Shiratori et al., 2021; Hartmann et al., 2010）来抵消重新初始化过程中界面的位移，这可以避免不自然的相变。𝜁是一个加权因子，在本研究中设置为0.5。𝐹是单元格𝑖处的强迫项，计算公式为

这里，𝑀𝑖是接口上相邻单元格的数量，其中满足以下关系。

●1.2 控制方程

气相和液相都被视为不可压缩的流体，并以欧拉方式求解。控制方程是以下质量、动量、能量和物质质量分数的守恒方程，这些方程使用称为FK3的内部热流分析代码应用于详细的数值模拟（Wen etal., 2020; Kurose, 2024）。

其中，𝜌是密度，𝒖是速度矢量，𝑃是压力，𝜇是粘度，𝑺是应变率张量𝑆𝑖𝑗=(𝜕𝑖𝑢𝑗+𝜕𝑗𝑢𝑖)/2，𝑭𝑠𝑡是表面张力项，𝑐𝑝是恒压下的比热容，𝑇是温度，𝜆是热导率。𝑌𝑘和𝐷𝑘分别是第𝑘个物质的质量分数和质量扩散系数。𝑇𝑠𝑎𝑡和𝑇𝛤分别是环境压力下的液体饱和温度和局部界面温度。下标𝑙和𝑔分别表示液相和气相的物理值。𝑚(时间导数)在方程中。（12）和（14）是蒸发速率，其计算方法详见第2.1.3节。方程（12）的RHS的第一项和方程（14）的RHS的第二项分别解释了由于速度和焓的相变引起的跳跃。焓跃包括潜热𝐿𝑉和液相和气相之间比热容的差异。𝐿𝑉由以下Watson方程计算（Watson，1943）：

其中𝐿𝑉,𝐵,𝑎𝑡𝑚为常沸点温度𝑇𝐵,𝑎𝑡𝑚下的蒸发潜热，𝑇𝐶为临界温度。这些由相变引起的跳跃被认为是使用式（1）中的三维狄拉克𝛿函数𝛿𝛤的源项，由式（5）估计。

动量守恒方程（13）在非保守平流公式中被离散化，并与VOF平流方程（4）弱耦合，已知该方程可能会引入每个相的动量守恒误差，如Vaudor等人（2017）和Huang等人（2019）所述。动量守恒误差在附录B中进行了讨论，并观察到了能量耗散，这是由平流项的非保守、不一致离散化引入的。该误差对目前对气液密度比相对较低的液体燃料射流的雾化-蒸发过程的模拟并无害处。然而，当应用于具有更高密度比的模拟时，对平流项实施保守、一致的离散化是必不可少的（Arrufat等，2021），并在未来的工作中设想。

使用VOF函数𝐶定义物理属性，如下所示。

这里，𝜉可以是密度𝜌、粘度𝜇、导热系数𝜆和热容𝐶𝑝=𝜌𝑐𝑝。

表面张力项𝑭𝑠𝑡由连续介质表面力（CSF）模型（Brackbill等，1992）考虑，计算公式为

其中𝜎是表面张力系数，𝜅是局部界面曲率。𝜅的值是使用水平集函数计算的，𝜙为

这里，𝒏是气液界面的法向量。由于交错网格中的速度分量是在像元面中心定义的，因此𝜅的值计算为基于VOF函数的体积加权平均值。在x方向上，使用单元格中心的𝐶和𝜅值计算单元格中心之间的面中心的𝜅𝑖+1/2值，如下所示。

在这里，𝜖是一个小常数，以避免除以零，在本研究中设置为1×10−10。该加权假设在VOF函数𝐶为0.5的界面上可以更准确地获得曲率的近似值，并将其性能与附录A中的几种插值进行比较。

𝑆𝜌、𝑆𝜌𝒖、𝑆𝜌ℎ和𝑆𝜌𝑌𝑘分别是质量、动量、能量和化学物质的源项，它们解释了气相和分散燃料液滴之间的相互作用，可以使用粒子源单元（PSI-Cell）方法计算如下（Crowe等，1977）。

其中，𝛥𝑉是计算单元的体积，𝑁是单元中的液滴数，𝑚𝑑、𝒖𝑑和ℎ𝑑分别是液滴的质量、速度和比焓。

控制方程中平流项的空间导数使用WENO方案（Jiang and Shu, 1996）近似，但使用Kawamura-Kuwahara方案评估的动量除外（Kawamura et al., 1986）。控制方程中扩散项的空间导数使用二阶精确中心差分方案进行近似。三阶显式TVD Runge-Kutta方案用于对流项的时间积分。隐式时间离散化用于温度和化学物质的控制方程的时间积分。气相和液相的热力学性质和输运系数分别根据 CHEMKIN（Kee et al., 1986, 1989）计算，并分别从NIST数据库（美国国家标准与技术研究所网络手册，2023年）中获得，两者都考虑了温度依赖性。

●1.3 蒸发模型

气液界面处的蒸发速率𝑚(时间导数)计算如下。

其中𝑌𝐹,𝑔和𝑌𝛤𝐹,𝑔分别是气相和气液界面处的燃料蒸气质量分数。𝑌𝛤𝐹,𝑔通过以下关系在界面处假设气液平衡态获得。

其中𝑋𝛤𝐹,𝑔为界面处的燃料蒸气摩尔分数，𝑊和𝑊𝐹分别为气相中燃料的平均分子量和分子量。𝑃𝑠𝑎𝑡是由佐藤方程（Sato，1954）得到的饱和蒸气压，它是非极性衬底的经验公式，写如下。

其中𝑇𝛤是界面温度，𝑐是常数。𝑃𝑠𝑎𝑡的单位是mmHg。𝑇𝛤的计算方法是将气态域中的液体温度和扩展液体温度𝑇𝑔ℎ𝑜𝑠𝑡平均为

假设单元𝑖处于气相（液级集函数𝜙𝑖<0），单元𝑖+1处于液相（𝜙𝑖+1>0）。气态域中的𝑇𝑖𝑔ℎ𝑜𝑠𝑡是使用 Aslam （2004）中提出的扩展技术获得的。

使用类似于Malan et al. （2021）中提出的方法计算燃料蒸气质量分数∇nΓYF，g）的界面法向梯度，使用PLIC重建确定界面位置，并采用插值技术获得界面值。该过程的示意图如图2所示。相比之下，我们确认本研究中应用的方法可以在静止液滴的界面处提供蒸发速率的均匀分布。蒸发速率的计算过程详细描述如下。

如图2（a）所示，每个单元的∇𝒏𝛤𝑌𝐹,𝑔值首先计算为：

其中𝑑𝛤是从纯单元节点到界面的法线距离，计算公式如下。

对于每个纯单元，存在几个相邻的界面单元，这意味着方程（33）中∇𝑝𝒏𝛤𝑌𝐹,𝑔值的几个候选单元，如图2（b）所示。第一和第二相邻的界面单元可以分别成为第一和第二纯单元的候选者。然后，根据以下因素，确定候选∇𝑝𝒏𝛤𝑌𝐹,𝑔的值，使𝒏𝛤和𝛥𝒙的共线性最大

图2 计算红方块内界面单元处蒸发速率 ̇m 的过程草图。蓝色单元格表示界面单元格。浅红色和深红色的细胞代表仅包含气相的细胞（称为纯细胞），它们分别是界面细胞的第一和第二邻居。蓝线是用PLIC方案重建的气液界面。（a）计算纯电池处的距离𝑑_𝛤。（b）根据共线性𝜉在纯单元上选择𝑑_𝛤。（c）通过对 5 × 5 模板中相邻纯单元中相邻纯单元中(𝑌_𝐹_,_𝑔−𝑌^𝛤_𝐹_,_𝑔)/𝑑_𝛤的值进行加权来计算界面单元处的传质速率。

最后，如图2（c）所示，通过对相邻纯细胞中∇𝑝𝒏𝛤𝑌𝐹,𝑔的值进行加权，计算界面单元处∇𝒏𝛤𝑌𝐹,𝑔的值，如下所示。

其中𝑞是相邻纯单元格的索引。𝜒𝑞是每个纯单元格的归一化加权因子，由下式给出

方程（2）中的VOF函数𝐶应与无发散液体速度𝒖𝑙平流，如第2.1.1节所述。由于相变，单流体速度𝒖不会在界面上保持无发散条件。因此，应将液速扩展到气态域，以获得整个域中的无发散速度场。在整个域中获取𝒖𝑙的过程详细描述如下。由相变引起的Stefan流速𝒖𝑠可以通过求解以下泊松方程来获得。

扩展的液体速度𝒖𝑙现在可以在整个域中计算为

通过这种结构，𝒖𝑙保留了如下无发散条件。

2. 欧拉/拉格朗日框架

●2.1拉格朗日颗粒的控制方程

在欧拉/拉格朗日框架中，通过求解液滴轨迹𝒙_𝑑、速度𝒖_𝑑、温度𝑇_𝑑和质量𝑚_𝑑的方程，以拉格朗日方式单独跟踪生成的分散液滴，如下所示（Wen等，2020）。

使用非平衡Langmuir-Knudsen蒸发模型计算蒸发速率（Miller et al.,1998; Miller and Bellan, 1999; Kitano et al., 2014）如下。

2.2 标记的欧拉液滴向拉格朗日粒子的转变

在本研究中，我们的目标是通过每个步骤的标记方法将欧拉液滴转化为拉格朗日粒子，其详细信息可以在Herrmann（2010）中找到。使用标记方法，需要检测要转化的液体夹杂物。对于本研究，VOF函数𝐶大于0.01的网格被视为用于检测的液体网格。

标记的欧拉液滴将转化为拉格朗日粒子。每个标记的欧拉液滴的体积、重心和温度由以下方程使用VOF函数𝐶计算。

对于拉格朗日粒子，速度是在重力中心给出的，其值是通过平均欧拉速度获得的。满足以下所有条件的欧拉液滴将转化为拉格朗日粒子。变换后的欧拉液滴的VOF函数𝐶的值设置为零。

欧拉/欧拉框架中蒸发模型的验证

通过将计算结果与理论解和湿度数据进行比较，验证了气液界面的欧拉蒸发模型。为理论解比较做出了假设，从而限制了可以验证的信息量。因此，进行了复杂度不断增加的不同数值模拟。将一种数值模拟的结果与理论解进行比较，然后将一种类型的数值模拟结果与湿度数据进行比较。

1.以恒定质量通量阵法液滴

在此基准案例中，考虑了蒸发的2D圆形液滴。液滴的初始直径为𝑑0，并且位于每个方向长度为L为4𝑑0的方形域的中心。在这种情况下，质量通量是以常数形式给出的，而不是由使用方程（28）的燃料蒸气浓度场计算的。因此，在考虑恒定质量通量的同时，可以验证方程（12）中的速度跳跃，方程（40）中构造扩展的液速𝒖𝑙，并更新方程（1）中的VOF函数。

在恒定质量通量条件下，液滴直径d的时间演化理论计算如下

根据式（53），液滴直径在无量纲时间t∗上的演变仅取决于密度比。

因此，在m0（时间导数）∕d0ρg固定在50的条件下，在10、20和50处进行不同𝜌𝑝的情况，并将结果与理论解进行比较。

图3 在256×256网格上计算的ρ_p=10的工况，在t^∗=0.42时（a）速度u和（b）扩展液体速度u_l的瞬时矢量场。气-液界面用C=0.5的白色等值线表示。

图3显示了在256×256网格上计算的ρp为10的情况的瞬时矢量场（a）速度u和（b）扩展液体速度ul。气-液界面由白色实线表示，白色实线是VOF函数C=0.5的等值线。如图3（a）所示，再现了界面上的速度跳跃。向外法向矢量到界面的速度场是在气相中形成的，液相中速度矢量的大小比气相中的速度矢量小得多。此外，如图3（b）所示，ul的向量场在整个域中形成，该向量场尊重界面上的无发散条件，可用于平流方程（2）中的VOF函数。虽然对于静止的蒸发液滴，ul预计为零，但可以在界面周围观察到杂散流动。

图4 （a）在256×256网格上计算的ρ_p=10、20和50的不同密度比；（b）在L∕Δ=64、128和256的不同网格上计算的ρ_p=10的情况下，液滴直径的时间演变与理论解的比较。

图4显示了液滴直径在非量纲化时间t∗上的时间演变，适用于（a）在256×256网格上计算的10、20和50处的不同密度比ρp，以及（b）在L∕Δ=64、128和256的不同网格上计算的10处的ρp。如图4（a）所示，结果与不同密度比的理论解吻合较好。在图4（b）中，即使网格相对粗糙，结果也与理论解一致。通过与理论解的比较，验证了式（1）中VOF函数C的更新方法，即式（2）中C的平流与扩展的液速ul，以及式（3）中考虑相变的源项的估计。

图5 通过与理论面积的比较对计算出的界面面积Sint的评估。（a）中心位于不同初始位置的液滴界面面积的相对误差取决于使用基于PLIC的方法计算的偏差δx值，该值在d₀∕Δ=2、4、8、16、32、48、64和80的不同分辨率的三维网格上计算。

由于上述情况模拟了二维液滴的蒸发，因此还验证了方程（3）中对三维网格中界面面积Sint的估计。直径为d₀的液滴最初设置在不同位置。δx是初始位置的偏差值，从零变为网格间距Δ，增量为Δ∕N。对于此测试，N设置为100。考虑了PLIC方案使用理论界面-法向量和单元中心给出的平面常数重建的气液界面，给出了每个单元的VOF函数的初始值。网格收敛测试是通过将液滴直径d₀∕Δ内的网格点数从2更改为80来执行的。图5（a）显示了中心位于不同初始位置的液滴的Sint的相对误差，这取决于具有不同网格分辨率的3D网格上δx的偏差。相对误差的计算公式为

这里，Stheoretical是液滴的界面面积，其体积与通过求和给定VOF函数计算得出的体积相同，Scomp是每个单元的界面面积的总和，通过第2.1.1节中描述的界面面积计算方法计算，其界面由PLIC方案重构。如图5（a）所示，相对误差随着网格分辨率的增加而减小。此外，在更高的网格分辨率下，无论偏差δx如何，都可以再现可比误差。图5（b）显示了在具有不同网格分辨率的2D和3D网格上计算的不同偏差δx的相对误差的L₂范数。L₂范数的计算公式为

如图5（b）所示，随着网格分辨率的增加，L₂范数停滞不前，网格分辨率指数d₀∕Δ超过8，可以再现小于1%的可比L₂范数。这种停滞可能是由于PLIC方案中界面正态向量的计算方法的低阶。对于本研究，使用第2.1.1节中描述的Youngs方法计算界面正态向量，并且可以通过实现高阶方法（例如有效的最小二乘体积流体界面重建算法（ELVIRA））来实现L₂范数的网格收敛（Pilliod and Puckett, 2004）。此外，2D和3D案例之间的L₂范数差异相对较小。因此，即使液滴的网格分辨率相对较低，基于PLIC的方法也可以很好地估计界面面积，并且可以以与2D情况相同的精度应用于3D情况。

2.浸入空气中的水滴蒸发

该基准案例还考虑了蒸发的2D圆形液滴。计算条件和配置与第3.1节相同。在这种情况下，蒸发速率由方程（28）计算，并求解了蒸气种类和能量传输方程。因此，蒸发是由界面处燃料蒸气浓度的梯度驱动的。因此，可以验证蒸发速率的估计方法，并将蒸气种类和能量传输耦合到纳维-斯托克斯方程。

模拟浸入空气中的水滴，同时将温度和相对湿度值固定在域边界处。最初，液滴和环境温度等于畴边界处的温度，即所谓的干球温度，环境相对湿度与畴边界处的相对湿度相同。随着蒸发的驱动，由于蒸发的潜热，液滴温度逐渐降低。当从环境到液滴的热通量抵消了由于潜热引起的吸热时，液滴温度会收敛到一定值。稳态液滴温度称为湿球温度，并与湿度数据进行比较。在这种情况下，由于方程（1）中的蒸发，不考虑源项而达到稳态。此外，域边界处的温度固定在所需温度。

表1 模拟浸入空气中的水滴蒸发的热物理性质

请注意，模拟中的表面张力系数σ和液体密度ρl的值设置小于水的值。这是因为这些修改可以有效地降低计算成本。较小的σ可以增加时间步长，较小的ρl可以减少达到稳态的时间。使用不同的Tdb和ψ进行模拟，并将稳态下计算的湿球温度Twb与湿度数据Twb,0的湿球温度进行比较。

图6 在128×128网格上计算的ψ=50%和T_db=313K的情况下，温度T和（b）水蒸气质量分数Y_W,_g在t=0.00011s、0.028s和1.1s时的时间变化。Y_W,g的色条的最小值对应于从方程（56）计算出的质量分数。气-液界面用C=0.5的白色等值线表示

图7 在128×128网格上计算的ψ=50%和T_db=313K的情况下，不同时间沿域水平中心线的温度T剖面

图6显示了在128×128网格上计算的温度和（b）水蒸气质量分数的时间变化，ψ和Tdb分别固定在50%和313K。如图6所示，界面周围的温度最初降低，然后液滴温度收敛到湿球温度，这可以从湿度数据中读取。此外，水从气液界面蒸发，水蒸气浓缩场达到稳态。温度场的收敛性也可以在图7中找到，图7显示了在128×128网格上计算的沿域水平中心线的温度曲线，在不同时间实例中，ψ和Tdb分别固定在50%和313K。

图8 在128×128网格上计算不同（a）干球温度T_db和（b）相对湿度的湿球温度T_wb与湿度数据的比较ψ

图8显示了在128×128网格上计算的不同（a）干球温度Tdb和（b）相对湿度ψ的湿球温度Twb与湿度数据的比较。结果表明，在很宽的Tdb和ψ范围内，与湿度数据具有良好的一致性。

错流液体射流的雾化-蒸发过程

1.计算设置

图9 错流雾化-蒸发过程的计算域和条件模拟示意图

图9显示了错流雾化-蒸发过程的计算域和模拟条件。底部安装高度为0.1毫米的边界墙。液体燃料从直径为0.1毫米的喷嘴喷射到空气横流中，液体燃料的抛物线速度分布如温等人（2020）所示。假设燃料为n-C10H22，温度最初为300K。环境压力设置为2.0MPa。具有均匀速度分布的空气横流从y-z平面进入域，温度最初为500K。计算网格是网格间距为5μm的均匀交错笛卡尔网格。计算域分别由x、y、z方向的360×320×240个网格点组成。无滑移和自由条件分别施加在底壁和流出边界的域边界上。

表2 错流雾化-蒸发过程的模拟参数

表3 错流雾化-蒸发过程模拟的计算条件。

表2和表3分别列出了错流雾化-蒸发过程模拟的参数和计算条件。在进行模拟时，空气动力学韦伯数W（例如）的差异以及蒸发过程的考虑。Weg给出如下。

调整喷射液体燃料的体积速度以固定液体与气体的动量比。案例1-3不考虑蒸发过程，而案例4-6考虑蒸发过程。案例1-3和案例4-6分别研究了韦伯数对雾化过程和雾化-蒸发过程的影响。此外，案例1-3和案例4-6之间的比较使我们能够研究蒸发过程本身的影响。其他无量纲参数，如射流韦伯数、错流雷诺数和射流雷诺数，也显示在表2中。通过在日本理化学研究所的超级计算机Fugaku上进行并行计算（使用具有6000个内核的MPI并行化），模拟错流中的雾化-蒸发过程所需的CPU时间约为426,000小时（计算所需的实际时间约为71小时）。

●2.1. 液体射流渗透

图10 （a）案例1、（b）案例2和（c）案例3中VOF函数C=0.5的瞬时等值面快照。每个拉格朗日粒子的颜色表示其大小。

图10显示了案例1-3中VOF函数C=0.5的瞬时等值面快照。每个拉格朗日粒子都按其大小着色。在图10（a）中，较低的韦伯数，发现在蓝色圆圈表示的液柱尖端剥离了相对较大的液滴，这与袋子破碎模式一致。相反，在图10（c）中，较高的韦伯数中，红色圆圈表示的液柱表面形成了许多液体螺纹，并产生了许多微小的液滴，这与剪切破碎模式一致。在图10（b）中，韦伯数为40，可以观察到袋子和剪切破碎模式的特性，这与多模式一致。

图11 （a）案例1、（b）案例2和（c）案例3中中心x-y平面（z=0mm）上的时间平均液体射流轨迹

图11显示了案例1-3中中心x-y平面（z=0mm）上的时间平均液体射流轨迹（红色区域）。黄线是使用Song等人（2015）提供的经验相关性计算的喷流轨迹，如下所示。

正如预期的那样，由于动量比q的固定值，喷流轨迹几乎相同。此外，观察到随着韦伯数的增加，射流穿透力略大，这与预测公式一致。

●2.2 雾化

图12 在案例1-3中，（a）上游区域和（b）下游区域的时间平均液滴体积分布。灰色和红色的部分分别表示在欧拉和拉格朗日框架中捕获的液滴体积

为了详细研究W对雾化过程的影响，图12显示了时间平均液滴体积分布，在案例1-3中，在靠近液柱的上游区域和下游区域获得。灰色和红色的部分分别表示在欧拉和拉格朗日框架中捕获的生成燃料液滴的体积。对于在欧拉框架中捕获的液滴，它们的大小是通过采用与第2.2.2节中描述的相同的标记方法来计算的，以检测具有相同VOF函数阈值的液体夹杂物。研究发现，欧拉框架中未分辨的微小液滴成功地转化为拉格朗日粒子。从这些液滴体积分布中，可以读取三个点。首先，在低韦伯数的情况下，整个区域都存在大液滴，并且在下游观察到更多的大液滴，如实心圆所示。相反，当韦伯数=40和150时，实心圈指示的上游区域存在大液滴，而在下游区域则无法观察到它们。此外，当韦伯数低至10时，下游的Sauter平均直径（SMD）大于上游区域的平均直径（SMD）。相反，当韦伯数较高时，40和150，下游的SMD小于上游区域的SMD。此外，SMD随着整个区域W的增加而降低。

图13 在案例1-3中，（a）上游区域（9<x∕d_inj<10.5）和（b）下游区域（15<x∕d_inj<16.5）垂直于壁的y方向上的时间平均液滴体积分布

为了从雾化过程的角度研究这三点，图13显示了上游和下游区域垂直于壁的y方向上的时间平均液滴体积分布。从这些分布中，可以读取两个点。首先，当韦伯数为10时，高体积比区从上游向下游向直径较大的区域移动，下游可以观察到一些直径大于50μm的大液滴，如实心圆所示。相反，当较高的韦伯数（40和150）时，高体积比区从上游向下游向直径较小的区域移动，上游区域实心圆指示的相对较大的液滴在下游消失。当韦伯数为40时，高体积比区域显示出与韦伯数为150的情况3相同的趋势，但即使在下游区域仍观察到一些大液滴。其次，参考图10，在与液柱尖端相同的高度处，体积比较高，韦伯数为低，而高体积比区在整个区域以高韦伯数为单位向壁扩展。从这些结果中，观察到的雾化过程取决于韦伯数，例如，可以用初级和次级雾化来解释，如下所示。在较低的韦伯数条件下，大液滴从液柱尖端撕裂，即使在下游区域，通过袋式初级雾化，也不会发生二次雾化，这导致SMD增加。相反，在较高的韦伯数条件下，相对较小的液滴通过剪切模式的初级雾化从液管中剥离出来，并由于其较大的局部韦伯数而进一步进行二次雾化，这导致SMD的降低。

表4 比较案例1-6中上游区域（9<x∕d_inj<10.5）和下游区域（15<x∕d_inj<16.5）的Sauter平均直径（SMD）

3. 韦伯数对雾化和蒸发的影响

为了研究韦伯数对雾化-蒸发过程的影响，讨论了案例4-6的结果，这些结果考虑了欧拉/拉格朗日框架中气液界面的蒸发和分散颗粒。

●3.1.蒸发对液体射流渗透和雾化的影响

图14 （a）案例4、（b）案例5和（c）案例6中VOF函数C=0.5和燃料蒸气质量分数=0.01的瞬时等值面快照。每个拉格朗日粒子的颜色表示其大小

图14显示了案例4-6的VOF函数C=0.5和燃料蒸气质量分数=0.01的瞬时等值面快照。每个拉格朗日粒子都按其大小着色。在案例4中，在蓝色圆圈指示的区域中可以观察到袋子破裂行为，韦伯数较低。相反，剪切破坏行为可以在案例6中红色圆圈指示的区域中找到，韦伯数为150，较高。这些主要雾化特性取决于韦伯数，例如案例4-6的雾化特性与案例1-3的雾化特性相似。此外，燃料从液柱和雾化液滴的气液界面蒸发，并被输送到下游区域。

图15 在案例4中，（a）上游区域（9<x∕d_inj<10.5）和（b）下游区域（15<x∕d_inj<16.5）的时间平均液滴体积分布。灰色和红色的部分分别表示在欧拉和拉格朗日框架中捕获的液滴体积

图16 在案例4中，（a）上游区域（9<x∕d_inj<10.5）和（b）下游区域（15<x∕d_inj<16.5）垂直于壁的y方向上的时间平均液滴体积分布

无论有没有蒸发，案例4-6中的时间平均液体射流轨迹与案例1-3中的相同。此外，当查看图15和图16中案例4的上游区域和下游区域的时间平均液滴体积分布时，可以发现与案例1相同的初级雾化趋势。在其他韦伯数条件下，例如在案例5和6中，结果分别与案例2和案例3的结果相同。此外，在比较案例1-6中4中上游和下游区域的SMD时，发现案例4中的SMD与案例1中的SMD不同。这是因为在韦伯数为10的低韦伯数的算例4中，在液柱尖端剥离的液滴相对较大，并且计算出的SMD受该区域存在的此类液滴的影响很大。对于韦伯数较高的案例5和6，例如，SMD分别与案例2和案例3的SMD几乎相同。最初，由于蒸发引起的液滴体积减少，那里的SMD可能会更小。然而，在本例中，计算域不够大，无法观察到由于蒸发而导致的SMD的这种降低。对于案例5和6，液滴体积分布和SMD下游分别与案例2和3相同。因此，从上面的这些结果可以看出，在当前条件下，蒸发过程不会直接影响雾化行为。

●3.2. W对蒸发的影响

图17 （a）算例4、（b）算例5和（c）算例6中燃料蒸气质量分数在中心x-y平面（z=0mm）上的时间平均分布。气-液界面用C=0.5的白色等值线表示

图18. 在算例4-6中，燃料蒸气质量分数在（a）x∕d_inj=10.5和（b）x∕d_inj=16.5的y-z平面上的时间平均分布

图19. 在算例4-6中，速度矢量和燃料蒸气质量分数在（a）y∕d_inj=2和（b）y∕d_inj=6的z-x平面上的时间平均分布。气-液界面用C=0.5的白色等值线表示

这些雾化过程取决于韦伯数，例如强烈影响蒸发过程。图17显示了案例4-6中z=0mm时中央x-y平面上燃料蒸气质量分数的时间平均分布。此外，图18显示了在案例4-6中x∕dinj=10.5和16.5时y-z平面上燃料蒸气质量分数的时间平均分布。此外，图19显示了案例4-6中速度矢量和燃料蒸气质量分数的时间平均分布，在y∕dinj=2与液柱尖端相交的y∕dinj=6时。白色等值线表示气液界面。结果表明，在所有韦伯数条件下，上游区域的燃料蒸气质量分数较高，下游区域的燃料蒸气质量分数较低。从图19中可以看出，在上游区域，来自气液界面（如液柱表面）和从液柱表面突出的液螺纹的燃料蒸气比拉格朗日粒子的燃料蒸气占主导地位，并且在所有韦伯数条件下都停滞在那里。为了进一步研究气-液界面的蒸发，图20显示了案例4-6中气-液界面上蒸发速率的瞬时分布。研究发现，沿从液柱表面突出的液体螺纹的蒸发速率最高。这种蒸发速率分布可以用图21来解释，图21显示了案例4-6中气液界面上界面温度的瞬时分布。由于来自环境的热传递，沿液体螺纹的界面温度最高，界面面积扩大提高了温度，这对应于图20中蒸发速率最高的区域。相反，在下游区域，当参考图19中的速度分布时，可以观察到空气横流夹带在液柱后面的反向旋转涡流对（CVP）中。在CVP的影响下，燃油蒸气与空气错流混合，导致下游燃油蒸气质量分数降低。

图20 VOF函数C=0.5的瞬时等值面快照，由（a）算例4、（b）算例5和（c）算例6中气液界面上的瞬时蒸发速率着色

图21 （a）算例4、（b）算例5和（c）算例6中，VOF函数C=0.5的瞬时等值面快照由气液界面上的瞬时温度

图22 在算例4-6中，燃料蒸气质量分数沿z-x平面上平行于x轴的中心线的剖面（a）y∕d_inj=2和（b）y∕d_inj=6

此外，随着韦伯数的增加，上游地区燃油蒸气浓度高而下游地区变低的趋势更加明显。在高韦伯数值下，例如在算例5和6中，上游区域的蒸发会增强。相反，尽管燃料蒸气浓度在下游变得低，韦伯数值高，但其倾向于在低韦伯数值下保持下游，特别是在液柱尖端的高度，如图19（b）所示。如第4.3.1节所述，在低韦伯数的情况下，即使在下游区域，大液滴也会从液柱尖端撕裂，由于界面面积扩大，这增强了蒸发过程。因此，在低韦伯数条件下，与空气错流混合的影响被袋式初级雾化的影响部分抵消。这种取决于韦伯数的趋势在图22中得到了定量验证，图22显示了在算例4-6中y∕dinj=2和y∕dinj=6的z-x平面上沿平行于x轴的中心线的燃料蒸气质量分数的剖面。

图23 （a）算例4、（b）算例5和（c）算例6中温度T在中心x-y平面（z=0mm）上的时间平均分布。气-液界面用C=0.5的白色等值线表示

图24 在算例4-6中，温度T沿平行于y∕d_inj=2的z-x平面上x轴的中心线的剖面

图25 在算例4-6中，温度T在（a）y∕d_inj=0.25和（b）y∕d_inj=2的z-x平面上的时间平均分布。气-液界面用C=0.5的白色等值线表示

图23显示了案例4-6中温度T在中心x-y平面（z=0mm）上的时间平均分布。白色等值线表示气液界面。在所有韦伯数条件下，由于低温液体燃料射流和高温空气错流之间的热传递，液柱后面的温度低于空气错流的温度。此外，还发现液柱后面的温度按情况4、情况6和情况5的顺序增加。这种温度随韦伯数变化的趋势也可以在图24中找到。这种温度对韦伯数的依赖性可以通过液柱后面的涡流结构和传热来解释。图25显示了在壁附近和在与液柱相交的z-x平面上温度的时间平均分布。在韦伯数为10的低韦伯数算例1中，逆时针循环流在液柱后面形成，如图23所示，能量主要从壁附近的高温区A输送，如图23和图25（a）所示。A区之所以出现这种高温，是因为液柱后面的CVP规模小于其他两种情况，并且与高温空气错流的能量混合更加增强。相反，在韦伯数为150的算例3中，CVP在液柱后面的逆时针循环流中占主导地位，如图23所示，能量主要从高温尾流区D传输，如图23和25（b）所示。D区的高温是因为高韦伯数，使液柱迎风侧的温度边界层变薄，并且被环境气体急剧加热的燃料蒸气被夹带入液柱后面形成的CVP，这可以在25（b）中观察到。在韦伯数为40的算例2中，观察到逆时针循环流和CVP，其中能量分别从B区和C区传输，如图23所示。与算例1和算例3相比，B和C区域的温度较低，这是由于液柱后面的CVP比图25（a）所示的算例1更大，并且液柱迎风侧的温度边界层比算例3更薄，如图25（b）所示。这导致液柱后面的最低温度，如图24所示。

在许多工作中忽略了气液界面的蒸发，例如（Wen等，2020）。然而，考虑到上述结果，可以说，在目前从欧拉向拉格朗日公式过渡的条件下，欧拉气液界面蒸发对总蒸发的贡献大于拉格朗日分散粒子的贡献。在欧拉/拉格朗日框架中考虑这两种蒸发对于精确预测燃料与氧化剂的混合和实际燃烧器中的点火至关重要。

结论

本研究通过考虑欧拉/拉格朗日框架中气液界面和分散粒子蒸发的详细数值模拟，研究了空气动力学韦伯数对液体燃料射流雾化-蒸发过程的影响。用CLSVOF方法捕获气液界面。特别是，通过对单个液滴蒸发进行二维数值模拟，并将计算结果与理论解和湿度数据进行比较，验证了欧拉方式处理的气液界面蒸发模型的准确性。获得的主要结果总结如下。

1.本详细的数值模拟合理地再现了不同破碎模式下错流中液体射流的一次和二次雾化特性，即W的影响。在韦伯数=10的低韦伯数条件下，在液柱尖端剥离出相对较大的液滴，这与袋式模式一致。另一方面，在韦伯数=150的高韦伯数条件下，液柱表面形成许多液线，产生许多微小的液滴，这与剪切模式一致。这些微小的液滴由于其较大的局部韦伯数而容易进行二次雾化，例如。此外，在目前从欧拉公式到拉格朗日公式过渡的条件下，欧拉气液界面蒸发对总蒸发的贡献大于拉格朗日分散粒子的贡献，这表明考虑这两种蒸发对于精确预测燃料与氧化剂的混合和实际燃烧器中的点火至关重要。

2.无论有没有蒸发，随着韦伯数的增加，生成的燃料液滴的Sauter平均直径（SMD）在整个区域内减小。这是由于雾化行为的差异。也就是说，在较低的韦伯数条件下，即使在下游区域，大液滴也会通过袋式初级雾化从液柱尖端撕下，并且不会发生二次雾化，这导致SMD增加。另一方面，在较高韦伯数的条件下，相对较小的液滴主要通过剪切模式的初级雾化从液螺纹中剥离出来，并进一步进行二次雾化，导致SMD的降低。

3.在所有韦伯数条件下，上游区域的燃料蒸气质量分数较高，下游区域变低。这是因为来自液柱表面的燃料蒸气和从液柱表面突出的液体螺纹比来自拉格朗日粒子的燃料蒸气占主导地位，并停滞在上游区域，而燃料蒸气与下游液柱后面形成的循环流中夹带的空气错流混合。这种趋势随着韦伯数的增加而更加明显，例如。这是因为，在较高的韦伯数条件下，由于从液柱表面突出的液体螺纹的界面面积扩大，蒸发过程得到增强。

4.随着韦伯数的增加，尽管液柱后面的燃料蒸气质量分数增加，如上文所述，即使蒸发引起的热损失增加，该区域的温度也不会降低。在较高的韦伯数条件下，液柱迎风侧的温度边界层变薄，被环境气体急剧加热的燃料蒸气被夹带入液柱后面形成的循环流中，导致那里的温度升高。

附录A.静止液滴

在这里，我们通过研究动量守恒方程以非保守形式离散化可能引入的非物理速度来验证本方法的性能。在此基准案例中，考虑了二维圆形静止液滴。基准工况的初始条件和热物理性质与错流中液体射流雾化-蒸发过程的模拟相同。液滴直径为0.1mm，与错流模拟中的喷嘴直径相同，以研究气液界面周围的这种非物理速度对雾化-蒸发过程的影响。

图26 （a）通过方程（21）的体积加权平均值将基于水平集的界面曲率插入到像元中心时，静止液滴在3ms时的速度大小分布；（b）基于水平集的界面曲率通过算术平均值插值到像元中心；（c）通过HF方法计算的界面曲率（Popinet，2009）通过体积加权平均值插值到网格中心

图26显示了静止液滴在3ms处的速度幅度分布，其中有三种界面曲率插值到细胞面中。当方程（19）基于水平集函数计算界面曲率时，液滴周围的速度大小最小，并使用方程通过体积加权平均值插值到像元中心。（21）（图26（a））。当通过方程（19）计算界面曲率并通过算术平均值插值到单元中心时，可以观察到液滴周围速度大小的不均匀分布，并使模拟不稳定[图26(b)]。相反，在[图26(c)]中，当通过HF方法计算界面曲率时，液滴周围的速度大小相当小。

图27 比较静止液滴的最大速度大小在各种界面曲率插值到单元中心（a）无和（b）有对数尺度的间值

图27显示了在三种界面曲率插值到单元中心的界面中，静止液滴最大速度大小的时间演变的比较。如图26（c）所示，在HF方法的情况下，速度幅度急剧降低。此外，体积加权平均值和算术平均值的速度大小分别逐渐减小和增加。就速度幅度而言，除HF方法外，体积加权平均值优于将界面曲率插值到像元中心的算术平均值。尽管如第3.1节所述，HF方法的速度幅度是三种情况中最小的，但基于水平集的方法比HF方法更适合估计横流模拟中大变形的界面的曲率。此外，与错流模拟中的速度场相比，图26（a）中观察到的速度大小可以忽略不计。因此，在本研究中，界面曲率由使用方程（19）的水平集函数估计，并使用方程（21）通过体积加权平均值插值到像元中心。

附录B.振荡液滴

图28 在L∕Δ=64、128和256的不同网格上计算的情况，动能的时间演变与理论总能量的比较

在这里，我们通过研究由于表面张力效应而振荡的二维粘性液滴的动能，确认了本方法的性能，其中动量守恒方程以非保守形式离散化。所呈现的案例与Vaudor et al.（2017）中的案例相同。半长轴和半短轴分别为0.15和0.1的椭圆液滴位于边L为1的方形域的中心。流体的物理性质与Vaudor等人（2017）中的物理性质完全相同。图28显示了使用不同网格分辨率计算的情况动能的时间演变与理论总能量的比较。能量值由初始能量E0无量纲化，计算为 E0=σ（Sellipse−Scircle），其中σ是表面张力系数，Sellipse和Scircle分别是具有相同面积的初始椭圆液滴和圆形液滴的总弧长。通过参考关于粘性液滴阻尼效应的理论工作（Rush和Nadim，2000），计算了总能量的时间演变如下。