复合喷管结构温度场流固耦合仿真分析①

傅学金，邓 恒，寇军强，宋 军，陈 刚

(1.第二炮兵工程大学，西安 710025;2.西安现代控制技术研究所，西安 710065;3.第二炮兵装备部驻航天科技集团型号办事处，北京 100076)

复合喷管结构温度场流固耦合仿真分析①

傅学金1，2，邓 恒2，寇军强2，宋 军2，陈 刚3

(1.第二炮兵工程大学，西安 710025;2.西安现代控制技术研究所，西安 710065;3.第二炮兵装备部驻航天科技集团型号办事处，北京 100076)

为了评估复合喷管热防护性能以及获取喷管烧蚀和结构应力分析的工况条件，运用Fluent流体动力学软件，对复合喷管的结构温度场进行了数值仿真。分析中，采用了两方程RNG k-ω湍流模型和增强型壁面函数，利用流固耦合的计算方法，获得了喷管结构瞬态温度场的计算结果，重点分析了结构温度场最终分布状态和初期传播特点，以及喉衬温度随时间的变化规律，估算了喉衬的烧蚀。分析结果表明，喷管结构热防护性能满足要求，温度最高区域位于喷管收敛段中后部，喉衬线烧蚀量约为2.1 mm，为喷管结构进一步优化设计提供了重要参考依据。

流固耦合;温度场;复合喷管;数值仿真

0 引言

在高压力、高燃温、长时间工作的中大型发动机设计中，多选用复合材料喷管的方案。复合喷管的结构设计是喷管设计中十分重要的内容，包括了耐烧蚀层、隔热层和喉衬的结构设计和各部件之间的接缝匹配设计。耐烧蚀层、隔热层和喉衬的结构设计主要是以对喷管的烧蚀和传热的评估作为依据，而各部件之间的接缝匹配设计则是以喷管的热应力分析作为依据的。无论是喷管的烧蚀、传热和热应力的分析，都要以喷管温度场的分析为前提条件。因此，在喷管结构设计中，喷管结构温度场分析是十分重要的工作。

目前，对喷管结构温度场的计算有3种计算方法:(1)对喷管内壁面施加温度边界条件，计算喷管结构的温度场［1］;(2)考虑燃气与壁面的对流换热作用，利用巴兹公式计算燃气与壁面之间的对流换热系数，以此作为边界条件计算结构的温度场［2－8］;(3)将喷管内流场与结构的交界面作为流体域和固体域的热耦合边界，进行流场和传热的耦合计算［9－11］。在第1种方式中，对喷管内壁面施加的是温度边界条件，用热传导的传热方式替代了喷管内壁面与燃气之间的对流换热方式，这改变了边界条件的性质。在第2种方式中，虽考虑到了喷管内壁面与燃气的对流换热边界条件，但巴兹公式所需的壁面温度、压力和马赫数都是从绝热壁面边界条件的喷管流场计算结果中取得的，没有考虑结构对温度场结果的耦合影响作用。第3种计算方式就将流场和温度场耦合起来进行计算，这符合喷管的实际状态，计算结果更为准确可信。文献［10］就明确指出，喷管内流场和结构温度场的解耦计算是不合适的;文献［11］则利用试验测试的喷管内壁面的温度数据验证了耦合计算的正确性，与解耦计算结果的对比分析表明，耦合计算比解耦计算结果更为准确。综上分析，本文将采用流场与结构耦合传热的方式计算某型喷管结构的温度场，分析初始阶段和工作结束时温度场，以及喉衬温度随时间的变化规律，为喷管结构优化设计提供重要参考依据。

1 计算方法

采用Fluent软件，对喷管的流场和流场对固体的传热进行整体分析。Fluent软件运用了有限体积法离散流体动力学微分控制方程，将其变换为一组代数方程来求解。流体动力学控制方程主要包括了质量守恒方程、动量守恒方程和能量守恒方程，可写为如下的通用形式:

式中 φ为通用变量，可代表速度u、v、w和温度T等变量;Γ为广义扩散系数;s为广义源项，它们在不同的控制方程中具有不同的表达式。

对于涉及到传热求解的问题，必须使用能量守恒方程来对温度场进行求解，此处特将能量方程展开如下:

式中 T为温度;k为传热系数;cp为比热容;sT为粘性耗散项。

从式(2)可看出，在流体域内，温度的变化与密度、速度和温度梯度都有关系，在固体域内，介质的速度为零，密度为常数，因而温度的变化仅与温度的梯度有关系，在流体域和固体域的耦合界面上，温度在要满足对流体域内温度求解的同时，还要满足对固体域的求解，即符合连续性的要求。从以上分析可看出，要获得喷管结构准确的温度场，必须将流场域的流动分析与结构的固体传热分析耦合起来。因此，设置流体域与固体域直接的边界为耦合边界。

在工业化和城市化快速推进过程中，欧美国家也经历了乡村发展滞后于城市、农村居民收入增长滞后于城市居民收入增长的不协调阶段。中共中央党校(国家行政学院)《行政管理改革》编辑部编辑张世贵梳理了发达国家破解城乡发展差距问题的经验。

由于喷管的结构温度场取决于喷管内的流动，所以对喷管内流场的准确分析显得十分必要。湍流模型对喷管内流场的计算具有很大的影响，在本分析中，选择RNG k-ε湍流模型。该模型应用的范围比标准k-ε模型更广泛，对低雷诺数和高雷诺数的流动都能进行很好的分析，能更好地反映出高切应变率和流线曲率对流动影响的效应。固体火箭发动机喷管内的燃气流动经历了收敛段的亚音速、喉部的音速和扩张段的超音速，雷诺数可能从低雷诺数变化到高雷诺数，并且先收敛、后扩张的构型，使得靠近喷管壁面的流线存在较大的曲率。因此，RNG k-ε湍流模型适合对喷管内流场的计算。

对喷管壁面函数的处理直接关系到流固界面的传热，从而直接影响到喷管结构的温度场，因而选择合适的壁面函数也是喷管结构温度场分析的一个重要方面。本分析中，将选择增强的壁面处理(Enhanced Wall Treatment)，该壁面函数结合了壁面的二层模型。当边界层网格细密到可满足直接求解粘性边界层时为二层模型，当临近边界层网格较稀疏时，则采用增强的壁面函数求解边界层。

2 复合喷管结构温度场分析

2.1 物理问题

为了给喷管热喷管的耐烧蚀层和背壁厚度设计提供依据，以及将来为喉衬热的应力分析和接缝设计提供输入条件，分析喷管结构的温度场分布。某发动机复合喷管结构如图1所示，收敛段和扩张段的耐烧蚀层为碳/酚醛材料，喉衬为炭/炭材料，喷管背壁为高硅氧/酚醛材料，最外层为钢壳体。发动机燃烧室平均压力为8 MPa，工作时间为35 s。

图1 喷管结构简图Fig.1 Nozzle structural schematic

对于该喷管的结构温度场分析，可忽略一些次要影响因素，对模型进行必要的简化，作以下假设和简化:

(1)燃气简化为理想气体

(2)不考虑喷管各部件之间缝隙中的粘接剂和填充物，各部件无缝隙接触;

(3)忽略壁面的炭化、烧蚀和辐射效应，只考虑燃气与喷管内壁面的对流换热效应;

(4)忽略各部件之间的接触热阻;

(5)真实的高硅氧/酚醛背壁是分段的，在该模型中合并为一段。

2.2 计算模型

针对喷管的轴对称特点，建立喷管的轴对称计算模型，将喷管的内流场区域设为流体域，将喷管的各部件设为固体域，如图2所示。对各区域分别采用四边形网格进行划分，如图3所示。喷管内流场区域与内壁面之间建立耦合边界条件，喷管外壁面边界设为与环境的对流换热边界，对流换热系数设为5 J/m2。喷管的入口设为压力入口边界，压力为燃烧室平均压力8 MPa，燃气温度3 200 K，出口为压力边界，为一个大气压，温度为300 K。材料参数如表1所示。

图2 喷管轴对称计算模型Fig.2 Nozzle axisymmetric model

图3 喷管网格模型Fig.3 Nozzle computational mesh

表1 材料参数Table 1 Parameters of nozzle materials

在复合喷管中，流场的发展速率远大于结构的热传导速率，因而在短时的点火增压过程，可忽略结构的热传导效应，喷管结构的热传导计算就可从发动机进入平稳工作阶段时开始。如此，在喷管结构的热传导计算过程中，喷管流场都是稳定状态。因此，在本文分析中，首先对喷管进行稳态内流场分析;然后，在喷管稳态内流场的结果基础上，对喷管的热传导进行瞬态分析。具体的分析过程如图4所示。

图4 喷管结构分析过程Fig.4 Nozzle fluid structrue analysis process schematic

2.3 计算结果及分析

2.3.1 工作35 s时温度场分布

喷管在35 s时温度场分布云见图5。从图5中可看出，热沿喷管内壁面近似平行向外壁面传递，喷管内壁面附近区域为高温区，逐渐向外壁降低，温度的传递仅通过了由炭/酚醛收敛段、喉衬和炭/酚醛扩张段组成的第一层结构，高硅氧/酚醛背壁还处于300～500 K的温度范围，钢壳体外壁温度为300 K。由此可见，该喷管的结构是满足热防护的要求的，可能还有较大的余量。

从图5中还可明显看出，最高温度位于喷管内壁在从收敛段中部到喉部的一段区间，在2 680～2 860 K之间，见图6。炭/酚醛收敛段与喉衬的接缝就位于该处，因此要注意此处烧蚀和接缝设计。

图5 35 s喷管温度场分布Fig.5 Nozzle temperature distribution at 35 s

根据喷管结构分布特点，分别在收敛段、喉部和扩张段取截面(见图7)，考察温度各部位沿径向截面的分布(见图8～图10)。由曲线图8～图10可看出，温度沿截面从内壁向外壁逐渐降低，在与高硅氧/酚醛背壁的交界面上曲线发生了转折，表明界面两侧材料参数的不同，导致了温度分布的一阶导数不连续，图中横坐标作了归一化处理。高硅氧/酚醛背壁的温度都在300～500 K范围内，其外壁面温度为300 K。显然，钢壳体的温度还属于常温300 K，没有发生温升。收敛段和扩张段截面的温度呈线性下降，而喉部的截面的温度下降逐渐趋缓，为非线性。喉部3个截面的温度分布趋势相同。其中，截面1温度最高，截面2次之，截面3最低。表明喉衬在入口部位是整个喉衬的高温区域，烧蚀可能最为严重。

图6 喷管最高温度区Fig.6 The highest temperature zone of nozzle

图7 提取温度的截面位置Fig.7 The positions where temperature is extracted

图8 35 s收敛段截面温度分布Fig.8 Temperature change with radial direction distance at nozzle convergent cross section at 35 s

图9 35 s喉衬截面上的温度分布Fig.9 Temperature change with radial direction distance at nozzle throat cross section at 35 s

图10 35 s扩张段截面温度分布Fig.10 Temperature change with radial direction distance at nozzle divergent cross section at 35 s

2.3.2 工作初期温度发展情况

考察了5.25 s内喉衬以及与收敛段和扩张段接缝附近的温度变化，见图11。由图11可看出，炭/炭材料喉衬的温度传导速度明显高于两侧炭/酚醛材料的收敛段和扩张段，热传导开始时，接缝两侧温差较大。随时间推移，接缝处喉衬一侧的温度向高硅氧/酚醛的收敛段和扩张段一侧传导，接缝两侧的温差逐渐缩小。

图11 5.25 s内温度变化情况Fig.11 Temperature evolution within 5.25 s

接缝应力是接缝宽度设计的主要依据，其最大值往往发生在工作初期阶段。接缝应力取决于接缝的接触力和接触面积。接触力会随温度由内壁向外传导得越充分而变大，因为收敛段、喉衬和扩张段在温度的作用下，将向轴向膨胀，相互之间发生挤压，温度越高，挤压越严重，接触力也就会越大。另一方面，温度逐渐向外层传导，具有一定宽度的接缝由开始不接触到接触，之后热传导越充分，接触面积会越大，有降低接触应力的趋势。最大接触应力发生在何时有多大与喷管的长度、材料的性质、缝隙宽度、结构的厚度都有关系，在接缝设计时要考虑到这些因素。

2.3.3 喉衬温度分布及变化

喉衬会因自身的热应力产生破坏，因此有必要考察喉衬的温度分布及变化情况。工作35 s结束时，喉衬的温度分布云图如图12所示。整体上，温度依据喉衬弧形内型面由内向外、由高到低的分布，喉衬入口段的温度最高，达到了2 833 K，外壁中部区域温度最低，在600～765 K之间，内外温差较大，到达了2 000 K左右。在这样的温度分布下，喉衬在两个方向上可能存在发生断裂的风险:一是轴向上，内壁温度高，外壁温度低，喉衬产生由内向外的弯曲变形，这使得喉衬有发生轴向断裂的风险;二是周向上，外壁约束了内壁周向膨胀，有发生周向断裂的风险。是否发生破坏，还需要应力的分析。

观察喉衬沿截面1(见图7)由内壁到外壁的温度分布随时间的变化，见图13。

图12 35 s喉衬温度分布Fig.12 Throat temperatue distribution at 35 s

图13 截面1温度随时间变化Fig.13 Temperatue change with time at the first cross section of throat in Fig.7

由图13可知，随着时间的推移，整个截面上的温度在逐渐增高，同时温度在逐步向外壁传播，靠近内壁的区域在开始阶段升温较快，然后趋缓，靠近外壁区域的升温则较为匀速。由图14可知，内壁面温度在5 s内急剧上升，之后趋于平缓。

图14 截面1内壁面温度随时间的变化Fig.14 Inner wall temperature change with time at the first cross section of throat in Fig.7

3 喉衬烧蚀估算

由于喉衬壁面温度较高，在湍流边界层内CO2和H2O与喉衬壁面的碳反应速率较快，大于CO2和H2O的向边界层的扩散速率。因此，炭的质量耗损计算采用扩散控制模型［12］，其炭的质量耗损率表示式为

本计算中，CCO2e=0.049 3，CH2Oe=0.426 7，cv=3 208 J/(kg·K)，=0.027 3 kg/mol，hg从仿真计算的结果中提取出来，其随时间的变化曲线如图15所示。经计算得到mw=3.766 kg/m2，折算为线烧蚀量r=2.1 mm(喉衬密度为 1 800 kg/m3)。

图15 喉径处对流传热系数时间曲线Fig.15 Time curve of throat heat transfer coefficient

4 结论

(1)工作35 s结束时，喷管温度仅传导过第一层结构，喷管的结构满足热防护要求。

(2)温度最高区域位于内壁收敛段中后部，在2 680～2 860 K之间，此处可能是烧蚀最为严重的地方。

(3)喉衬温度由内到外从高到低平行分布，最高温度为2 833 K，位于入口段，内外温差达到2 000 K。

(4)喉衬线烧蚀量估算为2.1 mm。

［1］沈勇.碳/碳固体火箭发动机柔性喷管热及热-结构耦合作用研究［D］.华中科技大学，2004.

［2］王天波，薛谈顺，周长省，等.复合结构喷管温度场及应力场数值模拟［J］.弹道学报，2012，24(2):88-105.

［3］王作良.喷管热防护的有限元数值分析［D］.哈尔滨工程大学，2006.

［4］藤亚莉，黄立威.复合材料喷管热应力分析［C］//第七届中国CAE工程分析技术年会论文集，2011:352-356.

［5］熊永亮，郜冶.喷管温度与应力场的数值研究［J］.哈尔滨工程大学学报，2007，28(8):852-855.

［6］胡江华，孟松鹤，常新龙，等.锥形套式连接C/C喷管扩张段温度场与应力场分析［J］.固体火箭技术，2012，35(1):64-68.

［7］田四朋，唐国金，李道奎，等.固体火箭发动机喷管结构缝隙设计［J］.推进技术，2005，26(5):448-451.

［8］田四朋，唐国金，李道奎，等.固体火箭发动机喷管结构完整性分析［J］.固体火箭技术，2005，28(3):180-183.

［9］胡海峰，鲍福廷，王艺杰，等.喷管分离流流动-热-结构顺序耦合数值模拟及试验研究［J］.宇航学报，2011，32(7):1534-1542.

［10］郑晓亚，陈凤明，蔡飞超.复合喷管结构耦合计算的一种策略［J］.航空动力学报，2011，26(1):223-227.

［11］黄宏艳，王强.轴对称喷管内外流场与结构温度场耦合计算［J］.推进技术，2008，29(2):195-200.

［12］陈汝训.固体火箭发动机设计与研究(下)［M］.北京:宇航出版社，1993.

(编辑:崔贤彬)

Coupled numerical simulation on flow field and thermo-structure of composite nozzle

FU Xue-jin1，2，DENG Heng2，KOU Jun-qiang2，SONG Jun2，CHEN Gang3
(1.Second Artillery Engineering University，Xi＇an 710025，China;2.Xi＇an Modern Control Technology Research Institute，Xi＇an 710065，China;3.Equipment Representive Office，Second Artillery Corps，Stationad in China Academy of Launch Vehicle Technology，Beijing 100076，China)

In order to estimate the heat protection performance and obtain the temperature load for erosion calculation and structural stress analysis，temperature numerical simulation of composite nozzle thermo-structure was implemented by the Fluent CFD software using the fluid-structure coupled method.In the computation，the RNG k-ω turbulence model and the enhanced wall treatment method were selected.The temperature distributions at the end of working time and its evolution in initial period were analyzed，and the feature of nozzle throat temperature field with time-varying was investigated.The throat erosion was also estimated by diffusion control model.The simulation results show that the thermo-structure meets the heat protection requirement.The highest temperature zone is located in the middle and posterior portion of nozzle convergent section.The linear erosion thickness of nozzle throat is about 2.1 mm.The results also offer technical support for further optimizing the design of the nozzle thermo-structure.

coupled fluid-structure;temperature field;composite nozzle;numerical simulation

V438

A

1006-2793(2014)06-0781-06

10.7673/j.issn.1006-2793.2014.06.008

2014-03-21;

2014-06-29。

傅学金(1978—)，男，博士，研究方向为固体火箭发动机结构数值仿真。E-mail:fxej_epgc@sina.com