1 引言

离心式喷嘴由于雾化性能好、点火性能好、稳定燃烧范围宽、结构简单、运行可靠等优点^[1]，被广泛应用于航空发动机（燃气轮机）、火箭发动机、内燃机及工业锅炉等热力设备上。特别是敞口型离心式喷嘴广泛应用于液氧/煤油液体火箭发动机中，如俄罗斯的液体火箭发动机RD120，RD170，RD180^[2]和我国新一代液氧煤油大型火箭发动机YF100，YF115。

目前离心式喷嘴的数值模拟研究主要集中在内流场，外喷雾场锥形液膜雾化破碎过程主要依靠实验研究。York对液膜破碎过程进行了理论分析，认为气液界面处波的形成与增长是导致液膜破碎的重要原因^[3]。岳明等对离心式喷嘴喷出的锥形液膜做了空间不稳定性分析，导出了无粘锥形液膜的色散方程，并进行了数值求解和分析研究，分析了扰动波增长率的变化规律^[4]；Youngbin等实验拍摄了不同喷注压降、不同反压下的液膜破碎图像^[5]；Paul等实验研究了高压离心式喷嘴液膜破碎机理，使用高速转鼓相机和铜蒸汽激光获得了雾化区域的细节图像，并使用“Void fraction”图像分析技术，清楚地看到连续雾锥上形成孔洞的细节^[6]；刘娟等运用高速阴影拍摄了锥形液膜破碎过程图像，得到液膜破碎长度、液膜锥角随喷注压降、喷孔直径以及几何特性参数的变化规律^[7]；Yue等应用PIV激光粒子图像测速仪对由离心式压力雾化喷嘴喷出的锥形液膜的破碎进行了实验研究^[8]；张新桥等使用高速摄影仪对单排/双排切向孔布局的离心式喷嘴雾化特性进行了实验研究，研究表明双排切向孔布局使得液膜雾化性能变差，雾锥角减小，液膜破碎长度增大^[9]。

由于涉及多相多尺度流动问题，外喷雾场的锥形液膜雾化过程数值模拟一直是难点。前人在这方面做了一些探索研究，但从搜集的资料来看，相关工作较少。徐让书等采用VOF和LES对离心式喷嘴内部流动和液膜破碎过程进行了耦合的CFD模拟，捕获了清晰的气液相界面以及油膜的不稳定波和初级破碎现象^[10]；刘娟采用VOF和LES方法，获得了离心式喷嘴锥形液膜一次破碎过程的结果^[11]；Olivier等基于DNS/LES和改进的守恒Level set方法，清楚显示了离心式喷嘴出口形成的锥形液膜和喷嘴里面的空气核^[12]。这些数值仿真都仅限于一次雾化过程，然而对于二次雾化过程仿真难度很大，目前缺少这方面的相关研究。借助仿真手段获得包含一次雾化和二次雾化在内的液膜破碎过程，可以获得液滴的空间分布和粒径分布，将为后续的蒸发、混合和燃烧提供较为有效的初边值条件。然而实现二次雾化过程中液丝破碎成尺寸大范围变化的液滴时精确捕捉，难点在于网格尺度的大范围变化和相应相界面的精确捕捉。基于自适应网格相界面捕捉技术的Gerris开源程序，通过建立合理的计算模型便能有效解决这一难题。

目前的雾化破碎过程研究主要靠实验，但实验测量有一定的局限性。本文针对特定结构的敞口型离心式喷嘴，对其锥形液膜破碎过程进行高分辨率的细致模拟，可视化揭示锥形液膜的形成、液膜破碎成液丝、液丝破碎成液滴的全过程，并获得更准确的雾化特性参数，与高速摄影拍摄的液膜破碎过程图像、PDPA测量的液滴SMD进行对比分析，验证数值方法的准确性，突出显示本文采用的雾化过程数值模拟方法的强大优势，为精细研究离心式喷嘴雾化机理及准确评估雾化特性提供一种新的手段。

2 数学物理模型 2.1 数值方法

Gerris求解的是不可压、变密度、带有表面张力的N-S方程，见文献[13]。

$ \left\{ \begin{array}{l} \rho \left( {\frac{{\partial \mathit{\boldsymbol{u}}}}{{\partial t}} + \mathit{\boldsymbol{u}} \cdot \nabla \mathit{\boldsymbol{u}}} \right) =-\nabla p + \nabla \cdot \left( {2\mu \mathit{\boldsymbol{D}}} \right) + \sigma \kappa {\delta _{\rm{s}}}\mathit{\boldsymbol{u}}\\ \frac{{\partial \rho }}{{\partial t}} + \nabla \cdot \left( {\rho \mathit{\boldsymbol{u}}} \right) = 0\\ \nabla \cdot \mathit{\boldsymbol{u }}= 0 \end{array} \right. $

(1)

式中u=(u, v, w)为流体速度；ρ≡ ρ(x, t)为流体密度；μ≡ μ(x, t)为动力粘度；D为应变张量，$ {D_{ij}} \equiv \frac{1}{2}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right) $；δ_s函数表示表面张力项集中在界面上；κ和n分别是界面的曲率和法向方向。

对于两相流动，引入第一相的体积分数c(x, t)来定义密度和粘度

$ \left\{ \begin{array}{l} \rho \left( c \right) \equiv c{\rho _1} + \left( {1-c} \right){\rho _2}\\ \mu \left( c \right) \equiv c{\mu _1} + \left( {1-c} \right){\mu _2} \end{array} \right. $

(2)

式中ρ₁, ρ₂和μ₁, μ₂分别是第一相和第二相的密度和粘度。

将密度的连续方程用体积分数的连续方程替换

$ \frac{{\partial c}}{{\partial t}} + \nabla \cdot \left( {c\mathit{\boldsymbol{u}}} \right) = 0 $

(3)

在Gerris中体积分数/密度和压力采用时间交错离散格式，使用经典的时间分裂投影法进行简化，达到时间二阶精度。使用四叉树/八叉树进行空间离散，使得自适应加密算法可简易灵活地实现，在不损失计算精度的情况下显著降低了计算量，非常适合处理多尺度流动问题。使用Godunov动量差分格式，基于四叉树/八叉树的多尺度求解器可以有效求解压力泊松方程，达到空间二阶精度。使用分段线性的VOF几何重构方法进行自适应网格界面捕捉，非常适合应用于包含破碎、聚合现象的雾化过程计算。通过将表面张力转化为某一区域连续的体积力并结合高度函数曲率估计实现表面张力的精确求解。采用单调集成大涡模拟MILES（又称隐式大涡模拟ILES）近似模拟亚格子SGS（SubGrid-Scale）的能量传递，这是由于数值计算不可避免地有数值耗散，这种数值耗散的效果与亚格子应力的作用类似，Boris在1992年系统提出并分析了MILES方法，见文献[14]。所以Gerris使用的网格自适应过程相当于是一个滤波的过程，对于大于最小网格尺度的涡结构进行直接数值模拟DNS，而小于最小网格尺度的涡结构由数值粘性模拟其耗散过程。

特别说明Gerris对于小液滴的处理方法：Gerris采用自适应网格技术和VOF相界面捕捉技术对两相界面界面进行捕捉，能够捕捉到小液滴，如图 1所示^[15~17]。这主要是由于使用自适应网格加密算法可以根据流场参数变化对局部网格进行动态加密或粗化，对两相界面处进行细化加密时，网格可以自动加密到液滴的尺度量级，再结合分段线性的VOF方法就可以比较精确地捕捉到相界面，自然可以捕捉到小液滴。

2.2 计算模型

计算所选用的基准喷嘴Chk-0结构参数如图 2所示，敞口型喷嘴总长度L=15mm，喷嘴直径D=5mm，切向孔直径d_t=0.7mm，切向孔个数n=4，切向孔轴线距离喷嘴底部l=2mm，流量为40g/s。

边界条件设置：进口流量为40g/s，使用速度入口边界，折算的切向孔入口速度为26.03m/s，出口为压力出口边界，壁面均为无滑移壁面边界；空气为第一相，水为第二相。使用Coupled Level Set+VOF方法计算喷嘴内部两相流动；压力离散方法为PRESTO；压力速度耦合方法采用SIMPLE；动量方程等均采用二阶迎风格式；表面捕捉采用CICSAM技术^[18]；湍流模型选用带旋流修正的k-ε Realizable模型，它在强流线弯曲、漩涡和旋转流动中有更好的表现；使用非均衡的壁面函数的方法来改进壁面部分的准确度；计算区域中的流体是不可压空气和水。

采用Fluent计算喷嘴内部流动，得到喷嘴出口截面的液膜厚度、轴向速度、切向速度和径向速度的分布，如图 3所示，其中u和w分别为轴向速度和切向速度，都用切向孔进液的速度U_in无量纲化，r为径向位置，用喷嘴半径R无量纲化，可以看出喷嘴出口截面轴向速度和切向速度具有明显的分区分相流动特征^[19]。依据这一分区分相流动的特征，用Matlab编写对应的接口程序，分区分相进行出口截面的参数传递。最终实现将Fluent采用贴体网格计算内流场的优势与Gerris采用自适应网格界面捕捉技术计算外喷雾场的优势相结合，更准确地将喷嘴出口的参数传递给Gerris作为外喷雾场的入口条件。

这里之所以采用Fluent和Gerris两种软件的结合分别计算内流场和外喷雾场，是为了解决单独使用Gerris计算喷嘴全流场时，用笛卡尔网格逼近轴对称固壁，存在网格非体贴效应。这种非贴体网格效应会使得从喷嘴流出的环形液膜周向参数分布不均，表面出现褶皱，离开喷嘴不远，在气动力和表面张力的作用下形成异常的沿流向不连续的纵向液束。经过分析这是由于喷嘴内部用笛卡尔网格逼近轴对称壁面，当网格不够细时会形成“台阶”圆柱壁。增加网格密度，一定程度上会减弱这种效应，但需要的网格量特别大。目前大多数计算两相流的LES和DNS方法均采用笛卡尔网格体系，这也是都仅仅局限于外流场计算而很少涉及带喷嘴整体计算的一个原因，实现离心式喷嘴内流场和外喷雾场一体化计算的文献也鲜有公开报道。离心式喷嘴不同于其他类型的喷嘴，计算离心式喷嘴的外喷雾场必须获得准确可靠的喷嘴出口液膜厚度和三个方向速度分布，因此本文将Fluent采用贴体网格的优势与Ger-ris采用自适应网格界面捕捉技术的优势结合，可以比较有效地解决这一问题，从而实现离心式喷嘴从内而外的一体化准确计算。

Gerris计算外喷雾场的边界条件：Gerris中计算域均由L×L×L的基本结构Box构成，如图 4所示，本算例由两个2个Box沿流向X方向串联拼接组成，即计算域为2L×L×L，如图 5所示。计算域左端面为入口，依据传递过来的喷嘴出口截面参数给定外喷雾场入口条件，采用空心环形液相速度入口边界，中间为空心气柱的速度入口边界，内环和外环半径分别为R_m和R，分别对应着中心气柱半径和喷嘴半径。出口截面参数沿周向近似均匀分布，因而出口截面的三个方向速度分量仅仅是半径r（y，z）的函数，其中轴向速度为u_x（y，z）=u（y，z），切向速度w（y，z）和径向速度v（y，z）可分解为互相垂直的u_y（y，z）和u_z（y，z），三个方向速度分量的具体设置详见式（4）。其余面为出口，采用Outflow边界，背压为大气环境。空气为第一相，水为第二相。计算域的L=40mm，最高网格等级采用Level=10加密，最小网格约39μm。网格自适应函数设置为体积分数梯度，即网格会实时根据流场中计算的体积分数梯度大小进行自适应加密或粗化。

$ \left\{ \begin{array}{l} {u_x}\left( {y, z} \right) = u\left( {y, z} \right)\\ {u_y}\left( {y, z} \right) = v\left( {y, z} \right)\cos \left[{\arctan \left( {z/y} \right)} \right] - w\left( {y, z} \right)\sin \left[{\arctan \left( {z/y} \right)} \right]\\ {u_z}\left( {y, z} \right) = v\left( {y, z} \right)\sin \left[{\arctan \left( {z/y} \right)} \right] + w\left( {y, z} \right)\cos \left[{\arctan \left( {z/y} \right)} \right] \end{array} \right. $

(4)

为了考察验证Fluent和Gerris之间参数传递的可行性与可靠性，故而对参数传递前后的计算结果进行对比，如图 6所示。图中红色部分为Fluent计算的全流场液膜形状，黑色部分为Gerris计算的外喷雾场锥形液膜形状，可以看出两者计算的液膜形状及锥角基本吻合，一致性很好，因此可以证明这种组合传递参数的方法是准确可靠的。

3 实验系统

实验系统简图如图 7所示，实验采用的相位多普勒粒子分析仪PDPA系统测量粒径的范围为1~ 10mm，速度测量范围为0~500m/s，测量精度为1%，单点测量取样点数设为5000，采样时间设为20s，两者满足其一，该点测量终止，转换至下一位置。实验中在喷嘴出口下方不同横截面位置进行测量，每个横截面采用在半圆形范围内，每隔15°设置一个取样点，即每个截面13个取样点。高速摄影系统可对雾化破碎过程进行高速瞬态的拍摄实现可视化研究。实验采用一款高度集成的黑白高速摄影机拍摄雾化过程图像。内置32G高速读/写硬盘。具备内/外触发功能，内置的CMOS传感器为宽幅结构，最大分辨率为1280×800，单像素8/12bit可选，分辨率128×8，最高拍摄频率10⁶fps，最小曝光时间1μs。

4 结果分析 4.1 锥形液膜的雾化破碎过程计算与实验结果对比

图 8为计算的不同视角下锥形液膜三维喷雾场结构形态，可以看到液滴随着轴向距离增加，液滴粒径在渐渐减小，径向同样，外围液滴较小，小液滴速度较小，说明在此过程中大液滴经过二次雾化形成了更小的液滴。图 9（a）和图 9（b）分别为数值计算和实验拍摄的锥形液膜破碎过程图像。从图中可以看到，本文的计算结果可以细致地捕捉到雾化过程中液丝、液滴等的各种细节结构，与高速摄影拍摄的锥形液膜雾化破碎过程基本吻合。可见，基于Gerris软件建立的计算方法很好地捕捉到了锥形液膜雾化破碎过程的细节。

Gerris中库函数GfsOutputDropletSums对全场所有的液滴数目、体积进行统计输出，近似认为每个液滴呈球形，通过Matlab编程对液滴进行体积换算得到液滴直径，再通过液滴粒径分段统计处理得到了概率密度分布曲线。图 10为计算捕获的全场液滴粒径概率密度分布函数曲线，可以看到液滴呈两端少，中间多的分布，与实际雾化过程结果相符。可见，基于Gerris建立的计算模型捕捉到的液滴数目众多且尺寸变化范围大，是研究雾化过程机理及雾化特性的强有力工具。

相比之下，图 11（a）为Fluent采用CLSVOF方法计算的雾化破碎过程细节，最小网格是20μm。从结果可以看出，计算只能获得喷雾场的宏观结构特征，捕捉到大块的液团和少量的大液滴。Fluent由于使用基本均匀一致的网格，即使局部网格加密，计算量也将异常大，难以适应雾化过程瞬时变化的多尺度问题求解，因而计算效果不佳，对雾化过程的细节捕捉远不如Gerris捕捉的细腻。如图 11（b）为文献[12]中Olivier Desjardins等基于LES和改进的守恒level set方法计算的液膜破碎过程，液丝、液滴捕捉效果也不佳。可见，使用基于Gerris的高分辨率的数值方法和自适应网格技术可以有效降低网格规模和计算量，获得细致的雾化场结构和液丝液滴捕捉效果。

4.2 锥形液膜雾化破碎过程分析

锥形液膜破碎形成液滴要经过复杂的过程，图 12揭示了离心式喷嘴出口形成的液膜从初始的环形液膜，到液膜锥角渐渐张开、液膜变薄，再到前缘开始破碎脱落液丝、大液滴，液丝、大液滴又开始二次雾化形成小液滴，最终形成整个雾化的细节过程。Gerris计算的雾化过程清楚地捕捉到了雾化破碎过程的结构特征，如图 13所示。

图 14更清楚地揭示了锥形液膜前缘液丝的脱落过程，可以看到液膜表面波动不断增大，连续液膜出现穿孔破裂，随着孔洞的增长接合逐渐形成液丝，接着液丝又断裂形成液滴。图 15为液膜前缘液丝脱落过程中两相界面自适应生成的网格，在相界面进行网格细化加密，正是由于这种网格的优势，才使得液丝脱落过程捕捉效果较好。

在速度矢量图 16（a）中用颜色表示体积分数，可以较为清楚地显示涡相对于两相界面的位置，看到气液界面相互作用的现象。随着雾锥张开，液膜表面的扰动不断增大，会在相界面附近形成许多大大小小的涡，这些涡又使得液膜表面的扰动进一步加剧，直至液膜破碎成液滴。在图 16（b）的涡量图中可以更清楚地看到，液膜、液滴的表面处涡量明显大于内部的，液膜破碎处涡量也显著增大，可见气液界面处由于剪切作用使得涡量增加，这些涡是使得连续液膜产生波动进而破碎的动力。

4.3 雾化特性计算

Gerris软件在计算三维喷雾场时，获得了全场所有的液滴，对这些液滴进行统计平均获得索泰尔平均直径SMD。实验中使用高速摄影相机拍摄锥形液膜图像，测得液膜锥角，使用PDPA测量液滴平均粒径。计算的工况如表 1所示，D为喷嘴直径，d_t为切向孔直径，n为切向孔数目，L为喷嘴总长度，A为几何特性参数$ A = \frac{{D\left( {D-{d_{\rm{t}}}} \right)}}{{nd_{\rm{t}}^2}} $。

流量或喷注压降是影响喷嘴雾化效果的重要工作参数。在喷嘴变工况工作过程中，需要对流量进行调节。如果流量在一定范围内变化，雾化特性也会跟着变化。通过30g/s，35g/s，40g/s，45g/s，50g/s五种不同流量工况下喷嘴的雾化过程，研究流量变化对雾化特性的影响。

对于编号Chk-0喷嘴，在基准工况40g/s流量下，实验测得的喷注压降为Δp = 0.6MPa，3D喷嘴模型数值计算的喷注压降为Δp = 0.57MPa，相对误差为5%。

对于给定的喷嘴结构，随着流量的变化，液膜锥角和液滴平均粒径的变化分别如图 17（a）和17（b）所示，数值计算结果与实验结果的变化趋势一致，二者吻合较好，液膜锥角的最大相对误差为3.19%，液滴平均粒径的最大相对误差为4.61%，因此实验结果验证了数值方法的准确性。另外与Lefebvre，Rizk，Ben-jamin的经验公式估算的液膜锥角以及与Jasuja，Jones的经验公式估算的液滴平均粒径对比，可以看出变化趋势一致，但经验公式估算结果与计算结果、实验结果在定量上存在一定的差别。

在给定流量的条件下，仅仅改变切向孔直径的大小得到如图 18（a）和18（b）所示的液膜锥角和液滴平均粒径的变化曲线。同样数值计算结果与实验结果吻合较好，液膜锥角的最大相对误差为4.18%，液滴平均粒径的最大相对误差为11.82%。经验关系式的估算值同样存在定量上的较大差异，经验公式的计算表达式见表 2和表 3。另外，还可以看到液膜锥角和液滴平均粒径随着切向孔直径的变化较显著，因此在设计加工离心式喷嘴时，切向孔直径是一个比较重要的结构参数。

5 结论

基于Gerris采用的自适应网格技术和VOF方法，建立了模拟离心式喷嘴锥形液膜雾化破碎过程的数值方法。本文的主要结论如下：

（1）本文所采用的计算方法，能够细致地揭示锥形液膜破碎成液丝、液丝进一步破碎成液滴的全过程，能捕捉到雾化过程中的细节特征，并能捕获到全场液滴空间分布，液滴捕捉效果逼真，与高速摄影拍摄的雾化破碎过程基本吻合，突显了新方法在精细可视化研究方面的巨大优势。

（2）在外喷雾场计算方面，Fluent采用CLSVOF方法只能获得喷雾场的宏观结构，捕捉到大块液团和少量大液滴，对雾化过程细节的捕捉效果远不如Gerris。Gerris捕捉到的液滴数目众多且尺寸变化范围大，液滴粒径分布与实际雾化过程相符，在液膜破碎过程中涡显著增多，它的形成对相界面的发展、破碎有明显作用。建议在初步计算外喷雾场时采用Fluent，精细计算液滴分布时采用Gerris。

（3）采用本文方法计算不同流量和不同切向孔直径下的液膜锥角、液滴平均粒径与实验获得结果吻合较好，液膜锥角的最大相对误差为4.18%，液滴平均粒径的最大相对误差为11.82%，表明计算方法相比较于经验关系式具有更高的精度，为精细研究离心式喷嘴雾化机理及准确评估雾化特性提供一种新的手段。同时通过研究发现，切向孔直径对液膜锥角和液滴平均粒径的影响较显著，在设计加工离心式喷嘴时，是一个比较重要的结构参数。

致谢 此项工作是在国家超级计算天津中心的“天河一号”超级计算机上完成，感谢“天河一号”的大力支持。