铝冰发动机内流场的数值计算①

刘平安，王 良，王 璐，王 革

(哈尔滨工程大学 航天与建筑工程学院，哈尔滨 150001)

铝冰发动机内流场的数值计算①

刘平安，王 良，王 璐，王 革

(哈尔滨工程大学 航天与建筑工程学院，哈尔滨 150001)

为了使用数值模拟的方法计算铝冰发动机的性能，用颗粒表面反应模型和气相反应模型模拟铝颗粒在铝冰发动机燃烧室中与水蒸气的燃烧过程，用欧拉-拉格朗日方法计算颗粒沿轨迹的参数，分析了数值模拟的结果，并进行了相同尺寸的铝冰发动机实验，把数值模拟结果与实验结果进行了比较。数值计算得到的燃烧室稳态工作压强约为9.38 MPa，与实验结果接近，燃烧室平均温度为2950.65 K，相比热力计算得到的推进剂燃烧温度略低。通过对铝冰发动机的内流场数值计算，得到了与实验相符合的结果，验证了数值计算模型的有效性。

铝冰发动机；内流场；数值计算

0 引言

近年来低温固体推进剂成为新的研究方向，而金属燃料特别是纳米金属燃料在推进剂中的应用也成为重要研究方向[1]。纳米金属燃料点火温度低、燃烧效率高，在推进剂中作为添加剂时能大大提高其燃速[2]，其超大的比表面积又使其能充当粘结剂以提高推进剂的力学性能，最终可作为推进剂基体和主要成分而不仅是添加剂。包覆后的纳米金属燃料用于低温固化的固体推进剂中能长期保持活性，燃烧产物无污染，有望被用于深空探测，目前已被实验验证[3]。铝冰发动机是低温固体推进剂和包覆纳米金属燃料结合的产物，配方中只有纳米铝颗粒和水，铝在推进剂中充当基体、燃料和粘合剂。实验表明，纳米铝水混合物的燃烧效率很高[4]，其理论性能参数相比某些复合推进剂要高很多[5]，虽然国外铝冰发动机(铝水比为0.71∶1)实验的性能效率很低[3]，但国内的铝冰发动机(铝水比为1∶1)实验结果很好，这证明了铝冰发动机的实际可行性。在铝冰发动机实验中发现按照传统理论(热力计算和零维内弹道理论)设计的燃烧室压强远高于实验压强，这使发动机实验屡屡失败，浪费人力物力成本。燃烧室压强理论与实验的误差主要是由推进剂燃烧和两相不平衡效应造成的，因此对铝冰发动机性能的预估，如燃烧室压强的理论估算需要采用考虑燃烧和两相不平衡效应的新方法。

与常规复合推进剂发动机相比，铝冰发动机主要不同在于推进剂中有大量的金属铝颗粒。铝颗粒在发动机中燃烧和运动时要经历相变、团聚、脱离燃面、点火、气相反应、表面反应、氧化铝凝结、颗粒成长、颗粒碰撞团聚和颗粒破碎等过程[6]。研究铝颗粒燃烧过程的理论模型主要有半经验模型、化学机理模型和CFD模型[7]。其中CFD模型可同时计算两相流动和燃烧过程，得到整个内流场的流动情况，用于发动机工作过程模拟。

为了计算铝冰发动机的理论性能，验证计算模型并得到铝冰发动机内流场参数，本文用CFD模型对铝冰发动机的燃烧和流动过程进行计算，并把得到理论燃烧室压强与实验比较，以验证数值模型在计算发动机性能(如燃烧室压强)时的可行性。

1 数值计算方法

1.1 物理模型

铝冰发动机装药为锥柱形，内部结构如图1所示。发动机的外径160 mm，内径100 mm，喷管喉径18 mm，装药锥角为19°，在工作初始时刻(燃层厚w=0 mm)，装药最大长度为300 mm。计算假定初始时刻燃面全部被点燃，工作过程中各处燃速相同，因此在不同燃层厚时燃面位置如图1所示燃层厚w=0、10、20 mm。

1.2 计算(CFD)模型

本文用Fluent软件中的欧拉拉格朗日模型和可实现k-ε模型来模拟铝冰发动机内流场中的湍流两相流，用组分输运和涡耗散概念(EDC)模型模拟气相混合物流动以及铝和水蒸气的气相反应过程，用颗粒表面反应模型模拟颗粒的表面燃烧过程。其中欧拉-拉格朗日模型假设颗粒不占体积，沿运动轨迹计算其参数，该模型的控制方程如下[8]：

(1)

(2)

其中各参数含义可参见文献[8]，该方程组通过有限体积法数值求解，计算格式为一阶迎风格式。湍流的参数也用相同方法求解k和ε的控制方程得到。

铝颗粒在燃烧时同时发生铝蒸气的气相反应及颗粒表面反应：

Alg+1.5H2O→0.5 Al2O3g+3H2↑

Als+1.5H2O→0.5 Al2O3s+3H2↑

从文献[8]可知，参加2种反应的铝质量比为4∶1。化学反应造成气相中铝、水蒸气等组分含量发生变化，而各组分含量通过气相组分输运方程求得：

(3)

在颗粒反应的过程中，首先发生气相反应，然后再发生颗粒表面反应。发生气相反应时，铝颗粒吸热蒸发为铝气，然后与水蒸气反应。反应的速率由颗粒蒸发速率和化学动力学速率控制，其中颗粒蒸发速率的计算式为

-dmp/dt=A0fv，0mp，0

(4)

式中mp为颗粒质量；mp，0为颗粒初始质量；A0为蒸发速率常数，计算假定A0=5×105s-1；fv，0为气相反应铝占铝颗粒初始质量的百分量，fv，0=0.8。

发生颗粒表面反应时，铝直接以凝相的形式与水蒸气反应，反应速率由气相组分扩散速率和化学动力学速率控制。气相反应铝已经蒸发，因此在计算铝颗粒表面反应的过程中，假定不存在被氧化铝包裹的未完全反应的铝，颗粒中铝最终全部变为氧化铝。这一过程中假设颗粒粒径不变，颗粒中铝/氧化铝的消耗/产生速率通过式(5)计算[9]：

(5)

式中Ap为颗粒表面积；Yj为颗粒组分j的质量含量；ηr为影响因子，ηr=1；Rj，r为由颗粒表面反应化学动力学速率、颗粒粒径、颗粒温度、气相温度、气相压强等共同确定的单元反应速率。

上述2种反应的化学动力学速率均由阿累尼乌斯公式计算，但是目前有关铝反应的动力学参数的研究很少且并不准确，本文计算暂先假定气相反应的活化能Ea，r1=100 J/mol，指前因子A1=1.14×107，不受温度影响；颗粒表面反应活化能Ea，r2=7.9 J/mol，指前因子A2=2×10-3。气相化学反应生成热全部被气相吸收，颗粒表面反应的生成热被气相吸收75%，剩余的全部被颗粒吸收。

在纳米铝冰发动机中，由于参与反应的铝颗粒粒径很小，燃烧时以单颗粒燃烧为主[10]。现有一发动机装药所用纳米铝粒径为以80 nm(0.08 μm)为中心的分布，下面计算假定该粒径分布如图2所示，粒径分布函数为Rosin-Rammler。

计算假设所有颗粒全部为球形颗粒，各颗粒内部温度均匀，颗粒表面的包覆层在颗粒进入燃烧室时已经被去除(包覆层通常是在高温下易分解的物质)。由于本文模拟的是发动机点火后装药稳态燃烧和内流场流动的情况，所以计算选定的初始条件为：pc=8 MPa，Tc=3000 K。

1.3 发动机燃烧室压强估算方法

在发动机实验前，由压强预估结果设计的装药和喷管决定了发动机能否正常工作。通常这一估算用热力计算和零维内弹道理论完成。热力计算实质上是用最小Gibbs自由能方法或其他方法计算推进剂的完全燃烧产物参数，并用“两相平衡流”假设计算喷管中的两相流[11]。热力计算得到铝冰推进剂的配方为Al∶H2O= 1∶1时，推进剂的特征速度C*=1373.4 m/s[12]。该特征速度值通常会直接用于零维内弹道(压强)计算。

发动机零维内弹道理论用质量守恒原则求解零维稳态燃烧室压强，由于计算简单而被广泛使用。零维平衡压强公式最终形式如下[13]：

pc=(ρpC*a·Ab/At)1/(1-n)

(6)

式中ρp、C*、a、n分别为推进剂密度、特征速度、燃速系数和压强指数；Ab和At分别为燃面和发动机喷管喉部面积。

除零维内弹道理论外，燃烧室压强的估算还可采用其他方法。这些方法通过流场计算求解燃烧室压强，包括考虑燃烧的CFD方法[14]、考虑两相流的数值方法[15]等。其中，CFD方法通过合适的简化数学模型计算推进剂反应物的燃烧过程和两相燃烧产物的流动过程，两相流数值计算则仅计算两相燃烧产物的流动过程，这2种方法相比零维内弹道理论考虑了两相不平衡效应(以及推进剂的燃烧)对燃烧室压强的影响，理论上结果更准确。

2 计算结果分析

计算选用Fluent软件中密度基稳态求解器，通过反复迭代得到发动机内流场稳态解(流量误差波动小于5%)，最终得到燃烧室压强、燃烧室温度、气相反应和颗粒表面反应区厚度、计算监测的颗粒数量及变化特征等参数如表1所示。

从表1可看出，发动机稳态工作时燃烧室压强在8.04～10.78 MPa变化，燃烧室压强平均值为9.38 MPa，燃烧室温度平均值为2950.6 K。

和所有固体火箭发动机一样，铝冰发动机工作过程中推进剂组元的分解/相变产物以较低温度进入燃烧室，在燃面附近迅速发生反应放热，由此建立起燃烧室的高温高压环境。在所选的发动机的3个工作状态下，计算得到发动机内流场中气相总压沿发动机轴线的变化情况如图3所示。

表1 不同发动机工作状态下数值计算结果

注：1)计算按燃层厚分为3个状态，已在图1中标出。

在铝冰发动机中，贡献压强的工质主要是氢气，随着发动机工作的进行，燃烧室压强呈下降趋势，压强值从10.78 MPa变化到8.04 MPa。燃气通过喷管时膨胀做功，压强迅速下降，在喷管出口下降为环境压力。铝和水在反应过程中会放出大量热量，建立燃烧室的高温环境，计算得到燃烧产物总温在发动机不同工作时刻的变化情况如图4所示。

由图4可见，推进剂反应物从燃面进入燃烧室时温度很低，但通过反应区后温度迅速升高，放出反应热并在燃烧室内形成了一个较均匀的高温环境。该高温环境维持了推进剂的持续燃烧，而随着燃面推移，燃烧室温度在2900～3000 K之间稍有变化，平均值为2950.65 K。用热力计算软件CEA[11]对相同配方的铝冰推进剂进行计算，在燃烧室压强为9 MPa的情况下得到其绝热燃烧温度为3054.76 K。可见数值模拟得到的发动机燃烧室温度相比推进剂绝热燃烧温度偏低，分析发现这是由于数值模拟忽略了氧化铝的凝结放热而造成的。铝的气相反应产物是氧化铝气，目前并没有研究证明这一物质实际存在，在简化分析和计算中，可用它来代替铝的气相燃烧产物[16]。实际上氧化铝在流动过程中会因温度下降而逐渐凝结为氧化铝颗粒，并最终随燃气流出。研究表明，这些新凝结的颗粒相比由铝颗粒直接变成的氧化铝颗粒尺寸更小，更容易与燃气达到两相平衡，所以这些颗粒往往被研究者们用欧拉方法处理成气相的一部分[8，17]。因此，本文为了简化计算而忽略了这一部分氧化铝的凝结过程。虽然计算得到的燃烧室温度与实际有偏差，但是由于氧化铝密度大，对压强影响小，计算表明即使氧化铝气全部变为氧化铝颗粒，压强的计算结果变化也并不明显。另外，从图4还可看出，在发动机头部某些位置存在温度局部变化较大的情况，这是由于一些颗粒在发动机头部运动不能向外排出，在这些区域聚集并放出大量热量导致的。在计算选定的三个状态下，内流场中气相马赫数沿发动机轴线的变化如图5所示。

从图5可见，燃烧室中气流速度很低，通过喷管时流速迅速增大，并在喷管出口达到超音速(Ma>2.2)，在发动机工作的不同时刻，出口气相马赫数变化不大。铝颗粒燃烧时，气相反应和颗粒表面反应同时发生，在图1中所选定的3条监测线上(这3条监测线分别从垂直于燃面方向远离燃面，长度均为10 mm)，可给出气相反应速率、颗粒表面反应速率和化学反应热沿监测线远离燃面方向的变化分别如图6、图7所示。

从图6可看出，2个反应在监测线长度范围内就已经完成，在监测线上，反应速率大于0的区域(反应区)很薄，其中气相反应区的平均厚度约为7.7 mm，颗粒表面反应的反应区平均厚度约为1.4 mm，具体的厚度值已在表1中列出。铝颗粒进入燃烧室后很快发生反应，反应速率最大值在贴近燃面的位置，最大气相反应速率为4×105～1.2×106kmol/(m3·s)，最大表面反应速率为50～110 kmol/(m3·s)。推进剂反应物经过反应区快速完成反应后，其产物把反应所产生的热量带入燃烧室，因此相比反应速率而言，反应的放热区向远离燃面的方向有一定程度的推移和放大。

从图7可看出，即使在反应速率为0的区域，反应热仍不为0，反应速率最大的区域反应热也并非最大，这正是由于反应物边流动边反应，把热量带出反应区，并向周围扩散造成的。假定反应放热量大于0的区域为反应放热区，在不同燃层厚下反应放热区厚度的近似值如表1所示，其平均厚度约为12 mm，可见放热区厚度比反应区厚度大，这是由于燃烧产物的流动造成的。通过放热区后混合物的温度迅速升高，在反应放热区附近形成很大的温度梯度，如图4所示。

颗粒的燃烧和运动过程对内流场影响较大，图8给出了不同燃层厚下数值计算所监测的颗粒中铝质量占颗粒总质量(铝+氧化铝)的百分数沿颗粒轨迹的变化过程。

由图8可见，由于颗粒粒径较小反应速率大，铝颗粒在燃面附近很快就反应完成变成氧化铝。大部分颗粒随燃气从喷管流出，但仍有少部分颗粒最终会沉积在燃烧室中，形成残渣。计算监测的颗粒变化的最终结果已在表1中列出，可以看出，不同燃层厚下计算监测的颗粒总数并不相同，而所有颗粒中包含气相反应颗粒、表面反应颗粒和沉积的颗粒，颗粒的这些变化过程可用1.2节所给列出的模型近似计算。在燃烧室中沉积的颗粒会造成沉积区域附近流场有较大温度变化(见图4)。

3 计算结果的验证与讨论

3.1 内流场数值结果的验证

为了验证数值计算结果的有效性，下面给出与图1相同尺寸的铝冰发动机在相同的工作条件(水燃比1∶1)下实验测量得到的燃烧室压强曲线，如图9所示。

从图9可看出，发动机稳态工作时间约为1 s，稳态工作时平均压强约为9.33 MPa，这一压强值与数值模拟得到的平均燃烧室压强(9.38 MPa，如表1)很接近。数值计算主要是对发动机的3个稳定工作状态进行的，假定在这3个状态下燃烧室实验压强值如压强曲线中3个“╋”点所示。为了比较流场数值计算得到的燃烧室压强、零维内弹道理论得到的燃烧室压强、铝冰完全燃烧产物两相流平衡流流场计算得到的燃烧室压强、完全燃烧产物两相非平衡流计算得到的燃烧室压强与实验测量得到的压强的差别，把各个燃层厚下压强结果列如表2所示(其中两相流数学模型参见文献[18])。

从表2可看出，数值模拟计算的燃烧室压强相比零维内弹道理论的计算结果更接近实验值，且考虑燃烧的CFD计算结果相比实验的平均误差仅为4.62%，小于工程误差(5%)。使用两相非平衡流方法计算铝冰推进剂两相燃烧产物的流动时，燃烧室压强的计算结果与实验误差为8.52%，相比零维内弹道的计算结果，计算精度已有很大的提高。在推进剂热力计算和零维内弹道理论中都使用了两相平衡流假设，从表2压强计算结果可看出，在两相平衡流假设下计算的燃烧室压强与实验值的误差为18.34%，相比零维内弹道计算结果的误差(50.74%)已经有很大的提高，但其计算精度仍比两相非平衡流模型和CFD模型低，这说明在计算中应该考虑两相非平衡、推进剂的燃烧过程等因素的影响。

表2 不同方法得到的铝冰发动机燃烧室压强

上述结果证明，目前常用的零维内弹道理论在铝冰发动机中使用时燃烧室压强的计算结果相比实验误差较大，而用CFD模型计算的铝冰发动机燃烧室压强结果与实验相比误差很低。这说明用CFD方法得到的铝冰发动机内流场参数是可信的，这一方法可推广到其他高金属含量发动机。

3.2 铝冰发动机燃烧室压强预测公式

为了得到简单的公式以在铝冰发动机实验前对发动机燃烧室压强进行准确估算，需要对传统压强计算方法进行修正。根据本文数值计算(CFD模型)得到的发动机不同工作状态下的燃烧室压强结果，与零维内弹道理论得到的燃烧室压强值进行比较，并对后者进行修正，得到铝冰发动机零维燃烧室压强pc，理的修正式如下：

pc，修正=0.7055pc，理

(7)

实验前可用该公式修正铝冰发动机的零维燃烧室压强计算结果，这样就在压强计算中近似考虑了推进剂燃烧和两相燃烧产物的流动过程对压强的影响。对于同类铝冰发动机，该公式有望提高压强计算的准确度并降低实验失败的机率。

4 结论

(1)从铝冰发动机的实验曲线和数值模拟结果可以看出，纳米铝的反应活性很强，铝和水蒸气在推进剂燃面附近就能迅速反应且放出大量热量，铝冰发动机能够正常工作。

(2)在含金属发动机中推进剂燃烧产物是两相不平衡流，使用传统零维内弹道理论计算燃烧室压强时，误差较大，而使用两相流数值计算的方法，求解两相流控制方程，能考虑两相不平衡效应。使用数值计算方法(CFD模型)可有效预估铝冰发动机的燃烧室压强，该方法可用于工程计算。

(3)用压强修正公式能够简单、快捷地在铝冰发动机燃烧室压强计算中考虑推进剂燃烧和两相不平衡效应对压强的影响，提高零维燃烧室压强计算结果的准确度，从而提高实验成功率，降低实验成本。

[1] 庞爱民,黎小平.固体推进剂技术的创新与发展规律[J].含能材料,2015,23(1):3-6.

[2] 胡润芝,张小平.高能推进剂燃烧效率研究和实测比冲预估[J].推进技术,2001,22(5):415-417+436.

[3] Risha G A,Connell Jr T L.Novel energetic materials for space propulsion[R].Nasa:DTIC Document:A546818,2011.

[4] Risha G A,Sabourin J L.Combustion and conversion efficiency of nanoaluminum-water mixtures[J].Combustion Science and Technology,2008,180(12):2127-2142.

[5] Ingenito A,Bruno C.Using aluminum for space propulsion[J].Journal of Propulsion and Power,2004,20(6):1056-1063.

[6] Geisler R L.A global view of the use of aluminum fuel in solid rocket motors[C]//38th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit,2002.

[7] 姜治深,梁晶晶.固体推进剂中铝燃烧理论研究展望[J].飞航导弹,2013(10):88-91,96.

[8] Sabnis J S.Numerical simulation of distributed combustion in solid rocket motors with metalized propellant[J].Journal of Propulsion and Power,2003,19(1):48-55.

[9] Ansys_Inc.Fluent 15.0 Help[R].2013.

[10] 李颖,宋武林.纳米铝粉在固体推进剂中的应用进展[J].兵工学报,2005,26(1):121-125.

[11] Gordon S,Mcbride B J.Computer program for calculation of complex chemical equilibrium compositions and applications[R].Nasa,Cleveland,OH,United States:NASA Lewis Research Center,1994:58.

[12] 韩秀杰.铝-冰固体推进剂燃烧性能研究[D].哈尔滨:哈尔滨工程大学,2013.

[13] 李宜敏,赵元修.固体火箭发动机原理[M]. 北京:国防工业出版社,1985.

[14] Puskulcu G.Analysis of 3-D grain burnback of solid propellant rocket motors and verification with rocket motor tests[D].MS.Thesis,Dept.of Mechanical Engineering,METU,2004.

[15] 陈军.固体火箭发动机零维两相内弹道研究[J].弹道学报,2013,25(2):39-43.

[16] Daniel E.Eulerian approach for unsteady two-phase solid rocket flows with aluminum particles[J].Journal of Propulsion and Power,2000,16(2):309-317.

[17] Najjar F M,Massa L.Effects of aluminum propellant loading and size distribution in BATES motors:A multiphysics computational analysis[C]//41st AIAA/ASME /SAE/ASEE Joint Propulsion Conference & Exhibit,2005.

[18] 方丁酉.两相流动力学[M].长沙:国防科技大学出版社,1988:187,545.

(编辑：吕耀辉)

Numerical simulation of aluminum-ice solid propellant motor

LIU Ping-an，WANG Liang，WANG Lu，WANG Ge

(College of Aerospace and Civil Engineering，Harbin Engineering University，Harbin 150001，China)

To calculate the theoretical performance of the Aluminum-Ice Solid Rocket Motor(ALICE SRM),the particle surface reaction model and the gas phase reaction model were used to simulate the burning process of aluminum particle and the steam in the combustion chamber of the ALICE SRM.The Euler-Lagrange model was used to calculate the particle parameters along their trajectories.The numerical results were analyzed and compared with the experimental results of the same sized ALICE SRM.The numerical chamber pressure of the ALICE SRM is about 9.38 MPa,which is very close to the experimental value.The calculated average chamber temperature is 2950.65 K,which is lower than that of the thermo-chemistry calculation result.By the numerical simulation of the internal flow-field of the ALICE SRM,a good agreement between the numerical results and the experimental data was obtained,which clearly verified the numerical models in predicting ALICE SRM performance.

aluminum-ice solid rocket motor；internal flow-field；numerical calculation

2016-07-12；

2016-12-15。

中央高校基本科研基金(HEUCFD1404；HEUCFD1502)。

刘平安(1982—)，男，副教授，研究方向为金属燃料发动机。E-mail：liupingan631@126.com

V435

A

1006-2793(2017)04-0425-07

10.7673/j.issn.1006-2793.2017.04.005