固液捆绑火箭上升过程喷流干扰特性及底部热环境研究

2022-12-03严立孙培杰赵一霖包轶颖秦川江泽敏刘宽

上海航天 2022年5期

严立，孙培杰，赵一霖，包轶颖，秦川，江泽敏，刘宽

（1.上海宇航系统工程研究所，上海 201109；2.上海航天技术研究院，上海 201109）

0 引言

固体推进系统和液体推进系统的动力组合是国外主要采用的动力组合方式，该方式拥有固体推力大、加速快，液体工作时间长、推力调节容易等技术优点，可实现运载火箭的系列化、模块化、通用化和组合化生产制造。目前，国外捆绑运载火箭中固体助推器约占81%，大型、重型运载火箭几乎都将固体火箭发动机作为助推级的首选动力［1］，如美国的德尔它系列（除重型德尔它4）、大力神系列、宇宙神系列以及航天飞机［2-4］，西欧的阿里安系列3 和5火箭，日本的M 系列、N 系列及H 系列火箭等［5］都是如此。

在固液捆绑火箭的发展历程中，底部热环境预示及防护一直作为关乎成败的问题受到关注。固体火箭发动机排出炽热高速气流含有的Al2O3固体颗粒，除了大大增加火箭底部热辐射外，还影响火箭底部气流的流场分布，在助推与主发动机火焰间造成排气交叉效应和附加干扰，加剧固体发动机羽流和液体芯级羽流间的相互作用。

美国在研制德尔它Ⅰ［6］和大力神Ⅲ［7-8］的过程中，开展了底部热环境地面风洞试验研究，探究底部不同区域加热热流呈现的变化趋势。荷兰的REIJASSE等［9］参照大力神III 的地面试验方法，对阿里安5 芯级及固体助推3 个喷管状态底部喷流流态进行了缩比风洞喷流的可视化试验，详细分析了不同喷流压力与外界压力比情况下，喷流与喷管喷流之间、来流与喷流之间的干扰情况，直观了解喷流状态变化的机理。文献［10］采用Fluent 软件，对H-IIA 固液捆绑火箭发动机底部加热问题进行了数值模拟研究，分析了芯级液体发动机与4 枚固体助推器喷流相互干扰、回流现象。文献［11］对双发动机喷流干扰进行了实验研究，拍摄了喷流干扰纹影图片，测量了2 股喷流之间中心轴线的皮托压力，研究了高空喷流干扰所引起的回流现象。

我国在固液捆绑火箭研究方面起步较晚，目前主要针对单一固体发动机喷流及热环境开展了研究。文献［12-16］采用Fluent 软件对单一固体发动机的两相喷流流场进行了数值仿真，其中：乌岳［12］分析了不同飞行马赫数、飞行高度和不同颗粒相粒径对流场的影响；黄颗［13］分析了后燃反应模型和颗粒相对流场的影响；武利敏［14］对颗粒轨道模型和拟流体模型进行了分析。文献［17-18］采用数值仿真和试验研究的方法研究了固体火箭发动机喷流反流对底部的影响，杨学军［17］利用线性化热学参数单一介质数值模拟方法分析了真实飞行轨道条件下的热环境预示结果，并与地面试验及飞行试验测量结果进行了比对；李国良［18］对冷喷与热喷方法进行对比分析，结果表明2 组分热喷方法获得的数据与试验值一致性较好，较冷喷方法精度提升。文献［19-20］分析了捆绑火箭底部喷流及热环境，史亚男［19］数值模拟的方法分析了固体芯级捆绑固体助推火箭，不同助推器个数、不同助推器与芯级间距离对底部热环境的影响；戴欣怡［20］采用数值模拟方法对助推器内部点火、燃烧和燃面推移等瞬态过程以及助推器尾部超声速气-固两相流场分布规律进行了研究。综上所述，至今未见含芯级液体火箭捆绑多枚固体助推器发动机底部喷流交互作用及底部热环境研究的分析报道。

本文用连续相模型模拟气相、离散相模型（Discrete Phase Model，DPM）模拟固体粒子相、离散坐标模型（Discrete-Ordinates Model，DOM）模拟含有固体粒子的介质辐射，对固液捆绑火箭在上升飞行过程中的气固两相喷流的演变过程进行了仿真研究，分析了不同飞行高度下固体发动机与液体发动机交会喷流干扰特性、高温固体颗粒影响作用以及箭体底部表面热环境，并将仿真结果与飞行结果进行了比较。研究成果对固液捆绑火箭底部热防护设计具有重要指导意义。

1 仿真计算模型

1.1 数学计算模型

分析过程考虑固体粒子对喷流流场的影响，模拟过程分为3 个步骤：①气相喷流场模拟；②加入固相粒子DPM 模型，计算气固两相喷流；③加入热辐射模型，计算底部热流。

1.1.1 连续相模型

连续性方程为

式中：ρ为密度；xi为i方向矢量；ui为i方向的速度。

动量方程为

式中：P为静压；μ为动力黏度；xj为j方向矢量。

能量方程为

式中：xk为k方向矢量。

采用标准k-ε方程模型，控制方程如下：

湍动能k方程为

湍动耗散率ε方程为

式中：K为湍流动能；ε为湍流能量耗散率，常数系数σμ=0.09；c1=1.44，c2=1.92。

1.1.2 离散相模型

采用拉格朗日离散相模型考虑颗粒相与连续相之间的相互作用力。首先是作用在颗粒物上的两相间的拖拽力，表达式为

其中

式中：FD为单位质量阻尼系数；CD为阻力系数；u为流体相对速度；up为颗粒速度；μ为流体动力黏度；ρp为颗粒密度；dp为颗粒直径；Re为相对雷诺数；a1、a2、a3为常数。

1.1.3 离散坐标方法

DOM 计算含有固体颗粒的喷流介质辐射作用。离散坐标法基于对辐射强度的方向变化进行离散，通过求解覆盖整个4π 空间立体角上一系列离散方向上的辐射传递方程而得到问题的解。它可以计算所有光学厚度的辐射问题，并且计算范围涵盖了从表面辐射、半透明介质辐射到燃烧问题中出现的介入辐射在内的各种辐射问题。基于谱带模型的吸收、发射、散射性介质内辐射传递方程的表达式为

在三维直角坐标系(x，y，z)下，采用离散坐标法，上式右端积分项近似由一数值积分代替，并在离散的方向上对辐射传递方程求解如下：

式中：ξm、ηm、μm为辐射传输方向的方向余弦；ωl为积分系数；l、m为空间方向离散的第l个和第m个立体角（l，m=1，2，…，NΩ），NΩ为4π 空间方向离散的立体角总数；=Φk(Ωm，Ωl)为离散后的辐射相函数。

介质的吸收系数采用基于网格的加权求和模型，散射相函数选择各项同性，散射吸收选择常数。

1.2 物理计算模型

某固液捆绑火箭芯级采用2 台液氧/煤油发动机，捆绑4 个固体助推器，如图1 所示。由于箭体外部结构轮廓和发动机喷管布置位置的对称性，仿真模型采用箭体周向1/4 进行建模，计算区域包含1/4芯级以及1 个助推器。

图1 固液捆绑火箭几何模型Fig.1 Geometrical model of solid-liquid bound rocket

模型网格采用结构化六面体网格，网格数量1.2×107。综合考虑计算资源和计算精度，针对不同的计算区域采用不同大小的网格，对发动机推力室内部及喷管出口附近进行网格加密处理。计算网格如图2 所示。

图2 火箭模型计算区域及网格Fig.2 Rocket model calculation area and grid

固体火箭发动机采用铝丁羟复合HTPB 推进剂（Hydroxyl-Terminated Polybutadiene，HTPB），铝粉质量分数17%，燃烧产物中Al2O3颗粒质量分数约32%，固体颗粒尺寸分布对固体发动机燃烧室效率和喷管效率影响较大，是发动机性能分析及仿真分析的重要输入参数。但推进剂燃烧过程复杂，影响因素众多，颗粒项直径分布在几微米到几百微米之间，目前还没有成熟的理论来对颗粒项直径分布定量描述。为了综合考虑不同直径粒子对喷流的影响，采用Rosin-Rammler 分布来描述高温颗粒分布，其函数关系为

式中：Yd为比指定粒径d大的颗粒的质量分数；dm为平均粒径；n为传播系数。

颗粒平均直径dm采用质量平均直径［21-22］，公式为

中国农业大学水利土木工程学院党委书记杨培岭以《节水灌溉技术的未来发展方向和趋势》为题进行了精彩演讲，他呼吁要深入基础理论研究，加快节水灌溉科研成果的转化，实现节水灌溉技术的创新。要推广自动化控制系统，加强节水灌溉设备质量的监管控制，加强水资源管理，合理确定水价，建立健全节水灌溉体系服务。

式中：D为喷管喉部直径；Cm为100 g 推进剂中铝粉的摩尔质量；Pm为燃烧室压强；τ为颗粒在燃烧室中的驻留时间。

边界条件设置如下：

1）入口边界，设定压力入口：固体发动机燃烧室总压6.3 MPa，喉部压力3.6 MPa，总温3 236 K；液体发动机燃烧室总压17.7 MPa，喉部压力10.2 MPa，总温3 810 K。

2）出口边界，在喷管的远场处设定为压力出口：地面常压环境设置为1 atm，温度300 K；高空设置为当地的环境压力和温度。

3）壁面边界，即设定壁温边界条件。

4）对称边界，即设定对称边界条件。

本文选取火箭飞行过程中6 个不同飞行高度下的工况进行分析，各工况来流参数见表1。

表1 火箭飞行计算工况来流参数Tab.1 Parameters of the coming flow under different calculation conditions

2 计算结果与分析

2.1 流场分析

提取2 个芯级轴线和芯级固发轴线截面上的参数，在不同飞行高度上外流场速度分布云如图3所示。由图可知，飞行高度10 km 以下时，火箭飞行马赫数较小，环境压力较高，发动机喷流扩张受到抑制，形成桶鼓状的激波，即膨胀波和压缩波交替的激波形式，发动机喷流集中在中心区域，芯级和助推喷流之间相互干扰作用较小，如图3（a）和图3（b）所示。随着飞行高度增高和速度增大，箭体逐渐进入超音速飞行段，外界气压降低，从工况3 开始，环境压力小于发动机出口压力，箭体周围逐渐产生不同类型的激波。箭体顶端形成弓形激波（Bow Shock），喷管喷流外缘与自由来流相互干扰，因喷流羽流对自由来的流阻碍形成羽流干扰激波（Plume Induced Shock），如图3（c）所示。且随飞行高度增加，压力比逐渐增大，喷管喷流扩张角逐渐扩大，自由来流与发动机羽流干扰效应逐渐加强，产生喷流交互作用，喷流向底部回流，对底部形成热、力冲刷作用，如图3（d）～图3（f）所示。

图3 不同工况下火箭飞行外流场速度分布云Fig.3 Contours of the velocity distribution under different working conditions

以上喷流流场分析可知，在低空状态下，发动机之间无直接喷流交互作用，而在高空状态下，喷流扩张，发生空气来流、芯级喷流和固发喷流之间复杂的交互现象，在交互过程中，产生喷流回流，造成对底部的力、热影响。喷流干扰回流高温区域如图4 所示。

图4 底部喷流交互作用高温区分布示意图Fig.4 Distribution diagram of high-temperature areas caused by the bottom jet interaction

工况2、工况3 和工况5 的喷流温度分布如图5～图7 所示。由图5 可知，工况2 时，由于外界气压较高，各发动机喷流集中在中心区域，芯级和助推喷流之间相互干扰较小，高温区出现在发动机的出口下方。在工况3 时，发动机之间有了明显的喷流交互作用，喷流干扰引起了局部高温区，如图6 所示，首先是芯级2 个液发喷流之间发生交互作用，在5区出现高温区。其次由于2 固发与2 液发共4 个发动机喷流之间的交互作用，在1 区和3 区也出现高温区，最高温达2 200 K，高于同截面喷管中心喷流温度。在工况5 时，随着喷流扩张角增大，喷流干扰高温区域向底部移动，在2 区和4 区出现明显的喷流交互高温区，该高温区由2 固发与1 液发共3 个发动机喷流之间的交互作用产生，如图7 所示。

图5 工况2 喷流温度分布Fig.5 Temperature distribution of the jet flow under working condition 2

图6 工况3 喷流温度分布Fig.6 Temperature distribution of the jet flow under working condition 3

图7 工况5 喷流温度分布Fig.7 Temperature distribution of the jet flow under working condition 5

从以上分析可以看出，喷流交互作用受发动机扩张角、大气环境、飞行速度和喷管相对位置等因素影响，发生交互作用的形态和时刻都有所不同。最高温区在1、3 和5 区，次高温度在2、4 区，且随着喷流的扩张，高温交互出现的轴向方向距离箭体底部越来越近。

2.2 固体颗粒影响分析

分别对小粒径（3 μm）颗粒和大粒径（100 μm）颗粒对固液捆绑火箭气固两相流流场的影响进行了数值仿真分析，数值仿真结果如图8 所示，由图可知，3 μm 直径的颗粒由于粒径较小，随流性较好，受气体膨胀波影响较大，分布在喷流与空气的边界混合层区域；而100 μm 的大粒径，随流性较差，受气体膨胀波影响较小，喷出喷管后，呈自由膨胀状态流动，对交互流场核心区的影响较大。但从2 种直径颗粒的交互流场来看，都不存在交互流场固体颗粒反流的情况，这与交互流场的气体反流的状态有本质区别。

图8 不同粒径的两相流流场分布Fig.8 Flow distribution of the two-phase flow with different particle sizes

2 种直径颗粒的温度场和速度场对比结果如图9 所示，2 个喷流流场的膨胀波外形尺寸相当，边界混合层的参数分布相似，固体颗粒大小对喷流流场的膨胀程度影响不大，但核心区温度与速度参数分布存在差异，100 μm 颗粒核心区的温度衰减滞后于3μm颗粒，且速度低于3μm颗粒，主要是由于100μm 颗粒的流场中，颗粒呈自由膨胀状态流动，4个固体助推器喷出的固体颗粒在喷流流场的中后部产生交互作用，对核心区的流场影响较大。

图9 不同粒径的温度场及速度场对比Fig.9 Comparison of the temperature field and velocity field with different particle sizes

固液捆绑火箭底部从芯级底部的防热板中心到流场计算域边界的轴线上的参数分布如图10 所示。由压强变化曲线可知，3 μm 颗粒的流场参数与100 μm 流场参数吻合较好，这与之前二维流场分析所获得的结果一致，固体颗粒对流场的压力影响可以忽略；由温度变化曲线可知，3 μm 流场的轴线温度与100 μm 流场的轴线温度变化趋势一致，但在轴向距离40 m 后，100 μm 颗粒流场的温度高于3 μm颗粒流场，主要是由于4 个固体助推器的固体颗粒在流场中心处产生交互作用，使温度进一步升高；由速度以及马赫数轴向曲线可知，在轴向距离40 m后，100 μm 颗粒流场的速度低于3 μm 颗粒流场，与温度场变化规律一致，这是由于固体颗粒交互作用对流场产生减速作用。

图10 流场沿轴线的参数分布Fig.10 Parameter distribution along the axis in the flow field

2.3 底部热环境分析

火箭底部热流分布云如图11～14 所示。由图可知，底部热流分布与发动机喷流状态相关，底部对流热流密度最大值最早出现在2 芯级发动机喷管之间，随着飞行高度的增加，回流冲击位置向径向移动，一直到工况5 时，移动到底板边缘，并对火箭底部侧壁产生热冲刷作用。回流冲击底板热对流大小受回流温度、速度和密度等影响，在工况4 时出现峰值，最大对流热流密度为90 kW/m2。喷流对底部的辐射作用与表面和喷流之间的相对可视距离密切相关，可以看出芯级发动机与助推发动机轴线之间的辐射热流较大，辐射热流变化范围为160～400 kW/m2。

图11 工况2 时底部热流分布Fig.11 Heat flux distribution at the bottom under working condition 2

在火箭飞行试验中，底部防热板上布置了3 个热流测点，如图15所示。

图15 底部测流点计算位置Fig.15 Locations for calculation

图12 工况3 时底部热流分布Fig.12 Heat flux distribution at the bottom under working condition 3

图13 工况4 时底部热流分布Fig.13 Heat flux distribution at the bottom under working condition 4

图14 工况5 时底部热流分布Fig.14 Heat flux distribution at the bottom under working condition 5

测流点的对流、辐射与总热流值如图16 所示。与测量值的比较如图17 所示。测流点位置热流在起飞10 s 达到最大值，最大值分别410 kW/m2和472 kW/m2，平均热流为191 kW/m2和222.7 kW/m2。仿真的热流最大值为434.1 kW/m2和452.4 kW/m2，平均热流为334 kW/m2和305 kW/m2。由图可知，仿真最大热流与测量值吻合较好，平均值是测量值的1.5 倍左右，变化趋势上略有不同。测量值在10 s 内快速上升，然后缓慢下降，在80～100 s 内增大，最后降低；仿真值是起飞后一直缓慢上升，到90 s 左右增大到最大值，然后开始下降。分析仿真值与测量值趋势不同的原因在于：1）仿真中未考虑起飞过程发射台的遮挡作用，起飞过程中的热流偏低；2）实际飞行过程中壁面温度会一直上升，而仿真中未考虑流场与壁面的耦合温升过程，热流值一直增大。