针栓式喷注单元膜束撞击雾化混合过程数值模拟

王凯，雷凡培，杨岸龙，杨宝娥，周立新

1. 西安航天动力研究所 液体火箭发动机技术重点实验室， 西安 710100 2. 中国船舶工业集团有限公司，北京 100044

与双组元液体火箭发动机上常用的典型的撞击式或同轴式喷嘴不同，针栓式喷注器具有独特的几何特性和喷注特性，采用轴向液膜与径向液膜或者径向孔液束垂直撞击形成喷雾扇雾化[1]，具体原理结构如图1所示。针栓式喷注器具有诸多优点：推力调节简单、变推力工况下能够产生很高的燃烧效率(典型值为96%～99%)、面关机、固有燃烧稳定性好、可在很宽的推力水平内按比例放大或缩小等工作特性[2]。不用设置诸如声腔或隔板之类的稳定装置，相比传统火箭发动机的喷注器由成百上千个精细的喷嘴组成，针栓式喷注器的结构大为简化，意味着可靠性提高，同时可以快速更换零组件实现改进升级[3]，减少质量和降低成本。这些工作特性给火箭发动机的设计、性能、稳定性和试验灵活性带来了极大的好处，也被看成是降低现有运载器大发动机成本的一种有效手段。

针栓式喷注器的概念最早出现于20世纪50年代喷气推进实验室(Jet Propulsion Laboratory, JPL)设计的一种简单而精致的试验装置[4]。后来经过TRW公司的发展与演变，形成了著名的针栓式发动机LMDE(阿波罗登月舱下降发动机)[5-6]和目前技术水平最高的SpaceX公司的Merlin 1D系列发动机[7-8]。关于针栓式喷注器的基础研究较少，前期主要在于TRW公司，依靠大量试验数据和丰富设计经验。他们的研究发现除了针栓式喷注器单元的几何尺寸参数外，重点关注的有4个无量纲参数，分别是动量比、阻塞比、越程比和直径比，2个喷注特性参数为喷注速度和雾化角。从TRW的研究成果可以发现，雾化角是其中相当重要的一个雾化特性参数，它对雾场结构和液雾空间分布起着决定性作用。于是以前研究者大多针对雾化角开展研究，以试验研究为主。代表性的研究主要有Heister基于横向射流雾化角模型提出针栓式喷注器液膜撞击的雾化角模型[9]；Boettcher等通过纹影法观测气针栓式喷注器的雾化结构，发现雾化角随动量比增大而增大，并提出了雾化角和动量比的计算关系式，并对比了三维环形膜与二维平面膜的理论结果，表明两者的雾化角接近，二维平面膜的雾化角略高于三维环形膜的，在计算针栓式喷注器雾化角时二维平面膜的假设是合理的[10]。于是Sakaki等率先通过平面针栓式喷注单元使用高速摄影相机拍摄了喷雾场结构，研究了动量比和喷雾场结构的关系，并将平面针栓与轴对称针栓特性进行了对比，结果表明两者特性的基本趋势一致[11-13]。这就为用平面针栓式喷注器替代实际的轴对称针栓式喷注器开展深入研究提供了有力证据。后来成鹏等开展了针栓式喷注单元高速摄影试验研究，拍摄了雾化过程图像，得到了雾化角和液滴粒径，并与理论推导的理论公式进行了对比[14-15]。

尽管光学测量试验手段很先进，然而由于针栓式喷注器独特的喷雾场结构，不仅范围大，而且液雾过于稠密，现有的光学测量设备如高速摄影、PDPA(Phase Doppler Particle Analyzer)等很难获得有效的喷雾场参数(液滴尺寸分布、混合比分布、流强分布等)，同时详细雾场结构也看不清楚。因此，数值仿真便成为研究针栓式喷注器喷雾过程机理和详细破碎混合特性的主要手段。相对于其他类型的火箭发动机喷嘴的数值仿真研究，目前采用CFD软件计算针栓式喷注器喷雾燃烧特性的研究成果极少。近期最新的研究成果是刘昌波为了深入认识影响针栓式喷注器雾化特性的因素及影响规律，基于虚网格LPT(Lagrangian Particle Tracking)法耦合VOF(Volume of Fluid)方法对液膜/液膜和液膜/液束撞击形式的针栓喷注器雾化过程进行了仿真研究[16]；郑刚使用Fluent中CLSVOF(Coupled Level Set and Volume of Fluid)方法对针栓式喷注器一次雾化过程进行了仿真计算[17]，采用了非常细致的网格划分加密，对雾化过程中液膜及大的液团液丝特征捕捉比较好，较之前研究者使用类似手段的研究结果更进了一步。然而以上2种方法虽有效捕捉了一次雾化过程，但在进一步捕捉液滴和二次雾化过程计算方面发展的潜力不大，同时他们仅对雾化特性的定性影响规律关注较多，未深入研究液膜/液束相互撞击作用机制及雾化过程细节。另外，2016年Son等为了克服仅仅靠试验难以深入分析雾化破碎和混合特性的问题，使用Fluent对简化的二维轴对称针栓式喷注器模型进行了喷雾过程仿真，与之前的试验形成验证补充，并分析了速度分布、雾化角及液体分布等[18]；2017年该团队的Radhakrishnan等基于DPM(Discrete Phase Model)方法开展了喷注条件对混合特性的影响，重点关注二维轴对称针栓式喷注器计算模型的雾化角和混合质量与动量比的关系，计算中采用模型化的K-H(Kelvin-Helmholtz)模型和Wave模型来引入一次破碎和二次破碎对雾化混合的影响[19]，虽然考虑了雾化破碎过程对混合的影响，但采用的雾化破碎模型有一定的缺陷，不能准确反映雾化破碎过程，计算的混合质量存在准确性问题。

图1 针栓式喷注器原理图Fig.1 Schematic diagram of pintle injector

近年来Popinet基于四叉树/八叉树的自适应网格(Adaptive Mesh Refinement，AMR)和VOF方法开发了一种新的开源软件Gerris[20-21]，在处理多相流问题中有着很大的优势，可以精确地描述雾化过程。国外一些研究人员已经将Gerris很好地运用于直射流和互击射流雾化过程研究，并得到了实验验证。如Salvador[22]和Mehravaran[23]等分别开展了低喷注压降和中等速度的直射流雾化过程仿真研究；Ma等使用Gerris对射流撞击的雾化形式进行了高精度的数值模拟研究[24]；Chen等采用Gerris对撞击式射流雾化特性及雾化破碎过程进行了高分辨率的数值模拟计算，得到了非常好的可视化效果[25-26]。国内的研究者主要有西安航天动力研究所的李佳楠和王凯、清华大学的张培玉分别采用Gerris对互击式喷嘴、离心式喷嘴和撞击雾化及旋流雾化做了系统性的仿真研究[27-29]，特别是李佳楠和王凯都通过喷雾实验从定性和定量两方面对雾化全过程和仿真计算结果进行了对比验证，雾化过程特征定性吻合很好，定量计算的雾锥角和液滴平均粒径最大相对误差分别在5%和10%左右[30]，为Gerris在雾化过程和雾化特性方面的推广应用提供了先验支撑。然而目前在国内外公开发表的资料中，还未看到将Gerris运用于针栓式喷注器雾化过程的仿真计算，更缺少对针栓式喷注器喷雾场结构特征的认识，这主要是由于针栓式喷注器靠2种不同推进剂撞击，之前数值模拟多将两路推进剂视为同相处理，与环境气体合起来仅有两相，即气相和液相，这对于模拟自击式喷嘴撞击雾化过程没有问题，然而对于2种推进剂互击式的喷嘴便无能为力，不能真实反映并获得2种推进剂撞击后各自的变形及雾化后的混合情况。

鉴于针栓式喷注器喷雾场基础研究认识的不足以及模拟2种不同推进剂撞击的难点，本文以平面针栓式喷注单元为研究对象，基于Gerris建立了采用三相VOF方法分别识别两路液体的计算模型(即将2种推进剂视为不相溶的液体，之间存在相界面)，对膜束撞击作用过程及雾化混合过程进行细致仿真，首先通过高速摄影试验对计算模型进行定性和定量对比验证，接着对液膜液束撞击喷雾场结构、撞击变形过程、流场涡结构和雾化破碎典型特征进行捕捉和分析，最后对破碎后的雾化混合分布特征进行识别分析。首次通过数值仿真从近场撞击变形到喷雾扇形成再到雾扇液膜的雾化混合开展系统研究，从而为认识和揭示液膜液束相互作用机制及雾化混合分布规律奠定重要基础。

1 数值方法

1.1 控制方程

数值仿真采用开源软件Gerris作为计算工具。Gerris采用的数值方法Popinet等在2009年发表的文章中作了详细的描述，下面对Gerris采用的数值方法做一个总结性的描述。Gerris求解的是不可压、变密度、带有表面张力的Navier-Stokes (N-S)方程[31]:

(1)

采用气相、液相1和液相2三相进行计算，分别定义液相1体积分数c1(x,t)和液相2体积分数c2(x,t)，得到的流体密度和黏性系数为

(2)

式中：ρl1、ρl2和μl1、μl2分别是液相1和液相2的密度和动力黏度；ρg和μg分别为气相的密度和动力黏度。

在Gerris中使用经典的时间分裂投影法进行简化，达到时间二阶精度。使用四叉树/八叉树进行空间离散，达到空间二阶精度，使得自适应加密算法(根据流场参数变化对局部网格进行动态加密或粗化)可简易灵活地实现，在不损失计算精度的情况下显著降低了计算量，非常适合处理多尺度流动问题。使用分段线性的VOF几何重构方法进行自适应网格界面捕捉，非常适合应用于包含破碎、聚合现象的雾化过程计算，图2所示为计算生成的自适应网格。通过将表面张力转化为某一区域连续的体积力并结合高度函数曲率估计实现表面张力的精确求解。采用单调集成大涡模拟(MILES，又称隐式大涡模拟ILES)[32]近似模拟亚格子SGS的能量传递，这是由于数值计算不可避免地有数值耗散。

图2 基于分段线性VOF几何重构方法的自适应网格加密Fig.2 Adaptive mesh refinement based on piecewise- linear geometrical VOF

1.2 气液两相界面捕捉

VOF方法通过定义体积分数函数c来描述界面(含有运动界面的网格满足0

分段线性的VOF方法(PLIC VOF)对界面进行重构，图3所示红色斜平面描述的是不同情况下的相界面，其中网格单元为单位长度的正方体，相界面的单位法向量n由液相一侧指向气相一侧，α为相界面斜平面和原点之间的最小距离，液相的体积分数c(体积)可由n和α唯一确定表达，H(x)为Heaviside阶跃函数。

图3 分段线性的VOF界面重构方法Fig.3 Piecewise-linear VOF interface reconstruction method

1.3 计算模型

计算所选用的结构参数及计算域如图4所示，轴向液膜厚度为h，轴向速度为u1，径向液束直径为d，径向速度为u2，流体密度分别为ρ1和ρ2。计算域由8个L×L×L基本结构Box构成，计算域左端面和底端面分别为轴向液膜和径向液束入口。图中标注的面为无滑移壁面，其余面为出口，采用Outflow边界，背压为大气环境。第一相为空气，第二相和第三相均为水，对应两路液体。计算域L=10 mm，最高网格等级采用Level=9加密，最小网格约19.5 μm，最低Level=6，且相邻网格等级之差不超过1级，这使得纯气相的网格尺寸绝大部分也均小于100 μm。网格自适应函数设置为体积分数梯度，即网格会实时根据流场中计算的体积分数梯度大小进行自适应加密或粗化。计算主要以追踪相界面运动及破碎过程为主，同时相界面附近也是各种参数变化较显著的位置，可以保证流场中速度、压力等都获得较高的计算精度。设定轴向液膜与径向液束撞击前的距离(即跳跃距离)为6 mm，对应模拟的针栓柱直径D=20 mm，液膜展向宽度取10 mm。

图4 计算域示意图Fig.4 Sketch map of computation zone

2 算例验证

将液膜与液束撞击后会形成的喷雾场数值仿真结果与高速摄影拍摄的试验结果进行对比，如图5所示，从不同视角定性看两者的喷雾扇形状非常相似。对试验获得的喷雾瞬态图像进行图像增强处理后，将1 000张增强后的瞬态图像进行平均处理，再取雾扇边界平均中心来测量获得雾化角。按照如图6所示的方式获取雾化角θ，将两者对比获得如图7所示的曲线图，可以看到数值仿与试验结果吻合较好。

然而液膜与液束撞击形成的喷雾扇是空间三维结构，仅通过定量对比试验和数值仿真的雾化角来验证数值仿真结果是不够充分的。因此，为了进一步对比试验结果和数值仿真结果，除了在喷雾扇侧视图中定义雾化角之外，又在喷雾扇正视图中补充定义了雾扇扩张角，如图8所示可测量获得雾扇扩张角。

试验和数值仿真中均选择轴向液膜厚度h=0.25 mm、径向液束直径d=1.0 mm，通过改变液膜速度和液束速度来改变工况，具体的工况参数如表2所示。

表1 液膜撞击液束的结构参数Table 1 Geometry parameters of different injectors with liquid sheet impinging on liquid jet

图5 液膜与液束撞击的喷雾扇试验结果Fig.5 Experimental results of spray fan of impingement between liquid sheet and liquid jet

图6 液膜与液束撞击的雾化角Fig.6 Spray angle of impingement between liquid sheet and liquid jet

对比9种不同工况下试验测量和数值仿真的雾扇扩张角，如图9所示，可以看出两者吻合较好，最大相对误差约为8.55%；两者均随有效动量比的增大而增大。因此，通过雾化角和雾扇扩张角2个定量参数从多个角度对比了试验结果与数值仿真结果，表明采用的数值方法计算的喷雾扇空间三维结构具有较好的准确性。

图7 不同动量比下的雾化角变化曲线Fig.7 Spray angles in various momentum ratios

图8 试验拍摄的喷雾场图像Fig.8 Experimental image of spray field

表2 液膜撞击液束的工况参数表Table 2 Different operating parameters with liquid sheet impinging on liquid jet

图9 雾扇扩张角试验结果与数值仿真结果对比Fig.9 Experiment and simulation results of divergence angle in comparison

3 数值仿真结果分析

3.1 液膜/液束撞击喷雾场结构

液膜与液束撞击后形成如图10所示的扇形喷雾场结构。通过数值仿真发现喷雾场结构会随着动量比的大小而发生变化。最明显的变化就是雾化角会随着动量比的增大而增大。随着动量比的增大，雾场结构基本相似。然而当动量比大到一定程度时，喷雾场会发生一个显著变化，液束会穿透液膜形成如图11(b)所示的喷雾扇。通过液膜与液束分相识别的仿真计算可以清楚地看到，小动量比和大动量比下形成液束未穿透液膜和液束穿透液膜2种不同的喷雾扇结构。从图11(a)可以看到，液膜与液束撞击后形成一股液膜“包裹”液束的合成雾扇，两者一起向下游运动破碎，产生大量液滴，一定距离内合成雾扇完全破碎；而从图11(b)可以看到，液膜与液束撞击后形成液束直接穿透液膜的雾扇，液束单独向下游运动，液束保持很长距离仍未完全破碎。总之，径向液束能够穿透轴向液膜并持续射向室壁，或者径向液束与轴向液膜合成“包裹”液束以较大角度射向室壁，撞壁后在室壁铺展形成有效的冷却液膜及在室壁边区形成合适的混合比分布，这2种喷雾场结构都可作为提高燃烧室身部可靠冷却的喷雾场备选构型。

图10 液膜与液束撞击扇形喷雾场结构Fig.10 Fan-shaped spray field structure of liquid sheet impinging on liquid jet

图11 液膜撞击液束的两种不同喷雾场结构Fig.11 Two different spray field structures of liquid sheet impinging on liquid jet

3.2 液膜/液束撞击变形过程

在撞击过程中液膜与液束均会发生弯曲变形和横截面变形，液膜也会发生绕液束的类圆柱绕流流动。仅有一定宽度的部分液膜与液束发生有效撞击，中心正撞滞止点达到最大雾化角，其余展向位置雾化角小于该撞击雾化角。

图12为不同垂直高度的横截面，图中，w为仅与径向液束发生有效撞击的液膜宽度，认为近似等于d；L为针栓喷注单元宽度，即相邻喷注单元间距，从图12的结果可以看出，液膜绕液束的流动与圆柱绕流有相似之处，同时也存在诸多不同之处。相似的是液膜绕圆柱射流后会在射流后缘两侧形成一系列对转的尾迹涡串结构，且涡自圆柱射流后缘产生，刚产生时涡较小，随着流向下游，涡逐渐增大。不同之处在于，液束受到撞击作用后会发生变形，横截面展向拉伸，致使喷雾扇的展向宽度随高度增加而增大。另外一个显著不同之处是液膜绕过液束后形成了分叉流动，底部液膜被分开成具有一定宽度的气腔区，自射流位置开始，沿着轴向气腔区的宽度逐渐增大。同时从三维空间角度来看，在射流后面形成了一个半开放的喇叭口(“Ω”)形喷雾扇腔，横截面呈“n”形。喷雾扇腔内的气体流动如图13所示，中心气体向内流动，贴液体附近的气体随着液体一起向外流动，这很容易理解，由于剪切层内黏性作用，液体表面附近的气体随液体一起向外流动，在射流根部处形成低压区，导致外界气体从中心倒灌，从而形成特殊的腔内流场结构。沿轴向液膜流动方向看喷雾扇，可以看到如图14所示的结构，液膜与液束撞击后形成的雾扇展向宽度随着垂直高度增加而逐渐增加，形成一定的展向角。这与液膜/液束撞击后两者的变形相互作用密切相关。

图12 液膜液束撞击变形流场三维结构Fig.12 3D flow field structure of impingement and deformation between liquid sheet and liquid jet

图13 液膜撞击液束形成的腔内速度场分布Fig.13 Velocity field distribution in hollow cavity formed by impingement between liquid sheet and liquid jet

从图13中可以清楚地看到,液束撞击液膜后轴向受压变形，展向变宽铺展成薄液膜，类似两股射流撞击和射流撞壁变形过程。另外，径向液束在轴向液膜的挤压和剪切双重作用下，液束横截面也发生变形，从初始的圆形逐渐发展形成如图15 所示的形状。可以看到液束迎风面受到正应力压缩作用，滞止点正应力较大，压力也较高，产生的压力梯度驱动由内而外的流动；同时两侧又受到切应力剪切作用，产生由前缘外向后缘内的流动，使得液束横截面两侧开始形成向背风面发展的表面褶皱，尾端出现拖尾现象，最后随着液束弯曲变形，其迎风面与背风面逐渐靠近，展向不断被拉伸呈膜状，厚度也逐渐减薄并趋于均匀。

图14 液膜撞击液束轴向视角的雾扇结构Fig.14 Spray fan structure impingement between liquid sheet and liquid jet in axial view

图15 液束横截面变形过程Fig.15 Cross-sectional deformation process of liquid jet

3.3 流场涡结构

当轴向液膜与径向液束相撞时，由于液体流动的速度差和黏性，在界面处会存在剪切层，产生复杂的涡流现象：对漩涡对、尾迹涡、前缘涡和马蹄形涡等。下面类比如图16[35-36]所示的横向射流流场中的涡结构进行分析。

当轴向液膜流过径向液束时，剪切层涡在轴向液膜与径向液束相互作用下涡旋面发生转向折叠，这助力形成了对漩涡对(Counter-rotating Vortex Pair，CVP)。CVP持续向下游发展，并主导了下游流场的流动特征，如图17(a)所示。这与图15所示的径向液束横截面发生的变形过程密切相关，液束横截面的变形进一步在流场中发展，诱导形成了这对由外向内卷、彼此方向相反的肾形涡流结构，此结构向下游发展最终导致了CVP的形成，这也与“Ω”形喷雾扇的形成有关。肾形结构及CVP会将包裹于液束外面的轴向液膜不断捲向径向液束，使得两者混合较为均匀。可以看到越往下流，由于CVP的增大，使得液膜与液束均向展向拉伸展开，使得两者的接触面增大，流场中相互作用区域增大，混合也越均匀。因此，经过分析可以表明由轴向液膜流与径向液束流相互作用产生的涡对混合有一定的促进作用，其中CVP的尺度最大，它的出现增强了针栓的雾化混合。

图16 横向射流流场中涡结构[35-36]Fig.16 Vortex structure in cross-flow field[35-36]

另外，液膜与液束撞击后，液膜绕液束流动会形成分叉流动，在液束根部后面形成类似圆柱绕流的尾迹涡区(Wake vortex)，如图17(b)所示，这些尾迹涡源于液膜与液束根部非稳态作用的分离涡，具有与圆柱绕流下游涡街类似的特性。径向液束在前缘驻点附近与平板壁面的液膜边界层发生相互作用产生马蹄形涡系(Horseshoe Vortex)。

图17 流场中的涡结构Fig.17 Vortex structure in flow field

液膜越厚，液束直径越大，剪切作用面越大，流速越大，剪切作用也越强。在某些相互作用强的工况可以看到，在液束的迎风面形成了一连串逆时针旋转的涡(前缘涡，Leading-edge Vortex)，如图17(c)所示，这些涡位于液膜与液束的相界面之间，源于液束根部的驻点涡发展演化而来，根部的驻点涡是撞击点形成的高压区强逆压梯度导致液膜边界层分离造成的。当形成的驻点涡使得边界层部分回流，回流后滞止点压力略降低，致使液膜与液束之间相互作用减弱，回流区减小，再次使得液膜与液束直接相互接触产生强的相互作用，进而再次形成高压区，致使驻点涡变大，该过程周而复始，使得液膜与液束相界面上的作用位置和形状发生周期性变化，使得驻点涡结构也会周期性脱落，从而加剧液膜液束的相互作用，致使两者过早破碎。图17(c)中对应的速度矢量场结构可以清楚地说明此现象的成因。受这些涡的影响，液束背风面也出现了明显的表面波动。这些涡流结构增大了液膜与液束作用面积，增强了液膜表面不稳定性，强化了雾化和混合。

喷雾场中还存在一些小尺度的涡，图18为使用Q准则[37]进行流场涡识别获得的结果与喷雾扇的叠加显示，可以看出液膜破碎产生的液丝空间发展方向似乎与这些涡有着密切关系[38]。

综上分析，可以发现与液膜/液膜撞击仅产生正向应力作用相比，液膜与液束撞击除了有正向应力的作用外，液束两侧也受到剪切应力作用，此剪应力导致了上述一系列复杂涡流现象，使得两路流体相互作用增强，雾化和混合均增强，这也是液膜撞击液束的喷注单元构型优于液膜撞击液膜的本质原因。

图18 Q准则识别的流场涡结构与喷雾扇叠加 (Q=0.04)Fig.18 Vortex structure identified by Q criterion and superimposed on spray fan(Q=0.04)

3.4 雾化破碎典型特征

液膜与液束撞击形成的喷雾扇破碎成液滴要经过复杂的过程，图19揭示了膜束撞击形成的喷雾扇液膜从初始铺展成膜，到液膜铺展张开、液膜变薄，再到前缘开始破碎脱落液丝、大液滴，液丝、大液滴又开始二次雾化形成小液滴，最终形成整个雾化的细节过程。Gerris计算的雾化过程清楚地捕捉到了雾化破碎过程的结构特征，如图19中液膜前缘液丝的脱落破碎、液膜表面波动及相应的波峰波谷。液膜液束各自铺展形成的膜完成破碎后主要分布在雾扇中心区液膜附近。轴向液膜受到径向液束撞击后被顶起形成雾扇，雾扇与轴向夹角区的液膜分布量最少，且最先减薄破碎形成少量小液滴，倾斜向下游运动。轴向液膜下游气腔区边缘受到液束的作用后也向上翻起，形成的翻边液膜变薄后产生向上的凸起液柱/液丝并破碎，破碎产生的少量大液滴向中间区域运动。

图19 数值仿真的喷雾场典型特征Fig.19 Typical characteristics of spray field in simulation

Hopper和Dombrowski[39]在1953—1963年间应用高速摄影技术将平面液膜破碎过程划分成4种典型的情况：边缘破碎(Rim Disintegration)、波动破碎(Wave Disintegration)、液膜穿孔破碎(Perforation Disintegration)和湍流破碎(Turbulence Disintegration)。参考上述平面液膜典型的破碎模式，依据液膜韦伯数Weh和液束韦伯数Wej从小到大的增长，将液膜撞击液束形成的喷雾扇破碎过程分为如图20所示的3种典型模式：边缘破碎(Rim Disintegration)、波动穿孔联合破碎(Wave and Perforation Disintegration)、和湍流破碎(Turbulence Disintegration)。低的液膜/液束韦伯数对应着边缘破碎模式，中等液膜/液束韦伯数对应着波动穿孔联合破碎模式，高的液膜/液束韦伯数对应着湍流破碎模式。

边缘破碎：由于液体表面张力的作用使液膜中间较薄的部分向两侧边缘处收缩，形成较厚的边缘，随后在气动力、表面张力的作用下开始破碎。当流速很低时，液体的黏性和表面张力作用相对显著。这种方式通常生成的液滴较大，主要分布在边缘下游附近，中间较薄的雾扇液膜下游少有液滴，典型的分布如图20(a)所示，Weh=300，Wej=280。

波动穿孔联合破碎：波动破碎模式与穿孔破碎模式往往相伴出现，较难单独区分开来，故此处将两者统一称为波动穿孔联合破碎模式。对于波动破碎，液膜上表面波不断增长，直至半个波长或一个波长的液膜被撕裂下来，形成液丝或液片，在表面张力作用一下收缩成液滴。对于穿孔破碎，在离开喷嘴一定距离处，液膜出现孔洞，孔洞尺寸不断变大，相邻孔洞间形成液带或液丝，接着液带和液丝相互分离，最后分离的液带和液丝再破碎成不同尺寸的液滴。该联合模式下在整个雾扇下游范围内均有液滴分布，雾扇脊部附近液膜主要以波动破碎模式为主，在K-H (Kelvin-Helmholtz)不稳定的作用下生成的液滴较大，尺寸也不均匀；中间较薄的雾扇液膜区主要以穿孔破碎模式为主，生成的液滴较小，尺寸也较均匀；雾扇根部液膜边缘受原轴向液膜下漏的影响，雾化较少，表面存在波动，在P-R (Plateau-Rayleigh)和R-T(Rayleigh-Taylor)不稳定的作用下产生垂直于液膜表面的凸起液柱或液丝，进一步发展从液膜表面脱落产生少量液滴，典型分布如图20(b)所示，两幅图对应的韦伯数分别为Weh=400，Wej=2 180和Weh=750，Wej=2 150。

湍流破碎：当液体的喷射速度很大时，液膜与液束在撞击后很快便破碎成小液滴，几乎看不到连续的液膜。雾扇根部液膜边缘，同样产生垂直于液膜表面的凸起液柱或液丝，进一步发展从液膜表面脱落产生少量液滴。在低动量比时，雾化主要由湍流液膜主导，液滴主要分布在轴向液膜附近的空间；在高动量比时，雾化主要由湍流液束主导，大量液滴主要由径向液束破碎产生，分布在液束附近的空间，雾扇中间液滴分布非常少；在中等动量比时，雾化由湍流液膜和湍流液束撞击共同主导，液滴在整个雾扇空间内均有分布。典型的分布如图20(c)所示，3幅图依次对应的韦伯数分别为Weh=4 610，Wej=2 140；Weh=2 840，Wej=9 950和Weh=1 265，Wej=10 600。

图20 喷雾场典型的破碎模式Fig.20 Typical fragmentation patterns of spray field

经综合分析雾化角(动量比)和破碎模式(韦伯数)2个特性参数，认为选取中等动量比下的波动穿孔联合破碎模式和湍流破碎模式，可以保证在有较大雾化角的前提下，喷雾场破碎产生的液滴空间分布也较均匀。

3.5 雾化混合分区结构特征

液膜与液束撞击后形成的喷雾扇在下游会破碎形成大量的液滴，喷雾场结构不同，液滴的空间分布区域也不相同。由于采用了分相计算技术，因此可对液膜和液束破碎产生的液滴进行分别捕捉，获得2种液滴的空间分布。

从图21可以看出，随着动量比的增大，液滴粒径显著减小。动量比越大，液膜与液束的撞击相互作用越强，使得撞击产生的液滴粒径也越小。同时也可以看出，液膜和液束运动的主路径上液滴粒径稍偏大，附近运动的液滴粒径较小。大液滴的运动速度自然也越大，附近小液滴的速度也越小，更容易受周围气体扰动，产生横向运动，离开主流运动方向。

通过数值仿真发现，液膜与液束撞击形成的喷雾场雾化和混合分布呈现分区结构特征，如图21 所示，分别是液束控制主导的上雾化区、液膜控制主导的下雾化区及夹在中间的混合区。从高中低3种不同动量比下的雾化混合分区结构可以看出，随着动量比的增加，由液束控制主导的上雾化区逐渐拉伸变长，区域变大，这很容易理解，动量比增大，液束的动量增大，雾化角也增大，液束越容易穿透液膜沿着垂直方向分布，破碎距离也越大，生成的液滴也会沿着垂直方向分布。随着动量比的增加，由液膜控制主导的下雾化区则变化不大，这主要是由于本文计算中动量比的增大主要是依靠增加径向液束的动量来改变的，轴向液膜的动量基本没变，因而下雾化区基本没变化。随着动量比的增加，中间混合区逐渐缩小，这主要是由于随着动量比增大，雾化角增大，液束越易穿透液膜，液束的雾化变差，从而使得液束控制主导的上雾化区与液膜控制主导的下雾化区分离，中间的雾化区也自然会减小。

在分区结构中，混合区的2种推进剂混合比较均匀，而上下雾化区中2种推进剂混合不足。当动量比变化时，液膜撞击液束的雾化混合比分布也在变化，可以看出动量比较大工况下，混合比分布不够均匀；动量比较小工况下，雾化后的液滴粒径不够小。在中等动量比(CMReff=1～2)，较高的液膜液束速度下，雾化混合较好。

图21 高中低3种不同动量比下雾化混合分区结构Fig.21 Spray and mixing regional structures in three different momentum ratios

总之，从两路液体雾化的角度来看，动量比越大，雾化区的液滴粒径越小，越有利于液滴的蒸发燃烧；从两路液体混合的角度来看，动量比越小，上下雾化区越靠近，中间混合区越大，越有利于2种液体推进剂混合比分布均匀。因此，实际中应兼顾雾化特性和混合特性，选取中等动量比的液膜与液束撞击。

4 结 论

为了全面认识针栓式喷注器喷雾场结构，本文采用自适应网格技术和分三相计算的PLIC VOF方法对针栓式喷注单元膜束撞击雾化混合过程进行仿真分析研究，可以得到以下结论：

1) 采用了自适应网格技术和分三相计算的PLIC VOF新方法，对针栓式喷注单元两路推进剂分别进行界面追踪，首次刻画了膜束撞击雾化混合过程的精细物理图画，获得了膜束各自撞击变形、液膜破碎、表面波及液丝液滴等雾化特征，这些特征通过常规的CFD软件(如Fluent、OpenFoam等)难以获得。另外与高速摄影试验结果定性定量对比均吻合较好，表明新的仿真方法在精细研究喷雾场方面具有较好的准确性，可为认识和揭示针栓式喷注器雾化机理提供一种有效的研究手段。

2) 通过分相识别的数值仿真发现膜束撞击形成液束未穿透液膜和液束穿透液膜2种不同类型的喷雾扇结构，即撞击后形成一股液膜“包裹”液束的合成雾扇和液束穿透液膜的雾扇。2种结构的液相空间分布不同，均可作为提高燃烧室身部可靠冷却的喷雾场备选构型。

3) 膜束撞击形成的喷雾扇呈半开放喇叭口(“Ω”)形，横截面呈“n”形，中间形成空心气腔，弧形雾扇展向呈一定角度展开，膜束同时发生弯曲变形和横截面变形。液膜撞击后被径向拉伸变形成薄液膜，同时绕液束形成分叉流动，类似圆柱绕流；液束撞击后轴向受压变形，展向变宽铺展成薄液膜，类似两股射流撞击和射流撞壁变形过程。

4) 与膜膜撞击仅产生正向应力作用相比，膜束撞击时液束两侧还受到剪切应力作用，这导致了喷雾场中一系列复杂涡流现象，使得两路流体相互作用增强，雾化和混合均增强。

5) 经综合分析雾化角(动量比)和破碎模式(韦伯数)2个特性参数，认为选取中等动量比下的波动穿孔联合破碎模式和湍流破碎模式，可以保证在有较大雾化角的前提下，喷雾场破碎产生的液滴位置的空间分布也较均匀。

6) 膜束撞击形成的喷雾场液滴分布呈现分区结构特征，分别是液束控制主导的上雾化区、液膜控制主导的下雾化区及夹在中间的混合区。实际中应兼顾雾化特性和混合特性，选取中等动量比膜束撞击。