基于King模型的含硼富燃燃气燃烧模拟①

冯喜平，刘 洋，任全彬，李进贤

(1.西北工业大学燃烧、热结构与内流场重点实验室，西安 710072 2.中国航天科技集团公司第四研究院，西安 710025)

0 引言

硼以高体积热值优势已成为理想的推进剂材料，在固冲发动机中表现出巨大的应用潜力[1]。但硼的高熔点(2 450 K)与高沸点(3 931 K)，使其点火困难。

国外自20世纪70～90年代就对硼粒子的点火机理进行了大量研究[2－6]，并建立半经验公式，其中硼粒子点火模型以 King模型[3]和 L-W 模型[4]为代表。King点火模型基于O扩散至B-B2O3表面与B反应生成B2O3(l)的假设，考虑了氧化层的产生、消除对B点火和燃烧的影响。L-W模型认为，在高温情况下，硼在氧化层界面上溶解，与B2O3发生反应生成聚合物(BO)n;该聚合物在氧化层中由内向外扩散，并在外界面上通过蒸发和化学反应消耗。Yeh和 Kuo[5－6]也研究了硼粒子点火过程，发现B颗粒在1 213 K就开始熔化，远低于其熔化温度2 450 K;在温度升高的过程中，B不断溶解于B2O3，并与B2O3反应生成了一种新的聚合物(BO)n。

国内霍东兴等[7]采用King点火模型和L-W燃烧模型研究了某发动机补燃室内硼粒子粒径对点火位置和点火效率的影响，得出2μm以下的硼粒子燃烧效率较高的结论。吴婉娥[8]建立了单个硼粒子的点火燃烧数学模型，通过数值方法研究环境温度、氧分压、水蒸气分压、粒径对硼颗粒点火的影响，得出高温环境、高水蒸气分压、小粒径有利于硼粒子点火，但过高氧分压会延长硼粒子点火时间的结论。李海波[9]对比了不同湍流模型、气相燃烧模型对模拟结果的影响，并给出冲压发动机补燃室中适用的一些模型[9]。

对于固冲发动机而言，含硼富燃燃气二次燃烧是其关键技术。对于二次燃烧过程的计算，燃烧模型是关键。尽管国内外研究者进行了大量相关研究，但主要集中于硼粒子的点火燃烧模型等细节性和局部性研究环节，许多研究采用的数值模拟方法过分依赖于商业软件自带模型，对燃烧过程过分简化，这些燃烧模拟普遍存在缺乏完整数值模型的问题(包括颗粒相燃烧模型、气相燃烧模型等)。

针对上述问题，本文以燃烧实验装置为物理模型，选用King模型作为硼粒子点火燃烧模型，选用有限速率/涡耗散作为气相燃烧模型，构建含硼富燃燃气流动与燃烧过程的数学模型，结合典型实验工况进行数值模拟，通过模拟结果与实验结果对比，验证模拟方法。

1 物理模型与计算方法

1.1 物理模型

燃烧室结构如图1所示。燃烧室由燃气发生器、燃烧室和尾喷管3部分组成。

图1 燃烧室结构Fig.1 Structure of the chamber

为方便实现燃烧流场的测量，燃烧室采用横截面为100 mm×100 mm、长635 mm的长方体结构;数值计算时，为减少计算量，将燃烧室简化为等流通面积的圆柱体，直径113 mm。其他参数见图1。整个燃烧室采用轴对称结构。

实验装置工作原理如下:含硼富燃燃气由燃气入口进入燃烧室，空气由位于燃烧室前部的环形空气入口喷入燃烧室，二者进入燃烧室后沿各自流线运动;在燃气和空气组分梯度的作用下相互扩散、燃烧，在燃烧室形成扩散火焰;燃烧后产物由喷管喷出。燃气流量通过燃气发生器喷管音速喉道控制，空气流量通过空气管道中的音速喉道控制。

1.2 数学模型

1.2.1 基本假设和控制方程

基本假设:

(1)燃烧室内燃气流动为准定常流动，与外界无热交换;

(2)忽略燃气各组分之间的辐射作用，忽略重力的影响;

(3)燃烧室燃气为理想气体，服从理想气体状态方程;

(4)不可燃气体假设为N2。

控制方程:

式中 各项表达式见文献[10]。

1.2.2 湍流模型

文献[9]通过对比不同湍流模型对模拟结果的影响，得出RNG k-ε湍流模型更适用于本物理模型内二次燃烧的结论。因此，湍流模型选取RNG k-ε模型。具体方程及说明见参考文献[11]。

1.2.3 燃烧模型

(1)气相燃烧模型

气相燃烧模型采用通用有限速率/涡耗散模型，分别计算Arrhenius和涡耗散反应速率，净反应速率取两个速率中较小的。Arrhenius反应速率作为一种动力学开关，阻止反应在火焰稳定之前发生。一旦火焰被点燃，涡耗散速率通常会小于Arrhenius反应速率，这样可延迟计算点火的开始，提供点火时间，更加符合实际[9]。

(2)碳颗粒燃烧模型

碳颗粒燃烧模型采用动力学/扩散控制反应速率模型。该模型假定表面反应速率同时受到扩散过程和反应动力学的影响。该模型的扩散速率常数为

动力学反应速率常数为

二者加权得反应速率为

式中 pox为颗粒周围的气相氧化剂分压;R为考虑了焦碳内表面的反应及其扩散的化学(动力学)反应速率常数。

(3)硼颗粒相燃烧模型

硼颗粒相点火燃烧模型选用King模型[3]。

硼粒子表面覆盖有致密的氧化层，阻止洁净硼颗粒与氧气直接接触，见图2。

硼粒子点火经历3个过程:

(1)表面B2O3蒸发;

(2)B2O3与H2O反应生成HBO2气体;

(3)O2扩散至B-B2O3边界，与洁净B发生反应，生成新的B2O3。

图2 King模型中硼粒子点火机理Fig.2 IgnitionofboroninKingmodel

King指出，只有B2O3的消耗速率大于生成速率时，才能认定B点火成功，进入燃烧阶段。

其数学描述[3]如下:

式中 dB、δ、f、TP分别为硼粒子直径、氧化层厚度、熔化分数、硼粒子温度;为点火过程中的能量变化，包括颗粒反应释放的热量及其与周围环境之间的对流换热和辐射换热;T、TRAD分别为环境温度和周围辐射温度;RB、RE、RH分别为 B消耗速率、B2O3蒸发速率、B2O3与水的反应速率;pO2和pH2O分别为环境气体中氧气和水蒸气的分压;ΔHM、p、XH2O、XO2分别为 B 的熔化热、环境压力、H2O的摩尔分数、O2的摩尔分数。

燃烧阶段，洁净硼粒子反应速率遵从有限速率模型:

式中 k为Arrhenius定律计算出的反应速率;dB为当前洁净硼粒子粒径。

为了完整、准确描述硼颗粒相点火燃烧过程，采用Fluent扩展应用UDF编写King模型程序。程序编制中，采用如下假设:

(1)假定硼粒子表面存在1%厚度的氧化层;

(2)认定H2O(g)存在能加快点火过程，但不考虑HBO2的生成反应;

(3)认为颗粒相中只存在一种成分B，且未点火前B颗粒质量、密度、粒径不发生变化，点火完成前B2O3的生成和消耗速率进行单独计算，不加入迭代方程中;

(4)点火过程完成后，硼粒子遵从Fluent自带颗粒表面反应模型。

1.3 网格划分

图3为计算网格。

图3 计算网格Fig.3 Meshofthechamber

采用四边形与三角形网格相结合，部分加密的方法对模型进行网格划分。其中，1区域为燃气通道，采用三角形网格并加密;2区域为燃烧室，采用四边形网格，中部区域加密;3区域为尾喷管，采用四边形网格，整体加密。总网格数量53854。

1.4 边界条件

边界条件见表1。

表1 边界条件Table1 Boundaryconditions

燃气流量为 0.03kg/s，空燃比为 10。

颗粒相采用离散相、面入射、随机轨道模型，初始粒径为10－5m，速度为100m/s，温度1800K。

燃气组分通过燃气发生器热力计算得到，结果见表2。燃气中，气相组分32.2%，颗粒相67.8%。

表2 富燃燃气主要组分Table 2 M ain component of fuel-rich gas

2 计算结果与讨论

2.1 燃烧室流场分布

图4为模拟得到的燃烧室流场分布。

图4(a)为燃烧室压力等值线分布。由图可见，整个燃烧室压力均匀，平均压力为0.48 MPa;燃烧室出口处，由于拉瓦尔喷管作用，压力逐渐降低。

图4(b)给出了燃烧室内温度分布。最低温度出现在空气入口处，为573 K;高温区域(2 800 K)分别出现在燃烧室前端燃气与空气接触面(210mm处)、燃烧室1/3(280～370 mm处)和燃烧室中后部靠近轴线区域(400 mm～尾部)。燃烧室前端燃气与空气接触面处，燃气与氧气发生剧烈的化学反应。因此，邻近微小空间内温度急速上升，形成第一个2 800 K高温区域;前部区域空气、燃气速度梯度与温度梯度较大，在接触面附近形成严重湍流，热量很快向周围传递，2 800 K高温区消失。

图4 燃烧室流场分布Fig.4 Distribution of flow field in the chamber

燃烧室280 mm处，气相燃料基本燃烧完毕，少量硼粒子点火成功，开始与氧气反应，形成第2个高温区域;由于仅有少量硼粒子在此区域点火成功，且很快完全燃烧，其他硼粒子尚未开始反应。因此，第2个2 800 K高温区持续90 mm后就消失。燃烧室400 mm处，氧气向中心轴线区域扩散，大量硼粒子点火成功，开始剧烈的燃烧反应，与前两个高温区域不同，第3个2 800 K高温区由空-燃接触面扩展至中心轴线处，表明氧气已扩散至中心区域，与大量处于轴线附近的硼粒子发生反应，释放巨大热量;高温区域一直延展到燃烧室尾部，说明硼粒子燃烧反应持续进行。

从整个温度场结构来看，随着流动进行，燃气向壁面方向扩散，与外部氧气发生反应，火焰宽度逐渐增大;尾部区域，温度场分层现象仍然严重，中心温度2 800 K，壁面温度仅为1 400 K，二者差值达1 400 K。温度场呈锥形扩散火焰形状。

图4(c)为燃烧室速度分布图。速度分布呈现射流特征。在头部区域靠近轴线处，中心射流速度高达140 m/s，而外部流场速度仅10 m/s左右。二者接触面速度梯度大，分层明显。由于气相的粘性作用、燃气与空气之间分子扩散产生的动量交换，中心气流速度逐渐降低，外部气流速度逐渐升高。燃烧室200mm左右位置，射流核心区消失，射流主体区域出现，此时轴线中心速度依然高达80 m/s。而在燃烧室中后部510 mm处，速度梯度已经很小，流场速度趋于均匀化。如果燃烧室足够长，流动到一定阶段后，整个流场速度将完全均匀，变为层流运动。

图4(d)给出了O2的分布。在燃烧室前部区域，氧化剂由空气入口处进入，O2质量分数最大;燃气中不含O2，故O2质量分数最小。随着流动进行，燃气向外扩散，二者在火焰面处进行燃烧，消耗氧气，高含氧区域缩小，氧组分等值线呈现向外弯曲形状;随着流动向下游进行，氧气进一步消耗。在燃烧室轴线附近区域，尽管存在氧气向中心区域的扩散，但由于燃烧的进行，氧气在火焰面位置大量消耗，因此燃烧室轴线附近区域氧含量非常小。

图4(e)给出了B2O3的分布。在燃烧室头部不存在B2O3。这是由于硼粒子刚进入燃烧室，尚未点火成功。210 mm处，B2O3开始生成，标志部分B粒子点火成功，并开始与O2发生反应。B2O3生成位置并非在中心轴线处，而是处于氧气与燃气的接触面附近区域。这是由于燃气中可燃气体成分在头部区域大量消耗O2，处于中心区域的硼粒子几乎很难接触到氧气。只有少量硼粒子偏离中心，射入燃气与氧气反应的区域中，在高温、高氧气分压、高水蒸气分压条件下，很快点火成功，生成B2O3。400 mm处，B2O3浓度明显提高，即大量硼粒子开始与O2反应。这一点在图4(b)温度分布中也有所反映:400mm处，燃烧室温度第三次达到峰值2 800 K。在这一阶段，中心区域也生成大量B2O3，这是由于燃气已经消耗完毕，O2扩散至中心，与轴线附近的B粒子发生反应。B2O3最高浓度出现在燃烧室尾部，说明B颗粒与O2的剧烈反应一直持续到燃烧室尾部;O2质量分数分布也支持这一结论，中后部区域中心轴线附近O2浓度仍为0，而靠近壁面处O2浓度基本达到一致，说明中心轴线区域持续发生燃烧反应。

2.2 硼粒子点火距离

硼粒子进入燃烧室到硼粒子剧烈燃烧的距离为硼粒子的点火距离，模拟结果见图5～图7。

图5 B颗粒群质量变化曲线Fig.5 M ass-length of boron particles

图6 B2O3质量分数分布图Fig.6 M ass fraction of B2 O 3

图7 硼粒子轨迹图Fig.7 Boron particle tracks

图5为燃烧室内不同入射位置的硼粒子质量变化曲线。图6为燃烧室头部B2O3的质量分数分布图(出现B2O3即认定点火成功)。图7为不同入射位置硼粒子轨迹。由于入射位置不同，硼颗粒在燃烧室中的轨迹也不同。靠近轴线方向的硼粒子周围氧气与水蒸气浓度较低，温度上升慢，点火距离较长。靠近边界层的硼粒子能更早地接触氧气与水蒸气，且温度较高，因此点火距离短。以图5中硼颗粒1、2为例，其轨迹如图7中所示。颗粒1最先点火成功，但燃速很慢，300 mm处突然燃速加快，并且很快完全燃烧，对应图7中颗粒1此前一直处于中心轴线处，温度、氧分压、水蒸气分压低下;在300 mm处开始偏离中心轴心，射入高温反应区(温度图见图4(b))，因此很快燃烧完毕，轨迹线消失。颗粒2也是由于偏离中心轴线，得以快速点火燃烧。

由图5和图6得出，硼粒子在210mm处质量发生变化，即成功点火，除去燃气入口及延长段147 mm距离，得出硼粒子点火距离为63 mm。大部分硼粒子在400 mm处才点火成功，开始剧烈燃烧，这种现象也可说明图4(b)中400 mm处2 800 K高温区产生的原因。由图5和图7得出，大量硼粒子没有完全燃烧，随气体一起流出尾喷管。

图8是实验者捕捉的稳定燃烧阶段燃烧室中硼与氧气反应时明亮火焰图像，即硼粒子点火位置。经过处理计算，可得出稳定阶段点火距离约为50 mm。数值模拟中硼粒子点火距离为63 mm，二者数值接近。

图8 硼粒子燃烧火焰图Fig.8 Flame of boron combustion

2.3 模拟结果与实验结果对比

为了验证King点火燃烧模型的准确性，将数值模拟结果与实验参数对比，结果见表3。

表3 模拟结果与实验结果对比Table 3 Com parison of simulation and experim ent result

由表3可得出，使用King模型进行数值模拟得出的燃烧室压力为 0.48MPa，实验结果为 0.43MPa，误差11.6%;点火距离为63 mm，实验结果为 50 mm，误差26%。可见，使用King模型的模拟结果较好地预测了实验工况。

3 结论

(1)进行含硼富燃燃气的二次燃烧时，粒子相采用King模型点火与燃烧，气相燃烧采用有限速率/涡耗散模型，湍流采用RNG k-ε模型，模拟结果与实验结果对比说明，该模拟方法较好预测了燃烧过程。

(2)采用Fluent软件UDF功能，依据粒子相点火燃烧模型自编程序，克服了目前数值模拟中燃烧过程过于简化问题，较准确地表述了硼粒子点火燃烧机理和燃烧过程，模拟结果更加可靠。

(3)数值模拟表明，粒子相在距燃气出口一定距离点火并开始燃烧，燃烧火焰头部呈锥形火焰结构。

[1] 汤中权，锐鑫.含硼固体推进剂[J].化学推进剂与高分子材料，1998，3:40-42.

[2] 王宁飞，关大林，范红杰.硼颗粒点火和燃烧研究进展[J].含能材料，2001，9:86-89.

[3] King M K.Ignition and combustion of boron particles and clouds[J].Journal of Spacecraft and Rockets，1982，19(4):294-306.

[4] Li SC，Williams F A.Ignition and combustion of boron particles[M].Combustion of Boron-Based Solid Propellants and Solid Fuels.Boca Raton，CRC Press，1993:248-271.

[5] Yeh C L，Kuo K K.Ignition and combustion ofboron particles[J].Progress in Energy and Combustion Science，1996，22(6):511-541.

[6] Kuo K K.Boron ignition model[M].NPU，1998.

[7] 霍东兴，陈林泉，刘霓生，等.硼粒子直径对点火位置及燃烧效率的影响研究[J].固体火箭技术，2004，27(4):272-275.

[8] 吴婉娥，裴明敬，郭耳铃，等.硼粒子在固冲环境中点火过程影响因素的数值模拟[J].火炸药学报，2008，6:79-82.

[9] 李海波.含硼燃气与空气扩散燃烧模型研究[D].西安:西北工业大学，2013.

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

[11] Yakhot，Orszag S A.Renormalization group analysis of turbulent[J].J.Sci.Comput.，1986，1:3-9.