燃气发生器绝热层烧蚀数值仿真*

周 源,齐 强,陈志刚,徐 明

(1 海军航空工程学院,山东烟台 264001;2 91515部队,海南三亚 572000)

燃气发生器绝热层烧蚀数值仿真*

周 源1,齐 强1,陈志刚1,徐 明2

(1 海军航空工程学院,山东烟台 264001;2 91515部队,海南三亚 572000)

为了研究燃气发生器燃烧室绝热层的传热烧蚀过程,建立了热化学烧蚀计算模型和绝热层传热计算模型。考虑到烧蚀过程中绝热层边界的移动,根据预测-校正格式对模型进行离散,并采用拟牛顿法对绝热层的烧蚀过程及传热过程进行了耦合计算。计算得到了绝热层表面温度、烧蚀率和烧蚀厚度等参数的变化规律。计算结果表明,预测-校正格式和拟牛顿法可以用于求解时动边界的瞬态传热模型。

绝热层;热烧蚀;时动边界

0 引言

燃气发生器是一种小型固体火箭发动机,是弹射发射系统的工质源和能源。药柱在燃烧室内点火燃烧,产生高温高压燃气使得燃烧室面临极其恶劣的内部热环境。因此,研究燃烧室内壁绝热层的传热烧蚀过程对确保燃气发生器安全工作具有重要意义。

对绝热层的传热烧蚀问题已有很多的研究,文献[1]和文献[2]对热解型绝热材料的烧蚀过程建立了物理结构模型,文献[3]和文献[4]采用有限差分法对绝热层的瞬态传热进行计算,文献[5]对于时动边界上的热传导问题的求解提出了人工边界的方法。

文中对燃烧室绝热层的烧蚀传热过程建立了碳化层—热解面—基体层的结构模型,采用预测-校正格式对时动边界模型进行离散求解,为绝热层的烧蚀传热研究提供了理论依据。

1 计算模型

1.1 热化学烧蚀计算模型

在燃气发生器内,燃气中参加烧蚀反应的气体组分主要是CO2、H2O和H2,将其余不参加烧蚀反应的惰性气体组分折合成N2,如表1所示。那么,热化学反应烧蚀率的计算采用以下3个化学反应方程:

表1 燃气中的主要反应气体组分

根据质量守恒定理,烧蚀反应气体组分来自于燃气主流扩散和绝热层热解气体,在烧蚀表面发生反应后被燃气主流带走。热化学烧蚀程序方程组[6]表示如下:

(1+F)K6,w=K6,g+Bf6

1.2 传热计算模型

绝热层在工作过程中形成碳化层—热解面—原始材料层的结构,如图1所示。将绝热层内部的导热过程简化为一维大平板的瞬态导热[1],不考虑绝热材料的热膨胀及其引起的热应力,且材料参数随温度的变化忽略不计[7],由此建立传热计算模型。

1)碳化层内的瞬态导热

(1)

2)热解面上能量守恒

(2)

式中:ΔHp表示绝热材料的热解潜热;λ1表示原始材料层的导热系数;Tp表示绝热材料的热解温度。

图1 传热计算模型

3)原始材料层内的瞬态导热

(3)

式中:ρ1表示原始材料层的密度;c1表示原始材料层的比热;λ1表示原始材料层的导热系数。

4)原始材料层与壳体层之间的界面上能量守恒

(4)

式中,λ2表示壳体层的导热系数。

5)壳体层内的瞬态导热

(5)

式中:ρ2表示壳体层的密度;c2表示壳体层的比热;λ2表示壳体层的导热系数。

6)边界条件

①根据假设,燃气发生器工作期间与外界是绝热的。壳体层与外部环境之间的界面上能量守恒,有

(6)

②当Tw

(7)

Qin=Qrad+Qcon

(8)

③当Tw≥Tp时,绝热层表面出现碳化层,根据边界条件方程,有

(9)

(10)

7)初始条件

绝热层内部初始温度是常数,取环境温度。

1.3 数值离散方法

在烧蚀过程中绝热层的厚度不断减少,对于具有时动边界的瞬态传热问题,采用预测-校正格式[8]进行离散求解。如图1所示,对每一层结构进行等分离散,则有:

(11)

通过与其它数值离散格式对比,该格式具有二阶精度且无条件稳定。

2 计算结果与分析

在计算过程中,使用变空间步长的差分方法,节点坐标不断更新。在每个时间步长上,根据烧蚀速率来确定移动边界节点的位置,然后重新进行离散,再计算温度场。

考虑到计算时间步长较小(取0.001 s),i时刻和i+1时刻初值变化很小,将i时刻程序的解作为i+1时刻程序初值[6],再采用拟牛顿法求解,以缩短计算时间。

图2给出了绝热层单位面积上热解气体质量流率随时间的变化规律。在燃气发生器开始工作的极短时间内,原始材料层不发生热解,表面未形成碳化层。随着表面温度的升高,绝热材料开始热解并产生极薄的碳化层。这时,通过热解面传入绝热材料内部的热流很大,热解速率较大。在0.009 s时,热解气体质量流率曲线出现了一个很大的峰值,热解气体质量流率达到6.697 1 kg/(m2·s);随着绝热材料的不断热解,碳化层厚度增加并逐渐稳定,热解速率减小,最后逐渐趋于定值。热解气体平均质量流率为0.519 4 kg/(m2·s)。

图2 热解气体质量流率

图3给出了绝热层内表面的温度变化规律。在开始阶段,温度迅速升高,温度越高,变化速率越小。在达到约2 500 K时,温度趋于平衡。

图3 绝热层内表面温度

图4给出了燃气发生器工作过程中烧蚀线与热解线的相对位置。在开始阶段,绝热层不发生热解和碳化。直到绝热层温度达到热解温度,碳化层厚度迅速增大,并随着温度的升高逐渐减缓。从图中可以看出,热解速率明显比烧蚀速率快很多。

图4 烧蚀线和热解线的相对位置

图5给出了绝热层内表面碳化层的线烧蚀率随时间变化曲线。在绝热层表面未形成碳化层的时候,烧蚀率为零。在绝热层内表面温度的升高,以及燃烧室压力的增加,根据Arrhenius定律,碳化层烧蚀率逐渐增大。在后效段,燃烧室压力急剧下降,线烧蚀率也急剧下降。根据曲线求得平均线烧蚀率为0.069 3 mm/s。

图5 碳化层表面线烧蚀率

3 结论

1)文中建立的热化学烧蚀模型和传热计算模型以及所采用的离散格式和计算方法实现了烧蚀变形与传热的双向耦合。

2)预测-校正格式可以用于求解时动边界的瞬态传热问题。

[1] 徐善玮, 侯晓, 张宏安. 固体火箭发动机内绝热层烧蚀质量损失计算 [J]. 固体火箭技术, 2003, 26(3): 28-31.

[2] 张涛, 孙冰. 三维烧蚀内部热响应数值计算研究 [J]. 宇航学报, 2012, 33(3): 298-304.

[3] 李宛珊, 王文洽. 二维热传导方程的有限差分区域分解算法 [J]. 山东大学学报: 理学版, 2011, 46(12): 1-5.

[4] 冯立伟. 热传导方程几种差分格式的MATLAB数值解法比较 [J]. 沈阳化工大学学报, 2011, 25(2): 179-182.

[5] 张雪艳, 周爱霞. 时动边界上热传导问题的求解方法 [J]. 廊坊师范学院学报: 自然科学版, 2010, 10(2): 9-12.

[6] 张斌, 刘宇, 王长辉，等. 长时间工作固体火箭发动机燃烧室热防护层烧蚀计算 [J]. 固体火箭技术, 2011, 34(2): 189-192.

[7] 杨春杰. 固体火箭发动机后效冲量研究 [D]. 长沙: 国防科技大学, 2011.

[8] 南京大学数学系计算数学专业. 偏微分方程数值解法 [M]. 北京: 科学出版社, 1979: 145-182.

Numerical Simulation for Ablation of Insulator in Gas-generator

ZHOU Yuan1, QI Qiang1, CHEN Zhigang1, XU Ming2

(1 Naval Aeronautical and Astronautical University, Shandong Yantai 264001, China; 2 No.91515 Unit, Hainan Sanya 572000, China)

In order to study thermal behavior and ablation of insulator in the combustor of gas-generator, the thermo-chemistry ablation model and heat transfer calculation model were established. Considering the fact that the boundary of insulator is moving during the ablation, the models were discredited with predicted-correction method. Coupling calculations were carried out for insulator ablation and temperature field with quasi-Newton method. Temperature of inner wall, ablation velocity and ablation thickness were calculated. The results show that the predicted-correction method and the quasi-Newton method can be applied in unsteady heat transfer and ablation problem with moving boundary.

insulator; ablation; moving boundary

2014-03-14

国家自然科学基金(51005242)资助

周源(1979-),男,安徽濉溪人,讲师,博士研究生,研究方向:兵器发射理论与技术。

V435

A