1 前言

燃气发生泄漏后，在受限空间内与空气混合会形成爆炸性混合物，遇点火源有燃烧爆炸危险。燃气管道泄漏周围空间总会存在家具、装置等障碍物，会扩大燃爆事故后果^[1]，而门窗、通风口等又会存在泄爆效应，降低事故危害。因此有障碍物和泄爆口同时存在的受限空间可燃气体燃爆机理^[2]一直是研究的热点。近年来国内外学者对可燃预混气体燃烧爆炸开展了大量的实验和数值模拟研究。Masri等^[3]采用大涡模拟与动态Germano火焰面模型耦合方式对丙烷/空气在障碍物形状和阻塞率变化空间内火焰传播过程进行模拟。Veynante^[4]，Nogenmyr^[5]，Knudsen^[6]等根据湍流尺度与火焰面厚度的相对尺寸大小对火焰面模型进行细化修正，以确保对火焰传播特性精确预测。余立新等^[7]研究发现在半开口管道中障碍物对火焰速度和压力提升作用显著。Wan等^[8-10]研究了通风口位置、大小及障碍物与通风口相对位置对管道内爆炸特性的影响。王世茂等^[11]研究了开口率和点火源类型对汽油蒸气泄压爆炸超压的影响。Chao等^[12]对带有排气口的小空间内可燃气云进行排气爆炸试验，探究压力变化的影响因素。上述工作主要集中在障碍物阻塞率、位置对火焰传播特性的影响，以及泄爆条件下火焰传播规律的研究，障碍物和泄爆耦合影响受限空间可燃气体燃爆机理问题尚未解决，分析不够全面，尚需进一步研究。而石化生产中障碍物与泄瀑口共存，且障碍物与泄瀑口分布和尺寸都是变化的，阻塞率和开口率对受限空间火焰传播影响较大，因此本文开展障碍物和开口率两因素耦合影响下的受限空间可燃气体爆燃机理研究，以期揭示障碍物和开口率耦合作用下的可燃气体燃爆机理，为相关事故预防提供理论支撑。

2 数值模型 2.1 控制方程

大涡模拟(large eddy simulation，LES)通过建立空间滤波函数，从湍流瞬时运动方程中过滤小尺度涡结构，分解出只含有大尺度涡的运动方程，而小涡对流场的影响是通过在控制方程中引入附加应力项建立亚格子尺度模型（subgrid-scale model）。大涡模拟引入滤波函数对纳维-斯托克斯方程(Navier-Stokes equations)方程进行滤波处理，使得高波数的波被截断，但能量传递过程仍被保留，即允许从大涡传递给小涡。大涡模拟控制方程^[9]为

$ \frac{{\partial \overline \rho }}{{\partial t}} + \frac{{\partial \left( {\overline \rho {{\widetilde u}_i}} \right)}}{{\partial {x_i}}} = 0 $

(1)

$ \frac{{\partial \left( {\overline \rho {{\widetilde u}_i}} \right)}}{{\partial t}} + \frac{{\partial \left( {\overline \rho {{\widetilde u}_i}{{\widetilde u}_j} + \overline p {\delta _{ij}} - {{\widetilde \tau }_{ij}} + \tau _{ij}^{\rm{sgs}}} \right)}}{{\partial {x_j}}} = 0 $

(2)

$ \frac{{\partial \rho {{\widetilde h}_s}}}{{\partial t}} + \frac{{\partial \rho {{\widetilde u}_i}{{\widetilde h}_s}}}{{\partial {x_i}}} - \frac{{\partial \overline p }}{{\partial t}} - \overline {{u_j}} \frac{{\partial \overline p }}{{\partial {x_i}}} - \frac{\partial }{{\partial {x_i}}}\left( {\lambda \frac{{\partial \tilde T}}{{\partial {x_i}}}} \right) = - \frac{\partial }{{\partial {x_j}}}\left[ {\rho \left( {{{\widetilde u}_i}{{\widetilde h}_s} - {{\widetilde u}_i} - {{\widetilde h}_s}} \right)} \right] $

(3)

式中：LES滤波的参量标注为横线上标，质量权重滤波的参量标注为波浪线上标；ρ为密度；p为压力；u_i、u_j为速度分量；t为时间；应力张量σ_ij由分子黏度μ决定，可表示为${\sigma _{ij}} = \mu \left( {\frac{{\partial {{\widetilde u}_i}}}{{\partial {x_j}}} + \frac{{\partial {{\widetilde u}_j}}}{{\partial {x_i}}}} \right) - \frac{2}{3}\mu \frac{{\partial {{\widetilde u}_k}}}{{\partial {x_k}}}{\delta _{ij}}$；ι_ij为亚网格尺度应力，将其定义为${\tau _{ij}} = \rho \widetilde {{u_i}{u_j}} - \rho \widetilde {{u_i}}\widetilde {{u_j}}$；h_s为显焓；λ为热导率；亚网格热焓通量可通过梯度假设近似为：$\rho \left( {\widetilde {{u_i}{h_s}} - {{\widetilde u}_i}\widetilde {{h_s}}} \right) = - \frac{{{\mu _{SGS}}{c_p}}}{{P{r_{SGS}}}}\frac{{\partial \widetilde T}}{{\partial {x_j}}}$，其中μ_SGS为亚网格黏度；Pr_SGS为亚网格的普朗特数；c_p为比定压热容；T为温度。

2.2 燃烧模型

对于预混燃烧，本文采用基于C方程的Zimont燃烧火焰面亚格子模型，c为反应过程变量，模拟过程中取c=0.1为火焰锋面。

$ \frac{\partial }{{\partial t}}\left( {\overline \rho \tilde c} \right) + \frac{\partial }{{\partial {x_j}}}(\overline \rho {u_j}\tilde c) = \frac{\partial }{{\partial {t_j}}}\left( { - \frac{{\mu + {\mu _t}}}{{S{c_t}}}\frac{{\partial \tilde c}}{{\partial {x_j}}}} \right) + {S_c} $

(4)

式中：u_j(j=1，2，3)为x、y和z方向的速度；S_c为过程变量源项；μ_t为亚格子黏性，采用动态Smagorinsky Lilly模型计算；Sc_t为湍流施密特数，取Sc_t=0.75。

$ c = \sum\limits_{i = 1}^n {{w_i}/\sum\limits_{i = 1}^n {{w_{i, \rm{eq}}}} } $

(5)

式中：n为燃烧产物数，w_i为产物组分i的质量分数，w_{i, eq}为平衡产物组分i的质量分数。根据式(5)的定义，未燃混合气体中c=0；燃烧产物中c=1；火焰面处0 < c < 1。

$ {\bar S_c} = A{\left( {{u^{'}}} \right)^{3/4}}{\rho _u}{U_1}^{1/2}{\alpha ^{ - 1/4}}\left| {\nabla \mathop c\limits^\sim } \right|{l_t}^{1/4} $

(6)

式中：模型常数A=0.52，u'为亚格子湍流脉动速度，ρ_u为未燃混合气体的密度，U₁为层流燃烧速度，α为反应物分子热传输系数，l_t为湍流特征尺度。

3 数值模拟 3.1 初始条件及边界条件

参考陈先锋^[13]的实验系统，建立截面为80 mm的正方形、长为500 mm的三维方形管道模型，空载管道内预混气体丙烷与空气以当量比1.0充分混合至静止状态，采用Spark点火，点火能为30 J，点火位置距管道左端面中心50 mm，时间步长0.000 1 s。模拟工况初始压力p₀=0.1 MPa(绝对压力)，初始温度T₀=293 K，反应进度变量c₀=0。压力速度耦合采用PISO算法，边界条件为刚性无滑移绝热壁面，对流项采用二阶逆风格式，扩散项采用中心差分格式。通过网格无关性验证，网格尺寸为2 mm ×2 mm ×2 mm，既能减少计算量又能满足所需模拟精度的要求。

3.2 模拟工况

管道内丙烷-空气预混气体中丙烷气体质量分数为4.5%，点火能大小为1 J，点火位置为管道左端面中心处，在距左端200 mm处分别设置阻塞率为0、0.5、0.7的障碍物(厚2 mm，1倍网格尺寸)构建3种置障管道模型，开口位置在管道右端面中心处，开口率分别为6.25%、14.1%、25%、39.1%、56.25%、64%(如图 1所示)的方形孔开口。通过改变各方形管道开口率模拟燃爆过程。

3.3 结果验证

实验采用高速纹影技术记录火焰阵面结构随时间的变化，大涡模拟得到火焰不同时刻发展过程，如图 2所示左侧和右侧分别为高速纹影技术记录的不同时刻火焰结构与大涡模拟得到的火焰结构随时间的变化，实验结果与数值模拟结果的火焰结构的变化规律一致，均在45 ms出现Tulip火焰；如图 3所示，火焰阵面传播速度随时间变化的实验结果与模拟结果对比，从图 3中可看出，火焰传播速度随时间变化的模拟结果与实验结果变化趋势吻合，实验峰值速度与模拟峰值速度的误差约7.0%，因此认为Zimont燃烧模型有效，后文采用该模型开展开口率对置障管道火焰传播特性影响的研究。

4 结果与讨论 4.1 火焰传播进程

如图 4(a)~(d)所示阻塞率为0时管道开口率分别为6.25%、25%、56.25%、64%时火焰发展进程，在火焰传播前期各开口率管道内火焰进程相差不大，火焰锋面极为相似，而当火焰传播经过管道中间位置时（即这4种开口率管道内火焰进程分别发展到25、24、23、22 ms时），火焰结构开始受到末端开口率影响，出现不同的“凸出”结构(图 4中方框标出)，并且这4种开口率管道火焰到达末端的时间为36、29、28、27 ms，可见随开口率增大，火焰传播加快，火焰结构凸起，出现更大形变(图 4中椭圆圈出)。分析认为，管道开口，未燃气体主要受到分布在开口轴线上的前驱冲击波压缩，不受壁面阻挡，拉伸轴线处火焰，并且随着开口率的增大，火焰拉伸作用越明显。当火焰传播到管道末端，火焰结构“凸出”增多，形变增大，分析认为，管道末端开口位置一部分未燃气体进入大气，另一部分未燃气体与到达末端壁面后产生的反射波相互作用，这两部分气体波动叠加对管内压力波产生扰动，从而影响火焰结构，开口率越大这种扰动越强，引起火焰结构变化越明显。

如图 5所示阻塞率为0.5时不同开口率管道内的火焰传播进程。从图 5中分析发现，点火初期(火焰传播到10 ms时)受障碍物及开口率影响不明显，4种开口管道内火焰前锋分别在20、19、19、19 ms时不同程度“跨越”障碍物(图 5中椭圆标出)，随后部分火焰受到湍流作用向障碍物方向逆向发展。以图 5(a)为例分析，20 ms火焰发生弯曲变形，火焰面积激增，燃烧更加剧烈，跨越障碍物后，火焰出现2种方向的改变：障碍物后方的回流区以及靠近末端开口的火焰湍流区。在阻塞率为0.7的管道内，由于阻塞率的增大，拉伸作用明显，与阻塞率为0.5的管道横向对比可知，开口率增大对火焰传播的促进作用增强。对比各开口率下火焰发展随时间变化发现，障碍物对火焰褶皱卷曲作用与末端开口对火焰的拉伸作用存在一种“竞争”关系，二者共同影响火焰进程。

4.2 流场分布规律

开口率接近的同一阻塞率管道流场分布大体一致，着重分析开口率为6.25%及56.25%管道内阻塞率为0、0.5、0.7的流场分布。如图 6、7所示开口率分别为6.25%、56.25%时置障管道内流场分布，从图 6和7中可以看出，在开口率相同时，随阻塞率增加，流场涡团规模增大，流线更加密集，对火焰的剪切作用加强，火焰褶皱破碎程度更显著，火焰面积增加。

对比图 6(a)和7(a)发现，图中流线向开口方向弯曲，越靠近开口端弯曲角度越大，流线越密集。管道末端开口位置一部分未燃气体进入大气，另一部分未燃气体与到达末端壁面后产生的反射波相互作用，这两部分气体波动叠加对管内压力波产生扰动，混合气体爆炸产生的压力波和火焰波传播至管道末端时，受末端开口影响，压力波发生反射并使气流扰动，加强了气流的湍流程度，而且推动已燃气体回传；同时火焰阵面发生扭曲和褶皱，增加了燃烧面积，引起火焰加速，加速的火焰又造成气流湍流程度的增强，从而影响火焰结构，开口率越大时这种扰动越强，引起火焰结构变化越明显。此时火焰锋面结构变化与密闭管道内出现Tulip火焰结构的情况不同^[14-15]。涡团的形成使火焰锋面两侧形成速度方向相反的剪切层，而在阻塞率为0的管道内，由流场图可知管道内没有形成规模涡团，随开口率增大，压力波和速度梯度减弱，其对涡团形成的促进作用逐渐减弱甚至消失。

对比图 6(b)和7(b)发现，当阻塞率为0.5时，由于障碍物的存在，10 ms时障碍物后方出现涡团，20 ms火焰跨越障碍物后，25 ms时火焰锋面和障碍物之间出现涡旋低压区，从图 6(b)中可以发现，当t=25 ms时，流线增多，湍流增强，火焰面结构剪切作用加强、褶皱破碎程度更显著；开口率的增大，加快了火焰传播，图 7(b)中24 ms时，火焰前锋较早波及管道末端。

对比图 6(c)和7(c)发现，当阻塞率为0.7时，较小开口率工况下，不仅在障碍物和火焰前锋之间出现涡团（20 ms），火焰前锋与管道末端之间也形成涡团(24 ms)，随着开口率增大，火焰前锋与管道末端之间涡团消失，末端开口的增大，促使障碍物与火焰前锋之间涡团规模增加，图中流线更加密集，湍流强度增加，开口率影响明显。流场内涡团波及范围较广，火焰面越过障碍物后横向水平拉伸更加明显，管道内火焰传播速度较大。

4.3 管道内火焰传播速度分布

如图 8~10所示阻塞率分别为0、0.5、0.7的置障管道内不同开口率条件下火焰传播速度随监测点位置变化的分布规律，发现在不同阻塞率管道内，随开口率的增加，火焰传播速度峰值整体呈现出增大趋势，当火焰锋面到达距点火位置100 mm处，开口率对火焰传播速度的促进作用开始显现。阻塞率为0、0.5、0.7时，开口率为64%管道内火焰速度峰值与6.25%管道相比，分别提升了192.61%，68.17%，98.19%，开口率在阻塞率为0时管道内加速效果最为明显。开口率为0、64%时，阻塞率为0.7的管道内火焰传播速度峰值相比于阻塞率为0管道增大了393.99%、233.92%。阻塞率与开口率都会增大火焰传播速度峰值，阻塞率对管道内火焰传播速度的加速效果更明显。

由图 8可知，阻塞率为0时管道开口率由6.25%持续增大时，火焰传播速度总体呈单调上升趋势，与密闭管道内火焰速度先升后降不同。密闭状态下，火焰锋面越逼近末端壁面，产生的与火焰传播速度方向相反的反射波越强，对火焰传播的抑制作用越强，火焰传播速度降低。

由图 9和10可知，开口管道末端处火焰传播速度也比密闭管道大得多，因为管道末端开口时，管内气体与大气之间对流作用，会促进管内气流向管道末端流动，在无障碍物管道内，随着开口率的增加，火焰传播速度会持续增加直到管道末端。

4.4 爆燃压力变化规律

开口率影响管内压力波动，对较大阻塞率管道内压力扰动作用更显著，选择开口率为6.25%、56.25%管道内压力变化进行分析，如图 11所示。当开口率增大为56.25%时，管道内压力波动较大，出现多个峰值，各阻塞率管道内峰值压力分别下降了96.83%、97.93%、98.58%，所以当阻塞率为0.7时，泄压百分比最大，开口率增大对阻塞率为0.7时管道内峰值压力削弱作用最强，变化趋势影响最大。

不同工况得到的峰值压力汇总如图 12、13所示。由图 12可知，管道中峰值压力与开口率呈负相关变化，开口率增大造成管道发生明显泄压。当开口率为6.25%时各管道内最大爆炸压力分别为0.069 4、0.236 7、0.210 7 MPa，当开口率为64%时分别为0.001 6、0.005 2、0.004 2 MPa，降幅为97.69%、97.80%、98.01%。根据爆炸条件下空气超压对人体伤害程度分类^[16]可知，当空气超压冲击波压力小于0.02 MPa时，对人体无伤害，3种不同阻塞率管道达到安全泄压的开口率分别为25%、39.1%、25%。

如图 13所示为不同开口率下峰值压力与阻塞率关系，从图中可以发现在开口率为6.25%时阻塞率对管道内峰值压力的影响更显著，随着开口率的增加，阻塞率对峰值压力的影响减小。在同一开口率状态下，阻塞率由0增加到0.7时，管内峰值压力先增大后减小，阻塞率为0时峰值压力最小，阻塞率为0.5时峰值压力达到最大，阻塞率为0.7时次之。阻塞率增大，湍流程度增加，管内燃烧反应加剧，压力会出现较大增加。而阻塞率为0.7的管道内，由于障碍物节流作用^[3]，通过障碍物时火焰厚度较阻塞率为0.5时薄，以致火焰传播速度增大而燃烧反应不如阻塞率为0.5的剧烈，而且流场扰动波及范围较广，造成未燃气体逸出管道，能量损失。这是阻塞率增大会引起管内压力发生较快加速，而阻塞率过大反而会削弱管内压力增大的原因。阻塞率为0.7的置障管道内，火焰面越过障碍物后与流场的相互作用加剧了管内气流扰动，末端开口时障碍物促进了泄压。

5 结论

采用大涡模拟对各置障管道末端开口率对火焰传播特性、管道内流场变化、压力、速度进行模拟计算，得到以下结论：

(1) 障碍物和末端开口对火焰传播进程有显著促进作用。开口率增加，管内气体与大气之间对流作用促进预混气体向管道末端流动，流线更加密集，涡团规模增加，火焰锋面褶皱增多，破碎程度增加，管道内燃烧进程加快，火焰传播速度增加；当障碍物存在时，受到壁面及障碍物的反射作用，流场内涡团波及范围较广，火焰面越过障碍物后水平拉伸更加明显，火焰受到扰动增加，湍流程度增强，加快火焰传播速度。

(2) 管道中峰值压力与开口率呈负相关变化，开口率大，管道内压力波动增大，管道泄压越显著；在同一开口率条件下，阻塞率影响峰值压力，阻塞率为0.7时泄压百分比最大，最大可达98.58%，泄压效果最明显。对于阻塞率为0的管道与阻塞率为0.7的管道，开口率为25%可达到安全泄压效果，而对于阻塞率为0.5的管道，在开口率为39.1%时才能实现安全泄压。

(3) 管道内设置障碍物和增大开口率均会加速火焰传播、改变火焰形态，但阻塞率对火焰传播速度发展趋势的影响占主导作用。各阻塞率下，开口率由0增加到64%时，管道内火焰传播速度峰值提升了68.17%~192.61%；而各开口率下，当阻塞率由0增加至0.7时，管道内火焰传播速度峰值提升了233.92%~393.99%。