不同压力工况下AP-HTPB推进剂微尺度燃烧的数值模拟

2019-04-11陈王琦余永刚

兵器装备工程学报 2019年3期

陈王琦，余永刚

(南京理工大学能源动力工程学院，南京 210094)

复合固体推进剂凭借其良好的力学特性、可加工性和储存性等诸多优点，广泛运用于火箭弹、导弹、和航天器发射发动机之中。研究固体推进剂燃烧过程以及其燃烧速率的相关影响因素一直是固体推进剂研究中的一大重点，这是因为其燃烧过程不仅受固体推进剂本身的影响，例如：各组分之比、是否使用催化剂以及固体推进剂颗粒大小等等，而且外界环境因素诸如燃烧初始点火温度、压力和辐射传热等等也影响其燃烧过程。高氯酸铵/端羟基聚丁二烯(AP/HTPB)是众多固体推进剂中最早运用于固体火箭发动机之中，目前依然在大量使用。Fitzgerald R P等[1]对AP/HTPB固体推进剂燃烧特性分别进行了实验研究和数值仿真，通过紫外线发射成像方法，获得在2到55个大气压下的燃烧火焰结构与燃烧表面的燃速。GaduPArthi T等[2]进一步采用72步化学基元反应和引入39种中间反应产物，对AP/HTPB复合固体推进剂燃烧过程中的气相火焰结构进行了实验测试与数值模拟。Price E W等[3]最先通过扫描电镜(SEM)对AP/HTPB的微观结构进行研究，发现其排列分布具有明显的周期性，将其简化为周期性单元三明治结构，这也是本文所采用模型的依据。压力是影响固体推进剂燃烧特性的一个主要因素，Taehokim I C等[4]设计了一套实验装置，并研究各种固体推进剂处于低压环境下的燃烧特性。研究表明，相比于更高的压力环境时，低压环境下燃烧速度有一定幅度的增加，并且氧化剂的含量变化在一定程度下又会影响燃烧压力。Ramakrishna P A等[5]也对AP/HTPB三明治模型进行了几个工况的超低压模拟，在不同的压力工况下其燃烧火焰结构明显不同，随着压力的不断增加，总体火焰结构经历了预混火焰、预混扩散火焰混合以及扩散火焰三种火焰结构。窦燕蒙、罗运军等[6]也对储氢合金AP/HTPB推进剂燃烧性能进行了研究，分析了在添入催化剂之后AP/HTPB燃速的变化。对于三明治燃烧模型，Jackson T L[7]进行了不同模型的调研并且总结了模型微观结构特征，徐冠宇等[8]考虑到传统的三明治模型不曾考虑相邻火焰的影响，采用了多层三明治模型结构进行数值仿真，孙迪等[9]对基于BDP模型的三明治模型的燃速参数敏感性进行了研究，分析了在此模型的基础上指前因子、活化能等参数对于仿真结果的影响。

为了进一步研究压力对于AP/HTPB复合固体推进剂的燃烧特性和火焰结构的影响，本文考虑了气固两相的反应，应用简化的2步总包反应动力学机理，建立了AP/HTPB二维三明治稳态燃烧模型，通过fluent软件数值模拟了AP/HTPB固体推进剂在0.3MPa到7.2MPa不同压力工况下的燃烧特性。

1 物理模型

根据AP/HTPB二维周期性三明治结构特点，采用如下基本假设：

1) 不考虑氧化剂AP和粘合剂HTPB的相互作用，即两种组分为独立单元，各自具有不同的且为定常的热物理参数；

2) 气相反应中各气体均为不可压缩理想气体，在整个气相的微观计算空间之内压力均匀分布，且热物理属性是关于温度的已知函数。

3) 固相热分解反应仅仅发生在气固耦合燃面的一层薄面内，固相内部仅考虑热传导换热；

4) 对于固相热分解反应的化学描述采用Arrhenius定律，而对于气相燃烧和火焰结构方面的假设则是基于BDP多火焰模型的包含四种组分两步化学反应的总包反应；

5) 不考虑气相反应生成的固相产物，忽略高温气体的热辐射作用；

6) 采用源项法描述固相热分解和气相扩散燃烧过程。

2 数学模型

基于以上物理模型，建立AP/HTPB固体推进剂二维稳态流动燃烧的基本控制方程：

2.1 气相控制方程

气相控制方程包括质量、动量、组分和能量守恒四个方程；

1) 质量守恒方程

(1)

2) 动量守恒方程

(2)

3) 组分守恒方程

(3)

4) 能量守恒方程

▽T=▽·(λg▽T)+SE

(4)

5) 状态方程

(5)

因为整个计算过程中进行的是稳态模拟，所以压力P看作常数，假设Pr数为常数,取0.79，Le数为1。气相导热系数随温度而变化，采用下面的关系式[10]:

λg=1.08×10-4T+0.013 3 W·m-1·K-1

(6)

2.2 固相控制方程

(λc,i/cc,iρc,i)▽2T+ST,i=0

(7)

式中:ST,i为固相能量源项，λc,i、ρc,i和cc,i分别为固相导热系数、密度和比热，而对于λc,i、ρc,i和cc,i的取值采用分段的常数形式，如果计算点位于AP颗粒内则取AP颗粒的导热系数、密度和比热，若是计算点落在HTPB内则使用HTPB的参数。如式(8)所示

(8)

式中:L为三明治模型的半宽度，α为AP的体积分数,x为横坐标。

2.3 气固耦合关系

采用燃面质量通量平衡方程、燃面温度连续性方程和燃面能量通量平衡方程描述气固耦合关系，对于燃面两侧的质量流动应该满足燃面质量通量平衡方程，即：

(9)

(10)

式中:AAP、AHTPB分别为分解速率常数；EAP、EHTPB为分解活化能；TAP,s、THTPB,s为AP和HTPB燃烧表面温度。

燃面两侧温度满足连续性方程，即：

Tg=Tc

(11)

燃面能量通量平衡方程，即：

(12)

2.4 化学动力学机理

基于BDP多火焰模型，AP/HTPB推进剂燃烧采用包含四种组分两步化学反应的总包反应[11]，即：

(13)

(14)

(15)

对于上述模型的化学反应速率计算，在此采用简化的Arrhenius定律：

(16)

(17)

式中:R1、R2为热分解反应速率和燃烧反应速率，C1、C2为化学反应速率常数，P为压力工况，n1、n2为压力指数，E1、E2为活化能，Ru为气体常数。

2.5 计算模型及初边条件

AP/HTPB复合推进剂燃烧模型如图1所示，整体计算区间2H=1 000 μm，2L=180 μm,Y=0处为燃面，Y>0处为气相燃烧区域，Y<0处为AP/HTPB固体燃料颗粒，燃面下方的|x|≤αL区域为粘结剂HTPB，αL<|x|≤L区域为氧化剂AP颗粒，AP体积分数为76%，AP颗粒粒径为110 μm，X=0 μm及X=90 μm的左右边界为对称边界，质量当量比β=6.3，Y=500 μm处为压力出口边界，Y=-500 μm处为绝热壁面。

图1 AP/HTPB推进剂燃烧模型示意图

网格采用矩形的结构性网格，网格总数为 101 199，并进行了网格无关性验证。对于网格的无关性验证方面主要分为每次设置全局网格尺寸为之前的1.4倍，即计算了在压力为7.2 MPa下10万、7.1万、5.1万和3.6万四种网格总数的工况，其计算结果温度最大值和燃气燃速最大值分别如表1、表2所示。

表1 不同网格总数时温度最大值

表2 不同网格总数时燃气燃速最大值

从表1、表2可以看出，对于温度最大值和燃气燃速最大值两个监测值的差值在4%左右，可以认为计算结果与网格疏密程度近似无关。

计算压力工况变化范围在0.3 MPa到7.2 MPa之间，气相和固相初始条件如式(18)、式(19)

P=Pi,Tg=Tg,i

(18)

Tc=Tc,i

(19)

而对于边界条件主要考虑气相远场、固相远场、燃面组分分布和对称边界四个部分，分别满足以下方程：

气相远场

(20)

固相远场

Tc=Tc,∞

(21)

燃面组分分布

(22)

(23)

对称边界

(24)

3 数值计算结果及分析

利用Fluent软件的有限体积法，采用二阶迎风格式，压力速度耦合问题采用SIMPLE(半隐式压力连接方程)求解。取AP体积分数为76%，对于0.3 MPa到7.2 MPa的不同压力工况下的燃烧进行仿真模拟，计算所需的参数如表3、表4所示。

表3 气相反应区参数[12-13]

表4 固相反应区参数[12-13]

3.1 气相温度分布

图2为AP/HTPB在不同压力工况下的气相温度分布云图，可以明显看出，随着压力的不断增加，火焰结构的形状发生了很大的变化，从最开始的AP热分解预混火焰结构模型渐变为扩散火焰模型。这是因为燃烧过程中的化学反应速率与压力成指数变化关系，压力的微小增加，都会导致化学反应速率的急剧增大。在压力为0.3～0.4 MPa时，火焰结构呈现为预混火焰，存在明显的温度分层现象。这是因为低压情况下AP的热分解反应较慢，所以热解的气相产物可以充分的和HTPB气相产物预混好，从而整体上以预混火焰为主。当压力大于0.4 MPa而小于2.4 MPa时，化学反应速率骤增，但与扩散混合作用还处于同一量级，所以整体上体现为预混火焰与扩散火焰并存，火焰形状呈现出了扁平的“W”特性，而且“W”形火焰突起部分中间高于两侧，由式(10)可知这是因为此时HTPB的燃烧速率明显高于AP燃速，导致在AP和HTPB接触面上形成了从AP面上向HTPB面上突的构造。而当压力大于2.4 MPa时，主要以扩散火焰结构为主，这时AP的燃烧速率要高于粘结剂HTPB，在两者的接触面上，HTPB面会有明显的向AP面上延伸的构造，所以火焰形状呈现出“瘦长的W”形状。并且在此初始压力工况下形成了长长的扩散反应带，温度的最大值点也出现在这扩散反应带区域。这是因为在压力增加的过程中热释放核心在往AP/HTPB重叠界面处移动。Chorpening等[14]对于燃烧过程中火焰结构跟随压力的变化进行了实验研究，发现压力在在较低的情况下，预混火焰结构模型占主导地位，随着压力不断升高，火焰呈现扩散结构，说明我们的仿真结果与实际结果是吻合的。

图2 不同压力工况下的温度分布云图

图3为0.3～7.2 MPa的压力工况下，燃面(y=0 μm)温度变化曲线，从图3中可以看出，在相同的压力工况下，燃面温度由边缘向中心处温度先缓慢升高，然后激增，并且在AP/HTPB交界面处达到温度极值，最后又下降。根据不同压力工况下相同位置的温度变化可知，随着压力的增加，燃面温度也持续增大。由图3还可知，AP燃面温度相比于HTPB燃面温度受压力影响要大的多。这是因为AP/HTPB不同的物性参数及燃面热通量的不同所导致的。

图3 不同压力工况下的燃面温度变化曲线

3.2 燃面速度分布

针对上述得出的燃面在不同压力工况下的温度分布，结合式(10)和表2可以进一步得出燃面燃速分布，如图4所示。

图4 不同压力工况下的燃面燃速变化曲线

从图4可以看出，整个燃面处的燃速随着压力的升高而增加，燃速的最大值出现在AP/HTPB的交界面附近。在整个0.3～7.2 MPa压力范围内，HTPB端的燃速变化很小，而在压力为0.3～0.4 MPa时，AP/HTPB两端的燃速对压力的敏感度不大，并且HTPB端燃烧速率与AP端燃烧速率相差不大，而当压力大于0.4 MPa时，AP.端燃烧速率随着压力的增加急剧增加，这与文献[3]的实验结果也十分吻合，说明模拟结果与实际结果是吻合的。

图5给出了AP体积分数为76%，AP颗粒直径为110 μm情况下数值模拟的平均燃面燃速与实验结果[14]：

由图5可知，两者吻合较好。同时，计算所得的燃速-压力指数为0.487，这也与文献[14]所给出的燃烧压力指数在0.4～0.6相吻合，这也验证了本文的AP/HTPB复合固体推进剂燃烧模型的有效性。

3.3 燃气速度分布

图6为不同压力下气相速度分布云图，图7为最大气相速度变化曲线。

图5 平均燃面燃烧与实验结果

图6 不同压力工况下的气相速度分布云图

图7 不同压力工况下的最大气相速度变化曲线

由图6、图7可以得出，在0.3～0.5 MPa时，燃气最大燃速随着压力增加而增加，而当压力大于0.5 MPa时，燃气最大流速随着压力增加而减小。这是因为在小于0.5 MPa时，随着压力增加，燃面温度和燃面速度也都在增加，根据式(9)可知，燃气流速也会增加。但当压力大于0.5 MPa之后，燃气密度也随之增大，所以导致燃气流速降低。