基于Level Set方法的液体燃料雾化DNS模拟
作为能源利用的主要形式之一，液体燃料的液雾燃烧广泛存在于能源动力、航空宇航推进、化工等重要工程领域。相比于固体燃料，液体燃料比冲高、可多次点火，在航天领域有着重要的应用价值。
燃料的雾化和后续蒸发以及燃烧过程密切联系，涉及装置燃烧的稳定性、燃料的燃烧效率以及污染物的排放。因此，雾化过程的深入研究具有重要意义。然而，雾化过程复杂，包含初次雾化和二次雾化过程，涉及到湍流、气液界面运动、液体的结构破碎和聚并等方面[1-2]。 图1 由液相特征长度与局部网格大小的比值及液相体积分数划分的雾化的不同区域[3]
在雾化的机理研究中，捕捉到精细的雾化结构，尤其是一次雾化及液体结构的复杂运动具有重要意义。目前，许多学者采用了界面捕捉方法，可以很好的捕捉到雾化发展过程。方法主要集中于Volume of Fluid（VOF）方法和Level Set（LS）方法，以及它们的一些衍生形式。 在数值计算中，直接数值模拟方法（DNS）区别于雷诺平均法（RANS）和大涡模拟法（LES），具有精确、计算数据信息详细等特点。通过结合研究气液两相流方法之一的LS方法，近些年学术界对液体雾化机理的研究取得了一定进展。
Level Set方法属于欧拉框架下的界面捕捉方法，最早由Osher和Sethian[5-6]提出，对处理复杂的拓扑结构具有优势。气液相界面运用较高维的LS函数的零等值轮廓表示，如图2所示。
图2 二维空间水平集表示的自由边界（蓝色实线），沿其法线方向移动，拓扑结构会自动处理 (红色为水平集方程) [7]
另外输运过程中需计算两个重要的量，即界面的法向和平均曲率： Level Set方法容易追踪界面，容易处理界面的断裂与合并，可以直接拓展到三维，并且曲率和法向矢量等易于计算。另外，距离函数本身是连续函数，对于连续函数的求解已经较为成熟，通过采用连续的距离函数对自由界面法向的求解精度高。
Kun Luo等人[8]、van der Pijl等人[9]结合了VOF方法的优点-质量守恒，对LS方法进行了改进，并通过benchmark验证了质量接近守恒。而针对符号距离函数的重新初始化问题，改进措施可归纳为4类[10，11]：
1)LS输运方程的高阶离散；
2)延展速度型LS方法；
3)双曲正切型LS方法；
4)重新初始化的改进方法。
Chang 等人[12]提出了一种求解摄动汉密尔顿-雅克比方程使其达到稳态。后期，McCaslin和Desjardins[13]进一步改进了初始化方程，考虑了 LS函数输运过程中较大的空间变化。同时，有些学者为了减小重新计算过程的计算量，提出了一种快速推进方法（Fast Marching Method）[6]。
Min Chai等人[14]提出了基于双S函数来改进初始化，通过对比发现计算量也得到了减小。
通过以上改进后的LS方法，可以使得燃料雾化模拟中，其质量守恒性得到极大的保障。
对多尺度（包括时间和空间）、多物理场、多相的雾化过程，早期对初次雾化的认识进展缓慢。近年来，计算机技术带动了计算流体力学的应用，关于初次雾化的直接数值模拟研究[1]得到不断发展，如图3所示，涉及旋流、横向射流、互击式、交叉射流等不同的雾化形式。
通过分析模拟计算结果，可以够清晰观察到液柱-液条-液丝-液滴的整个破碎演变过程；结合雾化的Kelvin-Helmholtz(KH)不稳定性和Rayleigh-Taylor(RT)不稳定性机理[15-16]，对不同雾化形式下液体相的发展产生进一步认识；同时可以得到雾化锥角、破碎长度等宏观雾化参数以及粒径分布等微观雾化参数的值。 在实际的航空发动机中，旋流雾化形式非常重要，下面以一例进行讨论。Shao等人[17]在质量守恒的LS方法[7]基础上，对比研究了层流入口和湍流入口条件下旋流雾化的结构和机理。其网格在十亿规模（1024×1024×1024），最小的解析尺寸为3.9µm。
如图6所示，初始时需要给定一个两相界面条件，此旋雾化属于气动雾化。中心圆孔喷动气体，外部环流喷动液体，内环直径为0.2mm，外环直径为0.4mm，整体旋流数S=1。
为了产生实际的湍流入口，采用液相的物性参数模拟了周期性的管道湍流流动，使其充分发展再作为湍流入口的初始条件。
图7 管道流速度场为湍流入口条件示意图（充分发展）
对旋流雾化的模拟计算罗列出几个时间的发展结果，可以看出从初始到一定时间过程，整个旋流雾化的演变(图8中从左往右依次是T = 0, 1.25, 2.50, 3.75, 5.00, and 6.25）。对比发现湍流会减弱液膜的刚度，并导致液体结构沿径向均匀分布。
在旋流射流中，主要的不稳定性包括KH不稳定性和离心不稳定性。KH不稳定性是由气液界面的剪切造成的，而离心不稳定性是由周向波动造成的。对于气液两相流而言，还应考虑由于密度差异引起的RT不稳定性。在模拟中，沿轴向的KH不稳定波并不明显，而沿周向的离心不稳定波和RT不稳定波在液膜端部出现，导致液桥的形成。研究还发现随着液膜向下游发展，液膜的厚度逐渐变薄，不稳定波逐渐增强。当不稳定波达到液膜厚度量级时，液膜发生破碎。与单相射流相比，两相射流在下游具有混沌速度，可以增强液滴和周围气体两相间的混合，并可以产生旋进涡流核（PVC）。
基于LS方法的雾化直接数值模拟对于认识初次雾化及雾化机理具有重要贡献，但也存在有问题需要去解决。 经验上，为了较好的在欧拉框架中解析一个近球形液滴的界面，需要在液滴直径上至少分布4-8个网格点。 在早期研究雾化及液雾燃烧时，常假定已经雾化完全，把液滴看作拉格朗日点源来进行追踪，与气相流动进行双向耦合。事实上，在雾化发展的后期成为了微小液滴，它的表面张力占主导地位，使它保持近球形。因此，如果把它看作球体，问题就会简化。这涉及到拉格朗日点源方法，小而近似球形的液滴被看作粒子。与欧拉描述的明显区别在于液滴是否有自己的体积，在拉格朗日框架中，没有考虑液滴的体积。只使用拉格朗日点粒子来表示它，传递它的主要信息，如质量或动量信息，不再需要大量的网格单元。 在整个雾化中用欧拉界面捕获方法来处理主要的液体结构，用拉格朗日点源描述来处理分辨率较低的液滴结构，计算成本也会得到降低。
采用这种不同框架的混合方法，结合了界面捕捉方法与点源方法各自的优势，相互弥补了不足，且能够适用整个雾化过程。Herrmann[20]提出了一套详细的欧拉/拉格朗日混合计算框架，并应用于直射流中进行了模拟，如图9所示。
图9 基于欧拉和拉格朗日混合框架对直射流的模拟结果
欧拉界面捕捉方法：大尺度雾化结构
拉格朗日点源：小尺度液滴
随着计算工况雷诺数、韦伯数等的增大，液滴破碎会更加剧烈，这种混合框架方法体现出优势。除了帮助我们深入认识雾化机理之外，还有助于发展LES下甚至RANS框架下的雾化模型，对解决工程实际中的雾化问题具有重要意义。
[1] Kun Luo, Changxiao Shao, Min Chai, Jianren Fan. Level set method for atomization and evaporation simulations[J]. Progress in Energy & Combustion Science, 2019, 73:65-94.
[2] Jenny P, Roekaerts D, Beishuizen N. Modeling of turbulent dilute spray combustion[J]. Progress in Energy & Combustion Science, 2012, 38(6). [3] Frank Ham, Sourabh Apte, et al. Unstructured LES of reacting multiphase flows in realistic gas turbine combustors[J]. Center for Turbulence Research, 2003. [4] Gueyffier D, Li J, Nadim A, et al. Volume-of-Fluid Interface Tracking with Smoothed Surface Stress Methods for Three-Dimensional Flows[J]. Journal of Computational Physics, 1999, 152(2):423-456. [5] Osher, S. and Sethian, J. A., Fronts propagating with curvature-dependent speed:Algorithms based on Hamilton-Jacobi formulations Journal of Computer Physics, 1988,79,12-49. [6] Osher S, Fedkiw R. Level Set Methods and Dynamic Implicit Surfaces[M]// Level sets methods and dynamic implicit surfaces /. Springer, 2003:B15 [7] F Gibou, R Fedkiw, S Osher. A review of level-set methods and some recent applications[J]. Journal of Computational Physics, 2018, 353:82-109. [8] Luo K , Shao C , Yang Y , et al. A mass conserving level set method for detailed numerical simulation of liquid atomization[J]. Journal of Computational Physics, 2015, 298:495-519. [9] van der Pijl, SP, Segal, A, Vuik, C, et al. A mass-conserving Level-Set method for modelling of multi-phase flows[J]. International Journal for Numerical Methods in Fluids, 2010, 47(4):339–361. [10] Nourgaliev R R, Wiri S, Dinh N T, et al. On improving mass conservation of level set by reducing spatial discretization errors[J]. International Journal of Multiphase Flow, 2005, 31(12):1329-1336. [11] Sussman M, Almgren A S, Bell J B, et al. An Adaptive Level Set Approach for Incompressible Two-Phase Flows [J]. Journal of Computational Physics, 1997, 148(1):81-124. [12] Changa Y C, Houa T Y, Merrimanb B, et al. A Level Set Formulation of Eulerian Interface Capturing Methods for Incompressible Fluid Flows[J]. Journal of Computational Physics, 1996, 124(2):449-464. [13] Mccaslin J O, Desjardins O. A localized re-initialization equation for the conservative level set method[M]. Academic Press Professional,Inc. 2014. [14] Chai M , Luo K , Shao C , et al. An Efficient Level Set Remedy Approach for Simulations of Two-Phase Flow Based on Sigmoid Function[J]. Chemical Engineering Science, 2017:S0009250917303986. [15] Gallaire F, Chomaz J M . Instability mechanisms in swirling flows[J]. Physics of Fluids, 2003, 15(9):2622-2639. [16] Gallaire F , Chomaz J M . Mode selection in swirling jet experiments: a linear stability analysis[J]. Journal of Fluid Mechanics, 2003, 494:223-253. [17] Shao C , Luo K , Yang Y , et al. Detailed numerical simulation of swirling primary atomization using a mass conservative level set method[J]. International Journal of Multiphase Flow, 2017, 89:57-68.D [18] Shinjo J , Umemura A . Simulation of liquid jet primary breakup: Dynamics of ligament and droplet formation[J]. International Journal of Multiphase Flow, 2010, 36(7):513-532. [19] Wu K J , Reitz R D , Bracco F V . Measurements of drop size at the spray edge near the nozzle in atomizing liquid jets[J]. Physics of Fluids, 1986, 29(4):941. [20] Herrmann M. A parallel Eulerian interface tracking/Lagrangian point particle multi-scale coupling procedure[J]. Journal of Computational Physics, 2010, 229(3):745-759. 
