近年来随着商业航天的迅猛发展，特别是SpaceX等商业航天公司的迅速崛起，国内外掀起了高性能、长寿命可重复使用、低成本、推力深度可调的液体火箭发动机技术研究的热潮。在变推力液体火箭发动机设计中，喷注器的设计尤为重要。与双组元液体火箭发动机上常用的典型撞击式或同轴式喷嘴比较，针栓式喷注器具有独特的几何特性和喷注特性，变推力工况下能够产生很高的燃烧效率(典型值为96%~99%)，并具有调节能力强、面关机、操作安全、固有燃烧稳定性好等工作特性^[1]。另外，针栓式喷注器可以在很宽的推力水平内按比例放大或缩小，并可适应不同的推进剂组合，不用设置诸如声腔或隔板之类的稳定装置。相比传统火箭发动机的喷注器由成百上千个精细的喷嘴组成，针栓式喷注器的结构大为简化，意味着可靠性提高，同时可以快速更换零组件实现改进升级^[2]，减少质量和降低成本。这些工作特性给火箭发动机的设计、性能、稳定性和试验灵活性带来了极大的好处，也被看成是降低现有运载器大发动机成本的一种有效手段。

针栓式喷注器典型的工作原理是通过伸入燃烧室内部的针栓结构，使一种推进剂流经针栓中心流道，并由针栓头端附近的一系列孔(或缝隙)呈放射状径向喷注；另一种推进剂通过针栓外侧的环形缝隙轴向喷注，两者呈90°交叉撞击，使推进剂雾化混合^[3]，具体原理结构如图 1所示。针栓式喷注器的概念最早出现于20世纪50年代喷气推进实验室(Jet Propulsion Laboratory, JPL)设计的一种简单而精致的试验装置。后来1960年TRW公司开始实践和研制针栓式喷注器，其后的40年间，至少有60多种方案进行了热试车^[4]。最著名的针栓式发动机为阿波罗登月舱下降发动机LMDE，其采用N₂O₄/A-50，推力变比10:1 (4.67~46.7 kN)^[5-6]。从20世纪80年代开始，TRW又对针栓式喷注器进行了一系列设计改进，实现了喷注面关机及凝胶推进剂的成功应用。从20世纪90年代开始，在进一步提升常规推进剂针栓式发动机性能的同时，重点开展了液氧/液氢^[7-8]、液氧/煤油^[9]、液氧/甲烷^[10-11]等针栓式发动机研究，在低温推进剂应用、低成本方案、按比例放大设计上取得了显著的成果，并将其应用于大推力针栓式发动机的研制。如2 900 kN的TR106 LOX/LH2发动机^[12-15]，SpaceX的Merlin 1D系列发动机^[16-17]。

中国对针栓式变推力发动机的研究主要从20世纪70年代开始，在2013年中国研制的5:1月球探测用7 500 N针栓式发动机成功将嫦娥CE-3送到了月球表面，2019年又将CE-4首次送到月球背面。为了实现载人登月，中国也开展了用于登月舱下降的液氧/甲烷针栓式变推力发动机方案研究^[18]。经过几十年的研究，中国已经基本掌握了针栓式发动机的设计技术，然而与工程研制的辉煌成就相比，针栓式喷注器工作特性的理论研究尚不足，许多机理尚不清楚。因此需要针对一些关键结构和工作参数进行系统深入的研究。

前期关于针栓式喷注器的研究主要集中在TRW公司，主要依靠大量试验数据和丰富的设计经验。他们的研究发现除了针栓式喷注器单元的几何尺寸参数外，重点关注的无量纲参数有动量比C_TMR和阻塞率C_BF，喷注特性参数有喷注速度和雾化角θ。

从TRW的研究成果可以发现，雾化角是其中相当重要的一个雾化特性参数，它对喷雾场结构和液雾空间分布起着决定性作用，然而对雾化角的准确预测仍然没有确切的计算公式。Heister基于横向射流雾化角模型提出针栓式喷注器液膜撞击液膜的雾化角模型为θ=180/π·C_TMR^0.5[19]。Boettcher等通过纹影法观测气气针栓式喷注器的雾化结构，发现雾化角随动量比增大而增大，并提出了雾化角和动量比的计算关系式，即$ {\rm{cos}}\;\theta = 1/\sqrt {1 + C_{{\rm{TMR}}}^2} ({\rm{tan}}\;\theta = {C_{{\rm{TMR}}}}) $，并对比了三维环形膜与二维平面膜的理论结果，表明两者的雾化角接近，二维平面膜的雾化角略高于三维环形膜的，在计算针栓式喷注器雾化角时二维平面膜的假设是合理的^[20]。Fang和Shen通过高速摄影冷态喷雾试验研究了主要的结构参数对气液针栓式喷注器雾化角的影响，结果表明液膜厚度和越程对雾化角影响显著^[21-22]。Santoro和Merkle使用LOX/RP-1研究表明液液针栓式喷注器的雾化角与动量比有重要关系，而与阻塞率几乎无关^[23]。Sakaki等也发现动量比是影响液液针栓式喷注器雾化角的关键主导性无量纲参数，但没有得出它们之间的关系式^[24]。Son等通过高速摄影试验建立了气液针栓式喷注器雾化角与动量比及韦伯数定量关系^[25]，并对简化的二维轴对称针栓式喷注器模型进行了数值仿真，与之前的试验形成验证补充^[26]，分析了雾化角随动量比的变化规律^[27-28]。Radhakrishnan等基于Eulerian-Lagrangian方法开展了喷注条件对混合特性的影响，关注了二维轴对称针栓式喷注器计算模型的雾化角与动量比之间的关系^[29]。Cheng等基于无黏不可压的动量方程及连续方程，忽略表面张力，推导出了液膜撞击液膜的雾化角理论公式，即θ=arccos[1/(1+C_TMR)]，并与仿真及高速摄影试验结果吻合较好^[30]。

针栓式喷注器结构形式独特，雾场范围大，过于稠密，光学观测困难，现有光学测量设备较难获得有效喷雾场参数，为了解决此问题，Sakaki等率先提出平面针栓式喷注器单元的设想，选取最小喷注单元为研究对象，研究了动量比和喷雾场结构的关系，结果表明平面针栓与轴对称针栓的特性基本一致^[31-32]。

以前研究者的结果均是液膜/液膜撞击，属于二维平面撞击，而实际的液液针栓式喷注器一般都是液膜/液束撞击，属于三维空间撞击，撞击过程中液膜液束均会发生显著变形，流场结构更加复杂，此时雾化角与动量比的关系存在较大差异，阻塞率及径向孔的形状也成为影响雾化角的重要参数。目前仍未见报道关于液膜撞击液束的基础理论分析和定量研究，其有待于进一步研究。因此，本文以平面针栓喷注单元为研究对象，首先通过理论推导引入了变形因子，将撞击的几何变形效应与雾化角关联，建立液膜/液膜撞击雾化角模型；接着通过引入阻塞率定义有效撞击动量比，同时将液膜/液束撞击变形的影响隐含考虑在变形因子中，建立液膜/液束撞击雾化角模型；并与常用雾化角模型进行对比分析，揭示雾化角模型之间的本质区别；然后通过数值仿真和试验结果对雾化角模型进行验证；最后针对圆孔液束与液膜相互作用，建立了对应的变形因子组合系数，获得适应性更广、准确性更高的雾化角预测模型。

1 雾化角理论模型分析 1.1 液膜/液膜撞击雾化角模型 1.1.1 基本定义

动量比的定义为：径向射流动量与轴向射流动量之比，即$ {C_{{\rm{TMR}}}} = {\dot m_2}{u_2}/{\dot m_1}{u_1} $，其中$ {\dot m_1} $和${\dot m_2} $分别为轴向和径向两路流量，u₁和u₂分别为轴向和径向速度，动量比反映了喷注器雾化混合的状态。

下面通过动量比的定义分析Boettcher等^[20]的研究结果中二维平面膜与三维环形膜雾化角差异的原因。两种结构的膜撞击示意图如图 2所示。

对于平面膜，动量比为

$ {C_{{\rm{TMR}}}} = \frac{{{{\dot m}_2}{u_2}}}{{{{\dot m}_1}{u_1}}} = \frac{{{\rho _2}u_2^2{A_{{\rm{p,2}}}}}}{{{\rho _1}u_1^2{A_{{\rm{p,1}}}}}} = \frac{{{\rho _2}u_2^2{\kern 1pt} {h_2}}}{{{\rho _1}u_1^2{\kern 1pt} {h_1}}} $	（1）

对于环形膜，动量比为

$ \begin{array}{*{20}{l}} {C_{{\rm{TMR}}}^* = \frac{{{{\dot m}_2}{\kern 1pt} {u_2}}}{{{{\dot m}_1}{\kern 1pt} {u_1}}} = \frac{{{\rho _2}{\kern 1pt} u_2^2{A_2}}}{{{\rho _1}{\kern 1pt} u_1^2{A_1}}} = \frac{{{\rho _2}u_2^2{\kern 1pt} {h_2}D}}{{{\rho _1}u_1^2{\kern 1pt} {h_1}({h_1} + D)}} = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{\rho _2}u_2^2{\kern 1pt} {h_2}}}{{{\rho _1}u_1^2{\kern 1pt} {h_1}}} \cdot \frac{1}{{1 + {h_1}/D}} < \frac{{{\rho _2}u_2^2{\kern 1pt} {h_2}}}{{{\rho _1}u_1^2{\kern 1pt} {h_1}}} = {C_{{\rm{TMR}}}}} \end{array} $

（2）

式(1)~式(2)中：A_{p, 1}和A_{p, 2}分别为二维平面膜的轴向路和径向路流通面积；ρ₁和ρ₂分别为轴向和径向流体密度；A₁和A₂分别为轴向和径向流通横截面积；h₁和h₂分别为轴向和径向液膜厚度；D为针栓柱直径。

从式(1)和式(2)可以看出，二维平面膜的雾化角θ=arctan C_TMR略大于三维环形膜的θ^*=arctan C_TMR^*。将三维环形膜简化为二维平面膜计算雾化角时，预测的雾化角略微偏大。当环形膜厚度远小于针栓柱直径(即h₁$ \ll $D)时，二者完全相等，因此可将三维环形膜撞击简化为二维平面膜撞击。

1.1.2 基本假设

根据撞击物理过程建立如图 3所示的控制体几何模型，做出以下假设：

1) 两种液体不相溶。

2) 流动过程是稳态，即流动与时间无关。

3) 入口和出口截面的速度方向垂直于该截面。

4) 假定控制体1的入口和2个控制体的出口均远离撞击点，因此这些面的压力均可视为环境压力。

5) 忽略体积力(如重力等)和表面张力。

6) 忽略壁面对液膜的摩擦阻力和环境气体对气液自由表面的剪切阻力。

7) 不考虑相间传热和相变(如蒸发)。

8) 轴向动量守恒，径向动量不守恒^[33]。

1.1.3 雾化角公式推导

由忽略非稳态项和体积力的积分型动量方程可得

$ \oint_A \mathit{\boldsymbol{V}} (\rho \mathit{\boldsymbol{V}} \cdot {\rm{d}}\mathit{\boldsymbol{A}}) + \oint_A p {\rm{d}}\mathit{\boldsymbol{A}} = {\mathit{\boldsymbol{F}}_{\rm{s}}} $	（3）

式中：V为速度矢量；A为面积；A为控制体表面矢量；ρ为流体密度；p为压力；F_s为相界面间剪切阻力。

根据假设，将式(3)应用于控制体1和控制体2，得到相应的动量方程在轴向(z方向)的投影方程为

$ \begin{array}{l} - \int_{{A_{{\rm{1in}}}}} {{\rho _1}} V_{{\rm{1in}}}^2{\rm{d}}A + {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta \cdot \int_{{A_{{\rm{1out}}}}} {{\rho _1}} V_{{\rm{1out}}}^2{\rm{d}}A + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta \cdot \int_{{A_{{\rm{ 2out }}}}} {{\rho _2}} V_{{\rm{2out }}}^2{\rm{d}}A = 0 \end{array} $

（4）

式中：V为速度大小；下标1和2分别表示控制体1和2；下标in和out分别表示入口和出口截面。

从式(4)可见，2种液相界面间的压力和剪切阻力属于内力(作用力与反作用力)，求和为零，故式(4)未出现此项。式(4)主要根据重要假设(8)获得，假设(8)在文献[33]中有详细的论证过程，并已得到了试验及数值仿真结果的充分验证。

从式(4)可获得雾化角为

$ {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta = \frac{{\int_{{A_{{\rm{1in }}}}} {{\rho _1}} V_{{\rm{1in }}}^2{\rm{d}}A}}{{\int_{{A_{{\rm{ 1out }}}}} {{\rho _1}} V_{{\rm{ 1out }}}^2{\rm{d}}A + \int_{{A_{{\rm{ 2out }}}}} {{\rho _2}} V_{{\rm{ 2out }}}^2{\rm{d}}A}} $	（5）

式(5)与Cheng等^[30]得到的方程相同，然而推导过程不同，接下来的推导过程则导致最终结果与其不同。

由于各个截面的参数具有分布特性，引入截面平均速度u和平均密度$ \bar \rho $(这里忽略平均速度和均方根速度的差异)，定义截面的平均动量为$ \bar \rho {u^2} = \frac{{\int_A {\rho {V^2}{\rm{d}}A} }}{A} $，则控制体1和2入口及出口截面的平均动量分别为为ρ₁u_1in²，ρ₁u_1out²，ρ₂u_2in²和ρ₂u_2out²。

将上述定义的参数代入式(5)可得

$ {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta = \frac{{{{\bar \rho }_1}u_{{\rm{1in}}}^2{A_{{\rm{1in}}}}}}{{{{\bar \rho }_1}u_{{\rm{1out }}}^2{A_{{\rm{1out}}}} + {{\bar \rho }_2}u_{{\rm{2out}}}^2{A_{{\rm{2out}}}}}} $	（6）

而$ {\dot m_1} = {\bar \rho _1}{u_{1{\rm{in}}}}{A_{{\rm{1in}}}} = {\bar \rho _1}{u_{{\rm{1out}}}}{A_{{\rm{1out}}}}, {\dot m_2} = {\bar \rho _2}{u_{{\rm{2out}}}}{A_{{\rm{2out}}}}$，同时

$ {C_{{\rm{TMR}}}} = \frac{{{{\dot m}_2}{u_{{\rm{2in}}}}}}{{{{\dot m}_1}{u_{{\rm{1in}}}}}} = \frac{{{{\dot m}_2}^2/({{\bar \rho }_2}{A_{{\rm{2in}}}})}}{{{{\dot m}_1}^2/({{\bar \rho }_1}{A_{{\rm{1in}}}})}} $，代入式(6)可得

$ \begin{array}{l} {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta = \frac{1}{{\frac{{{A_{{\rm{1in}}}}}}{{{A_{{\rm{1out}}}}}} + \frac{{{A_{{\rm{2in}}}}}}{{{A_{{\rm{2out}}}}}} \cdot \frac{{\dot m_2^2/\left( {{\bar \rho _2}{A_{{\rm{2in}}}}} \right)}}{{\dot m_1^2/\left( {{\bar \rho _1}{A_{{\rm{1in}}}}} \right)}}}} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{{\frac{{{A_{{\rm{1in}}}}}}{{{A_{{\rm{1out }}}}}} + \frac{{{A_{{\rm{2in}}}}}}{{{A_{{\rm{2out}}}}}} \cdot {C_{{\rm{TMR}}}}}} \end{array} $

（7）

定义$ {C_1} = \frac{{{A_{{\rm{1in}}}}}}{{{A_{{\rm{1out}}}}}} $和$ {C_2} = \frac{{{A_{{\rm{2in}}}}}}{{{A_{{\rm{2out}}}}}} $，则式(7)变为

$ {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta = \frac{1}{{{C_1} + {C_2} \cdot {C_{{\rm{TMR}}}}}} $	（8）

式(8)中C₁和C₂分别是2个控制体入口截面的面积与撞击后出口截面的面积之比，物理意义是表征由于撞击造成的截面面积变化。通过引入变形因子C₁和C₂将撞击的几何变形效应与雾化角关联，同时也可将撞击液膜的几何形状因素考虑其中，从而可以更加准确地描述雾化角与影响因素之间的关系。

1.2 液膜/液束撞击雾化角模型

液膜撞击液束属于三维空间撞击，流场结构更加复杂，此时雾化角不仅仅与动量比有关，阻塞率和径向孔的形状也是相当重要的参数。阻塞率是针栓头端部全部径向喷注孔的宽度之和与针栓周长之比。对于单个液膜/液束撞击单元，以圆柱径向孔为例，阻塞率为$ {C_{{\rm{BF}}}} = \frac{{nd}}{{{\rm{ \mathsf{ π} }}D}} = w/L $，它会影响与径向液束射流发生有效撞击的液膜动量，其中n为径向孔数目，d为径向孔直径，w为仅与径向液束发生有效撞击的液膜宽度；L为针栓喷注单元宽度，即相邻喷注单元间距。因而可引入阻塞率，定义有效撞击动量比(有效动量比)

$ \begin{array}{l} {C_{{\rm{MReff }}}} = \frac{{{{\bar \rho }_2}u_{{\rm{2 in }}}^2{A_{{\rm{2 in }}}}}}{{{{\bar \rho }_1}u_{{\rm{1in }}}^2{A_{{\rm{1in }}}}}} = \frac{{{{\bar \rho }_2}u_{{\rm{2in }}}^2{A_{{\rm{2in }}}}}}{{{{\bar \rho }_1}u_{{\rm{1in }}}^2{A_{{\rm{1in}}{\kern 1pt} {\kern 1pt} {\rm{total }}}} \cdot \frac{w}{L}}} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{C_{{\rm{TMR}}}}}}{{{C_{{\rm{BF}}}}}} \end{array} $

（9）

式中：A_{1in total}为轴向环缝液膜流动截面的总面积。式(9)中隐含假设认为w=d。

以单个液膜/液束撞击单元作为控制体研究对象，如图 4所示，基于液膜撞击液膜的雾化角分析模型，将式(9)代入式(8)中可推导出液膜撞击液束的雾化角模型为

$ \begin{array}{*{20}{c}} {{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta = \frac{1}{{{C_1} + {C_2}{C_{{\rm{ MReff }}}}}} = }\\ {\frac{1}{{{C_1} + {C_2}\frac{{{C_{{\rm{ TMR }}}}}}{{{C_{{\rm{ BF }}}}}}}}} \end{array} $	（10）

液膜/液束撞击过程中液束会同时发生弯曲变形和横截面变形，液膜也会发生绕液束流动，因此液膜/液束撞击后截面面积会发生较大变化，通过在式(10)中引入变形因子C₁和C₂，隐含考虑了液膜/液束撞击的显著影响，进而修正了雾化角模型。与液膜/液膜撞击相比，液膜/液束撞击的一个显著不同是液束的入口孔型有多种选择，可以是单一孔型(如圆孔、椭圆孔、方孔、长宽比不同的矩形槽等)，也可以是以上各种形状的组合孔型，这些不同形状的液束与液膜发生撞击后，雾化角变化规律类似，可通过相应的变形因子C₁和C₂的组合系数来表征，从而描述不同孔型液束与液膜撞击的相互作用，进而获得适应不同情况更加准确的液膜撞击液束的雾化角模型。本文中主要以简单的圆孔为例开展研究。

1.3 模型分析讨论

关于雾化角公式(8)的讨论：

1) 由于推导过程中以守恒量ρV²整体进行积分变换，后面又引入了平均密度，因此式(8)既适用于不可压的液体膜及气体膜撞击，也适用于可压缩的气体膜撞击。

2) 当认为两路液体撞击后出口截面面积与各自入口截面面积相等时，式(8)中2个变形因子均等于1，因此雾化角公式简化为

$ {\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta = \frac{1}{{1 + {C_{{\rm{TMR}}}}}} $	（11）

式(11)与Cheng等^[30]得到的结果一致，从而也可以看出Cheng等^[30]的结果仅适用于2个变形因子均为1的特殊情况，这也正是引入变形因子的新模型与Cheng模型的不同之处。对于二维平面膜撞击，各自的截面积变化相对小，变形因子均为1的假设近似成立。对于液膜/液束撞击，液膜绕流变形和液束变形效应更显著，变形因子更容易偏离1。具体的变形因子数值可根据不同的撞击形式及几何形状，通过试验和数值仿真获取，进而可获得适应性更广、准确性更高的雾化角预估模型。

3) 从式(5)或式(6)不难看出，该模型的雾化角余弦值等于轴向入口动量(即控制体1的入口截面动量)与撞击后出口总动量(即控制体1和2的出口截面的动量和)之比，是根据实际的出口轴向动量(等于入口轴向动量)和出口合成总动量来计算雾化角，隐含考虑了撞击造成的影响，也适合引入撞击液膜与液束的几何形状因素。该模型仅认为轴向动量守恒，径向动量不守恒。而相比于另一种常用的雾化角模型：tan θ=C_TMR，即

$ \begin{array}{*{20}{l}} {{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \theta = \frac{1}{{\sqrt {1 + C_{{\rm{TMR}}}^2} }} = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{{\bar \rho }_1}u_{{\rm{1in}}}^2{A_{{\rm{1in}}}}}}{{\sqrt {{{({{\bar \rho }_1}u_{{\rm{1in}}}^2{A_{{\rm{1in}}}})}^2} + {{({{\bar \rho }_2}u_{{\rm{2in}}}^2{A_{{\rm{2in}}}})}^2}} }}} \end{array} $

（12）

该常用模型的雾化角余弦值等于轴向入口动量与撞击前入口合成总动量(即轴向入口动量与径向入口动量平方和的算术平方根)之比，是根据撞击前入口的轴向动量和合成总动量来计算雾化角，未考虑撞击造成的影响，成立的前提条件是动量完全守恒(轴向动量守恒和径向动量守恒)，这也是式(12)和式(8)2个雾化角模型的本质区别。

由不等式1+x² < (1+x)²(x>0)易知，$ {\rm{cos}}\;{\theta ^*} = \frac{1}{{\sqrt {1 + C_{{\rm{TMR}}}^2} }} > \frac{1}{{1 + {C_{{\rm{TMR}}}}}} = {\rm{cos}}\;\theta $，即式(12)预测的雾化角总是小于式(6)或式(8)的预测值。而式(12)预测的雾化角总是小于试验测量的雾化角，这一结论在前人的研究文献中已得到证实，式(12)雾化角模型的计算值与试验测量值相比偏小，最大偏差20°以上^[33]。

2 数值物理模型 2.1 数值方法 2.1.1 控制方程

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

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

（13）

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

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

$ \left\{ {\begin{array}{*{20}{l}} {\rho ({c_1},{c_2}) = {c_1}{\rho _{{\rm{l1}}}} + {c_2}{\rho _{{\rm{l2}}}} + (1 - {c_1} - {c_2}){\rho _{\rm{g}}}}\\ {\mu ({c_1},{c_2}) = {c_1}{\mu _{{\rm{l1}}}} + {c_2}{\mu _{{\rm{l2}}}} + (1 - {c_1} - {c_2}){\mu _{\rm{g}}}} \end{array}} \right. $	（14）

式中：ρ_l1、ρ_l2和μ_l1、μ_l2分别是液相1和液相2的密度和黏度；ρ_g和μ_g分别为气相的密度和黏度。

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

2.1.2 气液相界面捕捉

VOF方法通过定义体积分数函数c来描述界面(含有运动界面的网格满足0 < c < 1，充满目标流体的网格c=1，不含目标流体的网格c=0)，通过求解函数c的输运方程跟踪界面运动^[40]。优点是可以方便地计算复杂的相界面变化过程，捕捉的相界面锐利程度高，对计算内存的要求较低，具有守恒的特性。

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

2.2 计算模型

计算所选用的结构参数及计算域如图 7所示，轴向液膜厚度为h₁，轴向速度为u₁，径向液膜厚度(或液束直径)为h₂(或d)，径向速度为u₂。图 7(a)为二维轴对称液膜撞击液膜的模型，计算域由8个L×L的基本结构Box构成，左端面和底端面为轴向液膜和径向液膜入口；图 7(b)为液膜撞击液束的模型，计算域由8个L×L×L的基本结构Box构成，计算域左端面和底端面分别为轴向液膜和径向液束入口，液膜展向宽度取10 mm。图中标注的面为无滑移壁面，其余面为出口，采用Outflow边界，背压为大气环境。第1相为空气，第2相和第3相均为水。计算域的L=10 mm，最高网格等级采用Level=9加密，最小网格约19.5 μm。网格自适应函数设置为体积分数梯度，即网格会实时根据流场中计算的体积分数梯度大小进行自适应加密或粗化。设定轴向液膜与径向液膜或液束撞击前的距离(即跳跃距离)为6 mm，对应模拟的针栓柱直径D=20 mm，其他具体的工况参数详见表 1~表 4。数值仿真获得雾化角的结果如图 8所示。

3 喷注单元试验装置及试验测量系统 3.1 喷注单元试验装置

为了深入研究不同几何结构参数、流动工况参数在大范围变化情况下对雾化角的影响规律，同时便于开展有效光学观测及将试验结果与数值仿真结果进行对比，冷态试验的喷注器设计借用文献[31-32]中平面针栓式喷注单元的设计理念，采用喷注器单元结构可更换的方案设计了平面针栓式喷注单元试验装置。平面针栓式喷注单元由一个可更换的液膜生成构件和一个可更换的液束生成构件组成，试验件如图 9所示。液膜生成构件可以产生一定厚度的平面液膜，平面液膜的宽度足够宽，设计中单孔喷注单元选取液膜宽度为10 mm；液束生成构件可以产生一定孔型型面和一定尺寸的射流。试验中液膜和液束两路单独供应，通过分别改变液膜和液束两路流量及更换不同结构尺寸的液膜和液束构件，从而形成不同动量比的工况。

为了充分研究无量纲参数动量比的影响，方案中工况设计选取的动量比从0.1到6.25，覆盖较宽的变化范围；同时动量比的变化包含了液膜、液束入口速度(流动参数)和液膜厚度、液束直径(结构参数)等的变化，可使研究结果更具一般性地说明动量比这一无量纲参数的影响规律。试验中选取的可更换液膜厚度h₁和液束直径d的结构参数如表 5所示。

3.2 试验测量系统

试验系统简图如图 10所示，试验中通过高压空气挤压贮箱实现试验用水供应，供应管路上设置科氏力质量流量计(型号：F050)，测量水的质量流量。试验中，喷前压力通过安装在靠近喷前的压阻型压力传感器测得，采样频率1 000 Hz^[41-42]。使用LED面光源照射喷注单元及其液雾场，通过Phantom V12.1型COMS黑白高速相机及其镜头拍摄喷注单元的撞击雾化过程图像，该相机最大分辨率为1 280 pixel×800 pixel，在最高1×10⁶ frame/s的帧频下能达到的最大分辨率为128 pixel×8 pixel，最小曝光时间1 μs。试验中观测方法如图 10所示，高速相机采样频率为3 000 Hz，曝光时间为10 μs，图像分辨率为1 024 pixel×768 pixel。

3.3 试验图像处理方法

试验中喷雾图像采用高速阴影法拍摄，然而雾化过程是瞬态过程，产生的雾扇会随时间波动，为了准确评估雾化角大小，对试验获得的如图 11(a)所示的喷雾瞬态图像进行图像增强处理，获得如图 11(b)所示的增强图像，然后将1 000张增强后的瞬态图像进行平均处理获得如图 11(c)所示的图像，再取雾扇边界平均中心来测量获得雾化角，如图 12所示。

4 结果讨论与分析 4.1 液膜/液膜撞击雾化角结果分析

针对1.1.3节建立的液膜/液膜撞击雾化角理论模型中引入的变形因子C₁和C₂，通过分相识别方法对2路液膜的撞击雾化过程进行数值仿真，分别提取轴向液膜和径向液膜撞击前后的截面面积，如图 13^[33]所示，图中M_1in、M_2in、M_1out和M_2out分别为轴向路和径向路撞击前后截面的动量，p_w为壁面压力，p_2in为径向路撞击前截面的压力，F_w为壁面摩擦阻力。出口截面的统计位置选在撞击点下游径向距离4 mm处。通过获得的2路液膜入口截面面积A_1in和A_2in及出口截面面积A_1out和A_2out计算得到的面积比C₁和C₂如图 14所示。可以看出随动量比的变化，面积比基本维持在0.9~1.1，特别是动量比小于2.5时，这也是实际应用中动量比的范围。因此，在计算雾化角时C₁和C₂在0.9~1.1范围内取值(推荐值如图 15所示，其是根据雾化角数据拟合获得，图中括号里还给出了C₁和C₂拟合获得的上限值和下限值)，采用推导的雾化角模型可准确预估雾化角。

根据上述获得的结果，将数值仿真计算的雾化角、根据面积比(变形因子)计算的雾化角、试验测得的雾化角、以前常用的tan模型预测的雾化角、Cheng的模型(即本文模型对应的C₁=1、C₂=1特例)预测的雾化角及根据图 15拟合得到的变形因子模型预测的雾化角进行对比，得到如图 16所示的雾化角随动量比变化曲线，可以看到雾化角随着动量比增大而增大；在动量比较小时，雾化角增长较快；在动量比较大时，雾化角增长开始变慢；在动量比很大时，雾化角增大很慢。从图 16还可以看出在动量比大范围变工况下，数值仿真结果、试验结果及变形因子模型预测结果均吻合很好，与以前常用的tan模型预测结果差异较大，最大相差20°以上，原因已在1.3节中讨论分析。另外，发现Cheng的模型也与计算结果及试验结果吻合较好，这是由于液膜撞击液膜属于正撞，仅有正向应力压缩作用，撞击过程液膜变形较小，撞击前后横截面变化很小，变形因子接近1，图 14和图 15的结果可以印证这点，因而该模型也成立。

4.2 液膜/液束撞击雾化角结果分析

相比于液膜撞击液膜，液膜撞击液束属于三维空间撞击，同时存在正向应力压缩和切向应力剪切的双重作用，撞击过程中液膜与液束均会发生较大变形，液束会同时发生弯曲变形和横截面变形，液膜也会发生绕液束的分叉流动，因此液膜/液束撞击的雾化角变化规律与液膜/液膜撞击的存在差异。根据1.2节建立的液膜/液束撞击雾化角模型，引入阻塞率定义有效动量比，获得了雾化角随有效动量比变化的曲线，如图 17所示。同样将数值仿真结果、试验结果、以前常用的tan模型预测结果及变形因子均为1(C₁=1、C₂=1)的模型预测结果进行对比可以发现，数值仿真结果与试验结果在有效动量比大范围变工况下均吻合很好，与变形因子均为1的模型预测结果趋势一致，与以前常用的tan模型预测结果差异较大，这类似于4.1节液膜与液膜撞击雾化角的结果。然而进一步对比图 17(a)可以发现，在中等有效动量比(1.2~2.0)下，变形因子均为1的模型(即Cheng模型)预测结果与仿真结果及试验结果吻合较好；在低的有效动量比(< 1.2)下，该模型明显高估了雾化角；在高的有效动量比(>2.0)下，该模型明显低估了雾化角。可见，对于膜撞束，由于变形显著增大，未考虑膜束变形影响的Cheng模型预测值误差增大，此时该模型适应的局限性凸显。

对于新模型，根据试验结果和数值仿真结果拟合变形因子C₁和C₂得到如图 17(b)所示的结果，在低的有效动量比时，得到的变形因子为C₁=0.95和C₂=1.05；在高的有效动量比时，得到的变形因子为C₁=0.55和C₂=1.45；综合高低有效动量比范围内的结果，得到的变形因子为C₁=0.75和C₂=1.25，可以看到在动量比大范围变工况下，拟合的变形因子模型预测结果与试验及数值仿真结果均吻合较好。以上结果是在动量比大范围变工况下获得的，包含了液膜、液束入口速度(流动参数)和液膜厚度、液束直径(结构参数)等的变化，因此可使研究结果更具一般性地说明动量比这一无量纲参数的影响规律。

对于液束为圆孔射流的情况，根据式(9)有效动量比为${C_{{\rm{MReff}}}} = \frac{{{\rho _2}u_{{\rm{2in}}}^2{\rm{ \mathsf{ π} }}d}}{{{\rho _1}u_{{\rm{1in}}}^24{h_1}}}\cdot\frac{d}{w} $。图 18所示为液膜与圆孔射流撞击发生绕流流动与变形过程的示意，由于液膜会发生绕圆射流的分叉流动，实际与圆射流发生有效撞击的液膜宽度w小于圆射流的迎风面宽度，因此，实际的有效动量比大于前文假设w=d计算的结果，这使得实际的雾化角大于变形因子为1的模型预测结果。通过变形因子C₂>1的修正可得到更准确的预测雾化角，而且有效动量比越大，圆射流的惯性越大，射流根部横截面越容易保持圆形不变，这使得发生绕流的液膜更多，即w比d小的更多，C₂更大；在中等/较小有效动量比下，根部横截面变形较大使得发生有效撞击的液膜宽度增大，此时w≈d近似成立，C₂接近于1，这也与图 17(b)中结果一致。同时，由于液膜与液束的相互作用，出口高压区使得液束进口压力p_2in大于出口压力p_2out(即环境压力p_a)，由伯努利方程可得u_2in < u_2out，对于不可压流可得到A_2in>A_2out，即C₂>1。对于C₁，液膜与液束撞击后液膜径向拉伸变形，形成包裹液束的半开放喇叭口形(“Ω”形)雾扇流动，使得液膜出口横截面积变大，如图 19所示，并产生了展向的流动速度，轴向速度便会减小，这都使得A_1in < A_1out，即C₁ < 1；有效动量比越大，雾化角越大，“Ω”形雾扇越大，C₁越小；有效动量比较小时，液束惯性较小，液束前缘更容易变形，且呈扁平状更明显，更接近液膜撞击液膜，因而C₁也更接近1。

至此，针对建立的雾化角模型，分析和验证了引入变形因子的合理性和准确性，也给出了定量预估雾化角建议的变形因子参考值，可为针栓式喷注器的理论研究和工程设计提供重要参考。

5 结论

为了研究动量比对针栓式喷注器撞击雾化角的影响规律，从动量守恒方程出发推导建立了液膜撞击液膜和液膜撞击液束雾化角的新理论模型，同时基于Gerris开源软件采用的自适应网格技术和PLIC VOF方法对液膜/液束撞击的雾化角进行了数值模拟，并结合试验结果对理论模型进行了验证。本文研究表明：

1) 基于轴向动量守恒和径向动量不守恒的假设，通过理论推导引入了2个变形因子，将撞击的几何变形效应与雾化角关联，推导建立了液膜撞击液膜的雾化角新模型。引入变形因子的理论模型预测值与试验及数值仿真结果吻合很好；变形因子基本维持在0.9~1.1，根据试验结果及仿真结果的推荐值为C₁=0.99和C₂=1.06。

2) 基于液膜撞击液膜的雾化角分析模型，通过引入阻塞率定义有效撞击动量比，同时将液膜液束变形的影响隐含考虑在变形因子中，推导建立了液膜撞击液束的雾化角新模型。引入变形因子和阻塞率的理论模型预测值与试验及数值仿真结果吻合很好，理论模型和数值仿真模型均得到了有效验证；根据试验结果及仿真结果获得的变形因子推荐值为C₁=0.75和C₂=1.25。

3) 建立的雾化角新模型根据实际出口的轴向动量(等于入口轴向动量)和合成总动量计算雾化角，隐含考虑了撞击变形作用的影响，较根据撞击前入口的轴向动量和合成总动量计算雾化角的常用tan模型预测值准确度显著提高，这是2种雾化角模型之间的本质区别，为针栓式喷注器的理论研究和工程设计提供了重要参考。