小推力液体火箭发动机沉降型液膜冷却液膜长度研究

兖立文　许坤梅　王慧洁

（北京航空航天大学宇航学院，北京　100191）

小推力液体火箭发动机沉降型液膜冷却液膜长度研究

兖立文许坤梅王慧洁

（北京航空航天大学宇航学院，北京100191）

沉降型液膜冷却是指从喷嘴向壁面以一定角度喷注某种液滴，液滴到达壁面沉降形成液膜，实现对壁面热防护的一种方法，常用于小推力液体火箭发动机热防护系统中。采用Stechaman半经验方法对决定液膜冷却效果的关键参数——液膜长度求解。采用k-w模型描述湍流流动、Eulerian-Lagrangian模型描述两相流，采用C/C++语言编写Bai-Gosman液滴撞壁模型程序，采用数值模拟推进剂液滴撞壁后复杂的状态变化过程，考虑壁面热辐射，对400N双组元MMH/NTO自然推进剂发动机推力室内的蒸发、流动、燃烧和传热过程进行了数值模拟。在采用半经验公式方法求解液膜长度的过程中，分析了液膜流量对液膜长度的影响，研究了该方法的实际应用情况。实验结果表明，该方法可以较好地计算液膜长度，计算结果与实验具有较高的一致性，对工程实践具有重要的指导意义。

液体火箭发动机，液膜冷却，液膜长度，液滴撞壁，数值模拟

引言

对于小推力液体火箭发动机来说，因为其推力和尺寸较小，所以，液态燃料喷入推力室后，不可避免地会出现液滴撞壁过程。液滴撞壁后，其运动过程非常复杂，在某些条件下会在壁面形成液膜。在对发动机推力室内喷雾湍流燃烧进行数值模拟的基础上，采用Bai-Gosman［1］液滴撞壁模型模拟了液滴撞壁后的复杂状态变化过程，并分析了液滴撞壁后的轨迹、推力室流场和发动机各项性能参数。

液体火箭发动机燃烧室温度高、热流密度大，可采用冷却方法保护发动机，目前常用的冷却方法有再生冷却、辐射冷却、热沉法、隔热法、液膜冷却或发汗冷却法［2］。其中，液膜冷却是一种应用较为广泛且有效的冷却方式。各国研究人员对液膜冷却开展了多年的深入研究。例如，文献［3～4］研究了推力室中心气流与液膜界面的结构，以及冷却液膜的不稳定性；文献［5～6］从实验和理论角度对液膜进行了研究，但未对液膜长度进行定量分析，且不针对小推力发动机。本文的研究对象为小推力发动机内沉降型液膜冷却，且需进行定量分析。分析发现，Stechaman半经验方法［7］针对小推力发动机液膜冷却进行了研究，提出了适用于小推力发动机液膜长度计算的经验公式，因此，可采用该方法分析液膜长度，计算中需确定液膜流率，但不同于喷入型液膜冷却，在沉降型液膜冷却中，液膜流率是未知的。本文针对这一点，采用Bai-Gosman液滴撞壁模型，通过自编程求解液滴撞壁后形成液膜时的流率，采用Stechaman半经验方法对400N小推力发动机［8］的液膜冷却过程进行数值模拟，分析液膜流量因素对液膜长度的影响。

1　物理数学模型

1.1气相湍流反应流的N-S方程

控制方程中的各种定律反映的都是单位时间单位体积内物理量的守恒性质，可以表示成以下通用形式［9］：

其中，Φ表示通用变量，可以代表u、v、w、T等求解变量；Г为广义扩散系数；S表示广义源项。式（1）中各项依次为瞬态项、对流项、扩散项和源项。

1.2液滴运动控制方程

液滴在空间的分布及其运动状态在很大程度上取决于其运动过程中所受到的力。在忽略液滴的旋转及流场中速度梯度产生的升力、Magnu力，以及重力等作用的情况下，液滴运动方程可以表示为：

通过对时间的积分，可得到液滴在各个时刻的速度和位置。CD是液滴阻力系数，rp是液滴的半径，和分别是介质气体和液滴速度矢量。阻力系数CD由下式求得：

其中，Re为液滴的雷诺数。

1.3Bai-Gosman液滴撞壁模型

Bai和Gosman［2］提出的模型将液滴撞壁过程详细划分成粘附、反弹、铺展、沸腾后破碎、反弹伴随破碎、破碎、飞溅等7种形式。在以临界韦伯数判定的湿壁情况下，不同状态的转变临界值情况为：反弹→铺展→飞溅。

式中，La-Laplace数据由壁面粗糙度决定，We为液滴的韦伯数，σ为液滴的弹性能，d为液滴的直径，v为液滴的速度，ρ为液滴的密度，μ为液滴的粘度。本文仅讨论液滴撞壁面后反弹和飞溅情况。

反弹后液滴的切向和法向速度如下：

其中，u和v分别是切向与法向分速度，下标e和r分别代表入射和反弹，ε为恢复系数。ε=0.993-1.76θ+1.56θ2-0.49θ3，θ是以弧度计算的液滴入射角。

飞溅时，Bai-Gosman模型假定液滴撞壁后、飞溅后液滴团的直径和速度分别为d1、d2和U1、U2。如果两个二次液滴团所含的液滴数分别为N1和N2，其质量守恒方程为：

借用气象学中的相关数据拟合确定总的二次液滴数

二次液滴的速度方程为：

参考雨滴尺寸与速度的关系式推导出液滴的速度比为：

上式表明，若破碎后的液滴尺寸与速度呈负相关关系，尺寸越大，速度就越小。

液滴的动量方程为：

式中，Cf是摩擦系数，取值为0.6～0.8，θ1和θ2是相应的角度。

1.4Stechaman半经验方法

在考虑推力室内热气流与液膜的热传递时，采用Stechaman［7］修正的Bartz［10］方法来计算热气流与液膜的传热系数［11］，传热系数为式（12）。

Stechaman半经验公式方法计算液膜长度为式（13）。

式中，*为喉部位置，μf为动力黏度，D为发动机燃烧室直径，A为面积，W为冷却剂流量，HW为壁面条件下的焓值，TW为壁面条件下的温度，Hr为恢复焓值，Tr为恢复温度，CPL为液膜冷却剂的等压比热容，Prf为液膜的普朗特数，η为与液膜雷诺数相关的无量纲数，σ为边界层处与密度和黏度有关的无量纲数，L为液膜长度，Tl为液滴初始喷注温度，Ts为饱和温度，P为发动机燃烧室周长，hg为热气流与液膜的传递系数，λ为冷却剂的潜热。

采用Stechaman半经验公式法计算液膜长度时，需确定形成液膜的流量，在喷入型液膜冷却中，可根据入口条件获得；而在沉降型液膜冷却中，由于燃料从喷嘴向壁面以一定角度喷注某种液滴，而液滴撞壁后状态复杂，不能直接确定其形成液膜的流量。因此，本文采用Bai-Gosman液滴撞壁模型，通过自编程求解液滴撞壁后形成液膜时的流率。

2　数值模拟

2.1计算参数

发动机的喷注器采用双组元离心式，向其内喷嘴注入燃料甲肼（MMH），外喷嘴喷入氧化剂硝基-三唑-酮（NTO），氧燃推进剂流率的混合比为1.667，全部流量为130g/s，初始温度为300K。液滴的其余参数由前面计算得出。MMH的喷注速度Vinj，MMH为16.2m/s，NTO的喷注速度Vinj，NTO为2.6m/s，与对称轴的夹角为55°。液滴进入燃烧室时，MMH的入口位置为：X=3.5mm，Y=5mm；NTO的入口位置为：X=3.5mm，Y=8mm。假定喷雾液滴的最小粒径dmin为5μm，最大粒径dmax为200μm或110μm；假设喷注器所产生的液滴尺寸分布服从Rosin-Rammler分布。按Rosin-Rammler分布将液滴分组，分组数目为50，则质量分布函数公式为：

n为液滴分布指数，表示液滴直径分布的均匀性，一般n=2～4。n越大，液滴直径分布越均匀，本文取n=4。其液滴参数详见表1。

表1　喷雾喷注参数

2.2计算网格

图1是某400N轨控发动机的计算网格，共包括24600个网格单元。根据湍流对壁面附近网格的要求，对壁面附近的网格进行合理加密。由于发动机燃烧室喷嘴入口附近流场的温度、浓度等参数梯度较大，所以，需对喷注面附近的网格加密。

图1　计算网格

2.3壁面边界条件

气相壁面为热辐射无滑移边界；液相壁面为Bai-Gosman液滴撞壁模型。

2.4数值解法

控制方程采用有限体积法进行离散，离散方程的求解采用SIMPLE算法。图2为计算方法流程：首先，对推力室内流场进行数值仿真，收敛时计算内壁面温度Tw， in［i-1］、获取Stechaman经验公式法所需参数、计算液膜长度，更新设置壁面条件；其次，对推力室内流场再次进行数值模拟至收敛，计算内壁面温度Tw， in［i］。重复以上两个步骤，直至满足以下公式：

式（16）中，Tw， in［i］是当前步第i次时推力室内壁面温度值，而Tw， in［i-1］是上一步第i-1次时推力室内壁面温度值，下标“max”为发动机推力室计算域中相应参数的最大值。

图2　计算方法流程图

3　计算结果与讨论

本文采用不同的液滴参数对400N轨控发动机的工作过程进行了数值模拟。

3.1Case 1数值模拟

3.1.1静压与静温

图3、图4显示了当前工况下的静压与静温，由图可以看出，压力在燃烧室内最大，从喷管收敛段处开始减小，在喷管扩张段继续减小直至喷管出口。在拉瓦尔喷管中，在收敛段压力开始降低直至喷管出口，燃烧室压强为1.1MPa，比实验所测（1.02MPa）略高，约为7.8%。发动机中心的温度较高，发动机头部的温度相对较低，这是由于头部有MMH液膜保护。

图3　压力分布

图4　温度分布

3.1.2MMH轨迹分布

图5为液滴撞壁时采用反弹模型，即撞壁液滴全部反弹，且不考虑液膜时的MMH液滴轨迹分布，液滴全部反弹且全部蒸发，壁面无液膜；图6为液滴撞壁时采用Bai-Gosman液膜撞壁模型且考虑液膜时，MMH液滴轨迹分布，MMH撞壁后一部分反弹，一部分形成液膜，但没有飞溅，这是由于撞壁液滴的韦伯数不满足Bai-Gosman液膜撞壁模型中液滴飞溅的条件。MMH撞壁形成的液膜可以对壁面起到热保护作用。

图5　反弹模型MMH轨迹

图6　Bai-Gosman模型MMH轨迹

3.1.3壁面温度分布

图7为采用Stechaman半经验法获得的壁面温度分布，在当前液滴参数下，计算的液膜长度小于实验数据；发动机头部温度较低，未出现明显的高温区，这是由于MMH撞壁后在壁面形成了液膜，对头部壁面起到了热防护作用，使得温度较低。

图8为不考虑液膜冷却时的壁面温度分布，发动机头部出现了局部高温。这是由于MMH/NTO射流以相应角度喷入推力室后，在壁面头部存在MMH/NTO混合较均匀的位置，MMH/NTO发生化学反应，而数值计算中未考虑液膜冷却作用，导致该处壁面温度较高。图7、图8对比分析表明，在发动机工作过程中，需考虑沉降型液膜冷却，形成的液膜可对发动机壁面起到热防护作用。

图7　壁面温度分布

图8　壁面温度分布

3.2增加MMH的质量中径（Case 2）

Case 2增加了MMH的质量中径。图9为采用Stechaman半经验法获得的壁面温度分布，计算的液膜长度大于Case1，这可能是由于MMH质量中径增大，MMH蒸发变慢、撞壁概率增大，导致液膜长度增大；计算液膜长度稍大于实验数据。图10为液滴撞壁时采用Bai-Gosman液膜撞壁模型且考虑液膜时的MMH液滴轨迹分布情况，图示MMH撞壁后一部分反弹，一部分形成液膜，这部分液膜可以对壁面起到热防护作用。

图9　壁面温度分布

图10　Bai-Gosman模型MMH轨迹

3.3增加NTO的质量中径（Case 3）

Case 3增加了NTO的质量中径。图11为采用Stechaman半经验法时壁面温度分布，在当前液滴直径下，其液膜长度比Case 2中低，这可能由于NTO质量中径增大，蒸发变慢，增加了NTO撞壁面的概率，增大NTO与壁面MMH液膜接触的概率，接触后，在附近发生化学反应，消耗了形成液膜的MMH的流量，导致液膜长度比Case 2低，但计算的液膜长度与实验中数据非常吻合。图12为液滴撞壁面时采用Bai-Gosman液膜撞壁模型且考虑液膜时，MMH液滴轨迹分布，MMH撞壁后一部分反弹，一部分形成液膜，这部分液膜可以对壁面起到热防护的作用。3.4液膜流率对液膜长度的影响

图11　壁面温度分布

图12　Bai-Gosman模型MMH轨迹

Case 1、Case 2相比，仅有液滴的平均质量中径不同，其余参数均相同，而Case 3中液滴的最大直径与Case 1和Case 2不同。采用Bai-Gosman液滴撞壁模型，采用自编程求解液滴撞壁后形成液膜时的流率并统计分析，分析结果如表2所示。

表2　液膜长度与液膜流率的分布

由表2可知，Case 1所得的液膜流率最小，Case 3所得的液膜流率其次，Case 2所得的液膜流率最大，Stechaman所得的液膜长度的顺序与液膜流率成正比，即Case 2的液膜长度最长，Case 3的次之，Case 1的最小。这表明，液膜长度与液膜流量成正比，即液膜流量越大，液膜长度越大。

4　结束语

本文对小推力液体火箭发动机燃烧室内的沉降型液膜冷却进行了数值模拟，得到以下主要结论：

（1）在Bai-Gosman液膜撞壁模型中，根据韦伯数将撞壁后液滴分为反弹、形成液膜和飞溅3种主要过程，可以较合理地数值模拟液滴撞壁后复杂状态的变化过程。

（2）在沉降型液膜冷却中，根据Bai-Gosman液膜撞壁模型，撞壁液滴只有在一定的韦伯数范围内，才形成液膜，而这部分流率是未知的，本文采用自编程求解液滴撞壁后形成液膜时的流率，并代入半经验公式方法中，液膜长度结果与实验数据较吻合。

（3）本文在采用半经验公式方法分析液膜时，形成的液膜可以较好地对燃烧室的头部起到热防护作用，使头部不出现明显高温。

（4）当工况条件确定时，在沉降型液膜冷却中，液膜长度主要由形成液膜的流率确定，液膜流率越大，液膜长度越大，对壁面的保护作用越好。

1Bai C X， Gosman A D. Development of methodology for spray impingement simulation ［J］. Sae Transactions， 1995，（104）: 21

2周红玲， 杨成虎， 刘犇. 液体火箭发动机液膜冷却研究综述［J］. 载人航天， 2012， 18（4）: 8～13

3Ueda T， Tanaka H. Measurements of velocity， temperature and fluctuation distributions in liquid films［J］. International Journal of Multiphase Flow， 1975， 2（3）: 261～272

4Wasden F K， Dukler A E. Insights into the hydrodynamics of free falling wavy films［J］. AIChE Journal， 1989， 35（2）: 187～195

5Gater R A， L' Euyer M R， Warner C F. Liquid-film cooling， it's physical nature and theoretical analysis， AD627028［R］. Lafayette，Indiana: Jet Propulsion Center， Purdue University， 1965

6William M G. Liquid film cooling in rocket engines［Z］. 1991

7Stechman R C， Joelee O， Howel J C. Design criteria for film cooling for small liquid - propellant rocket engines［J］. Journal of Spacecraft and Rockets， 1969， 6（2）: 97～102

8Knab O， Preclik D， Estublier D. Flow field prediction within liquid film cooled combustion chambers of storable bipropellant rocket engines［R］. Aiaa Journal， 1998， 1998～3370

9Gaffney R L， White J A， Girimaji S S， et al. Modeling turbulent / chemistry interactions using assumed PDF method［Z］. Aiaa， 1992

10Seamans T F， Vanpee M， Agosta V D. Development of a fundamental model of hypergolic ignition in space-ambient engines ［J］. Aiaa Journal， 1967， 5（9）: 1616～1624

11林庆国， 杨成虎， 刘彝. 射流角度和壁面曲率对撞壁液膜的影响［J］. 国防科技大学学报， 2013， 35（2）: 17～21

1009-8119（2016）09（1）-0059-04