固体火箭发动机喷管化学烧蚀的动态数值模拟

于勇，张航维

（北京理工大学 宇航学院，北京 100081）

固体火箭发动机通过喷管将高压高温燃气中的热能转换为动能而产生推力[1].因为具有低密度、高比强度、耐热冲击和耐烧蚀等一系列优异性能[2]，碳基材料(石墨和C/C)常用作固体火箭发动机的喷管材料.在发动机工作过程中，喷管长时间处在高温高压的流场环境下，同时受到化学反应烧蚀和凝相粒子的侵蚀作用，引起碳基材料的损失，导致喷管型面衰退[3].其中化学烧蚀是导致喷管烧蚀的重要原因，主要表现为燃气中的H2O、CO2、OH 等氧化性组分和喷管中碳基材料的化学反应[4].喉部区域型面的变化直接影响着发动机性能，因此喉部的烧蚀性能是喷管考核的重要指标[5-6].

在针对化学烧蚀的仿真研究中，BIANCHI 和NASUTI[7]建立了固体火箭发动机碳基喷管的化学烧蚀模型，证明喷管烧蚀率与燃烧室压力成线性关系.目前研究认为喷管化学烧蚀主要受化学动力学和组分扩散控制，并在此基础上发展出了多种控制机制模型.1985 年KOU 和KESWANI[8]率先提出了计算喷管化学烧蚀速率的双控制模型.考虑到双控制模型的局限性，THAKRE 和YANG[9]于2008 年提出了最小机制控制模型，该模型综合考虑了化学反应速率和燃气扩散速率的影响，分析了不同Al 质量分数工况下的烧蚀率，结果与试验吻合较好.张斌等[10]利用多种化学烧蚀机制分析喷管烧蚀过程，发现计算H2O 引起的烧蚀率时, 应综合考虑壁面温度和燃气组分中H2O 浓度；而计算CO2引起的烧蚀率时,可直接选用化学反应速率控制机制，其研究工作同时表明温度升高导致烧蚀率增加主要是通过影响化学反应速率实现.针对燃气中参与反应的氧化性组分，张晓光等[11]的研究结果表明H2O 在喷管化学烧蚀过程中占主要地位.鲁文涛[12]在模拟喷管化学烧蚀的过程中考虑了H2O 和CO 之间的可逆气相反应，发现在预测C/C 喉衬退移规律时，可以忽略气相反应带来的影响.巩伦坤等[13]以半潜入喷管的收敛段为研究对象，分析碳/酚醛材料的烧蚀规律与气流方向的关系，认为当气流与材料表面平行时，气相对流引起的化学烧蚀起主要作用.

此外，国内外的研究利用试验分析了影响喷管烧蚀过程的因素.EVANS 等[14]使用固体火箭发动机模拟装置研究了金属推进剂和非金属推进剂两种情况下石墨喷管的喉部瞬态烧蚀特性，发现采用非金属推进剂测试的喷管喉部表面更粗糙，被认为是由更显著的非均相反应造成的.陈博等[15]进行了燃气发生器条件下C/C 复合材料喷管的烧蚀性能试验，结果表明压强和流速影响了氧化性组分通过边界层的质量传输，氧化性组分浓度和温度的升高均会导致喷管烧蚀率上升.苏君明等[16]对工程应用的5 种C/C 复合材料和2 种石墨喉衬进行了地面点火试验，证明喷管热环境下的燃气氧化性组分浓度和燃烧室压强是影响碳基材料喉衬烧蚀率的主要原因.

在材料烧蚀导致边界移动方面的研究，国内外也取得了一系列成果.CHEN 和MILOS[17]利用TITAN程序和动网格技术成功模拟了流固热耦合作用下热防护材料的型面衰退.MENG 等[18]、ZHANG 等[19]利用Abaqus 中的任意拉格朗日−欧拉(ALE)自适应网格分析了多物理场作用下的表面烧蚀现象.薛海峰等[20]以碳/酚醛燃气舵为研究对象，建立了化学烧蚀模型，并利用弹性光顺法和网格重构法模拟了由热化学烧蚀引起的材料表面退移过程.

综上所述，目前国内外对喷管化学烧蚀的研究多集中在喷管内壁面上的化学反应机理，以流场影响喷管内部传热的单向流固耦合为主，忽略了二者的双向耦合过程，即喷管型面变化和流场变化在固体火箭发动机工作期间的相互作用.本文建立了动态烧蚀耦合模型，综合考虑了传热，化学烧蚀和型面衰退对喷管烧蚀率的影响，实现了喷管型面衰退和喷管流场的双向流固耦合，着重分析了动态烧蚀过程中喷管内壁面和流场的变化情况.

1 数学物理模型

1.1 基本假设

本文研究基于以下假设：

①燃气组分为理想气体；

②不考虑燃气中各组分气体之间的化学反应；

③不考虑凝相颗粒导致的机械侵蚀的影响；

④忽略辐射换热，忽略重力等体积力的影响.

1.2 控制方程

①气相控制方程与湍流模型.

喷管中多组分气体湍流流动的Reynolds 时均方程组通用形式[21]如下：

式中：ϕ为通用变量；Tϕ为输运系数；Sϕ为各方程源项.对于连续、组分、动量和能量方程，上述变量具有特定的形式.

理想气体状态方程为

式中：Yi为某气体组分i的质量分数；Wi为某气体组分i的摩尔质量.

湍流模型采用Realizable k-epsilon 两方程模型，采用增强型壁面函数处理近壁面的流动.

② 固相控制方程.

固相材料瞬时导热方程为

式 中：ρs为 固 相 材 料的 密 度；Cs为 固 相 材 料 的 比 热；Ts为固相材料的温度；λs为固相材料的导热系数.

③ 气−固界面边界条件.

质量平衡：

组分平衡：

能量平衡：

式中：ρg为气体混合物的密度；r˙为固相材料消耗率；Di,m为 某 气 体组 分i的 扩 散 系 数；ω˙i为某 气 体 组 分i的生成或消耗速率；λg为气体混合物的导热系数；qi为某气体组分i的热流密度；下标 w表示气−固交界面.

④ 表面异相反应.

喷管固相材料和燃气在喷管内壁面上发生反应，该反应的反应速率遵循Arrhenius 定律，见式(7).反应中固相材料的消耗率见式(8).

式中：kj为反应速率；Aj为反应指前因子；Tw为喷管内壁面温度；βj为温度指数；Ej为反应活化能 ；R为通用气体常数；pj,w为某氧化性组分j在喷管内壁面处的分压；nj为反应压强指数；ρs为喷管固体材料密度.

根据式(7)利用用户自定义函数UDF 和DEFINE_SR_RATE 宏定义化学反应速率，建立化学反应模型，完成喷管内壁面处的化学反应计算.

⑤化学烧蚀模型.

化学烧蚀模型表征喷管烧蚀率的表现形式，本文采用化学动力学控制模型，见式(9)，该公式表示喷管由于化学反应导致的固相材料消耗率即为喷管烧蚀率.

式中r˙erosion为喷管烧蚀率.

1.3 物理模型和参数

本文以70-lb BATES 发动机[22]为原型，模拟计算其喷管内流场以及喷管内壁面的化学烧蚀情况.发动机几何结构如图1 所示，喷管喉径为50.8 mm，收敛半角为45°，扩张半角为15°，扩张比为9.5，燃烧时间为5 s.喷管材料为石墨，其物性参数 ρs=1 830 kg/m3，Cs=1 050 J/(kg·K)，λs=70 W/(m·K).

图1 70-lb 发动机几何结构Fig.1 70-lb BATES motor configuration

为精确模拟喷管内壁面处的传热传质现象，在网格划分时需保证近壁面处网格的y+<1，近壁面第一层网格高度设定为5×10−7m.

结构网格在保证近壁面处y+值大小的同时，其动网格更新往往会导致远场处的网格变形较大，难以满足后续的计算需求.为适应网格变形并减少网格数量，本文采用混合网格进行计算，即在流体区域近喷管内壁面处采用结构网格，其余流体区域和固体区域采用非结构网格，如图2 所示.

图2 喷管计算网格Fig.2 Computational grid for nozzle

动网格更新方式为弹性光顺法和网格重构法，利用用户自定义函数UDF 和DEFINE_GRID_MOTION宏定义动网格，控制喷管内壁面网格节点的移动，节点的移动量由计算得到的烧蚀量给出.

图3 为采用动网格和不采用动网格两种计算方式下15% Al 质量分数喷管内壁面处的y+值分布对比，可以看出混合网格划分方案在保证网格运动的同时基本满足了y+值的需求.

图3 喷管壁面y+Fig.3 y+ value of nozzle wall

喷管入口设为压力入口、喷管出口设为压力出口，喷管内壁面设为传热耦合壁面及表面化学反应发生面，喷管侧壁、外壁设为绝热，喷管固体域初始温度设定为300 K.

本文化学烧蚀模型主要考虑C 与氧化性组分H2O、CO2的反应，其反应方程式及动力学参数见表1[23].

表1 表面反应动力学参数Tab.1 Kinetic parameters of chemical reaction

1.4 耦合机理和方法

在喷管的化学烧蚀过程中，喷管固体域中的石墨和燃气发生化学反应，喷管型面发生衰退，该过程中流场和喷管型面产生双向流固耦合作用.耦合计算流程如图4 所示:

图4 耦合计算流程Fig.4 Couping calculation scheme

①在初始时刻，给定流场计算的边界条件；

②求解流场区域，获得喷管内壁面附近的温度、组分浓度和组分分压等参数；

③根据步骤②中的参数计算化学反应速率和烧蚀率，得到喷管内壁面上的烧蚀率分布，当前计算时间为t，计算时间步长为 ∆t；

④利用动网格技术，将喷管内壁面上的网格节点以对应的烧蚀率进行移动；

⑤重复步骤②～④，直至完成所有计算.

2 计算结果及分析

2.1 Al 质量分数对烧蚀的影响

表2 为不同Al 质量分数时的喷管入口条件，包括燃气压强、温度及组分质量分数[9].图5 为计算得到的t=5 s 时不同Al 质量分数工况下喷管内壁面烧蚀率沿轴向的分布结果，其中喉部的烧蚀率与试验值[22]吻合较好，具体结果见表3.

表2 不同Al 质量分数工况下喷管入口条件Tab.2 Nozzle inlet conditions with different aluminum mass fraction

表3 不同Al 质量分数工况下喉部烧蚀率与试验值对比Tab.3 Comparison between calculated and measured throat erosion rates

图5 不同Al 质量分数工况下喷管壁面烧蚀率分布Fig.5 Erosion rate of under different aluminum mass fraction

结果表明，喷管烧蚀率随推进剂中Al 质量分数的增加而降低，其主要原因是燃气中参与化学反应的氧化性组分H2O 和CO2浓度的降低.喷管烧蚀主要分为化学烧蚀和颗粒侵蚀两部分，推进剂中Al 质量分数的增加导致燃气中Al2O3颗粒浓度增加，进而增强了Al2O3颗粒对喷管内壁面的机械侵蚀强度，而实际的喷管烧蚀率却不升反降，由此证明了喷管烧蚀中颗粒侵蚀所占比重较小，由此证明了化学烧蚀才是导致喷管烧蚀的主要原因[9].

图6 为15%的Al 质量分数工况下烧蚀率随时间变化的分布情况.随着烧蚀过程的进行，烧蚀率的上升变慢，这主要是由于喷管内壁面温度上升，使得喷管内壁面和燃气之间的热交换趋于平衡，导致燃气向固相材料的对流换热减弱.

图6 不同时刻下15%Al 质量分数的喷管烧蚀率分布Fig.6 Erosion rate with 15% aluminum mass fraction at different times

图7 和图8 分别显示了15%的Al 质量分数工况下t=5 s 时的喷管压力分布和速度分布.图9 对比了15%的Al 质量分数工况下不同时刻的喷管温度分布，显示了喷管内部的传热过程.

图7 15%Al 质量分数的喷管压力分布 (单位：MPa)Fig.7 Pressure distribution of nozzle with 15% aluminum mass fraction(unit: MPa)

图8 15%Al 质量分数的喷管马赫数分布Fig.8 Mach number distribution of nozzle with 15% aluminum mass fraction

图9 15%Al 质量分数的喷管温度分布(单位：K)Fig.9 Temperature distribution of nozzle with 15% aluminum mass fraction (unit: K)

2.2 动态烧蚀过程

在喷管工作过程中，由于内壁面发生化学烧蚀，造成了固体域材料损失和喷管型面变化，研究表明喷管喉部面积增加5%就会对发动机的性能造成巨大影响[24]，因此想要深入了解喷管的化学烧蚀过程，需结合喷管型面衰退进一步分析.非定常工况下喷管受到化学烧蚀导致流场和内壁面相互耦合作用，在该作用下流场和内壁面发生变化的过程即为喷管的动态烧蚀过程.本文通过耦合化学烧蚀模型和动网格模型，实现了喷管的动态烧蚀.

2.2.1 喷管烧蚀率对比

选用表2 中的工况分别模拟喷管的动态烧蚀过程.图10～12 分别是15% Al、18% Al、21% Al 质量分数工况下t=5 s 时动态烧蚀对喷管内壁面烧蚀率的影响，并与未采用动网格的计算结果进行对比，具体结果见表4.图13 是计算得到的t=5 s 时不同Al 质量分数工况下动态烧蚀喷管内壁面烧蚀率沿轴向的分布结果.

表4 动态烧蚀工况下喉部烧蚀率与试验值对比Tab.4 Comparison between calculated and measured throat erosion rates under dynamic erosion

图10 15% Al 质量分数的喷管烧蚀率分布对比Fig.10 Comparison of erosion rate with 15% aluminum mass fraction

图11 18% Al 质量分数的喷管烧蚀率分布对比Fig.11 Comparison of erosion rate with 18% aluminum mass fraction

图12 21% Al 质量分数的喷管烧蚀率分布对比Fig.12 Comparison of erosion rate with 21% aluminum mass fraction

图13 不同Al 质量分数工况下动态烧蚀喷管壁面烧蚀率分布Fig.13 Dynamic erosion rate of under different aluminum mass fraction

可以看出，3 种工况下采用动网格的喉部烧蚀率都更贴近试验值，验证了本文动态烧蚀计算方法的有效性.同时3 种工况下采用动网格的喷管烧蚀率峰值都较高.图14 为15%Al 质量分数工况下2 种计算方式的喷管内壁面温度分布.对比不采用动网格的计算方式，可以看出采用动网格的喷管由于喷管型面和流场的变化，导致喷管喉部及喉部上游区域温度较高，两种计算方式的温度差最高可达40 K，温度的升高增强了化学反应速率，导致了较高的峰值烧蚀率.

图14 15% Al 质量分数的喷管内壁面温度分布对比Fig.14 Comparison of temperature distribution of nozzle with 15% aluminum mass fraction

2.2.2 喷管型面变化和温度分布对比

图15(a)和图15(b)分别为截取喉部上游区域和下游区域的喷管内壁面型面变化趋势.由于烧蚀率在喷管上的分布，使得喷管上游型面变化最大.可以看出t=1 s 的初始时刻喷管型面变化较小，随后逐渐加大.这是由于初始时刻喷管内壁面温度较低，化学反应速率和烧蚀率较小，随着喷管内壁面温度的上升，化学反应速率加快，进而增大了单位时间内的烧蚀量.喷管内壁面温度变化见图16.

图15 15% Al 质量分数的喷管喉部上游型面及下游型面变化Fig.15 Variation of upstream and downstream profile of throat with 15%aluminum mass fraction

图16 不同时刻下15% Al 质量分数的喷管内壁面温度分布Fig.16 Temperature distribution of nozzle with 15% aluminum mass fraction at different times

2.2.3 喷管压力变化

图17 和图18 分别是15% Al 质量分数工况下沿喷管内壁面和沿喉部径向的压力分布.两种计算方式下喷管收敛段的压力差异较小，而相比于未采用动网格的计算方式，动态烧蚀的喷管在喉部及下游区域的压力更高，沿喉部径向的最大压力差约为0.3 MPa.这种压力变化的差异可以根据图10 的烧蚀率分布分析得出，图10 中15% Al 质量分数工况下两种计算方式在收敛段的烧蚀率曲线较为吻合，采用动网格的喷管烧蚀率在喉部及下游区域较高.喷管型面衰退和流场变化的耦合作用导致了两种计算方式下烧蚀率的差异，进而产生了不同的压力分布.

图17 15% Al 质量分数的喷管内壁面压力分布对比Fig.17 Comparison of pressure distribution of nozzle with 15% aluminum mass fraction

图18 15% Al 质量分数的喷管喉部径向压力分布对比Fig.18 Radial distribution of pressure at the nozzle throat with 15% aluminum mass fraction

2.2.4 喷管喉部组分变化

图14 的温度分布表明，考虑型面衰退影响的喷管在喉部的温度更高，导致喉部化学反应速率加快，同时在喉部壁面附近形成了组分浓度梯度，组分的浓度分布同时受到化学反应速率和扩散速率的影响.图19 是15% Al 质量分数工况下沿喷管喉部径向的组分分布，包括参与反应的H2O 和CO2.可以看出，动态烧蚀的喷管喉部壁面附近H2O 和CO2的浓度偏低，质量分数分别为0.001 7 和0.003 3，而不采用动网格计算得到的H2O 和CO2的质量分数分别为0.001 6 和0.002 5.这是由于高温状态下H2O 和CO2的扩散速率大于化学反应的消耗速率，2 种反应物的浓度分布主要受扩散速率的影响.

图19 15% Al 质量分数的喷管喉部径向组分分布对比Fig.19 Radial distribution of species at the nozzle throat with 15% aluminum mass fraction

3 结 论

通过本文的研究，可以得出以下结论：

①随着推进剂中Al 质量分数的增加，喷管燃气中参与化学反应的氧化性组分H2O 和CO2的浓度减小，喷管烧蚀率随之降低；

②建立的动态烧蚀模型综合考虑了传热、化学反应和喷管衰退等现象.相比于未采用动网格的计算方式，采用动态烧蚀模型得到的计算结果更贴近试验值，能够有效模拟喷管的动态化学烧蚀过程；

③喷管内壁面的衰退使得喷管流场发生变化，导致采用动网格技术计算的喷管峰值烧蚀率更高，喉部附近的温度和压力更大.喉部附近的温度差最高可达40 K，沿喉部径向的最大压力差约为0.3 MPa；

④受高温状态下化学反应速率和扩散速率的影响，采用动态烧蚀模型计算的喷管在喉部壁面附近的H2O 和CO2的浓度更高.