吴奉亮，李智胜，常心坦

（1.西安科技大学安全科学与工程学院，陕西 西安710054；2.西安科技大学西部矿井开采与灾害防治教育重点实验室，陕西西安710054）

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］。尽管采空区中瓦斯与氧气的运移过程都是“对流-扩散”问题，但在数学模型描述中具有不同的源项表现，采用普通对流扩散模型的控制方程进行求解迭代量大，并不高效。文中基于迎风有限元法思想，导出采空区“瓦斯-氧”对流扩散的统一模型，通过自编程将其用于煤自燃三带与可爆瓦斯浓度覆盖区域的判定。

1 采空区渗流特性参数分布规律

采空区流场的渗透系数E，孔隙率ε及冒落高度H等参数主要受采空区的冒落煤岩碎胀系数影响，研究表明碎胀系数呈“O”形圈分布，以图1所示U形工作面通风系统为例，各参数可用以下公式计算［16］

图1 U形通风工作面采空区布置图Fig.1 Gob layout of U-shape working face ventilation system

式中 Kp（x，y）为给定点（x，y）的碎胀系数；Kp，max和Kp，min分别为初始和压实后的碎胀系数；x0为开切眼到Y轴的距离；u为回采速度；t为回采时间；a0和a1为倾斜与走向方向冒落煤岩碎胀衰减系数；d0=y-｜y-0. 5W｜，其中 W为采空区宽度；ξ为“O”形圈的几何调整系数；h为采高；υ为空气运动粘性系数；g为重力加速度；d为粒子平均直径。

2 采空区漏风问题的定解模型

2. 1 采空区流场非线性渗流定解模型

采空区中的风流在冒落煤岩体中的流动是各向同性非均匀介质渗流，由于采空区宽度与深度远大于高度，因此三维采空区流场完全可简化为如图1所示的平面流场。若先假定其流动符合达西定律，则流速V，渗透系数E和压力P之间的关系为

进一步考虑连续性方程，得到平面采空区流场的定解模型为（6）式微分方程［9-10］

式中 G为整个采空区区域，采空区与风网连接处的边界L1为压力边界，p0为L1上的已知压力值；L2为不漏风的边界；nx，ny为边界外法线上单位向量的分量。实际上采空区内，特别是邻近工作面的范围，一般为非线性的渗流，其流动符合Bachmat方程［17］

式中 β为粒子形状系数。采用对渗透系数进行迭代修正的方法，可基于（6）式的线性模型求解以上非线性渗流，详细的迭代原理及过程见文献［12］。

2. 2 采空区中“瓦斯-氧”运移的定解模型

2. 2. 1 采空区瓦斯涌出与耗氧强度

采空区中单位厚度上的瓦斯涌出强度（mol/（m2h））计算式为［16］

式中 w′0和w0分别为采空区局部地点恒定的瓦斯涌出强度和下部煤层沿底板裂隙持续涌入瓦斯的强度，mol/（m2h）；w1和 w2是采空区遗煤和邻近煤层瓦斯初始涌出强度，mol/（m2h）；λi是瓦斯涌出量的衰减系数，d-1.根据 Arrhenius公式［4］，煤的耗氧速率WO2可用下式描述

式中 A为指前系数；CO2为氧气摩尔浓度，mol·m-3；Ea为表观活化能，kJ/mol；R是气体常数，J/（mol·Ｋ）；T绝对温度，Ｋ.由于式中exp（-Ea／RT）项在煤自燃的早期阶段变化很小，因此上式可简化为

即煤的耗氧速度与氧浓度成线性关系（A′是线性系数），与文献［2］中通过实验得到的结论相同。考虑遗煤厚度H1的影响，则采空区中煤的耗氧速率应为

2. 2. 2 采空区中“瓦斯-氧”运移的数学模型

由于工作面回采速度远小于采空区内气体对流扩散的速度，因此采空区中的瓦斯与氧气运移可视为稳态“对流-扩散”过程［18］，描述这一过程的偏微分方程及定解边界条件为

Θ为 CH4或 O2中的一种；CΘ为气体浓度；V′Θ为采空区气体流速，V′Θ=V／ε；DΘ为气体扩散系数；IΘ为源项中与CΘ无关的部分，sΘ为源项与CΘ存在线性关系的系数，IΘ与 sΘ的取值分别见表1；ГC为采空区的入风边界，Гq表示回风边界及不漏风的边界；D为气体扩散系数张量。D的计算式为［9］

式中 i，j取值 1或 2，分别为 x，y方向；D′为分子扩散系数，m2/s；τ为介质曲率；当 i=j时，δij=1，否则δij=0；aT，aL分别为纵向与横向弥散度。

表1 公式（12）中IΘ和 sΘ的不同表现形式Table 1 Different forms of IΘand sΘ in Eq. （12）

3 “瓦斯-氧”对流扩散迎风有限元解算法

3. 1 用伽辽金法求采空区流场

若将待求采空区问题（如（6）式中的压力P或（12）式中的浓度 CΘ）的真解记为 U（x，y），则它的有限元解是指将采空区离散成单元（文中采用三角形单元，设节点总数为n）后，用每个单元上的插值函数 UΔ（x，y）来近似 U（x，y），这个插值函数为

式中 ui，uj，uk为三角形单元的3个节点（按逆时针方向编号为 i，j，k）对应问题的待求量；Nη（η=i，j，k）是基函数，记 AΔ为单元面积，Ni=（ai+bix+ciy）／（2AΔ），Nj=（aj+bjx+cjy）／（2AΔ），Nk=（ak+bkx+cky）／（2AΔ），ai=xjyk-xkyj，bi=yj-yk，ci=xk-xj，aj=xkyi-xiyk，bj=yk-yi，cj=xi-xk，ak=xiyj-xjyi，bk=yi-yj，ck=xj-xi. 因此有限元方法的关键步骤是导出关于待求问题 ui（i=1，2，…，n）为未知数的控制方程组。加权余量法是最方便实施的一种方法，当在单元内为每个节点取定对应的权函数 ωη（η=i，j，k）时，对（6）或（12）式对应的问题导出的控制方程组的矩阵表示形式均可记为

伽辽金法是将权函数取为与基函数相同的一种加权余量法，对（16）式取ωη=Nη，积分得到［0，0，0］T，单元矩阵为一对称阵，其计算式为

3. 2 采空区“瓦斯 -氧”运移的迎风有限元法

在（18）中取权函数ωη为基函数Nη（即采用伽辽金法）时，导出方程组的解是振荡的［19］，解决办法是对权函数进行修正以增大来流方向的效应，修正式为［14］

式中 hΔ为单元的特征长度，可取三角形单元的最长边；α为由佩克莱特数Pe决定的系数，它们的计算式分别为

将（19）式代入（18）式积分后得到

式中 Sij项是源项中与浓度存在线型关系的项的积分结果，它的计算式为

综上，在求得采空区流场后，将（23）式代入（15）式即可求得采空区气体的浓度场。

4 应用实例

4. 1 模型实现与算例概况

文中基于 VC++与 ObjectARX［20］技术将以上模型及前后处理过程集成于AutoCAD中，方程组的求解采用PARDISO［21-22］并行求解器。应用实例为某矿U形工作面通风系统采空区，宽180 m，煤层开采厚度6m，工作面设计风量为25m3/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%的区域划定了氧化升温带，如图中黄色区域。由于风速过大将使热量无法积聚，因此更合理的方法是用风速与氧浓度双指标进行划分，即氧化升温带应满足

图2 采空区渗透系数Fig.2 Distribution of hydraulic conductivity in gob area

表2 模拟参数Table 2 Parameters for the numerical simulations

式中 CO2-min=8%，Vmax=0.004 m/s.按（25）式将图3与图4（a）叠加得到图4（b）所示的三带（仍用绿、黄、蓝区分）划分结果。模拟结果中三带的分布及自燃带的形状与文献［23］采用FLUENT模拟的结果相似：采空区自燃带位于工作面后一定位置，其宽度沿进风侧向回风侧逐渐减小。

图3 采空区流网Fig.3 Flow net in gob area

图4 采空区氧浓度场及自燃三带Fig.4 Distribution of oxygen concentration and“three zones”for analyzing spontaneous combustion in gob

图5为模拟得到的采空区瓦斯浓度场。蓝色显示采空区深部是瓦斯富积的区域，浓度大于15%，绿色是浓度低于5%的区域，黄色区域瓦斯浓度介于5%～15%、在正常空气氧浓度下具有爆炸危险。从图3，图5可知漏风主要从2个隅角流入采空区与流回工作面，工作面与采空区边界内侧、近上隅角的49.5 m长范围内（回风隅角）瓦斯浓度大于5%（5%～23%）。由于遗煤耗氧，使得采空区中氧浓度降低，不一定在任何地方都满足瓦斯爆炸的需要。图6左图的瓦斯爆炸三角形给出了瓦斯与氧气形成可爆气体的浓度范围［24］，当氧浓度低于12%时，将不满足瓦斯爆炸的需要。据此将图4（a）与图5叠加，得到图6右图红色区域所示的具有瓦斯爆炸危阶性区域，其中绿色、黄色分别表示瓦斯含量过低和过高的区域。因此在图6右图所示的红色区域内，任何点火源（包括煤自燃）均可引发瓦斯爆炸。因此针对模拟结果，对图中标定的危险区域进行惰化，以及保证必要的推进度来防止煤自燃对避免采空区火灾或爆炸事故的发生是十分必要的。

图5 采空区瓦斯浓度场Fig.5 Field ofmethane concentration in gob area

图6 瓦斯爆炸三角形（左）及采空区具有瓦斯爆炸危阶的区域（右）Fig.6 Methane explosibility triangle and presumed location of the explosive gas zones in the gob

5 结 论

1）基于伽辽金法建立了采空区流场的求解模型。给出了工作面动态回采下采空区渗流参数的计算方法。针对采空区流场的定解模型，基于三角形单元离散的采空区，采用标准伽辽金法导出了以所有节点压力为未知数的有限元控制方程组，为自编程研究采空区漏风问题提供了基础；

2）基于迎风有限元法建立了采空区“瓦斯-氧”输运过程的统一模型。分析了采空区瓦斯涌出与遗煤耗氧的特点，建立了采空区中“瓦斯-氧”输运过程的对流扩散模型。针对标准伽辽金法求解对流扩散模型时解的振荡问题，给出了迎风有限元法求解“瓦斯-氧”输运过程的控制方程组；

3）开发了采空区漏风问题模拟软件。基于以上模型，以及VC++，ObjectARX等技术开发了可视化模拟软件，应用实例表明使用此软件可方便地进行采空区模型的前后处理，并根据模拟结果给出采空区煤自燃三带及具有瓦斯爆炸危险性的区域。