固体火箭发动机喉衬流场及热结构耦合分析①

张晓光，王长辉，刘 宇，任军学

(北京航空航天大学宇航学院，北京 100191)

0 引言

喉衬构成燃气流动的最小通道，控制燃烧室压强，决定发动机性能，是固体火箭发动机的关键部件;而其所处燃气环境又最为恶劣。因此，喉衬也是发动机设计工作的薄弱环节。

流场及热结构耦合分析是评估喉衬气动特性、热防护有效性及结构完整性的重要手段［1］。传统的分析过程通常采用流场及热结构顺序计算的方法:先由一维等熵气体动力学公式获取燃气温度、压强等流场参数，由半经验公式(如Bartz公式)估算设定内壁温度下的内壁对流换热系数;再以燃气温度及内壁对流换热系数为热边界条件，进行固相材料导热计算，得到温度分布;最后，以流场及传热分析所得喉衬内壁压强与温度场为边界条件和载荷条件进行结构分析［2－5］。计算中一般未考虑内壁对流换热系数随内壁温度的变化，未实现流动计算与传热计算的耦合，这与实际情况不符;另外，流场简化计算模型及内壁对流换热系数半经验估算公式的应用，也对分析精度有较大影响。

针对上述问题，本文基于FLUENT流体计算软件与ANSYS结构分析软件，开展了喉衬流场及热结构耦合分析。计算分为2部分:首先，由FLUENT完成流固耦合烧蚀传热计算;然后，由ANSYS基于FLUENT的流场及热分析结果进行结构计算。计算实现了流场、热的双向耦合及流场、热到结构的单向耦合，兼顾了分析精度与计算资源两方面的因素，并充分发挥了2个软件各自优势，有助于准确预示喉衬在燃气环境中的烧蚀传热及热结构行为。

1 计算模型

本文的计算模型为Swaminathan［6］的传热试验发动机石墨喉衬，如图1所示。在喉部位置，4个热电偶沿圆周方向间隔90°埋置于不同深度:T1(81，20)、T2(81，24)、T3(81，29)、T4(81，34)。发动机采用聚氨酯推进剂，燃烧室压强3.43 MPa，燃烧温度3 203 K，燃烧时间7 s。喉衬石墨材料密度1 810 kg/m3，热导率37.7 W/(m·K)，比定压热容 711.3 J/(kg·K)。

图1 喉衬结构及热电偶布置Fig.1 Schematic of throat insert and thermocouple locations

2 流场与烧蚀传热耦合计算

由FLUENT流体计算软件采用整场离散、整场求解的方法，模拟燃气与固相材料间的非稳态流固耦合烧蚀传热。计算区域包括燃气流动区与固相材料区，各区域采用通用的控制方程，区别仅在于广义扩散系数及广义源项的不同;流固耦合界面成为计算区域的内部边界，其对流换热系数由软件本身的计算得出，而不再是分区求解时，需要人为设定的外边界条件;控制方程由有限体积法离散，可保证流固耦合界面满足温度连续、热流密度连续条件［7］。湍流模型采用Realizable k-ε两方程模型，近壁处理采用增强型壁面函数。

喉衬的烧蚀计算则依据热化学烧蚀理论，由FLUENT流体计算软件中的Wall Surface Reaction模型完成。热化学烧蚀理论认为，喉衬烧蚀由内壁面C与燃气氧化性组分H2O、CO2、OH等发生异相化学反应所致，反应速率遵循Arrhenius定律，而由反应所致的壁面退移速率就是烧蚀率［8］。

由NASA-CEA［9］程序进行热力计算，结果表明，燃气组分浓度沿喉衬轴向变化较小。因此，可认为喉衬流动为冻结流。氧化性组分 H2O、CO2质量分数为0.176、0.09，OH因浓度较小而不予考虑。烧蚀反应方程式及其动力参数见表1［8］。

表1 反应模型及动力参数Table 1 Reaction schemes and kinetics parameters

建立喉衬的二维轴对称计算模型，并划分四边形结构化网格，见图2。在喉部附近沿轴向加密，在内壁附近沿径向加密，且使流场近壁面第一层网格y+＜1，以满足增强型壁面函数法的要求，并更好地识别湍流边界层。由于喉衬入口就是喷管入口［6］，且为亚声速，因此给定总压3.43 MPa、总温3 203 K及进口气流角0°，静压则由内场外推;喉衬出口为超声速，全部参数外推;喉衬入口、出口的湍流条件则通过给定湍流强度和水力直径来确定。

图2 喉衬网格划分Fig.2 Computational grid of throat insert

本文喉衬非稳态流固耦合烧蚀传热计算只针对发动机稳态工作段，而不涉及达到稳态工作之前的压强上升段。计算时，首先设置内壁为绝热，外壁及端面为环境温度，进行稳态计算获得喉衬流场分布;然后，设置外壁及端面为绝热，内壁为流固耦合壁面及烧蚀反应发生表面，以稳态流场计算结果为初场，对流场区域与固体区域作整场求解，进行流场与烧蚀传热的非稳态耦合计算，计算时间7 s，时间步长取为0.000 5 s。

各热电偶测点温度计算值与试验值的对比见图3。二者基本相符，但计算值普遍高于试验值，且距离内壁越近，该现象越明显。分析其原因:在计算方面，外壁及端面设置为绝热，而实际是导热的;在试验方面，主要是热电偶动态热响应特性的影响。计算与试验均表明，距离内壁越近，温度及温升速率也越高。

图3 温度计算值与试验值对比Fig.3 Comparison between calculated and measured temperatures

图4为t=1、4、7 s时刻喉衬流场与固相材料的温度分布。由图4可见，随发动机工作进行，热响应逐渐深入到固相材料内部，而燃气与内壁之间强烈的对流换热及固相材料内部的导热，使固相材料形成明显的径向温度梯度，且以喉部区域最为显著。

图4 不同时刻喉衬流场与固相材料的温度分布Fig.4 Temperature distribution of flow field and solid material at various time

图5为t=1、4、7 s时刻喉衬内壁温度及烧蚀率的轴向分布。由图5可知，喉衬内壁温度及烧蚀率在喉部区域高，离开喉部区域则迅速下降，且峰值不在几何喉部截面，而是位于喉部上游。这基本遵循内壁热流密度的分布规律，见图6。文献［10－11］中，喷管烧蚀试验结果也证实，烧蚀率峰值位于喉部上游，而不是在喉部截面。

图5 不同时刻喉衬内壁温度及烧蚀率轴向分布Fig.5 Axial distribution of temperature and erosion rate of throat inser wall at various time

图6 不同时刻喉衬内壁热流密度轴向分布Fig.6 Axial distribution of heat flux at various time

图7为以燃气总温为定性温度时的内壁瞬时对流换热系数轴向分布。与图7不同，由Bartz公式预示的内壁对流换热系数，其峰值总在喉部截面。这是因为Bartz公式是基于管内旺盛湍流换热试验关联式及一维等熵气体动力学公式推导而来，未详细考虑湍流边界层发展及喉衬流道二维效应的影响［12－13］。实际上，从收敛段开始，边界层将逐渐减薄，传热速率不断增大;而在喉部上游某处，边界层外缘将达到声速，该处边界层最薄，燃气质量流率ρeue、热流密度均达到峰值，材料的热化学烧蚀率也最大;之后在扩张段，边界层又将迅速增厚，传热速率也急剧下降［12－14］。

另外，由图7也可知，随传热过程的进行，燃气与喉衬内壁的温差逐渐减小，对流换热逐渐减弱，对流换热系数也逐渐降低，从中可见将流动与传热耦合计算的必要性。

图7 不同时刻喉衬内壁对流换热系数轴向分布Fig.7 Axial distribution of convection heat transfer coefficient at various time

3 结构分析

在发动机工作过程中，喉衬承受燃气压强及由传热烧蚀所致变温载荷的作用，为了分析喉衬的结构完整性，需进行基于流场、热分析结果的结构计算。由于结构位移、应变对流场、温度场影响较小，因此忽略结构到流场、热的耦合。

将FLUENT计算的网格及流场温度场数据导入ANSYS结构分析软件。删去燃气区域网格，将固相材料区域网格单元类型设置为PLANE 42;喉衬瞬态温度场及内壁燃气压强作为载荷及边界条件加载;施加位移约束条件，使喉衬外壁径向位移为0，端面轴向位移为0。取弹性模量E=10.68 GPa，泊松比 μ =0.3，线膨胀系数 α =1.7×10－6/K。

由图4所示温度载荷及图8所示压强载荷计算了t=1、4、7 s时刻喉衬的应力分布。表2列出了各应力分量的最大、最小值，符号为正代表拉应力，符号为负代表压应力。由于为空间轴对称问题，因此仅有4个应力分量:径向正应力σρ、轴向正应力σz、环向正应力σφ及作用在圆柱面而沿轴向的切应力τρz。由表2可知，各应力分量中最主要的是环向压应力。

图9为t=1、4、7 s时刻喉衬的环向应力分布。热响应的深入，使喉衬温度分布不均加剧，喉衬整体应力水平随之升高，而喉部区域温度梯度最大，相应地其应力值也最大。

图8 喉衬内壁压强载荷分布Fig.8 Pressure loading distribution on the inner wall

表2 不同时刻喉衬应力最大值Table 2 Maximum/minimum stress at various time

图9 不同时刻喉衬环向应力分布Fig.9 Hoop stress distribution at various time

4 结论

(1)基于FLUENT流体计算软件及ANSYS结构分析软件，建立了固体火箭发动机喉衬流场及热结构耦合分析模型，实现了流场与烧蚀传热的双向耦合及流场、热到结构的单向耦合。

(2)算例分析结果表明，喉衬温度计算值与试验值吻合较好;内壁对流换热系数随壁温升高而降低;内壁热流密度的峰值在喉部上游入口处，内壁温度及烧蚀率也遵循相同的分布规律;喉部区域热流密度大，烧蚀率大，温度及其梯度高，应力水平也较高，环向压应力是主要应力分量。

(3)该方法兼顾了分析精度与计算资源两方面的因素，并充分发挥了2个软件各自的优势，有助于准确预示喉衬在燃气环境中的烧蚀传热及热结构行为。

［1］Lemoine L.Solid rocket nozzle thermostructural behavior［R］.AIAA 75-1339.

［2］吴景福.喉衬热结构数值分析［J］.推进技术，1985，6(4):38-52.

［3］李耿，周为民，张钢锤.喷管喉衬温度场及热应力计算分析［C］//中国航空学会航空动力分会火箭发动机专业委员会2006年学术年会文集，2006:183-188.

［4］Cozart A B，Shivakumar K N.Stress analysis of a 3-D braided composite ablative nozzle［R］.AIAA 99-1277.

［5］付鹏，蹇泽群，张钢锤，等.发动机喷管喉衬烧蚀及热结构工程计算［J］.固体火箭技术，2005，28(1):15-19.

［6］Swaminathan V，Setty S，Ramamadhav B S，et al.A method for the evaluation of temperature profiles and surface temperature in rocket motor nozzles with varying heat transfer coefficients［R］.AIAA 75-1347.

［7］陶文铨.数值传热学［M］.西安:西安交通大学出版社，2001:485-486.

［8］Thakre P，Yang V.Chemical erosion of carbon-carbon/graphite nozzles in solid-propellant rocket motors［J］.Journal of Propulsion and Power，2008，24(4):822-833.

［9］Gordon S，McBride B J.Computer program for calculation of complex chemical equilibrium compositions and applications［R］.NASA RP-1311，1994.

［10］陈博，张立同，成来飞，等.3D C/SiC复合材料喷管在小型固体火箭发动机中的烧蚀规律研究［J］.无机材料学报，2008，23(5):938-944.

［11］陈博，张立同，成来飞，等.燃气发生器条件下穿刺C/C复合材料喷管的烧蚀性能研究［J］.无机材料学报，2008，23(6):1159-1164.

［12］Back L H，Massier P F，Gier H L.Convective heat transfer in a convergent-divergent nozzle［J］.International Journal of Heat and Mass Transfer，1964，7:549-568.

［13］Keizers H L J，Veraar R G.Analysis of heat transfer in nozzles［R］.AIAA 96-3290.

［14］Williams F A，Huang N C，Barrere M.固体推进剂火箭发动机的基本问题(上册)［M］.京固群，译.北京:国防工业出版社，1976:122-124.