基于FDS与CFD组合的快速烤燃数值模拟

肖 游，智小琦，王 琦

(中北大学 机电工程学院，山西 太原 030051)

引 言

快速烤燃试验是评估弹药热易损性的重要指标之一，国内外学者已经对弹药的快速烤燃进行了许多研究。戴湘晖等[1]对大尺寸侵彻弹进行了快速烤燃试验，检验在大火中的热敏烈度；美国犹他大学的Ciro W等[2]进行了一系列在钢壳限制下、以HMX为主的混合炸药的快速烤燃试验，结果表明，试验可以直观评估弹药热安全性。

由于快速烤燃试验的高成本和对环境的污染性较重，大部分学者通过数值模拟的方法预测弹药的点火时间与点火温度。Yang H W等[3]采用FLUENT软件对试件施加不同升温速率和不同热流密度等边界条件进行数值模拟，再与以酒精为燃料的快速烤燃试验对比，发现点火位置为靠近外壁的环形区域，热流密度增加、点火时间缩短、点火温度提高。Gross M L等[4]以试验所测的平均热通量为边界条件，采用一维瞬态模拟方法，对以HMX为主的混合炸药进行快速烤燃研究，提出了对于小尺寸弹药使用经验数据预测点火时间的方法。Sahin H等[5]通过将试验测得平均温度作为边界条件，利用FLUENT软件进行仿真，研究火箭发动机的快速烤燃特性，并用MATLAB软件计算了快速烤燃过程中燃烧室内部的压力。

从上述情况可见，目前国内外对快速烤燃仿真研究，主要是通过FLUENT软件设定较快的升温速率以控制壳体的升温，或者对试件表面采用完全相同的热通量作边界条件进行模拟。这与实际快速烤燃状态不完全吻合。因为快速烤燃试验时，试件表面温度不可能完全相同，而且火焰的特征量是随时间变化的，这样试件表面的温度也是变化的。

本研究采用火灾动力学软件(FDS)模拟装填熔铸B炸药(60%RDX/40%TNT)试件的快速烤燃试验，通过池火火焰产生的对流与辐射热通量对试件的传热监测试件的温度变化，并将试件不同位置的真实温度作为边界条件用于FLUENT软件中，研究火焰不同阶段的特征及油池尺寸对其的影响，利用B炸药多步分解机制计算试件内部炸药的熔化情况，以期为弹药快速烤燃模拟方法的研究与试验装置设计提供参考。

1 理论模型

快速烤燃过程中，火焰燃烧环境为开放空间，液体燃料燃烧产生的热对流和热辐射共同作用于试件，试件通过热传导使内部装药温度升高，直到炸药达到点火温度发生点火反应。

1.1 热通量模型

(1)

辐射热通量是试件吸收的辐射热通量和试件表面发射的辐射热通量之差，忽略试件表面向内传递的热通量，不考虑各种波长的影响。因此，当试件表面吸收率和发射率相等时，表面接收的辐射净热通量可以表示为：

(2)

火灾的特点是温度分布不均匀，理论上入射辐射热通量应该包括附近火焰和热气等所有辐射源的贡献，即入射辐射热通量可以写成所有辐射源贡献的总和，设Fi为视图因子，则[6]：

(3)

式中：εi为不同物质的发射率；Ti为不同物质温度，K。

对流热通量取决于周围热气流温度和试件表面温度之差。由牛顿冷却公式[8]：

(4)

式中：h为换热系数，W·m-2·K-1；Tg为试件表面的气体温度，K。

由此可知，通量与温度差值成正比。则传入试件表面单位面积的总净热通量可以表示为：

(5)

1.2 绝热表面温度

将暴露在火焰中的试件表面视作完美绝热体，其温度应称为绝热表面温度TAST。故壳体表面的总净热通量为零。

(6)

根据Wickström U[9-10]的研究，用平板温度计测量火灾中物体的表面温度，这个温度接近绝热表面温度，以绝热表面温度作为计算火灾作用下物体温度变化的边界条件，并应用于仿真计算中。

1.3 接触热阻

接触热阻广泛存在于相互接触的两个固体之间，实际上的接触仅仅发生在一些离散的面积元上，在未接触的界面间隙之间常常充满了空气，增加了热量通过间隙传递的阻力。在炸药装药中，由于壳体与炸药之间并非完全贴合，而是通过离散的点接触，这就造成了壳体与炸药之间接触热阻的存在。接触热阻的计算公式为[11]：

(7)

式中：R为接触面间接触热阻，m2·K·W-1；A为垂直于热通量方向的面积，m2；k为热导率，W·m-1·K-1；Δx为接触面之间的距离，m。丁洋等[12]的研究与本研究参数相近，故炸药与壳体之间的接触热阻取值为0.003m2·K·W-1。

1.4 炸药分解动力学

对B炸药的仿真计算做如下假设：(1) 忽略炸药体积变化；(2) 炸药的自热反应遵循Arrhenius方程；(3) 忽略气体产物对传热的影响。

炸药在烤燃过程中的质量、动量、能量的连续方程可用以下通用形式来表示[13]：

(8)

式中：φ为通用变量，可表示质量、动量、能量等；ρ为密度，g/cm3；Γ为通用的扩散系数；μ为黏度，Pa·s；S代表炸药自热反应源项。

为精确计算炸药的分解过程，并监测中间产物质量分数的变化，采用多步反应动力学模型描述RDX和TNT炸药的分解过程。炸药中RDX的分解机理为[14-15]：

A→B,r1=Z1exp(-E1/RT)φA

(9)

B→D,r2=Z2exp(-E2/RT)φB

(10)

D→F,r3=Z3exp(-E3/RT)φD

(11)

式中：r为反应速率，mol/(L·s)；Z为指前因子，s-1；E为活化能，J/mol；R为普适气体常数，8.314J/(mol·K)；A为RDX炸药；B为固体中间产物；D为气体中间产物；F为气体终产物；φA、φB、φD分别为其对应的质量分数。

炸药中TNT的三步分解机理为[16]：

G→H,r4=Z4exp(-E4/RT)φG

(12)

H→I,r5=Z5exp(-E5/RT)φH

(13)

I→J,r6=Z6exp(-E6/RT)φI

(14)

式中：G为TNT炸药，H和I均为中间产物，J为气体最终产物；φG、φB和φI分别为对应的质量分数。

对于每一步反应，热量生成速率可表示为：

S=riφiQi

(15)

式中：Q为反应热，J/kg；i表示多步反应动力学的第i步，i=1,2,3。

故RDX和TNT在多步反应中生成热量为：

SRDX=S1+S2+S3

(16)

STNT=S4+S5+S6

(17)

将计算网格设为混合热分解单元[17]，RDX和TNT反应吸收或释放的热量为各自多步反应吸收或释放的热量，单元总热量为RDX和TNT吸收或释放热量的总和。这样，炸药热分解过程中生成的总热量为：

SB=0.6SRDX+0.4STNT

(18)

炸药与壳体的物性参数及反应动力学参数分别见表1和表2。

表1 炸药与壳体的物性参数Table 1 Physical parameters of explosives and shells

表2 炸药反应动力学参数Table 2 Kinetic parameters of explosive reactions

为了验证快速烤燃数值模拟的可行性，基于快速烤燃试验进行模拟验证。

2 B炸药快速烤燃试验

烤燃装置由油池、支架、航空煤油、热电偶、数据采集系统、烤燃试件等组成。试件尺寸为Φ76mm×256mm，壳体和端盖厚度均为7.5mm。B炸药药柱尺寸为Φ61mm×241mm，装药量为1190g，充满空间。油池的长×宽×高=1000mm×800mm×100mm，油池与支架材料均为45号钢，燃料为JP-8航空煤油。热电偶为WRN-130的K型热电偶，量程为0～1200℃，精度0.1K。试件水平放置，在试件周围(试件几何中心径向的X轴负方向下部、X轴负方向上部和X轴正方向下部距试件表面均10cm的位置)设置3个火焰温度监测点，使用FLUKE多通路测温仪(测量精度0.01K)采集火焰温度。采用起爆器击发电点火头点燃航空煤油。烤燃试验现场布置如图1所示。

图1 快速烤燃试验现场布置图Fig.1 Site layout of the fast cook-off test

快速烤燃过程中各监测点的温度—时间变化曲线如图2所示。由图2可知，监测点3在点火后15s内温度达到500℃以上，其余测点温度由于环境因素影响而稍低，监测点1、2、3火焰达到稳定后火焰平均温度分别为633、538和679℃，距离油面较高的测点2温度最低。从点火开始计时，49s试件发生响应，伴有一声巨响，各监测点断路。根据现场情况可见，油池里及外部有不同程度的燃烧火光，油池底部被破片击穿，煤油泄漏，放试件的支架扭曲并损坏。由于试验环境复杂，只回收到部分破片残骸，壳体沿轴向撕裂，破片尺寸较大。

图2 温度—时间变化曲线Fig.2 Temperature-time curves

油池里的火光是没有完全反应的继续燃烧的洒落小药块，火光较亮。旁边的火光是溅飞的燃料点燃易燃物着火所致。综合判定，快速烤燃试验响应等级为爆炸反应。

3 数值模拟

3.1 快速烤燃试验仿真

用FDS软件计算航空煤油燃烧对试件的传热，建立的快速烤燃模型，包含油池、燃料和试件。由于FDS采用大涡模拟(LES)方法进行火灾计算，只计算湍流场中大尺度的涡流，将小尺度的涡流简化或者忽略。因此，网格对计算结果准确性有很大影响，为尽可能提高计算精度，计算域为5m×5m×4m。计算网格尺寸为0.02m×0.02m×0.02m。快速烤燃模型如图3所示。由于实际快速烤燃情况为开放条件，将计算域壁面和顶部设置为开放类型，将试件表面设置为绝热表面，忽略试件支架对传热的影响。

图3 快速烤燃模型Fig.3 Fast cook-off model

在FDS中采用集总组分法对燃料、空气、燃烧产物和烟气4项进行定义，空气中的氧气和氮气为反应物，燃烧产物由水、二氧化碳和氮气组成，烟气的生成量取燃料的6%[21]。设置点火点，直至火焰将试件完全包裹，快速烤燃模拟场景如图4所示，其截面云图如图5所示。火焰将试件完全包裹，试件侧表面温度保持均匀，火焰温度峰值出现在试件下方。

图4 火焰结构示意图Fig.4 Schematic diagram of flame structure

图5 不同方向的火焰结构云图Fig.5 Flame structure nephogram in different directions

设置与快速烤燃试验位置相同的3个温度监测点，通过池火数值计算所得到的试件周围火焰温度变化曲线如图6所示，火焰在点火18s后升温至500℃，测点1、2、3平均温度分别为684、567和693℃，温度误差分别为8.1%、5.4%和2.1%，与试验曲线基本吻合。

图6 火焰中测点温度变化曲线Fig.6 Temperature curves of measuring points in flame

由图6可以看出，在靠近试件壁面处的火焰温度较外层的低，且试件上方火焰温度略低于试件下方的火焰温度。这是因为热流绕试件形成湍流所致，越靠近壁面，氧气含量越低，燃烧越不充分；此外，试件上方火焰结构较为稀薄，故火焰温度较低。

对比实际烤燃试验，试件不同位置所吸收的热量不同，在模型中将试件分为6个部分(将试件侧表面等分为均匀的4份，最靠近池底编号为1，2～4顺时针排列，5、6为左右两个端面)，监测点位于每个面中心，如图7所示。同时监测试件6个位置的绝热温度平均变化，试件表面的温度变化如图8所示。

图7 试件模型与测点示意图Fig.7 Schematic diagram of the test specimen model and measuring points

图8 壳体不同位置温度变化曲线Fig.8 Temperature curves at different positions of the shell

比较不同位置处的温度曲线，并将燃料点火后火焰的变化分为3个阶段：点火阶段、发展阶段和稳定阶段。点火阶段，火苗由点火点迅速蔓延至整个油面并开始逐渐上升，在火焰的作用下试件下部温度首先升高，试件侧壁和端面温度上升较后，而试件上部温度几乎没有升高。

发展阶段，火焰开始波动，液体燃料蒸发，燃烧速率加快，火焰逐渐将试件包裹其中，此时试件侧表面温度均迅速增加。由于试件端面与气流流动方向平行，所接受的热量不稳定，且试件轴向与池宽度方向平行，两端火焰结构较稀薄，故温度波动较大。而上表面温度较稳定地持续上升。

稳定阶段，火焰将试件完全包裹，火焰的变化仅有自身脉动，各表面温度趋于动态稳定，但试件上测点温度较下面两测点温度低。这主要是由于试件放置位置距离油面较高所致。由此可见，试件放置方向与位置高低是影响表面温度的重要因素。

油池火焰中未充分燃烧的C原子等形成的碳烟颗粒是辐射热的主要释放体。试件表面各个部位的总热通量以及所占的辐射热通量与对流热通量比例分别见图9和图10。

图9 试件表面不同位置处的总热通量变化曲线Fig.9 Variation curves of total heat flux at different positions on the test specimen surface

图10 试件表面不同位置热通量占比Fig.10 Percentage of the heat flux at different locations on the test specimen surface

从图9和图10可知，在点火阶段，试件下表面总热通量与辐射热通量均为最高，试件上表面还未受火焰作用，其总热通量与辐射热通量均最低，但相对而言对流热通量占比较高；当火焰引起空气绕试件流动后，上表面的对流作用强于辐射的作用；火焰达到稳定阶段前，火焰的高度在变化，试件表面各个部位的热通量与温度均在升高，且上表面仍是对流热通量占主要作用；随着烟气上升，辐射热通量占比逐渐增加。当火焰逐渐包裹试件进入稳定阶段后，烟气继续上升，表面各个部位的温度与热通量均趋于稳定，总的辐射热通占比为91%左右，与Faghri M[22]以JP-8为燃料的快速烤燃研究结果、辐射热通量占比90%相近。因此，辐射热通量是快速烤燃的主要热源。

3.2 快速烤燃试件的传热模拟

基于上述第二部分的试验，在FLUENT软件中建立试件有限元模型，模型尺寸与试验尺寸一致，网格尺寸为0.5mm×0.5mm×0.5mm，将烤燃试件绝热表面温度作为CFD计算的边界条件，分别赋予不同温度，将试件与炸药接触面设置为耦合边界条件，并设置接触热阻。B炸药参数按照炸药组份比例组合所得，其密度ρ为1690g/cm3，比热容C为1322.4J/(kg·K)，导热系数λ为0.23W/(m·K)。B炸药达到熔化点时开始熔化，使用Boussinesq近似简化Navier-Stokes方程来模拟药柱内部的自然对流，进而对试件内部炸药进行快速烤燃模拟，直到炸药装药发生点火。

仿真结果显示，油池点火51.4s后试件发生点火响应，与试验的49s点火相比，误差为4.9%。药柱表面温度变化云图如图11所示。

图11 药柱表面温度变化云图Fig.11 Nephogram of the temperature variations on the charge surface

从图11可知，12s时火焰处于发展阶段的前期，药柱上表面温度明显低于其余部位温度，表面最高温度低于81.3℃，没有相变发生。随着火焰的发展，24s时，药柱表面温度升高并达到熔化温度，即表层炸药熔化，最高温度区域在底部及两端面棱角处，最高温度低于177.3℃，最低温度仍然在药柱上表面较窄的轴向区域；36s时药柱表面温度持续提升，温度分布情况依然如上，但是最高温度出现在药柱两端面的棱角处，达到201.8℃以上，RDX发生分解反应；51.4s时发生点火，点火区域发生在两端面棱角处。

点火时，药柱纵剖面温度分布云图如图12所示。从图12可知，药柱只有表层很少的部分发生熔化，由于熔化量极少，熔化部分几乎监测不到流动速度；通过多步反应机制计算，点火时B炸药中RDX与TNT分解的质量百分比分别为0.01%和0.0001%。与徐瑞[15]所做的同等状态的慢速烤燃相比差距很大，说明快速烤燃点火时炸药分解质量远远小于慢速烤燃。如此少的炸药分解产生的气体压力不足以使壳体破裂，可见炸药发生响应主要是表层棱角处炸药达到爆发点所致。点火时，药柱内部绝大部分区域的温度低于33.1℃，这是由于炸药热导系数很低，热传导能力较差的缘故。可见快速烤燃响应时，药柱表层较薄的范围内温差较大，会产生较大温度应力，这或许也是快速烤燃时表层炸药发生点火的因素之一。

图12 点火时刻药柱温度分布Fig.12 Temperature distribution of the charge column at the ignition time

由于B炸药的熔点在81℃左右，通过监测计算，发现药柱表面不同位置处的熔化层厚度不同，下表面熔化最多，厚度最厚；其次是两侧表面，熔化层最薄的是药柱的两端面。药柱各温度边界处的熔化层厚度分别为：1号界面2.19mm；2号界面2.01mm；3号界面1.82mm；4号界面1.93mm；5号界面1.28mm；6号界面1.30mm。

3.3 油池尺寸增加对快速烤燃的影响

为研究池火对快速烤燃的影响，使用尺寸更大的相同长宽比的油池(1500mm×1200mm×130mm)进行Φ76mm试件的快速烤燃数值模拟，试件高度距油面距离与试验相同，3个测点位置也与试验一致。仿真所用的物性参数同上。图13为3个测点的火焰温度—时间曲线。

图13 3个测点火焰温度—时间曲线Fig.13 Temperature—time curves of three measuring points

由图13可见，当使用较大的油池，火焰升温速率更快，火焰平均温度比小油池的更高，且3个测点所测温度值更加接近，均在830℃左右，比小油池温度提升近140℃以上。

图14为试件6个面(划分与前面相同)的温度—时间曲线。由图14可见，大尺寸油池，火焰发展期的时间缩短，稳定期的时间延长，试件6个表面的温度也更接近，且温度脉动幅度减小。可见，快速烤燃时，采用尺寸较大的与烤燃试件相匹配的油池更合理，且烤燃试件的安放位置距油面的高度也可以提升。

图14 大油池快速烤燃时试件表面温度Fig.14 Surface temperature of the test specimen during the fast cook-off in a large oil pool

图15为总热通量、辐射热通量随时间变化的曲线。由图15可见，油池尺寸增大，总热通量提高，且辐射热通量的占比也随之提高。

图15 试件不同位置热通量变化曲线Fig.15 Curves of heat flux variations at different positions of the test specimen

通过模拟计算可知，Φ76mm烤燃试件使用较大的油池快速烤燃，火焰发展期所用的时间比小油池缩短10s左右，响应时间为46.2s，比小油池缩短5.2s，点火区域仍位于两端面的棱角区域。

4 结 论

(1)将FDS与CFD组合模拟计算的快速烤燃更能真实地体现火焰结构的变化特性与边界条件的差异，计算的温度—时间曲线与试验曲线吻合，且能获得快速烤燃过程中辐射热通量与对流热通量的变化特性；装填B炸药的试件进行快速烤燃，响应时刻药柱只有表层很少部分炸药熔化，绝大部分药温仍然是常温，响应烈度为爆炸反应。

(2)快速烤燃过程中，辐射热通量的贡献起主导作用，占比为91%左右；火焰温度随着油池尺寸的变化而变化，油池尺寸越大，火焰温度越高，火焰达到稳定阶段的时间越短，辐射热通量的作用也越大，对于确定的试件而言，响应时间也越短。