温小萍,余明高,邓浩鑫,陈俊杰,王发辉,刘志超
(1河南理工大学机械与动力工程学院,河南 焦作 454003;2河南理工大学安全科学与工程学院,河南 焦作 454003)
小尺度受限空间内瓦斯湍流爆燃大涡模拟
温小萍1,余明高2,邓浩鑫1,陈俊杰1,王发辉1,刘志超1
(1河南理工大学机械与动力工程学院,河南 焦作 454003;2河南理工大学安全科学与工程学院,河南 焦作 454003)
构建了150 mm × 150 mm × 500 mm小尺度受限空间三维模型,基于火焰表面密度模型和Charlette湍流燃烧模型,对两侧连续障碍物条件下瓦斯爆燃火焰与湍流耦合过程进行了大涡模拟(LES)。模拟结果均与实验结果进行了比较。结果表明:大涡模拟可以很好预测瓦斯爆燃过程中的火焰结构、火焰锋面位置、火焰传播速度及超压,验证了大涡模拟及湍流燃烧模型对于瓦斯爆燃的适用性。此外,通过Karlovitz数定量描述了瓦斯爆燃火焰与湍流之间的相互作用及其变化规律,并对不同时刻的火焰模态进行了判别,在两侧连续障碍物条件下瓦斯湍流爆燃火焰先后经历波纹小火焰和薄反应区两种模态。
爆燃;湍流;模型;数值模拟;超压;火焰
DOI:10.11949/j.issn.0438-1157.20151219
Lee等[5]通过实验研究指出,大尺度障碍物形成的湍流是爆燃火焰加速机制的最重要因素,并将这种在大尺度湍流激励下的爆燃过程称为湍流爆燃。国内外不少学者研究过瓦斯爆燃火焰加速机理,并进行了相关实验[6-10]或数值模拟[11-18]的研究,但由于现有实验方法和测试手段的限制,均未能很好捕捉详实的火焰传播过程和湍流流场。特别在数值模拟方面,对于瓦斯爆燃过程的研究大都采用传统的雷诺平均模拟(RANS)方法,其主要缺点在于它只能提供湍流的平均信息,这对于研究瓦斯爆燃的复杂湍流过程远远不够,因此无法有效预测瓦斯爆燃火焰-湍流耦合的行为细节。
大涡模拟(LES)将耗散尺度的脉动进行过滤,只对大尺度脉动进行求解,与雷诺平均模拟相比计算量更大,但其空间分辨率更高,因此可以有效捕捉瓦斯爆燃过程中的湍流特征。目前国内外关于瓦斯爆燃的有效大涡模拟仍然较少,未能揭示瓦斯爆燃火焰与湍流的耦合特性,更缺少这方面的实验佐证。本文对两侧连续障碍物条件下瓦斯爆燃火焰与湍流耦合过程进行大涡模拟,将数值模拟得到的火焰形状、火焰锋面位置、火焰传播速度及超压等数据与笔者实验结果[19]进行比较,进而揭示瓦斯爆燃过程中火焰-湍流耦合机理。
1.1大涡模拟控制方程
实验[19]中涉及小尺度受限空间内的瓦斯爆燃气体流速不超过100 m·s-1,即Ma<0.3,故能量方程中的压力梯度项可以忽略。但由于反应区的温度变化较大,爆燃气流应视为可压缩理想气体。一般情况下,可压缩湍流的大涡模拟控制方程可以通过N-S及其他方程的过滤导出,然而直接过滤可压缩N-S方程将得到十分复杂的可解尺度湍流的运动方程,因此本文采用Favre过滤方法,得到封闭的可压缩湍流大涡模拟方程
式中,ρ为密度;p为压力;T为温度;t为时间;i、j为坐标方向;ui是i方向速度;xi为直角坐标参量;σij为黏性应力张量;h为焓;λ为热导率;c为反应进程变量;Tω.为化学反应热;cω.为归一化的化学反应速率;D为扩散系数;上标“—”及“~”分别表示物理空间过滤量和Favre过滤量;亚网格应力、亚网格物质流、亚网格热流和滤波后的反应速率采用亚网格湍流模型和亚网格燃烧模型进行封闭,其中亚网格Prandtl数PrSGS和亚网格Schmidt数ScSGS均取0.7[11]。
1.2亚网格燃烧模型
过滤后的反应进程变量分子扩散及燃烧速率可写为
式中,ρu为未燃混合气体密度;Sl为层流火焰速度。将式(8)、式(11)代入式(4)中可得到
为了考虑火焰锋面与亚网格湍流之间的相互作用,Charlette等[20]给出了湍流燃烧模型,即式(12)中亚网格火焰褶皱系数ΞΔ表达式为
式中,u'Δ为亚网格湍流脉动速度;Δ为过滤尺度;亚网格湍流Reynolds数ReΔ=u'ΔΔ/ν;层流火焰厚度δf通过δfSl/ν=4估算得到,ν为运动黏度。
图1 实验装置及物理模型(单位:mm)Fig.1 Experimental setup and physical model (unit: mm)
为了提高管道出口边界条件的准确性,减小出口回流对管道内部压力的影响,物理模型在管道出口往x、y及z方向均延伸至1000 mm,即增加了一个尺寸为1000 mm × 1000 mm × 500 mm的扩展区域。这个扩展区域与原管道之间相通,在模拟过程中允许压力波在扩展区域继续传播,因此可以更加真实地反映管道出口条件(扩展区出口设为压力出口)。计算网格采用六面体非均匀网格,即在障碍物和点火源的附近选择较细的网格,而在远离障碍物和点火源的区域(特别是扩展区域)采用较粗的网格,其三维计算网格如图2所示。
图2 三维计算网格(单位:m)Fig.2 3D computational grids (unit: m)
在大涡模拟中,亚网格尺度的物理扩散是随着网格的加密而减小的,因此,对于大涡模拟而言,并不存在通常数值计算中所谓“网格无关性”的概念,较好的做法是采用十分精细的网格和很小的时间步长,以确保其数值扩散远小于雷诺平均方法。本文在障碍物和点火源附近区域的网格大小为δx = δy = δz = 1.25 mm,远离障碍物和点火源的区域的网格逐渐增大。计算单元总数约为180万,其中管道区的单元数约为150万,约占计算单元总数的83.3%。数值求解是利用CFD求解器ANSYS FLUENT实现的,整个计算区域采用有限容积法离散控制方程,空间离散采用二阶中心差分格式,选取二阶隐式时间推进法来提高计算精度,利用SIMPLE算法对速度和压力耦合方程组进行解耦。质量方程、动量方程、能量方程、反应进程变量方程残差分别小于1×10-6、2.5×10-5、1×10-6及1×10-5。
3.1火焰结构
图3为采用大涡模拟和Charlette湍流燃烧模型时不同时刻火焰结构与实验结果的比较,其中模拟结果中的火焰结构采用反应进程变量c=0.5等值面表示[14]。可以清楚地看到,Charlette湍流燃烧模型不仅很好地预测了火焰结构演变过程(包括火焰结构和火焰锋面达到各个位置的时间),而且准确捕捉到了瓦斯湍流爆燃火焰的结构特征。在爆燃初期(如t=20 ms),火焰结构主要为半球形,火焰传播速度较慢,致使火焰锋面到达第1个障碍物的时间约为25 ms,约占用火焰在整个管道内传播时间的59%。至t=32 ms时,火焰锋面从第1对障碍物之间的间隙喷出,并逐渐向两侧卷曲,说明火焰明显受到两侧障碍物后的湍流影响。火焰绕过第1对障碍物之后,继续向第2对和第3对障碍物传播。至t=42 ms时,火焰锋面不断发生褶皱,出现了明显的火焰湍流现象,这种现象是由于受到多个障碍物持续作用后,火焰被不断褶皱、变形,导致未燃气体与已燃气体能够快速相互掺混,由此出现明显的火焰湍流特征。在整个过程中,大涡模拟和实验结果非常相似。
图3 大涡模拟的不同时刻火焰结构与实验结果的比较Fig.3 Comparison between LES flame structures at different instants and experimental data
3.2火焰锋面位置与火焰传播速度
图4、图5分别表示利用Charlette燃烧模型模拟的不同时刻火焰锋面位置及火焰锋面传播速度与实验结果的对比图。其中,火焰锋面位置是指每张火焰图片中火焰锋面与点火点的最大距离,火焰速度则是通过前后连续两个时间点的火焰锋面位置差与时间间隔之比计算而得。从图4可以看出,模拟与实验在爆燃初期的火焰锋面位置存有偏差,实验数值略高于模拟结果,这是由于实验中高速摄像的角度偏差造成的。但是,在瓦斯湍流爆燃的整个过程中,模拟预测的火焰锋面位置与实验吻合较好。图5则反映了模拟的火焰传播速度与实验也非常接近,尤其是很好再现了火焰绕过障碍物时的加速、减速特征。稍有区别的是数值模拟中火焰加速程度略低于实验数据,这是由于数值模拟的计算网格精度有限,使得大涡模拟不能捕捉到极小尺度的湍流涡团对火焰加速的激励作用。显然,这种影响对模拟的精确度影响较小。
图4 大涡模拟的不同时刻火焰锋面位置与实验结果的比较Fig.4 Comparison between LES time histories of flame front position and experimental data
图5 大涡模拟的不同时刻火焰锋面位置与实验结果的比较Fig.5 Comparison between LES simulated time histories offlame front position and experimental data
3.3爆燃超压
图6为大涡模拟得出的瓦斯爆燃超压度与实验结果的对比。从实验结果来看,爆燃超压在15 ms ≤t ≤ 25 ms时间段出现了一个小幅度峰值,该峰值是由于管道顶端薄膜破裂形成的,而数值模拟并未考虑薄膜的影响,因此模拟超压并未出现这一小幅度峰值。还可以看到,数值模拟预测的最大超压为10.1 kPa,与实验最大超压11.6 kPa的误差为12.9%。这里值得注意的是,数值模拟和实验的最大超压的时间非常吻合,均出现在t=41 ms,再对照图3中的火焰结构可以看出,此时火焰锋面已经越过最后一对障碍物并迅速往两侧发展,火焰在障碍物后卷曲变形呈现蘑菇状,附近的未燃气体则以更快的速度掺混其中并燃烧,导致超压的急速上升。在41 ms ≤t ≤ 42 ms时间段内,虽然第3对障碍物下游的火焰继续膨胀,但管道底部的火焰在抵达壁面后开始出现局部燃尽熄灭现象,这可能使管道内火焰面密度减小,并导致超压开始下降。此后,由于火焰锋面已经越过容器出口,管道超压迅速回落。
图6 模拟的不同时刻爆燃超压与实验结果的比较Fig.6 Comparison between LES time histories of deflagration overpressure and experimental data
图7为火焰锋面分别经过3对障碍物不同时刻的速度矢量分布。当火焰经过第1对障碍物时(t=29 ms),3对障碍物之后均形成了一对左右对称的大尺度旋涡,左侧旋涡逆时针旋转,而右侧旋涡顺时针旋转。这是由于燃烧区的高温气流不断膨胀,迫使下游未燃气体向管道出口方向流动,而这些未燃气体经过障碍物时,由于主流流速与障碍物后气体流速存在较大的速度差,从而在障碍物后形成旋涡。这些旋涡的尺度大小与障碍物尺度、障碍物位置、流体流速、流体黏度有关。当火焰继续经过第2对和第3对障碍物时,随着气流流速的增大,旋涡尺度和速度大小都在增加。
为了揭示爆燃过程中火焰与湍流的耦合关系,本文引入Karlovitz数(Ka)来表征湍流与燃烧相互作用的程度,从而定性分析瓦斯爆燃火焰在连续障碍物条件下的火焰模态。Ka定义为火焰面的时间尺度τf与流动中最小涡(Kolomogrov涡)的时间尺度τη之比
图7 不同时刻速度矢量分布Fig. 7 Velocity vector maps when flame passes sequentially over three obstacles
对于预混湍流燃烧的大涡模拟,Pitsch等[22]提出Ka按式(15)近似计算
图8表示3个不同时刻(29、37及41 ms)火焰面上Ka的空间分布,火焰面取反应进程变量c = 0.5的等值面。在这3个时刻,火焰锋面正在分别绕过3个障碍物。由图8可以看出,由于位于点火源附近的湍流强度很小,其火焰面比较光滑,Ka始终处于较低水平。但随着火焰锋面先后绕过3个障碍物时,位于障碍物附近的火焰锋面上Ka出现明显增长,其中35 ms时刻位于第3个障碍物附近火焰面上的Ka较高,最大值可达到10~11。其原因是在障碍物湍流激励作用下,障碍物附近流体的应变率和亚网格黏度较高,使其亚网格湍流脉动处于较高水平造成的,根据式(15)可知,较高的湍流脉动所对应的Ka也就越大。因此,可进一步利用整个流场的Ka最大值(在障碍物附近)来反映火焰与湍流之间的作用程度。
图8 火焰面(c=0.5)上Ka在不同时刻的空间分布Fig.8 Spatial distribution of Ka over flame surface (c=0.5) at different time instants
Pitsch等[22]将大涡模拟过滤尺度和Ka引入之后,对燃烧模态区间作了进一步改进,改进后燃烧模态区间的划分如图9所示。图中的方形点表示瓦斯爆燃火焰在中间连续障碍物条件下不同时刻的Ka-Δ/δl数据点,对应时刻分别为2、14、34、37、41及42 ms,其中Ka均取对应时刻整个流场的最大值。由图9可以看出,在瓦斯爆燃初期Ka略低于1,还处于波纹火焰面模式;随着火焰向下游传播,Ka不断升高,火焰开始向薄反应区转变;当t = 42 ms时(此时火焰锋面抵达管道出口),Ka达到9.8,但仍然属于薄反应区。由此可知,在连续障碍物条件下,瓦斯爆燃火焰先后经历了波纹小火焰和薄反应区两种火焰模态,但对于火焰-湍流耦合作用的时间段(34 ms < t < 42 ms)而言,则主要处于薄反应区火焰模态。
图9 大涡模拟结果在湍流火焰模态区间上的分布Fig.9 Distribution of LES data on turbulent flame regimes
(1)利用大涡模拟及Charlette湍流燃烧模型对两侧障碍物条件下的瓦斯湍流爆燃过程进行数值模拟,并将模拟获得的火焰结构、火焰锋面位置、火焰速度及超压等参数特征与实验进行比较,验证了大涡模拟及Charlette湍流燃烧模型对于瓦斯爆燃的适用性。
(2)通过Karlovitz数定量描述了瓦斯爆燃火焰与湍流之间的相互作用及其变化规律,对中间连续障碍物条件下不同时刻的火焰模态进行了判别。在两侧连续障碍物条件下,瓦斯湍流爆燃火焰先后经历了波纹小火焰和薄反应区两种模态。
References
[1] 周心权. 煤矿井下瓦斯爆炸的基本特性[J]. 中国煤炭, 2002, 28(9):8-11. DOI: 10.3969/j.issn.1006-530X.200209002. ZHOU X Q. Basic characteristics of gas explosion in coal mines [J]. China Coal, 2002, 28(9): 8-11.DOI: 10.3969/j.issn.1006-530X. 200209002.
[2] ZHU C J, LU Z G, LIN B Q, et al. Effect of variation in gas distribution on explosion propagation characteristics in coal mines [J]. Mining Science and Technology, 2010, 4(20): 516-519.
[3] JOHANSEN C T, GICCARELLI G. Visualization of the unburned gas flow field ahead of an accelerating flame in an obstructed square channel [J]. Combustion and Flame, 2009, 156(2): 405-416.
[4] KIM W K, MOGI T, DOBASHI R. Flame acceleration in unconfined hydrogen/air deflagrations using infrared photography [J]. Journal of Loss Prevention in the Process Industries, 2013, 26(6): 1501-1505.
[5] LEE J H S, KNYSTAUTAS R, FREIMAN A. High speed turbulent deflagrations and transition to detonation in H2/air mixtures [J]. Combustion and Flame, 1984, 56(2): 227-239.
[6] BYCHKOV V, VALIEV D, AKKERMAN V, et al. Gas compression moderates flame acceleration in deflagration-to-detonation transition[J]. Combustion Science and Technology, 2012, 184(7/8): 1066-1079.
[7] IBRAHIM S, HAEGRAVE G, WILLIAMS T. Experimental investigation of flame/solid interactions in turbulent premixed combustion [J]. Experimental Thermal and Fluid Science, 2001, 24:99-106.
[8] CICCARELLI G, JOHANSEN C, PARRAVANI M. The role of shock-flame interactions on flame acceleration in an obstacle laden channel [J]. Combustion and Flame, 2010, 157(11): 2125-2136.
[9] WEN X P, YU M G, JI W T, et al. Methane-air explosion characteristics with different obstacle configurations [J]. International Journal of Mining Science and Technology, 2015, 25(2): 213-218.
[10] 郑立刚, 吕先舒, 郑凯, 等. 点火源位置对甲烷-空气爆燃超压特征的影响[J]. 化工学报, 2015, 66(7): 2749-2756. DOI:10.11949/j.issn.0438-1157.20141789. ZHENG L G, LÜ X S, ZHENG K, et al. Influence of ignition position on overpressure of premixed methane-air deflagration [J]. CIESC Journal, 2015, 66(7): 2749-2756. DOI: 10.11949/j.issn.0438-1157. 20141789.
[11] SARLI V, Benedetto A, RUSSO G, et al. Large eddy simulation and PIV measurements of unsteady premixed flames accelerated by obstacles [J]. Flow Turbulence and Combustion. 2009, 83: 227-250.
[12] SARLI V, Benedetto A. Effects of non-equidiffusion on unsteady propagation of hydrogen-enriched methane/air premixed flames [J]. International Journal of Hydrogen Energy, 2013, 38(18): 7510-7518.
[13] BI M S, DONG C, ZHOU Y. Numerical simulation of premixed methane-air deflagration in large L/D closed pipes [J]. Applied Thermal Engineering, 2012, 40: 337-342.
[14] SARLI V, Benedetto A, RUSSO G. Large eddy simulation of transient premixed flame-vortex interactions in gas explosions [J]. Chemical Engineering Science, 2012, 71: 539-551.
[15] JOHANSEN C, CICCARELLI G. Modeling the initial flameacceleration in an obstructed channel using large eddy simulation [J]. Journal of Loss Prevention in the Process Industries, 2013, 26(4):571-585.
[16] 朱传杰, 林柏泉, 江丙友. 受限空间内爆燃波瞬态流速与超压的耦合关系[J]. 燃烧科学与技术, 2012, 18(4): 326-330. ZHU C J, LIN B Q, JIANG B Y. Coupled relationship between gas velocity and peak overpressure of deflagration wave in confined space [J]. Journal of Combustion Science and Technology (China),2012, 18(4): 326-330.
[17] WEN X P, YU M G, LIU Z C, et al. Effects of cross-wise obstacle position on methane-air deflagration characteristics [J]. Journal of Loss Prevention in the Process Industries, 2013, 26(6): 1335-1340.
[18] XIAO H H, MAKAROV D, SUN J H, et al. Experimental and numerical investigation of premixed flame propagation with distorted tulip shape in a closed duct [J]. Combustion and Flame, 2012, 159(4):1523-1538.
[19] 温小萍, 武建军, 解茂昭. 瓦斯爆炸火焰结构与压力波的耦合规律[J]. 化工学报, 2013, 64(10): 3871-3877. DOI: 10.3969/j.issn. 0438-1157.2013.10.052. WEN X P, WU J J, XIE M Z. Coupled laws between flame structure and pressure wave of gas explosion [J]. CIESC Journal, 2013, 64(10):3871-3877. DOI: 10.3969/j.issn.0438-1157.2013.10.052.
[20] CHARLETTE F, MENEVEAU C, VEYNANTE D. A power-law flame wrinkling model for LES of premixed turbulent combustion(Ⅰ): Non-dynamic formulation and initial tests [J]. Combustion and Flame, 2002, 131(1): 159-180.
[21] IBRAHIM S S, GUBBA S R, MALALASEKERA W. A parametric study on large eddy simulations of turbulent premixed flames [C]// The 8th Asia-Pacific Conference on Combustion. 2010: 87-94.
[22] PITSCH H, DUCHAMP L. Large-eddy simulation of premixed turbulent combustion using a level-set approach [J]. Proceedings of the Combustion Institute, 2002, 29(2): 2001-2008.
Large eddy simulation of gas turbulent deflagration in small-scale confined space
WEN Xiaoping1, YU Minggao2, DENG Haoxin1, CHEN Junjie1, WANG Fahui1, LIU Zhichao1
(1School of Mechanical and Power Engineering, Henan Polytechnic University, Jiaozuo 454003, Henan, China;
2School of Safety Science and Engineering, Henan Polytechnic University, Jiaozuo 454003, Henan, China)
A 3D model of small-scale confined space with an inner size of 150 mm × 150 mm × 500 mm was set up. Based on the flame surface density model and the turbulent combustion model by Charlette et al., a large wddy simulation (LES) had been carried out on the process of gas deflagration flame-turbulence interaction with continuous obstacles at two sides of the chamber. All numerical results have been compared to experimental data. It showed that the LES is capable to predict the flame structure, position, speed, and overpressure in the process of gas deflagration, and the applicability of the LES and turbulent combustion model on gas deflagration was verified. In addition, the interaction between gas deflagration and turbulence and the relationship were described quantitatively by the Karlovitz number, and the transient flame regimes also identified. Under condition of continuous obstacles at double sides of the chamber, the gas turbulent deflagration flame experienced in the subsequent states of corrugated flamelets zone and thin reaction zone.
deflagration; turbulent flow; model; numerical simulation; overpressure; flame
发生在煤矿井下的瓦斯爆炸绝大多数属于爆燃[1]。当发生瓦斯爆燃时,火焰在传播过程中往往会遇到各种障碍物,如矿车、采掘机电设备、液压支架、管道、巷道风门等,火焰前方未燃气体不断受到热膨胀波的压缩推动作用,使得在障碍物下游流场形成大尺度涡团,其湍流尺度远大于层流火焰厚度,湍流脉动将诱导紧随而至的火焰产生明显褶皱、卷吸,从而使火焰锋面附近的已燃高温气体和未燃气体快速掺混,燃烧速率和火焰表面积增大,进而使火焰传播获得加速,并伴有爆燃超压的增加,引起更强烈的湍流[2-4]。
date: 2015-07-29.
WEN Xiaoping, wenxiaoping666@163.com
supported by the National Natural Science Foundation of China (51176021).
TD 712
A
0438—1157(2016)05—1837—07
2015-07-29收到初稿,2016-01-25收到修改稿。
联系人及第一作者:温小萍(1977—),男,博士,副教授。
国家自然科学基金项目(51176021);河南省教育厅科学技术研究重点项目(14A410007)。