基于RNG k-ε湍流模型的3D采空区瓦斯上浮贮移

李宗翔,顾润红,张晓明,毕 强,温永宇

(1.辽宁工程技术大学安全科学与工程学院,辽宁阜新 123000;2.辽宁工程技术大学矿山热动力灾害与防治教育部重点实验室,辽宁阜新123000;3.中交一航局第五工程有限公司,河北秦皇岛 066002;4.铁法能源公司,辽宁调兵山 112700)

基于RNG k-ε湍流模型的3D采空区瓦斯上浮贮移

李宗翔1,2,顾润红3,张晓明1,2,毕 强4,温永宇4

(1.辽宁工程技术大学安全科学与工程学院,辽宁阜新 123000;2.辽宁工程技术大学矿山热动力灾害与防治教育部重点实验室,辽宁阜新123000;3.中交一航局第五工程有限公司,河北秦皇岛 066002;4.铁法能源公司,辽宁调兵山 112700)

为研究工作面后方采空区3D空间瓦斯上浮特征,根据质量守恒、动量守恒和Fick定律,建立综放采空区风流-瓦斯变密度混合气体非线性渗流-扩散控制方程。结合大兴矿N2-706工作面实例,运用CFD软件模拟三维采空区瓦斯运移及分布的状态;流场冒落非均质按照“O”型分布描述,计算采用RNG k-ε湍流模型,并考虑重力场条件。理论计算调整与现场实际条件和瓦斯监测记录数据相拟合,指出三维采空区形成的瓦斯分布及上浮态势,是风流移动、瓦斯组分扩散-弥散和含瓦斯风流密度差引起上浮迁移的结果,也是流场瓦斯不断解吸涌出与漏风不断流入所形成的一种动态平衡结果。研究结果表明,采空区高位抽采流量与抽采获得的瓦斯体积分数近似呈反比例函数关系、与回风巷瓦斯体积分数呈负指数函数关系。

3D采空区;瓦斯上浮;CFD仿真;k-ε湍流模型;瓦斯抽采

工作面后方采空区具有用常规方法处理较难、人员又无法进入的特殊性,应用数值模拟手段来研究采空区瓦斯与自燃规律问题已经是不可或缺的主流方法[1-6]。近年来国内外学者关于采空区瓦斯贮移规律的建模和应用Fluent等CFD软件进行模拟求解的文献很多,但在理论计算与实际贴近度,以及处理方法方面还有很多不足,尚需在以下几方面加以完善:①采空区内风流流动是非线性的,全区域内层流、过渡流和紊流同时存在,仅仅把采空区内风流流动视为Darcy流是不够的;②紊流区域应选择合适的湍流模型;③CFD模拟3D采空区应使用如“O”型圈分布的连续性非均质各向同性介质模型;④ 采高较大的综放采空区3D模型计算,如笔者编制的G3程序[7-9]会存在较大的误差;⑤ 应考虑重力场作用下的瓦斯上浮问题[6]。基于上述理由,有必要就应用CFD软件对采空区场数值模拟再作研究,使综放采空区瓦斯涌出运移及分布规律特征的讨论更趋完备。

1 采空区3D模型及风流移动-瓦斯运移基本方程

1.1 采空区“O”型冒落特征及计算模型

综放工作面3D冒落采空区如图1所示。冒落空间内部为非均质冒落体,按“O”型圈模型描述。

图1 采空区三维模型及冒落碎胀系数变化Fig.1 3-D numerical model of the goaf and change of broken expand coefficient

按“O”型圈压实分布的碎胀系数表达式[9]为

其中,KP,0,KP,1分别为初始冒落和压实后的碎胀系数,无因次;a0,a1,ξ1为控制“O”型圈模型分布形态的调整系数,ξ1∈[0,1],根据矿压观测通过试算达到符合实际情况。图1(b)中 KP,0=1.5,KP,1=1.05, a1=0.036 8,a0=0.156,ξ1=0.233。

1.2 采空区变密度混合气体组分输运方程

采空区混合气体的组分输运方程为

其中,Θ为采空区多孔介质内气体组分;ρ为Θ组分密度,kg/m3;c(Θ)为Θ组分气体体积分数;DΘ为Θ组分气体的动力弥散系数张量(各向同性介质取标量值,以下同),m2/s,DΘ=Dt+Dm,其中Dm为气体分子扩散张量(在大雷诺数下,可将Dm忽略),Dt为多孔介质紊流机械扩散系数张量,Dt=eiejDt,ij,ei,ej为沿i,j坐标方向的单位矢量,Dt,ij=αTVδij+(αL-αT)× uiuj/V,其中V为渗流流速,m/s;δij为Kroneckey数; αL,αT分别为速度V纵、横向弥散度,在三维Oxyz直角坐标系下流速V沿x,y,z轴的分量为ui(i,j=x,y, z),在方程源汇项中c′为场析出的瓦斯浓度,kg/m3; ω为采空区单位时间、单位体积内冒落煤体解吸出的瓦斯体积量,m3/(m3·s)。

设采空区混合气体为理想气体,满足气体状态方程。

1.3 RNG k-ε湍流模型方程的引入

关于采空区与工作面风流交换关系,有的文献将工作面通风作为采空区边界条件考虑[10-11]。这里建立的采场物理模型是把工作面及其部分巷道与采空区合成一体来考虑。工作面及其进、回风巷内的风流为大雷诺数湍流,工作面附近采空区漏风流速大,漏风在大块矸石空隙间流动也为湍流流动,需建立湍流计算模型。在众多的湍流模型中,RNG k-ε模型能更好地处理低雷诺数和流线弯曲程度较大的流动。采空区冒落岩石,风流流线弯曲度较大,将这些看作符合RNG k-ε模型的多孔介质域,能更好地模拟采空区内气体在的真实流动[12-13]。

RNG k-ε模型来源于严格的统计技术,由暂态N-S方程推出。在ε方程中引入附加生成项,考虑了湍流旋涡,考虑了低雷诺数流动黏性处理近壁面区域,可有效地改善精度,使得 RNG k-ε模型比标准k-ε模型在更广泛的流动中有更高的可信度和精度。k和ε是两个基本的未知量。对于稳态、不可压缩气体、无源项,RNG k-ε模型可以写成

式中,k为湍流脉动动能,J;ε为k的耗散率,无因次; C1ε,C2ε,C3ε为经验常数,这里取 C1ε=1.42,C2ε= 1.68,C3ε=1.72,αk=αε=1.393;μeff=μ+μt,μ为风流黏度,μt为湍流黏度,μt=ρcμ(k2/ε);Gk为由于平均速度引起的湍动能k的产生项;Gb为由于浮动引起的湍流动能k的产生项;v^=μeff/μ;经验常数Cν=100; Rε代表平均应变率对ε影响的附加项,其中,η=Sk/ ε为无量纲应变或者平均流时间尺度与湍流时间尺度之比,经验常数取 η0=4.38,β=0.015,Cμ= 0.084 5;S为多孔介质风流运移过程中附加的动量损失源项,对各向同性多孔介质

式中,B为黏性阻力损失系数(为渗透率的倒数), 1/m2;C2为局部阻力因子。

此处以Blake-Kozeny推导思想,在非线性流态条件下定义适用的经验公式为

其中,dm为采空区介质几何平均粒径,m;n为采空区孔隙度。采空区孔隙度与碎胀性系数有

根据现场测试的数据以及CFD模拟经验[2],通过反复模拟试验确定采空区多孔介质模型所对应的识别参数C0的值。

2 不同因素下采空区瓦斯分布解的特征

算例是铁法矿区大兴煤矿N2-706工作面采空区,工作面长度为188 m,煤层厚度为5.62～8.53 m,煤层平均厚度7.6 m,工作面推进度为1.2 m/d。现场实测得工作面两端压差为32 Pa,工作面风量为1 308 m3/min。瓦斯涌出强度按距离工作面呈负指数衰减变化,衰减率为0.037 6。风流温度和采空区初始温度为21℃。

当前,模拟计算工具可供选择的仿真平台很多,近年来流行的如 Fluent或 COMSOL流体力学软件[1-3],笔者采用 Fluent软件模拟计算。在利用Gambit建立采空区3D模型后,通过Fluent软件提供的UDF编程来实现更改采空区内部的参数和一些功能的设置,使之尽可能贴近采空区实际情况。这里取C0=1.041×10-5,dm=0.15 m。

图2是结合大兴煤矿N2-706采空区流场条件采取几种不同处理方法时的模拟结果。通过模拟应用发现,Fluent不同处理方法会对采空区瓦斯分布解产生很大影响,图2(a)是没有考虑体积力(重力)作用时的层状瓦斯分布结果,瓦斯未出现上浮;图2(b)是均质多孔介质流场的结果,采空区漏风流衰减大、瓦斯分布变化范围过小的不准确结果;图2(c)是采空区渗透率控制不当的失真结果,渗透率过大,瓦斯分布变化范围过大;图2(d)是因素全面考量但计算收敛较差的一种结果。以上均为因处理不当导致的常见的典型错误结果,可见正确处理和设置软件平台的模型是计算成败的关键,否则会对瓦斯分布结果产生很大的偏差,尤其必须考虑重力场的效应。

3 采空区瓦斯分布及瓦斯上浮机理

瓦斯密度为空气密度的0.554倍,瓦斯的上浮现象是采空区内瓦斯分布的一大特征,这是已被瓦斯抽采实践所证实的结论。图3是综合考虑“非线性-非均质-变渗透率-重力场”等多因素完备条件下的计算结果。图3中模拟抽采时的流量为1.6 kg/s,即80 m3/min,此时抽采小口的负压值为-2 190 Pa,在抽采口局部周围出现了漏斗形的大梯度风压变化区,如图3(b)所示。模拟产生的采空区瓦斯上浮效果非常明显。采空区3D模拟现出的瓦斯上浮效应,是真实采空区瓦斯分布的写照,对瓦斯抽采研究具有重要的指导意义。

Fluent模拟出现的采空区瓦斯上浮分布现象是基于流体力学原理得到的客观理论结果,但在现实中有时人们还存在理解上的混淆。采空区瓦斯的上浮不能简单理解为是甲烷分子(组分)的上浮,而主要是因风流中瓦斯体积分数不同引起风流密度变化导致的上浮运动。含瓦斯风流一旦与漏入的新风交汇时,因两风流的密度不同,必然会出现运动的变化,出现瓦斯上浮运动以及瓦斯上浮的分布。所谓瓦斯上浮分布是指采空区沿着高度方向上高位处要比低处瓦斯浓度高。瓦斯上浮分布是瓦斯解吸涌出、组分扩散-弥散与密度差引起的自然对流迁移的动态平衡结果,也是瓦斯不断涌出和漏风流流入所形成的一种流体运动扩散的平衡状态。一旦瓦斯涌出或风流消失,这种平衡即会打破,在分子的扩散-弥散作用下,瓦斯分布将趋于均匀,形成新的平衡。煤矿现场实际中常表现出采空区解吸释放瓦斯源一般在顶板上层的高位处,使得在采空区高处的瓦斯积存较多的现象,这既与瓦斯源位位置有关,也与上浮运动有关。

图2 软件不同计算设置时采空区瓦斯分布的几种典型特征解(常见的错误解)Fig.2 Characteristics of gas distribution simulation in goaf with different factors

尽管甲烷分子在空气中扩散速度很快,且采空区冒落多孔介质对瓦斯机械弥散起到的加速作用,但组分扩散与弥散速度低于对流速度,这可从巷道中高浓度瓦斯出现分层分布现象得到证实。

图3(a)中在未进行抽采时出现的采空区瓦斯整体上浮,在采空区的回风隅角处,瓦斯随着漏回风流下扎涌入工作面回风隅角,汇入回风巷;在未进行抽采或抽采流量不足时,回风道上靠顶帮处出现瓦斯层积流动现象,模拟结果与综放工作面回风道瓦斯超限时观测到的现象一致,即当抽采不足时,在工作面出口及回风道的高冒顶板处出现了瓦斯积聚层,如图3(a)所示,一度被迫用草垫子隔离出二层棚,用风筒抽排瓦斯;当加大抽采流量后有效控制了工作面上隅角和回风巷顶层的瓦斯超限问题,如图4所示。

图4 N2-706采煤工作面上隅角瓦斯体积分数变化曲线Fig.4 Change curves of gas concentration at upper corner of N2-706 working face

图5是采空区瓦斯抽采时对应的氧气体积分数分布的数值结果,氧气体积分数与瓦斯分布数值相反。工作面附近采空区氧气体积分数变化主要原因还是被瓦斯涌出流稀释,深部氧气体积分数变化则还要与煤的自燃氧化耗氧有关[8]。

图5 采空区抽采时对应氧气体积分数分布Fig.5 O2concentration distribution during drainage in goaf

4 高位抽采采空区空间瓦斯分布及抽采参数量化分析

采空区瓦斯高位抽采是包括高位瓦斯巷道、高位水平钻孔和回风道斜上钻孔(一定高度位置)等几种,不论哪种抽采方式,客观上都在抽采口周围一定范围形成的高位负压区,形成抽采流[14-15]。因此抽采边界条件基本相同,这里抽采边界条件设为2类流量边界条件。

为获得采空区抽采流量与抽采效果的关系,笔者进行了多次有限模拟实验,结果如图6,7所示。实际3D采空区RNG k-ε模型的计算量非常大,使用AMDIIX6-1100T六核心3.3G的CPU电脑服务器做运算,每种抽采流量计算点运行时间需2 d左右,收敛很慢,不同流量也导致收敛性存在差异,结果出现一定的震荡现象。这里对试验结果用回归分析加以调整,以获得3D调节下抽采一般规律性的认识。

图6 抽采流量与抽采获得瓦斯体积分数的关系Fig.6 Relationship of gas drainage flow and gas volume concentration

图7 抽采流量与回风巷瓦斯体积分数的关系Fig.7 Relationship of gas drainage flow and gas volume concentration in return airway

式中,回归系数f1=0.102 238,f2=0.041 105,相关系数为0.999 1,如图6所示。

同样,抽采流量q与回风巷中瓦斯体积分数c1两者也严格地符合负指数函数关系,即

式中,回归系数b1=1.187 06,b2=1.052 37,相关系数为0.990 92,如图7所示。

从图6,7可看出,随抽采流量的增加回风巷瓦斯体积分数逐渐降低,回风巷瓦斯体积分数达到0.01% 以下,临界抽采流量为 1.13 kg/s(即57 m3/min)。如果抽采流量增大到 1.6 kg/s(即80 m3/min)时,回风巷瓦斯体积分数仅为0.67%,此时抽采获得的瓦斯体积分数也随之降低,仅为35%。再结合前面的图4,可以明显看出,当抽采流量为80 m3/min时,抽采的氧气体积分数为15.6%,说明此时抽采瓦斯的效率明显降低,即抽放流量越大,采

随着抽采流量q的增加,抽采瓦斯体积分数c逐渐降低,得到抽采流量q与抽采瓦斯体积分数c近似符合反比例函数关系,即空区氧气体积分数分布范围越大。经验公式(9)和(10)说明,随着抽采流量的增加,工作面向采空区漏风量也随之加大,这样加大了采空区遗煤自燃的可能性。大兴矿是煤层属于易自燃等级,采空区自然发火十分严重,因此,大兴矿N2-706工作面对瓦斯治理的同时,还要考虑防止由于抽采所导致遗煤自燃灾害的发生,这样就要求N2-706工作面的抽采流量要控制在一个合理的范围内,实际工作中应切实避免盲目过度抽采。根据数值模拟结果,大兴矿N2-706工作面的抽采流量应控制在60～65 m3/min内较合理(与实际情况相符),既能控制瓦斯不超限,又确保不能引起采空区煤自燃(配合注氮气)。

5 结 论

(1)不同因素下采空区瓦斯分布模拟解的特征差别很大,对采空区瓦斯问题产生很大偏差,因此多因素完备的条件和重力因素是准确分析高抽巷采空区抽采的关键问题。对大采高尤其是放顶煤开采的高瓦斯采空区,重力影响数值模拟结果的正确性。

(2)高位瓦斯道(或水平钻孔或斜上钻孔)瓦斯抽采是大采高采空区瓦斯抽采的最有效方法。通过模拟实验获得了抽采流量与工作面瓦斯治理结果和抽采流量与抽采瓦斯浓度的关联关系。

(3)避免抽排力度的盲目图大图强(从排放瓦斯的角度上看,抽排瓦斯浓度和抽排效果与抽排流量成衰减关系),同时,增大瓦斯抽排力度的前提也必须考虑防治自燃措施和效果。

[1] 金龙哲,姚 伟,张 君.采空区瓦斯渗流规律的CFD模拟[J].煤炭学报,2010,35(9):1476-1480.

Jin Longzhe,Yao Wei,Zhang Jun.CFD simulation of gas seepage regularity in goaf[J].Journal of China Coal Society,2010,35(9): 1476-1480.

[2] 胡千庭,梁运培,刘见中.采空区瓦斯流动规律的CFD模拟[J].煤炭学报,2007,32(7):719-723.

Hu Qianting,Liang Yunpei,Liu Jianzhong.CFD simulation of goaf gas flow patterns[J].Journal of China Coal Society,2007,32(7): 719-723.

[3] Zhang Xinhai,Xi Guang.Study on partition of spontaneous combustion danger zone and prediction of self-ignition in coalmine based on numeric simulation[J].Journal of Coal Science&Engineering(China),2006,12(1):56-59.

[4] 梁运涛,张腾飞,王树刚,等.采空区孔隙率非均质模型及流场分布模拟[J].煤炭学报,2009,34(9):1203-1207.

Liang Yuntao,Zhang Tengfei,Wang Shugang,et al.Heterogeneousm of porosity in gobs and its airflow field distribution[J].Journal of China Coal Society,2009,34(9):1203-1207.

[5] Sensogut C,Kaufmann M,Petit E.An approach to the modeling of spontaneous combustion in the goaf[J].The Journal of the South African Institute of Mining and Metallurgy,2002,43(4):311-313.

[6] 顾润红.综放采空区3D空间非线性渗流及瓦斯运移规律数值模拟研究[D].阜新:辽宁工程技术大学,2011.

[7] 李宗翔,王晓冬,王 波.采空区场流数值模拟程序(G3)实现与应用[J].湖南科技大学学报(自然科学版),2005,20(3):16-20.

Li Zongxiang,Wang Xiaodong,Wang Bo.Realization and application of numerical simulating program(G3)for field flow of goaf[J].Journal of Hunan University of Science&Technology(Natural Science Edition),2005,20(3):16-20.

[8] 李宗翔,吴 强,肖亚宁.采空区瓦斯涌出与自燃耦合基础研究[J].中国矿业大学学报,2008,37(1):38-42.

Li Zongxiang,Wu Qiang,Xiao Yaning.Numerical simulation of the coupling action mechanism of spontaneous combustion and gas effusion in goaf[J].Journal of China University of Mining&Technology,2008,37(1):38-42.

[9] 李宗翔,衣 刚,武建国,等.基于“O”型冒落及耗氧非均匀采空区自燃分布特征[J].煤炭学报,2012,37(3):484-489.Li Zongxiang,Yi Gang,Wu Jianguo,et al.Study on spontaneous combustion distribution of goaf based on the“O”type risked falling and non-uniform oxygen[J].Journal of China Coal Society,2012, 37(3):484-489.

[10] 李宗翔.自燃采空区耗氧-升温的区域分布特征[J].煤炭学报,2009,34(5):667-672.

Li Zongxiang.Study on distribution characteristic of remaining coal oxygen consumption and spontaneous combustion heating-up in goaf [J].Journal of China Coal Society,2009,34(5):667-672.

[11] 李宗翔,贾进章,周志林.通风换向采空区场量分布变动过程数值模拟研究[J].北京科技大学学报,2010,32(6):691-696.

Li Zongxiang,Jia Jinzhang,Zhou Zhilin.Numerical simulation of the change process of field variables distribution in goaf caused by air reversing in working faces[J].Journal of University of Science and Technology Beijing,2010,32(6):691-696.

[12] Speziale C G,Gatski T B,Fitzmaurice N.An analysis of RNG-based turbulence models for homogeneous shear flow[J].Physics Fluids A,1991,3(9):2278-2281.

[13] Yakhot V,Orszag S A,Thangam S.Development of turbulent models for shear flows by a double expansion technique[J].Physics Fluids A,1992,4(7):1510-1520.

[14] Ren Tingxiang.CFD modelling of Iongwall goaf gas flow to improve gas capture and prevent goaf self-heating[J].Journal of Coal Science&Engineering(China),2009,15(3):15-19.

[15] 张西斌,张 勇,刘传安,等.基于采空区瓦斯运移规律的抽采钻场设计[J].煤炭科学技术,2012,43(3):56-61.

Zhang Xibin,Zhang Yong,Liu Chuanan,et al.Based on the extraction of the law of goaf gas migration from drilling design[J].Coal Science and Technology,2012,43(3):56-61.

Simulation of gas migration in 3D goaf based on RNG k-ε turbulence model

LI Zong-xiang1,2,GU Run-hong3,ZHANG Xiao-ming1,2,BI Qiang4,WEN Yong-yu4

(1.College of Safety Science and Engineering,Liaoning Technical University,Fuxin 123000,China;2.Key Laboratory of Mine Thermodynamic Disasters and Control of Ministry of Education,Liaoning Technical University,Fuxin 123000,China;3.No.5 Engineering Company Ltd.of CCCC First Harbor Engineering Company Ltd.,Qinhuangdao 066002,China;4.Ventilation Department of Tiefa(Group)Limited Liability Corporation,Diaobingshan 112700,China)

For the study of gas floating characteristics and distribution in 3D goaf,according to mass conservation,conservation of momentum,and Fick’s law,suitable for fully mechanized goaf,describing the variable density of the airgas mixture,nonlinear flow-diffusion equations were built.Combined N2-706 coal face in Daxing Mine,using CFD software,the author simulated the transport and distribution of gas in 3D goaf.Heterogeneity of caving in flow field was described as the“O”type.RNG k-ε turbulence model was used,and taking into account the gravity field conditions.The results of theoretical calculations were consistent with the actual conditions and gas monitoring data.Authors point out that gas floating characteristics and distribution in 3D goaf,is due to the air flow moves,gas diffusion-dispersion, and the floating migration of gas-flow which caused by its density difference,and it is also a dynamic balance,which depends on constant desorption of gas in flow field and continuous inflow of air leakage.The results show that the flux of extracting gas from high drill holes is approximately inversely proportional to the volume percent concentration of extracting gas,and has a negative exponential relationship with volume percent concentration of gas in return air tunnel.Key words:3D goaf;gas floating;CFD simulation;k-ε turbulence model;gas drainage

TD752.2;TD712

A

0253-9993(2014)05-0880-06

李宗翔,顾润红,张晓明,等.基于RNG k-ε湍流模型的3D采空区瓦斯上浮贮移[J].煤炭学报,2014,39(5):880-885.

10.13225/j.cnki.jccs.2013.0640

Li Zongxiang,Gu Runhong,Zhang Xiaoming,et al.Simulation of gas migration in 3D goaf based on RNG k-ε turbulence model[J].Journal of China Coal Society,2014,39(5):880-885.doi:10.13225/j.cnki.jccs.2013.0640

2013-05-16 责任编辑:毕永华

国家自然科学基金资助项目(51074086,51174109)

李宗翔(1962—),男,黑龙江绥化人,教授,博士生导师。E-mail:lzx6211@126.com