偏二甲肼池火灾热动力学特性数值模拟

史宇航, 胡宇鹏, 李明海, 郭历伦

(中国工程物理研究院, 总体工程研究所, 绵阳 621999)

液体推进剂具有高比冲、推力大、能多次点火和停运等特点,在大型运载火箭的推进系统和各种航天器的姿态和轨道控制系统中得到广泛应用[1],目前较多使用的液体推进剂燃料主要为偏二甲肼(UDMH)。UDMH在生产、运输、贮存、加注等任务剖面中可能发生意外泄漏,泄漏到地面遇到防火堤或低洼边界时容易积聚形成液池,遭遇火源或氧化剂后极易燃烧形成池火灾,往往会造成人员和经济的重大损失。已有统计,近60年来国内外共发生推进剂突发燃爆事故300多次,伤亡400多人[2]。因此,研究UDMH池火灾热动力学特性对于液体推进剂在全寿命周期内的安全应用具有重要意义。

目前,关于池火灾的研究主要集中在化工等领域涉及的如煤油、乙醇等油料池火灾。实验研究方面,程海涛等[3]开展了直径为1 m、宽度为10～60 cm 的初始线性火源的秸秆燃烧实验,发现火焰蔓延速率与秸秆宽度成正相关关系,火焰温度呈现“快速升高-缓慢降低”的趋势,秸秆表面温度超过500 ℃。Zhu等[4]研究不同压力下乙醇汽油混合燃料的特性,发现随着环境压力的降低,火焰高度显著增加。Einar等[5]开展了10、20、30 cm三种尺寸的庚烷池火燃烧试验,发现边缘高度对庚烷质量燃烧速率的影响超过池火直径的影响。Tian等[6]研究了通风条件下甲醇池火灾的燃烧速率变化,发现对流换热对燃烧速率的提高起着重要作用,甲醇初始温度越高,燃烧速率越高。针对油池火辐射特性,May等[7]提出了点源模型,认为液化天然气池火的热辐射和距离的平方成反比,总的热辐射大概占燃烧总释放能量的15%。之后Mcgrattan等[8]研究提出了更适用于大尺寸油池火辐射特性的正圆柱固体壁面模型。数值模拟方面,Stewart等[9]和Ahmadi等[10]利用火灾动力学模拟工具(fire dynamics simulator,FDS)模拟乙醇和煤油的池火灾,在质量燃烧速率和热辐射方面的模型计算结果与池火灾的实验数据之间的误差在8% 以内。吕辰等[11]通过数值模拟技术准确的得出凹型外立面火灾竖向蔓延烟囱效应以及火焰前锋高度随时间符合指数函数变化关系,据此提出防火挑檐的合理设置可有效阻止火势沿凹型外立面向高处相邻住户蔓延。可见,关于UDMH池火灾特性的研究尚未见相关报道。

近年,已有关于UDMH燃烧的研究主要集中在液体发动机内部的燃烧特性及机理等方面,其与池火在燃烧方式、体系边界等方面都有较大不同。Moghaddam等[12]开展了特征长度对凝胶UDMH/IRFNA双推进剂燃烧性能影响的试验,结果表明凝胶双推进剂性能随特征长度增加的提升效果较低,最佳性能特征长度为162 cm。尹婷等[13]开展了压力振荡环境下UDMH液滴燃烧过程实验,发现增加燃烧室温度能够显著促进液滴燃烧,使得液滴燃烧速率和表面燃烧流率迅速增加。Ma等[14]对亚临界和超临界压力下的UDMH液滴的燃烧特性进行研究,发现亚临界条件下液滴燃烧时间随环境压力的增加而迅速缩短,而超临界条件下液滴燃烧时间没有减少。

综上,目前对于UDMH池火灾缺乏深入认识,需对其热动力学特性开展研究。然而,偏二甲肼毒性较大(表1),难以通过实验这一火灾动力学研究惯用方法对UDMH池火灾进行研究。因此,拟针对UDMH池火开展系列数值模拟,分析油池尺寸和环境风速对UDMH池火灾热动力学特性的影响规律,探索火焰结构、温度分布及辐射热流密度等热特性,为液体推进剂火灾事故应急响应提供理论参考。

表1 偏二甲肼毒性

1 数值模型

1.1 物理数学模型

UDMH遭遇意外情况发生泄露,泄露到地面之后向四周流淌,在低洼地形中容易积聚形成一定厚度的液池,若遭遇火源便会形成池火灾,池火灾结构示意图如图1所示。

图1 池火灾示意图Fig.1 Pool fire diagram

当泄露UDMH向四周流淌扩展形成具有一定厚度的液池时,其面积是由泄露量和地面性质决定的,液池面积计算公式为

(1)

式(1)中:S为油池的最大可能面积,m2;W为 UDMH 泄露质量,kg;Hm为最小液层厚度,m;ρ0为 UDMH 密度,kg/m3。

把形状不规则的UDMH液池等效为边长为L的方形液池,换算公式为

(2)

质量燃烧速率的计算公式为

(3)

式(3)中:mf为燃烧速率,kg/(m2·s);c为常数,c=0.001 kg/(m2·s);Hc为液体的燃烧热,kJ/kg;Hvap为液体汽化热,kJ/kg;Cp为液体的定压比热容,kJ/(kg·K);Tb为液体的沸点,K;Ta为环境温度,K。

由式(3)得出mf=0.049 75 kg/(m2·s)。表2给出了UDMH理化参数。

表2 偏二甲肼理化参数

池火灾燃烧是典型的非预混燃烧现象,计算流体动力学(computational fluid dynamics,CFD)通过控制方程对燃烧现象进行描述,控制方程包含连续性方程、动量守恒方程、能量守恒方程、气体组分方程及理想气体状态方程,这些方程构成的数学模型是数值模拟计算的理论支撑。具体如下。

连续性方程为

(4)

动量守恒方程为

(5)

能量守恒方程为

(6)

气体组分方程为

(7)

理想气体状态方程为

(8)

式中:ρ为气体密度,kg/m3;V为速度矢量,m/s;p为压力,Pa;g为重力加速度,kg/m2;τ为黏性应力张量,Pa;h为比焓,J/kg;Di为组分i的扩散系数,m2·s;q″为体积热源,W/m3;q为辐射热通量矢量,W/m3;φ为耗散函数,W/m3;R为通用气体常数8.314 J/(mol·K);m”为生成率或消耗率,kg/m3;M为气体摩尔质量,kg/mol;Y为质量分数,%。

1.2 计算方法

采用专业火灾CFD软件FDS开展计算工作,计算中燃烧模型选用适于工程应用和大规模火灾场景的混合分数控制模型,湍流计算选用大涡模拟(large eddy simulation,LES)。计算区域为3L×5L×3L(L为油池边长)的长方体区域,建模时所有油池均采用正方形油池,油池中心设置在原点,考虑到液体燃料泄漏时无边界,因此油池直接贴附在地面,底部为混凝土地面。模型的顶部和四周在无风环境下均设置为开放边界条件,有风环境下模型左侧面设置为风速入口。

为确保计算模型的合理性,对已有航空煤油池火模型[15]进行计算对比,油池边长为5 m,池火轴向温度分布如图2(a)所示,竖直方向辐射热流密度对比如图2(b)所示,本模型与已有文献在温度与辐射热流密度的误差均在10%以内,模拟结果与文献结果基本一致,即本数值模型适用于池火灾的计算。

图2 模型验证Fig.2 Model validation

进一步开展了网格无关性验证工作,对计算网格进行了优选。以边长5 m的UDMH池火计算为例,计算结果如图3所示,网格数量分别为70万和140万时得到的火焰轴向温度相差在8%以内,可认为70万网格满足计算精度要求。

图3 网格数量对计算结果的影响Fig.3 The influence of the number of grids on the calculation results

2 结果与分析

2.1 油池尺寸对池火灾热动力学特性的影响

2.1.1 火焰形态

图4给出了L=1、6 m油池的温度场分布及氧气质量流量。如图4(a)温度切片所示,UDMH油池火羽流整体呈锥形,可以划分为三个不同的区域:连续火焰区,此处时火焰持续存在的区段,中心温度是整体火焰最高温度;间歇火焰区,此处是火焰间歇存在的区段,该区温度随轴向高度增加急剧下降;浮力羽流区,没有火焰,燃烧气体席卷周围空气而上升的区段,该区温度继续下降直至接近环境温度。在整个燃烧过程中,火焰会不断地卷吸周围空气,从图4(b)中可以看出,着火3 s后火焰底部左右两侧氧气均向油池中心流动,同时油池上方由于燃料蒸汽的扩散作用使氧气向两侧退开,着火40 s后,火焰两侧氧气从底部到顶部受火焰卷吸作用均向中间流动,连续火焰区为支持充分燃烧卷吸的氧气速率最高为0.53 kg/(m2·s)。

图4 偏二甲肼池火燃烧状态Fig.4 Combustion condition for UDMH pool fire

2.1.2 火焰温度

图5为L=1～9 m UDMH油池火轴向温度分布曲线。轴向温度曲线可划分为温度上升区、温度不变区和温度下降区。其中温度上升区对应于火焰根部,靠近油池液面的位置由于氧气不足,燃料燃烧不充分,燃料蒸发的能量主要来源于外部火焰热反馈。随着高度增加,火焰卷吸周围空气使燃料充分燃烧,温度快速上升达到最高温度并在一定范围内基本不变是为温度不变区,如9 m边长油池火焰中轴线1～5 m范围内温度都在1 100 ℃以上。油池边长L从1～9 m逐渐增大,油池火焰的最高温度从1 006.23 ℃逐渐增大到1 160.92 ℃,温度最高点的位置随油池直径增加逐渐远离油池液面。之后进入温度下降区,燃料燃烧结束且外界大量冷空气卷入,温度逐渐下降,在火焰顶端,温度下降到 600 ℃ 以下,温度下降的速率随油池边长增大而逐渐降低。

图5 不同边长油池火轴向温度分布Fig.5 Axial temperature distribution of oil pools with different side lengths

图6给出了L=1～9 m油池中距离液面 0.6、2 m 的横向温度分布情况。0.6 m高处横向温度分布呈现先升高后降低的趋势,原因是火焰底部中心位置氧气不足,燃烧不充分,之后在距离到达火焰边缘之前随着燃料充分燃烧温度会逐渐升高直至最高温度,此处各油池横向最高温度为 1 021.00、1 039.40、1 067.40、1 109.70、1 157.10 ℃,离开火焰范围后随距离增大温度逐渐降低。

图6 不同高度处横向温度分布Fig.6 Radial temperature distribution at different heights

2 m高处各边长油池的横向温度最高点均在轴线上,分别为789.90、1 005.30、1 063.60、1 098.60、1 155.30 ℃,均低于接近液面的0.6 m处。同时可知,在超出油池边界之后,火焰温度均快速下降至200 ℃以下。以800 ℃为基准能得到0.6 m高处不同边长油池火轴心到火焰边缘的距离分别为:0.19、0.38、0.64、1.32、2.29 m,基本与油池边长呈正比关系,随高度增大火焰宽度逐渐减小,符合火焰整体上窄下宽的形态。

2.1.3 辐射热流密度

开放空间池火对周围人员及设备的破坏机理以热辐射为主。图7给出了L=1～9 m油池火的轴向及横向辐射热流密度分布规律。轴向热辐射在接近油池液面的位置较小,之后呈现先增大后减小的规律,同时可以看出其中各条曲线增长和减小的斜率较接近,说明不同边长油池火轴向辐射变化趋势一致。在同一高度上,辐射热流密度随油池边长的增大而增大。辐射热流密度的最大值随油池边长的增大从109.4 kW/m2增大至266.8 kW/m2,位置逐渐远离油池液面。图7(b)给出了距液面高0.6 m, 横向距中心3、5、10、15、20 m五个不同距离处辐射热流密度变化规律,可以看出UDMH油池火辐射热流密度随着距离增加而逐渐衰减,衰减速率逐渐变慢。同一距离的辐射热流密度随油池边长的增大而增大。

图7 不同边长油池火辐射热流密度分布Fig.7 Distribution of radiant heat flux of oil pools with different side lengths

2.2 风速对池火灾热动力学特性的影响

2.2.1 火焰形态

图8为边长6 m的UDMH油池在不同风速下火焰中心剖面的温度分布。在环境风作用下,火焰会倾斜并拉长,风速较大时火焰会部分贴附地面。轴向温度由于脱离火焰而逐渐降低,下风侧温度逐渐升高。如图8所示,风速为1 m/s时,火焰没有倾斜,温度场关于轴线呈对称分布;风速达到1.5 m/s后,接近油池液面的底部火焰未受影响,但顶部火焰已经发生倾斜,火焰高度降低为8.01 m,火焰倾角为24.5°;风速增大至2 m/s时,火焰部分贴附地面,火焰高度为6.16 m,火焰倾角达到46.3°;风速进一步增大至4 m/s时,火焰大部分贴近地面,只有火羽流顶端抬离地面且向下风侧倾斜,火焰高度降低至4.23 m,火焰倾角增大至82.3°。

图8 不同风速下油池火温度场分布Fig.8 Temperature distribution of pool fire under different wind speed

2.2.2 火焰温度

图9给出了不同风速下UDMH油池火轴向温度变化规律。当风速较小为1 m/s时,轴向温度分布规律与无风时较接近,最高温度能达到1 100 ℃,但有风环境加快了热量散失,导致较高位置温度降低。风速增大到1.5 m/s后,轴向温度变化趋势发生显著变化,呈单调下降趋势,表明油池火焰已经发生倾斜,偏离了油池中心,与图8中火焰中心剖面温度分布情况一致。风速超过1.5 m/s后,轴向温度单调降低的速率随风速增大而增大。

图9 不同风速下油池火轴向温度分布Fig.9 Axial temperature distribution of pool fire under different wind speeds

有风环境下UDMH油池火横向温度也会发生显著改变。图10为距油池液面h=0.6 m、2 m处横向温度变化曲线。在距液面较近的0.6 m高处,当风速较低时,火焰还未发生大幅度倾斜,横向温度变化与无风时较一致;风速达到2 m/s后,火焰出现明显倾斜,轴心处的温度显著降低,0.6 m高处轴心温度由919.8 ℃下降至605.9 ℃,下风侧温度则逐渐升高,横向温度分布变为单调升高趋势;风速达到4 m/s后,火焰大部分脱离油池范围,不同位置处温度明显降低。在距离液面较远的2 m高处,在风速作用下比0.6 m处温度变化更明显,风速为 1.5 m/s 风速时中心温度由1 125.6 ℃下降至485.7 ℃,温度分布为先升高后降低;风速增大后,横向温度同样变为单调升高,但各点温度值低于0.6 m处;风速达到4 m/s后,2 m高处横向温度维持在100 ℃不变,此时火焰已经大部分贴附地面。

图10 不同高度处油池火横向温度分布Fig.10 Radial temperature distribution of oil pool fire at different heights

2.2.3 辐射热流密度

与温度不同,图11(a)给出的轴向辐射热流密度在1～4 m/s风速作用下呈现先增大后减小的趋势,与无风环境一致,轴向辐射热流密度最大值随风速增大逐渐降低。横向辐射热流密度分布受风速影响较大,如图11(b)所示,风速较低为1～1.5 m/s 时距油池中心不同距离的辐射热流密度较低,且随着距离增加而逐渐降低;风速达到2～3 m/s 时,距油池中心3 m和5 m的位置由于火焰的倾斜辐射热流密度增大,最高达到168.33 kW/m2;风速增大到4 m/s时,温度场分布显示火焰倾角进一步增大,火焰贴近地面,此时横向辐射热流密度分布趋势变为先升高后降低,5 m处接受的热通量升高到119.08 kW/m2,通过图8中温度场分布可以看出5 m位置位于火焰内部。

图11 不同风速环境下辐射热流密度变化规律Fig.11 Distribution of radiant heat flux under different wind speeds

3 工程关联式

针对推进剂燃料在运输、发射等任务剖面可能遭受的火灾事故环境,需建立相应工程关联式对UDMH火灾事故的危险性进行理论预测,为火灾蔓延的控制和救援提供理论指导[16]。基于上述数值模拟结果拟合得到适用于UDMH池火灾的工程关联式,如表3所示。

表3 UDMH池火灾工程关联式

表3所示火焰高度为采用间歇率的方法定义的平均火焰高度,将火焰间歇率I=0.5时对应的火焰高度值作为平均火焰高度,即在所有观察次数中,火焰高度超过该高度值的次数占总观察次数的一半。无风环境下,不同尺寸油池火的平均火焰高度分别为2.67、4.43、6.11、10.03、13.22 m。图12给出了无量纲火焰高度的工程关联式与模拟数据的对比结果,偏差为1.09%。

图12 无量纲火焰高度与Fc2/3的关系Fig.12 Relation of dimensionless flame height to Fc2/3

图13给出了有风环境下火焰高度和火焰倾角的工程关联式与模拟数据的对比结果,偏差分别为4.23%和11.37%。不同环境风速下6 m UDMH池火的平均火焰高度分别为9.85、8.01、6.16、5.61、4.23 m。图14给出了无风环境下横向辐射热流密度的工程关联式与模拟数据的对比结果,偏差为14.75%。由此可知本文所建立的工程关联式可以对UDMH池火灾的热动力学特性进行较好预测。

图13 无量纲火焰高度及火焰倾角与U的关系Fig.13 Relationship between flame height and flame inclination angle and U

图14 6 m油池火横向辐射热流密度变化Fig.14 Distribution of transverse radiant heat flux in 6 m pool fire

4 结论

(1)无风环境下1～9 m偏二甲肼油池火焰高度及温度随油池尺寸增大而增大;火焰轴向最高温度从1 006.23 ℃逐渐增大到1 160.92 ℃,温度最高点的位置逐渐远离液面;横向温度在0.6 m高处呈先升高后降低趋势,2 m高处则呈单调下降趋势。

(2)有风环境下6 m偏二甲肼油池火焰高度随风速的增大而逐渐减小,火焰倾角从0°逐渐增至82.3°,火焰整体逐渐贴近地面;轴向温度分布在风速达到1.5 m/s后由先升高后降低变为单调下降,且风速越大下降速率越快,横向温度受风速影响中心温度降低,下风侧油池边缘温度升高。

(3)无风环境下,偏二甲肼轴向辐射热流密度随油池边长增大由109.32 kW/m2增至266.85 kW/m2,横向辐射热流密度随距离增大单调降低;有风环境下,随风速增加同一位置轴向辐射热流密度逐渐降低,但整体变化趋势没有改变,下风侧辐射热流密度由于火焰倾斜显著增加。

(3)基于数值模拟结果拟合获取偏二甲肼油池火的工程关联式,可较好地预测无/有风环境下偏二甲肼池火灾的火焰高度(偏差在4.23%以内)、火焰倾角(偏差在11.37%以内)以及距火焰不同距离处的辐射热流密度(偏差在14.75%以内),可为液体推进剂火灾事故应急响应提供理论参考。