摘要：为模拟矿井采空区漏风引起的煤自燃与瓦斯问题，建立了采空区“瓦斯-氧”对流扩散过程的迎风有限元统一求解模型。分析了基于伽辽金法的采空区流场有限元求解技术。给出了采空区“瓦斯-氧”对流扩散过程的统一模型，针对标准伽辽金有限元法求解统一模型时计算结果振荡的问题，采用迎风有限元方法得到模型的稳定解。基于三角形单元离散的采空区，给出了以所有节点的压力以及气体浓度为未知数的控制方程组。使用VC++开发了模拟软件，采用典型算例对软件的稳定性与实用性进行了验证，模拟结果表明：煤自燃耗氧导致采空区氧浓度降低，可爆瓦斯浓度范围5%～15%不再适用，采空区内具有瓦斯爆炸危险性覆盖的区域随遗煤耗氧与瓦斯涌出情况而动态变化;开发的软件可以直接划定采空区内煤自燃三带及具有爆炸危险性的采空区区域。
关键词：矿业工程;采空区漏风;煤自燃;瓦斯;迎风有限元
中图分类号：TD 75+2.2文献标志码：A
DOI：10.13800/j.cnki.xakjdxxb.2019.0208文章编号：1672-9315（2019）02-0234-07
0引言
煤矿采空区漏风容易引发煤自燃，形成可爆瓦斯区域[1]甚至引发瓦斯爆炸。由于对采空区进行直接观测十分困难，数值模拟成为国内外学者研究采空区漏风问题的有效手段，如判定采空区内的煤自燃三带[2-3]，模拟煤自燃过程[4-5]及瓦斯运移规律[6-7]。目前模拟采空区漏风问题的途径主要有采用商业CFD与自编程2种，相关的数值方法包括有限元、有限体积法等。尽管商业CFD（如Fluent，COMSOL等）已被大量学者使用，但它对理论知识要求高，目前还很难被现场人员所接受。自编程研究采空区漏风问题很可能是将CFD技术推向现场应用的有效途径，因为通过自编程可将矿井通风问题的相关概念作为软件的界面元素进行设计，将CFD算法封装到后台[8]，以此来降低CFD技术的使用门槛。有限元法相对其它数值方法具有对网格划分要求低、易编程的特点，许多偏微分方程都可通过标准的伽辽金法拆分成代数方程组，这些特点使其成为自编程研究采空区漏风问题的主要数值方法：章梦涛、赵阳升等给出了采空区线性渗流模型的有限元求解方法[9-10];丁广嚷、李宗翔等研究了采空區非线性渗流的有限元求解模型[11-12]，基于Matlab开发了可以模拟采空区瓦斯涌出与耗氧情况[13]的G3程序。采空区漏风引起的瓦斯运移与遗煤耗氧驱动的氧浓度场变化都属于“对流-扩散”过程。在对流起主导作用的采空区流体输运过程中，基于标准伽辽金法导出的有限元控制方程给出的是一个不正确的振荡解，Zienkiewicz，章本照等的研究表明这个问题可以通过迎风有限元法来解决[14-15]。尽管采空区中瓦斯与氧气的运移过程都是“对流-扩散”问题，但在数学模型描述中具有不同的源项表现，采用普通对流扩散模型的控制方程进行求解迭代量大，并不高效。文中基于迎风有限元法思想，导出采空区“瓦斯-氧”对流扩散的统一模型，通过自编程将其用于煤自燃三带与可爆瓦斯浓度覆盖区域的判定。
4应用实例
4.1模型实现与算例概况
文中基于VC++与ObjectARX[20]技术将以上模型及前后处理过程集成于AutoCAD中，方程组的求解采用PARDISO[21-22]并行求解器。应用实例为某矿U形工作面通风系统采空区，宽180 m，煤层开采厚度6 m，工作面设计风量为25 m3/s，入风隅角与回风流隅角的相对压力分别为308.1和268.2 Pa;工作面推进度为3 m/d，图2为回采90 d后形成的270 m长采空区（按渗透系数渲染）。模拟中涉及的其它参数取值见表2.
4.2计算结果
图3采用流网显示采空区流场的模拟结果，其中流线的流函数值间距为1.9×10-3，等压线间距为0.8 Pa，流网颜色按流速大小渲染，绿、黄、蓝分别表示按流速划定的采空区自燃三带：散热带（流速大于0.004 m/s），氧化升温带（流速位于0.001 7至0.004 m/s）和窒息带（流速小于0.001 7 m/s）。图4（a）为模拟得到的采空区的氧浓度，采用氧浓度划分三带时，通常将浓度介于8%至18%的区域划定了氧化升温带，如图中黄色区域。由于风速过大将使热量无法积聚，因此更合理的方法是用风速与氧浓度双指标进行划分，即氧化升温带应满足
模拟结果中三带的分布及自燃带的形状与文献[23]采用FLUENT模拟的结果相似：采空区自燃带位于工作面后一定位置，其宽度沿进风侧向回风侧逐渐减小。
图5为模拟得到的采空区瓦斯浓度场。蓝色显示采空区深部是瓦斯富积的区域，浓度大于15%，绿色是浓度低于5%的区域，黄色区域瓦斯浓度介于5%～15%、在正常空气氧浓度下具有爆炸危险。从图3，图5可知漏风主要从2个隅角流入采空区与流回工作面，工作面与采空区边界内侧、近上隅角的49.5 m长范围内（回风隅角）瓦斯浓度大于5%（5%～23%）。由于遗煤耗氧，使得采空区中氧浓度降低，不一定在任何地方都满足瓦斯爆炸的需要。图6左图的瓦斯爆炸三角形给出了瓦斯与氧气形成可爆气体的浓度范围[24]，当氧浓度低于12%时，将不满足瓦斯爆炸的需要。据此将图4（a）与图5叠加，得到图6右图红色区域所示的具有瓦斯爆炸危阶性区域，其中绿色、黄色分别表示瓦斯含量过低和过高的区域。因此在图6右图所示的红色区域内，任何点火源（包括煤自燃）均可引发瓦斯爆炸。因此针对模拟结果，对图中标定的危险区域进行惰化，以及保证必要的推进度来防止煤自燃对避免采空区火灾或爆炸事故的发生是十分必要的。
1）基于伽辽金法建立了采空区流场的求解模型。给出了工作面动态回采下采空区渗流参数的计算方法。针对采空区流场的定解模型，基于三角形单元离散的采空区，采用标准伽辽金法导出了以所有节点压力为未知数的有限元控制方程组，为自编程研究采空区漏风问题提供了基础;
2）基于迎风有限元法建立了采空区“瓦斯-氧”输运过程的统一模型。分析了采空区瓦斯涌出与遗煤耗氧的特点，建立了采空区中“瓦斯-氧”输运过程的对流扩散模型。针对标准伽辽金法求解对流扩散模型时解的振荡问题，给出了迎风有限元法求解“瓦斯-氧”输运过程的控制方程组;
3）开发了采空区漏风问题模拟软件。基于以上模型，以及VC++，ObjectARX等技术开发了可视化模拟软件，应用实例表明使用此软件可方便地进行采空区模型的前后处理，并根据模拟结果给出采空区煤自燃三带及具有瓦斯爆炸危险性的区域。
参考文献（References）：
[1]Brune J F.Methaneair explosion hazard within coal mine gobs[J].Transactions of the Society for Mining，Metallurgy，and Exploration，2013，334：376-390.
[2]文虎，张泽，赵庆伟，等.煤层分层前后采空区自燃“三带”的数值模拟[J].煤矿安全，2017，48（3）：178-191.
WEN Hu，ZHANG Ze，ZHAO Qingwei，et al..Numerical simulation on spontaneous combustion “Three Zones” of goaf before and after coal seam layering[J].Safety in Coal Mines，2017，48（3）：178-191.
[3]时国庆，胡方坤，王德明，等.采空区自燃“三带”分布规律的四维动态模拟[J].中国矿业大学学报，2014，43（2）：191-194.
SHI Guoqing，HU Fangkun，WANG Deming，et al.Unsteady simulation on distribution of three zones for spontaneous combustion in goaf areas[J].Journal of China University of Mining & Technology，2014，43（2）：191-194.
[4]YUAN Liming，Smith Alex C.Numerical study on effects of coal properties on spontaneous heating in longwall gob areas[J].Fuel，2008，87（15-16）：3409-3419.
[5]朱紅青，刘星魁.尾巷瓦斯抽采下采空区煤自燃升温的数值模拟[J].西安科技大学学报，2012，32（1）：1-7.
ZHU Hongqing，LIU Xingkui.Numerical simulation of heating up of coal spontaneous combustion in gob area under tail roadway methane drainage pattern[J].Journal of Xi’an University of Science and Technology，2012，32（1）：1-7.
[6]XIA Tongqiang，ZHOU Fubao，WANG Xinxin，et al.Controlling factors of symbiotic disaster between coal gas and spontaneous combustion in longwall mining gobs[J].Fuel，2016，182：886-896.
[7]XIA Tongqiang，XU Mingjing，WANG Yonglong.Simulation investigation on flow behavior of gob gas by applying a newly developed FE software[J].Environmental Earth Sciences，2017，76（14）：485.
[8]高科，刘剑，邓立军.集于AutoCAD的采空区风火瓦斯耦合软件开发研究[J].中国安全生产科学技术，2015，11（5）：97-103.
GAO Ke，LIU Jian，DENG Lijun，et al.Development of airfiregas coupling software for goaf based on AutoCAD[J].Journal of Safety Science and Technology，2015，11（5）：97-103.
[9]章梦涛，王景琰，梁栋.采场大气中沼气运移过程的数值模拟[J].煤炭学报，1987，11（3）：23-30.
ZHANG Mengtao，WANG Jingyan，LIANG Dong.Numerical modelling of methane migration in air in underground workings[J].Journal of China Coal Society，1987，11（3）：23-30.
[10]赵阳升.有限元法及其在采矿工程中的应用[M].北京：煤炭工业出版社，1994.
[11]丁广嚷，柏发松.采空区混合气体运动基本方程及其有限元解法[J].中国矿业大学学报，1996，25（3）：21-26.
DING Guangxiang，BAI Fasong.Basic equations of mixedgas motion in gob and their finite element solution[J].Joual of China University of Mining & Technology，1996，25（3）：21-26.
[12]李宗翔，刘玉洲，吴强.采空区流场非线性渗流的改进迭代算法[J].重庆大学学报，2008，31（2）：186-190.
LI Zongxiang，LIU Yuzhou，WU Qiang.Improved iteration algorithm for nonlinear seepage in flow field of caving goaf[J].Journal of Chongqing University，2008，31（2）：186-190.
[13]李宗翔，王晓冬，王波.采空区场流数值模拟程序（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，2005，20（3）：16-20.
[14]Zienkiewicz O C，Taylor R L，Nithiarasu P.The finite element method for fluid dynamics（7th edition）[M].Waltham：ButterworthHeinemann，2014.
[15]章本照.流体力学中的有限元方法[M].北京：机械工业出版社，1986.
[16]李宗翔.高瓦斯易自燃采空区瓦斯与自燃耦合研究[D].阜新：辽宁工程技术大学，2007.
LI Zongxiang.Study of coupling of gas and spontaneous combustion in highly gassy and spontaneous combustion goafs[D].Fuxin：Liaoning Technical University，2007.
[17]Bear J.Dynamics of fluids in porous media[M].New York，USA：America Elsevier Publishing Company Inc.，1972.
[18]秦跃平，刘伟，杨小彬.基于非达西渗流的采空区自然发火数值模拟[J].煤炭学报，2012，37（7）：1177-1183.
QIN Yueping，LIU Wei，YANG Xiaobin.Numerical simulation of impact of nonDarcy seepage on spontaneous combustion in goaf[J].Journal of China Coal Society，2012，37（7）：1177-1183.
[19]Heiich J C，Huyakorn P S，Zienkiewicz O C，et al.An upwind finite element scheme for twodimensional convective transport equation[J].International Journal for Numerical Methods in Engineering，1977，11（1）：131-143.
[20]吳奉亮，周澎，李晖，等.基于智能对象的通风CAD模型研究[J].煤炭科学技术，2009，37（5）：54-57.
WU Fengliang，ZHOU Peng，LI Hui，et al.Research on the model of mine ventilation CAD based on the intelligent object[J].Coal Science and Technology，2009，37（5）：54-57.
[21]Schenk，Grtner K.Solving unsymmetric sparse systems of linear equations with PARDISO[J].Journal of Future Generation Computer Systems，2004，20（3）：475-487.
[22]曹大志，强洪夫，任革学.基于OpenMP和Pardiso的柔性多体系统动力学并行计算[J].清华大学学报（自然科学版），2012，52（11）：1643-1649.
CAO Dazhi，QIANG Hongfu，REN Gexue.Parallel computing studies of flexible multibody system dynamics using OpenMP and Pardiso[J].Journal of Tsinghua University（Sci. & Tech.），2012，52（11）：1643-1649.
[23]张辛亥，席光，徐精彩，等.基于流场模拟的综放面自燃危险区域划分及预测[J].北京科技大学学报，2005，27（6）：641-644.
ZHANG Xinhai，XI Guang，XU Jingcai，et al.Partition of spontaneous combustion dangerous zone and prediction of spontaneous combustion based on numerically modeling the flow field at fully mechanized caving face[J].Journal of University of Science and Technology Beijing，2005，27（6）：641-644.
[24]Coward H F，Jones G W.Limits of flammability of gases and vapors[J].Journal of the tranklin Institute，1928，203（1）：161.
[25]CHENG Jianwei，YI Luo.Modeling atmosphere composition and determining explosibility in a sealed coal mine volume[J].Achieves of Mining Science，2014，59（1）：25-40.