双组元推力器雾化燃烧过程数值模拟及液膜冷却作用分析

曹 顺，汪凤山，陈 健

（北京控制工程研究所，北京 100190）

双组元推力器雾化燃烧过程数值模拟及液膜冷却作用分析

曹 顺，汪凤山，陈 健

（北京控制工程研究所，北京 100190）

在带有液膜冷却双组元推力器的设计过程中，对燃烧室的数值仿真是一项非常重要的工作.充分考虑了常被忽略或简化处理的推进剂雾化、液滴破碎、液膜形成等重要过程，选取EDC燃烧模型和rea lizable k-ε湍流模型，采用有限体积计算方法得到燃烧室内部液滴分布、液膜分布、静温分布等重要数据.最后，根据所得结果阐述了燃烧室内的详细工作过程，并着重分析了液膜的分布特点以及对壁面的冷却作用，得到了两个与液膜冷却相关的结论.

双组元；液膜冷却；燃烧室；数值模拟

双组元推力器是主要用于航天器的姿态控制、轨道控制、姿态机动和位置保持等功能的执行机构，具有非常重要的作用.为了设计性能优化的推力器，需要对燃烧室内的雾化、液膜、燃烧、流动等过程有清楚的认识.但燃烧室内部工作过程极其复杂，工况非常恶劣，目前很难完全通过试验的方法得到温度场、压力场、速度分布等数据.随着计算流体力学、计算燃烧学等数值方法的逐渐成熟，仿真已经成为双组元推力器设计过程必要的一个环节.

然而到目前为止，国内对燃烧室的仿真绝大多数都忽略了液滴破碎、液膜形成等重要物理过程.总装备部指挥学院的聂万胜虽然对双组元推力器燃烧室内做了大量的仿真工作［1-3］，但均没有考虑推进剂雾化过程，同时也忽略了液膜对燃烧室的影响；许坤梅博士对双组元推力器的燃烧、流动过程进行仿真［4］，但没涉及液滴二次破碎和液膜形成等重要过程；张锋等［5］虽然分析了液膜对燃烧室壁面的冷却作用，但是对推进剂由液滴到液膜形成过程缺乏足够的研究，也没有得到液膜在壁面分布的详细情况.本文总结了前人的工作经验，充分考虑了推进剂喷注、液滴破碎和液膜形成等因素，建立合理的燃烧模型和湍流模型，对国内新近研制的25N双组元推力器燃烧室的工作过程进行了数值仿真，着重分析了液膜冷却对推力器工作过程的重要影响，为推力器的设计提供了重要的参考.

1 数学模型与仿真计算

1.1 推力器物理过程描述

双组元推力器内部工作原理极其复杂，涉及到多个复杂过程的耦合.图1简单介绍了双组元推力器的基本结构和工作原理.在忽略壁面传热的前提下，本文将推力室内部工作过程合理地简化为以下几步：

1）推进剂由喷注器喷出，在各种力的作用下连续液体破碎为细小液滴；

2）液滴与燃烧室内壁相互作用，形成液膜；

3）液滴在运动过程中与气相作用，蒸发为气体；

4）NTO（N2O4）遇高温快速分解为NO2；

5）MMH（甲基肼）与NO2反应产生高温燃气；

6）高温燃气由燃烧室经扩张喷管喷出.

图1 双组元推力器内部工作原理示意图Fig.1 The working principle scheme of bipropellant thruster

1.2 连续相控制方程

双组元推力器内部的流动可以通过N-S方程来描述.由于推力器采用双漩涡的喷注方式，燃烧室内燃气带有很强的旋流.因此，本文采用带旋流修正的realizable k-ε湍流模型.燃烧模型采用工程上应用比较成熟的EDC模型.离散方式采用二阶迎风格式.连续相控制方程包括连续方程、动量方程（忽略体积力）、能量方程、组元输运方程（组分 k），分别如式（1）～（4）所示：

式中，ρ为密度，t为时间，xi、xj为位移变量，ui、uj为速度变量，p为压力，h为混合物静焓值，hj、hk分别为组元j和k的静焓，T为温度，λeff为有效传热系数，jki是组分k的扩散流量，τij为粘性应力，Yk是组分k的质量分数，D为组分Yk的扩散系数.上述方程中的源相来源于离散相（液滴）和气相燃烧.其中质量源相Sm是由液滴蒸发产生的；动量源相SM，i等于在一个网格内所有液滴动量的损失值；Sh包括由燃烧产生的能量和液滴蒸发的能量交换两项；SYk是化学反应的反应产物.可假设气相为理想气体，其状态方程为

式中，Ru为通用气体常数，Yi为i组元的质量分数，Mi为i组元的分子量.

MMH与NTO之间的化学反应采用一步总包反应模型.其化学反应方程式为

由于N2O4在常温下就非常容易离解为NO2，在燃烧室内的高温条件下，其自身离解速度将远远高于N2O4与MMH的反应速度.因此可以认为N2O4在由液态转化为气态的一瞬间便全部离解为NO2.所以，燃烧室内实际发生的化学反应可近似为

1.3 离散相模型

离散相模型主要包括推进剂雾化、液滴蒸发、液膜模型等.考虑到双组元推力器的结构和喷注方式，本文采用压力漩涡喷注模型，液滴的破碎过程采用TAB液滴破碎模型.

液膜的形成过程非常复杂，一般可以这样描述［6］：推进剂液滴由喷注器喷出以一定的速度撞到燃烧室内壁，与壁面相互作用，部分液滴被粘滞，部分铺开附着，部分反弹，还有一部分液滴破碎并溅射，其中附在壁面的液滴形成壁面液膜，液膜在各种力的作用下在边缘处与壁面发生分离.在液膜形成过程中，液膜表面与气相进行相互作用，包括液膜蒸发、对流换热、受到气相的拖拽力，液膜底部与固体壁面通过传热进行能量交换.其原理参见图2.

图2 液膜形成机理示意图Fig.2 Scheme of liquid film formation mechanism

本文综合考虑了双组元推力器双漩涡喷注方式、燃烧室内工作特点以及工程应用成熟度等因素，采用了Stanton和O’Rourke的壁面液膜模型［7-9］.此液膜模型充分考虑了液滴与壁面相互作用区域内的主要信息，如液滴速度、物性、温度、壁面温度、壁面粗糙度等.模型描述的液滴与壁面的相互作用方式与双组元推力器液膜形成机理相符.因而，采取Stanton和O’Rourke的壁面液膜模型是合理的.

基于液滴的撞击能量和液滴的沸点区分四种不同作用方式.其中液滴撞击能量定义为：

ρ是液滴密度，Vr是液滴速度，D是液滴直径，σ是液滴表面张力，δbl是边界层厚度，H0是液膜厚度.

当液滴撞击能量小于16J且液滴温度低于其沸点，液滴会粘滞到壁面；当液滴温度高于其沸点时，并以一特定的角度和速度与壁面相互作用，液滴会铺开附着于壁面；当液滴温度和壁面温度均高于液滴沸点，同时撞击能量小于临界撞击能量Ecr（根据相关文献Ecr取为57.7J［7］）时，液滴将发生反弹，其反弹能量恢复系数为：

其中：θi为液滴对壁面的撞击角度；当液滴温度与壁面温度均高于液滴沸点，且撞击能量大于临界撞击能量Ecr时，液滴与壁面相互作用后，将溅射出多个新的液滴.在液膜的边缘处，当液膜边缘应力超过壁面对液膜的粘滞力时，液壁将发生分离.二者受力关系非常复杂，主要取决于液膜边缘的液体粘度、温度、速度以及表面情况等当地条件.关于液膜控制方程等其他详细信息参见文献［7-9］.

1.4 计算几何模型

本文的计算模型选用国内新近研制的25N双组元推力器.推力器采用MMH/NTO推进剂，设计比冲为295N·s/kg，喷管膨胀比为150，壁面粗糙度为5μm.

本文采用了二维轴对称的几何模型，计算网格总数为281340.推力器初始条件见表1.燃烧室壁面设定为绝热无滑移的边界条件；燃气经过推力室喉部后成为超音速流，超音速流下游不影响上游，在确定出口边界条件时只需将内点值直接赋于出口边界上的流场参数即可.

表1 推力器稳态工作初始条件Tab.1 Initial condition of thrust in stable working state

2 计算结果与分析

经过数值计算，本文得到了双组元推力器稳态工况下燃烧室内重要物理参数的分布数据.如图3所示，液态推进剂由喷注器喷出后，在各种力的作用下破碎成大小不一的液滴（最大液滴直径为408μm，最小直径为1.27μm）.液滴表面在燃烧室高温的作用下不断的蒸发，由MMH和NTO摩尔分数和化学反应图可以看出，未反应的气态推进剂主要集中在燃烧室靠近喷注器的一侧，二者在这里通过充分的掺混，发生化学反应，放出大量的热.其反应速率最快的位置出现在喷注器出口附近，达到了337mol/（m3·s）.

图3 燃气室雾化与化学反应Fig.3 Particle atomization and reaction in combustion chamber

图4给出了在绝热壁面条件下燃烧室静温分布情况.在靠近喷注器的前半部分为化学反应不完全区域，同时由于液滴蒸发过程也需要吸收大量的热，因而其温度相对后部较低；后半部分到喉部以前属于化学反应完全区域，绝大多数化学能转化为燃气的热能，其温度相对较高，最高温度达到2880K.从图中也可以看出越是靠近燃烧室内壁面的区域温度越低，这是由于液膜的冷却效应引起的.

图4中标记了两点，A点和B点.A点位于燃烧壁面，B点与A的横向距离相同，纵向位置处燃烧室内壁与轴线之间，此处温度达到2880K.从图5中看出A点的温度约为2290K.A、B点的温差接近600K，这说明液膜对于燃烧室壁面的热防护作用非常明显.此外，从图6中还可以看出，燃烧室壁面温度是沿着轴线方向逐渐升高的，在喉部位置温度最高，超过2400K.

图4 燃烧室静温分布情况Fig.4 Static temperature distribution in combustion chamber

图5 燃烧室内壁温度（绝热壁面条件）Fig.5 Inwall temperature distribution（thermal insutated inwall）

图6 放大了燃烧室壁面局部区域的温度分布，进一步说明了液膜的冷却效果.从图中可以看到，越靠近壁面的区域温度越低，由于液膜的冷却作用，在紧贴壁面处会形成一层很薄的低温层.

壁面温度分布具有很强的波动性，主要是由液膜分布的不均匀性和间断性所致.而液膜这种不均匀和间断性是因为燃烧室内的燃烧和流动进程非常复杂，使得液滴的初始条件（速度、方向）各不相同，其与壁面的相互作用方式极其不规律造成的.

如果以喷注位置为轴线零点，液膜在壁面的总体分布区域从轴线方向0.35mm一直延伸至45mm左右.本文给出了4.6mm到8.6mm之间燃烧室内壁液膜厚度的分布图，如图7所示.从图中可以看出，A、B、C三点是明显的液膜间断处，不同位置的液膜厚度最大值为18.2μm.正是由于液膜分布的这种不均匀性，造成燃烧室壁面温度分布的波动.

图6 壁面区域静温分布放大图Fig.6 Static temperature distribution in zone adjacent inwall（zoomed picture）

图7 燃烧室内壁面液膜分布（局部）Fig.7 Distribution of liquid film on the combustion chamber inwall（local）

3 结 论

本文对双组元推力器燃烧室内的工作过程进行了数值仿真，考虑了在以往仿真工作中常被忽略的液滴破碎及液膜形成等重要进程，得到了液滴雾化、化学反应、静温、液膜分布等重要结果.根据数值仿真的结果，阐述了双组元推力器燃烧室内详细的工作过程，着重分析了液膜的分布特点以及对燃烧室的冷却效果，最后得到了两个关于液膜冷却的结论：

1）在双组元推力器的稳态工作过程中，液膜冷却的效果明显，对燃烧室壁面具有非常好的保护作用；

2）液膜在燃烧室内壁的分布是间断的，厚度不均匀（0～18.2μm），但这并不影响液膜的冷却效果.

［1］ 聂万胜，庄逢辰.自燃推进剂火箭发动机稳态燃烧过程的数值模拟［J］.推进技术，1998，19（5）：7-9 Nie W S，Zhuang F C.Numerical simulation of stable combustion process in earth storable bipropellant rocket engine［J］.Journal of Propulsion Technology，1998，19（5）：7-9

［2］ 聂万胜，庄逢辰.喷雾特性对液体火箭发动机燃烧稳定性的影响［J］.推进技术，2000，21（3）：56-59 Nie W S，Zhuang F C.Effect of spray characteristic on the liquid rocket combustion stability［J］.Journal of Propulsion Technology，2000，21（3）：56-59

［3］ 聂万胜，丰松，解庆纷，庄逢辰.液体火箭发动机燃烧不稳定性主动控制数值仿真［J］.工程热物理学报，2008，29（7）：1253-1255 Nie W S，Feng S J，Xie Q F，Zhuang F C.Numerical simulation of liquid rocket combustion instability active control［J］.Journal of Engineering Thermophysics，2008，29（7）：1253-1255

［4］ 许坤梅.双组元小推力液体火箭发动机推力室工作过程数值模拟［D］.北京：北京航空航天大学，2005 Xu K M.Numerical simulation of the working process in thruster chanber of bipropellant thruster［D］.Beijing： Beihang University，2005

［5］ 张锋，仲伟聪.膜冷却推力器室传热计算研究［J］.火箭推进，2009，35（4）：35-37 Zhang F，Zhong W F.Computational investigation of heat transfer for film cooling thrust chamber［J］.Journal of Rocket Propulsion，2009，35（4）：35-37

［6］ Naber JD，Reitz R D.Modeling engine spray/wall impingement［J］.General Motors Research Laboratories，SAE Transactions，1989，97（6）：118-140

［7］ O’Rourke P J，Amsden A A.A particle numerical model for wall film dynamics in port-fuel injected engines［C］.International Fuels and Lubricants Meeting and Exposition，San Antonio，TX，Oct 14-17，1996

［8］ O’Rourke P J，Amsden A A.A spray/wall interaction submodel for the kiva-3 wall film model［R］.Detroit，SAE Technical Paper，2000-01-0271，2000

［9］ Stanton DW，Rutland C J.Modeling fuel film formation and wall interaction in diesel engines［R］.Detroit，SAE Technical Paper 960628，1996

Num erical Sim ulation of W orking Process in Cham ber of Bipropellant Thruster and Analysis of Film-Cooling Effect

CAO Shun，WANG Fengshan，CHEN Jian
（Beijing Institute of Control Engineering，Beijing 100190，China）

In the design of film-cooled bipropellant thruster，the numerical simulation of the combustion chamber is a necessary step.In this paper，process of atomizing of the liquid propellant，breaking up of drop lets，liquid film formation are taken into account，and the EDC combustion model and the rea lizable k-εturbulent flow model are selected.Distributions of droplets，liquid film and static temperature are obtained by using the numerical simulation.At last，based on the analysis of simulation results，the working process in combustion chamber are interpreted，and two conclusions on the liquid film cooling are derived.

bipropellant；liquid film cooling；combustion chamber；numerical simulation

V434

A

1674-1579（2012）02-0058-05

曹 顺（1983—），男，硕士研究生，研究方向为航天器推进技术；汪凤山（1981—），男，工程师，研究方向为航天器推进技术；陈 健（1969—），男，研究员，研究方向为航天器推进技术。

2011-11-15

DO I：10.3969/j.issn.1674-1579.2012.02.012