气体泄爆压力分步计算模型及其湍流修正*

孙 松，高康华，邱艳宇，王明洋

（陆军工程大学爆炸冲击防灾减灾国家重点实验室，江苏 南京 210007）

可燃性气体因其来源广、热能高、污染小等特点已经成为人们生产、生活中最常用的能源之一；与此同时，频繁发生的气体爆炸事故也给人们的生命财产安全带来了巨大威胁。泄爆是工程实际中最常用的防护措施，研究者们对于泄爆条件下空间内的爆炸压力进行了大量的实验研究[1-5]，对气体爆炸的物理过程及各影响因素有了较全面的了解。

同时，出于安全性、经济性和实用性的考虑，研究者们以大量实验为基础结合气体爆炸物理过程进行了理论模型的研究，以期能够预测空间内气体爆炸压力的变化情况。Bradley 等[6]根据守恒方程、理想气体状态方程和等熵压缩方程推导出关于升压速率的偏微分方程组，成为其后理论模型发展的基础。Han 等[7]通过对比模型计算值与实验结果发现，压力计算模型可以较好地描述建筑空间内的压力变化，随后通过模型计算探讨了泄放面积、泄放压力、屋室尺寸、气体浓度和点火位置对爆燃压力的影响。Ugarte 等[8]根据不同点火位置、容器形状在一维等熵泄放模型的基础上对火焰面积、火焰传播速度和外部爆燃压力进行修正，对比发现修正后的模型与实验结果吻合更好。美国防火规范NFPA68-2013[9]中总结大量实验数据得到不同工况下泄爆压力及泄爆面积的简化确定方法，为实际工程防护提供了重要指导。Sustek 等[10]使用常用经验与半经验公式对不同实验工况进行计算，对比发现，由于各经验公式中的参数很多都是基于实验数据拟合得到，因此各公式适用范围较窄，使用经验公式计算不同工况可能会产生较大误差。

现有的气体爆燃压力计算方法多是以大量实验数据为基础评估空间内的最大峰值超压，无法准确描述空间内的压力变化过程；数值模拟方法计算成本高，耗时长，不利于工程快速评估应用；而常用的压力计算模型[11-13]多为关于升压速率的偏微分方程组，其求解过程较复杂，同时这些模型大多假设泄放构件开启后先泄放已燃气体，再泄放未反应气体，无法计算同时泄放已燃气体和未反应气体的工况。本文中，通过将爆燃过程进行微分，假设每个微小时间段内爆燃泄放过程均按照先燃烧、再泄放、最后压力平衡的过程独立分步进行，得到了泄爆压力分步计算模型；同时在尺寸为2 m×1.2 m×0.6 m 的爆炸腔体一端安装不同面积的泄压结构进行泄爆实验，对分步压力计算模型进行验证，并根据实验结果对模型进行湍流修正。该模型清晰地反映了气体爆燃泄放机理，通过对微分时段内的分步计算过程进行迭代，可以较好地描述泄爆条件下空间内的压力及各状态参数的变化过程，计算过程简单，可为工程实际中快速评估内爆炸压力提供参考；经过湍流加速修正的压力计算模型可以较好地描述因泄放截面改变引起的湍流区域内的压力变化情况。

1 理论模型

1.1 基本假设

由于爆燃过程非常复杂，涉及到燃烧学、热动力学、流体动力学等多方面内容，因此在模型计算时假设可燃气体密度均匀充满密闭空间且符合理想气体状态方程；同时假设密闭容器为绝热壁面，忽略热交换及内部反应的能量耗散过程。

泄放条件下气体爆燃的超压变化过程受燃烧产生的压力和泄放降低的压力共同影响，变化过程较复杂，为简化计算将燃烧过程和泄放过程单独考虑，将气体爆燃泄放过程分成3 个独立阶段：（1）只燃烧不泄放(confined combustion stage)。将此阶段简化为密闭空间内的爆燃过程，不考虑泄放对于爆燃压力的影响。（2）泄放过程(venting stage)。第1 阶段燃烧结束后室内压力上升，高于外部环境压力，此时燃烧产物和未燃气体将分别发生泄放。（3）泄放后压力平衡阶段(pressure equilibrium stage)。该阶段泄放过程结束，此时空间内部燃烧产物与未燃气体压力不同，空间内部会发生气体流动最终达到压力平衡。密闭空间内气体爆燃泄放过程如图1 所示，图1(a)～(c)为燃烧反应阶段，图1(d)～(e)为泄放阶段，图1(f)为压力平衡阶段。

图 1 爆燃泄放分步计算模型示意图Fig.1 Diagrams of the sub-steps model for combustion venting

1.2 控制方程

1.2.1 燃烧阶段

通常假设预混气体以规则球形燃烧并向四周扩散，则在 Δt时间内所燃烧的体积为：

式中： VB,j为第j 个时间段内燃烧的气体体积；S 为火焰传播速度； Vh,j-1为j-1 时刻已燃气体体积。

现定义燃尽率 Δµj为 Δt时间内燃烧的气体质量与可燃气体总质量之比，则：

式中： ρc,j-1为 j-1 时刻未燃气体密度， V0为 爆炸空间总体积， ρ0为初始气体密度。

又因为根据等熵压缩方程，在绝热过程中有：

式中： γu为未燃气体绝热指数， pj-1为j-1 时刻空间内的压力，p0为空间内的初始压力。

则 Δµj可以写成：

由于空间内的压力随气体燃烧分数的增大而升高，假设室内压力的升高与燃烧产物生成量成线性关系[14]，则图1 中(c)阶段室内压力为：

式中： pb为 密闭空间内气体爆燃产生的最大压力，该值与气体种类和浓度有关； pn,j为j 时刻密闭燃烧阶段结束时内部压力值。

由于整个反应过程为绝热过程且压力 pn,j已知，则根据等熵压缩方程：

可得图1 中j 时刻(c)阶段未燃气体体积为：

已燃气体体积为：

式中： Vc,j-1为 j-1 时刻未燃气体体积； Vh,j-1为j-1 时刻已燃气体的体积。

由质量守恒可得图1 中(c)阶段未燃气体密度为：

式中： ρn,c,j为j 时刻燃烧阶段结束时未燃气体密度, ρc,j-1为j-1 时刻未燃气体密度。

已燃气体密度为：

式中： ρn,h,j为j 时刻燃烧阶段结束时燃烧产物密度, ρh,j-1为j-1 时刻燃烧产物密度。

此时第j 个时间段内燃烧阶段结束，开始进行泄放阶段的压力计算。

1.2.2 泄放阶段

考虑实际情况，假设该阶段未燃气体与燃烧产物分别发生泄放且外界气压恒为标准大气压，则以次临界速度（ p0/ pj-1≥βk）泄放时[6,14]，单位时间内泄放的燃烧产物和未燃气体体积分别为：

以临界速度（ p0/ pj-1≤βk）泄放时，单位时间内泄放的燃烧产物和未燃气体体积分别为：

式中： Ah,j和 Ac,j分别为 Δt 时间内燃烧产物和未燃产物的泄放流动面积； γb为 燃烧产物的绝热指数； βk为临界比值，其计算公式为 βk=[2/(γu+1)]γu/(γu-1)。

由于发生气体泄放，未燃气体和燃烧产物的密度及压力均发生变化，根据等熵压缩方程可得图1 中(e)阶段j 时刻泄放阶段结束时空间内未燃气体和燃烧产物的压力分别为：

j 时刻泄放阶段结束时空间内未燃气体和燃烧产物的密度分别为：

在第2 阶段中假设只发生气体的泄放流动，但空间内未燃气体与已燃气体所占体积不发生改变，此时第2 阶段结束，开始进行空间内压力平衡阶段的计算。

1.2.3 压力平衡阶段

由等熵压缩方程可知:

将式(19)代入式(20)并变形后得完整反应过程结束时空间内的压力：

式中： β= pr,j/pj。

求得图1 中(f)阶段的压力pj后，根据状态方程与等熵压缩方程可以求得(f)阶段的气体物理参数：

式中： ρh,j和 ρc,j分别为j 时刻整个爆燃泄放过程结束后空间内未燃气体和燃烧产物的密度； Vh,j和 Vc,j分别为此时空间内未燃气体和燃烧产物的体积。

此时第3 阶段结束，一个微分时间段内腔体中气体压力、密度、体积的完整变化过程结束，开始迭代进行下一个时间步的计算。

反应后期若空间内混合气体已全部燃尽，则燃烧过程结束，此时气体继续由内部流向外部空间，直至内外压力平衡。该阶段在迭代计算过程中省略图1 中(b)、(c)过程，即该情况下完整反应流程图变为图1(a)、(d)、(e)、(f)。

2 实验验证

2.1 实验装置

实验装置由配气系统、长方体爆炸容器、点火系统、泄爆系统与数据采集系统构成，如图2 所示。

图 2 实验系统示意图Fig.2 Schematic of the experimental system

爆炸腔体为自行设计的长方体容器，长、宽、高分别为2、1.2、0.6 m，最大可承受1 MPa 的压力。爆炸容器一端为开敞的用于安装泄爆结构的卡箍式法兰盘，其尺寸为1.2 m×0.6 m，泄爆构件示意图如图3所示；另一端壁面中部通过法兰结构与电阻丝点火器连接，构成点火系统，其最大点火温度可达500 ℃。实验中选用乙烯作为可燃气体，通过在爆炸腔体侧面等距离设置进气孔和抽气孔运用压力分配法配气，气体浓度由流量计和压力表控制。进气孔内配置有细长进气管道，管道上均匀分布有小孔，这种均布直管式进气系统可以保证腔体内预混气体混合较均匀。在爆炸腔体顶部装有PCB113B26 系列压电式高频压力传感器，分别距泄放端口1.2、0.1 m。数据采集系统将采集到的电信号转化为压力信号并输出，其数据采样频率为200 kHz。高速摄影仪置于距泄爆口侧面8 m 处，拍照频率为1 000 s-1。

图 3 泄爆构件示意图(单位为mm)Fig.3 Schematic diagrams of vent covers (unit in mm)

如图3 所示，实验中泄爆构件采用击穿压力均为7.5 kPa 的泄爆板和泄爆膜。泄爆板泄放面积为1.1 m×0.5 m，泄爆膜泄放面积为0.3 m×0.6 m。采用泄爆膜作为泄爆构件开展实验时，在腔体泄爆端安装尺寸为1.2 m×0.6 m 的隔板，并在隔板中部开设面积为0.6 m×0.3 m 的泄爆口，使用泄爆膜对其进行封闭，如图3(a)所示；采用泄爆板作为泄爆构件开展实验时，直接在腔体泄爆端安装尺寸为1.2 m×0.6 m 的泄爆板，如图3(b)所示。构件击穿压力通过在泄爆结构背面粘贴应变片以应变信号消失时刻腔体内对应的压力确定，实验时选取的可燃气体的体积分数为4%～8%。

2.2 结果对比

使用分步压力计算模型对实验工况进行计算，由于实验工况为矩形容器端部中心点火，因此需针对火焰的传播形状对计算模型进行修正。由高速摄影观察到火焰由泄爆端流出时其阵面可近似认为是平面，因此为计算简便假设火焰未接触壁面时为球形传播，接触壁面时为平面传播，则端部点火时火焰的体积表达式为：

式中：2a、2c 为爆炸腔体的长和高，x 为火焰半径。

泄爆板和泄爆膜工况的计算值及实验测量结果分别如图4～5 所示。由图4 可知：采用泄爆板作为泄爆构件时，2 个传感器测得的压力时程曲线基本重合，符合准静态压力形式[15]；此时腔体内部为单峰值压力曲线，分析认为当构件击穿后由于泄放面积较大使得短时间内大量气体外泄，引起腔体内部压力下降至环境压力。该情况下模型计算值与实验结果吻合较好。

图 4 可燃气体的体积分数不同、采用泄爆板作为泄爆构件时，泄爆压力的计算值与实验测量结果的比较Fig.4 Comparison of gas explosion venting pressure between calculated results by the model and experimental ones in the experiments with vent plates at different volume fractions of combustible gas

图 5 可燃气体的体积分数不同、采用泄爆膜作为泄爆构件时，泄爆压力的计算值与实验测量结果的比较Fig.5 Comparison of gas explosion venting pressures between the model and the tests with vent films at different volume fractions of combustible gas

观察图5 发现，当采用泄爆膜作为泄爆构件时腔体内压力曲线均为双峰值曲线，分析认为由于泄爆膜面积较小，构件击穿后腔体内部气体外流引起压力下降产生第1 个压力峰值；气体外泄同时又会加剧腔体内部湍流扰动，扭曲火焰阵面，提高燃烧速率与升压速率，引起腔体内压力再次上升产生第2 个峰值[1]。同时由图5 发现当乙烯体积分数较高时靠近泄爆口的压力传感器T1 测得的第2 个峰值压力要大于内部传感器T2 相应的测量值，分析原因认为由于泄爆膜开口面积较小，外流气体流动面积的突然改变、结合容器壁的限制作用使得近泄爆口处湍流效应明显，可燃气体得到充分燃烧，腔体内部产生压力梯度，泄爆口处压力大于腔体内部压力，且越接近最佳体积分数该现象越明显[16-17]。此时理论模型可以较好地描述腔体内部的压力变化情况，但与近泄爆口处压力变化情况存在误差。

为了进一步验证计算模型的适用性，运用该模型计算其他研究者的实验工况[5,7,18]，其结果如图6所示。通过对比发现计算结果与实验结果吻合较好，该方法具有较好的普适性。

图 6 运用本文计算模型计算其他研究者的实验工况[5, 7, 18]下的泻爆压力Fig.6 Gas explosion venting pressures calculated by the model proposed in this paper for the experimental conditions[5, 7, 18] of other researchers

3 模型修正

由于传感器T1 测得的第2 个峰值压力大于腔体内部相应的压力值是因泄放截面改变引起的火焰湍流传播引起，因此在此引入湍流加速因数λ 对火焰传播速度进行修正，以计算靠近泄放口位置的压力情况。湍流加速作用因泄放截面的改变引起，因此将传播路径中的障碍物截面与火焰传播截面的比值，即阻塞率, 作为其增强因子的影响因素。研究者们对障碍物对火焰传播的加速作用进行了大量的实验研究[16,19-20]，虽然实验环境不同，但是障碍物阻塞率对于火焰的湍流加速效果表现出一定的规律性。现根据美国防火规范NFPA68 确定因截面改变引起的湍流加速因数[9]：

式中：Aobs为障碍物截面面积，As为火焰传播路径截面面积， λ0为因自湍流等因素引起的湍流因数。在本文的泄爆膜实验工况中Aobs取0.37， As取0.55， λ0取2。

则未因传播截面改变引起湍流加速的火焰传播速度为[21]:

近泄爆口处因截面改变产生湍流加速的火焰传播速度为：

式中： S0为基本火焰传播速度；p 和p0分别为空间内的压力和初始压力；T0和Tu分别为初始温度和未燃气体温度； α 和 η为 根据实验确定的参数，此处 α取2，η取-0.15。

现将修正后的火焰传播速度代入压力计算模型，则计算结果如图7 所示。由图7 可知，经湍流修正后的压力计算模型得到的峰值压力及升压时间与传感器T1 测得的实验结果吻合较好，使用湍流修正因数预测因传播截面改变引起的湍流区域的压力值具有可行性。

图 7 针对体积分数不同的可燃气体，修正后的压力计算模型得到的泻爆压力与实验结果的比较Fig.7 Comparison of gas explosion venting pressures calculated by the modified model and the test data at different volume fractions of combustible gas

4 结 论

（1）通过将爆燃泄放过程进行微分，假设每个微小时间段内爆燃泄放过程均按照先燃烧、再泄放、最后压力平衡的过程独立分步进行，得到了泄爆压力分步计算模型。该模型计算简单，不仅可以计算爆燃泄放过程中的压力变化情况，还可以得到各时刻空间内的气体密度、体积等物理状态参数。

（2）在尺寸为2 m×1.2 m×0.6 m 的长方体容器中充满乙烯-空气混合气体进行实验验证，发现当采用大面积泄爆板作为泄爆构件时爆燃压力曲线均为单峰值曲线，该条件下压力模型计算结果与实验结果吻合较好。

（3）当采用小面积泄爆膜作为泄爆构件时会因外流气体流动截面改变引起火焰湍流加速，使得腔体内产生压力梯度，近泄爆口处压力值大于腔体内部。通过湍流加速因子对火焰速度进行修正，则修正后的压力计算模型可以较好地描述近泄爆口处的压力变化情况。