采空区注氮软件模拟及开发研究

王月红，吴 怡，张九零，关雅洁

(1.华北理工大学 矿业工程学院，河北 唐山 063009；2.河北省矿山开发与安全技术重点实验室，河北 唐山 063009)

由于采空区中存在的大量遗煤和适宜漏风环境为煤氧化创造了条件，增大了自燃的危险性[1]。而注氮作为防治遗煤自燃的有效措施之一，它具有稀释采空区氧浓度、冷却遗煤、增大采空区压力、减少漏风量、延长自然发火期等优点，近年来被广泛应用[2-5]。

由于数值模拟成本低廉，只需在计算机上进行模拟和数据处理，有时可达到理想化克服现场的困难。国内外学者关于模拟采空区注氮技术的成果很多。在模拟漏风流场和温度场中，李宗祥[6-7]针对Y形通风的工作面，综合各种传热方程，并对模型参数及边界设计后进行注氮问题的数值模拟，得到了各种气体(瓦斯、CO、氧气等)、温度和漏风流的分布情况及变化规律，并且计算出了注氮与工作面推进速度的配比关系。美国学者Yuan等[8]根据CFD模型研究采空区漏风流场变化规律，并通过研究结果完善现场通风技术，控制灾害的发生。罗新荣[9]采用Fluent软件建立物理模型和数学模型，设置相应的参数，并根据采空区渗透率和遗煤氧化范围编制自定义函数，模拟不同抽采量和注氮形式对采空区温度场的影响规律。在关于注氮前后氧浓度的变化规律上，刘玉良[10-11]将理论和数值模拟结合对浅埋深煤层注氮技术进行研究，比较了注氮前后“氧化带”的范围。澳大利亚学者Ren[12]利用软件模拟详细研究了瓦斯的运移规律和氧浓度分布情况，并结合理论知识控制自然发火问题。唐冠楚[13-15]在CFD模型下，模拟不同的抽采方式与注氮技术相结合自燃危险区域的分布情况，在瓦斯抽采效率提高的同时防止采空区遗煤自燃的发生。

但是，国内外均用流体动力学软件进行模拟，不能更针对性地研究采空区注氮防灭火技术。笔者基于有限体积法，利用VB语言编制注氮防灭火模拟系统，通过Tecplot软件作为图显，进行温度场和氧浓度场的模拟。利用编制的软件模拟采空区注氮，减少了现场的大量实验，应用于各种类型下的采空区注氮问题。

1 理论基础

在采煤工作持续进行中，工作面向前推进，相对区域也持续扩大，因此选取一部分采空区作为研究对象，引入动态坐标。在注入氮气后，采空区气体压力、氧气浓度、温度均会产生影响。基于有限体积法建立数学模型，利用数学语言对实际问题进行定量化，使指标体系可计算。

1.1 有限体积法

有限体积法(FVM)是将计算区域划分为一系列不重复的控制体积，使每个网格点周围有一个控制体积；将待解的微分方程对每一个控制体积积分，然后对积分式进行离散化处理，再导出离散化方程。

有限体积法在流体流动和传热数值计算领域被广泛应用，与其他有限元和边界方法比较，其基本思想简洁，能够在物理意义上表示控制体的通量平衡，建立的模型方程意义明确。因此，将有限体积法应用到采空区注氮防灭火的模型构建上，为解决矿井火灾易复燃、难扑灭等问题提供了有效支撑。

1.2 采空区注氮氧浓度场数学模型

根据质量守恒定律：在单位时间内，漏风风流通过任意封闭区域内氧气质量的变化量(WC)主要由以下几个方面引出：由于渗透风流中的氧气流出与流入任意封闭曲面之差(W1)；氧气消耗量(W2)；浓度差异引起的弥散进出任意封闭曲面的氧气之差(W3)，如公式(1)，(2)，(3)所示。

WC=W1+W2+W3

(1)

(2)

(3)

因为在采空区存在气体的动力弥散现象，所以氧气的扩散系数与风流速度大小成正比，故W3：

(4)

(5)

即将式(2)，(3)，(4)以及(5)代入(1)得：

(6)

在实际情况中，氧气流动速率要大于采空区工作面的回采速度，即单位时间内选取的微元ΔS内的孔隙率和氧气密度变化可以忽略不计，即式(6)可变成：

(7)

为了和Fick动律保持一致，将(7)式中质量浓度改为摩尔浓度。

(8)

式中，n为有效孔隙度；ρO2为氧气的密度，kg/m3；CO2为氧气的浓度，mol/m3；kO2为氧气的扩散系数常数；vx为x方向气体渗流速度，m/s；vy为y方向气体渗流速度，m/s；v为气体的渗流速度，m/s；u(t)为松散煤体内的氧气消耗速度，mol/(s·m3)。

确定边界条件Γ1，Γ2，Γ3，Γ4后，得到氧浓度场的数学模型，如(9)所示。

(9)

式中，CO2(x，y)为回采工作面连接的边界上的氧气浓度函数，mol/m3。

根据模型可以分析出氧浓度在注氮整个过程中的变化，从而研究出注氮对采空区火灾的影响作用。

1.3 采空区注氮温度场数学模型

在采空区中，除了岩石间导热的过程外，还有气体与气体间互相导热，它们导热的方式是有区别的，除此之外，固体间的导热还会受到气体的影响，因此在采空区注氮模型中的温度场分为固体温度场和气体温度场。

1.3.1 固体温度场

在一定的时间内，导入和导出冒落岩石的热量差Qsc、冒落岩石氧化放热的热量差Qf之和，也就是进入控制体的热量，应该与导出控制体的热量相等，也就是冒落岩石及孔隙中气体的对流换热量Qd和冒落岩石移动带出控制体的热量差Qv之和。即为公式(10)所示。

Qsc+Qf=Qd+Qv

(10)

由傅里叶定律得知，由边界Γ流入曲面F全部热量Qsc：

(11)

式中，k为渗透率，m2。

单位时间内控制体的氧化放热量为：

(12)

在松散煤体内同时存在对流换热和导热两种过程，这两种热量传递方式也相互影响，同时还受气体运动规律的影响，因此这是一个十分复杂的过程，利用牛顿冷却公式进行计算，如公式(13)所示。

(13)

(14)

将式(11)，(12)，(13)，(14)代入式(10)得：

(15)

根据傅里叶定律可知：

(16)

将式(16)代入(15)得：

(17)

式中，λy为采空区冒落岩石的导热系数，W/(m·℃)；ts为冒落岩石温度，K；qx为沿x轴方向的热流通量；qy为沿y轴方向的热流通量；Ke为对流换热系数，J/(m2·s·K)；Sn为单元体内冒落岩石与气体对流换热的表面积，m2；tg为气体温度，K；ts为冒落岩石温度，K；q(t)为单位时间单位体积内冒落岩石的放热量，kJ/(m3·s)；ρs为固体的密度，kg/m3；cs为固体的比热容kJ/(kg·K)；a为单元体的比表面积,a=Sn/ΔxΔyΔz，L/m；v0为工作面推进速度，m/s。

即式(17)是移动坐标下采空区内冒落岩石能量方程。

1.3.2 气体温度场

设有温度场Qg(Q,t)=Qg(x,y,t)，由能量守恒定律可知，单位时间内，汇入封闭曲线范围的气体热量Qgc，同控制体体内固体对流散热量Qd之和，应等于单元体内孔隙中气体内能的增加值Qh，即：

Qgc+Qd=Qh

(18)

由傅里叶定律得知，单位时间内曲面边界Γ流入曲面F全部热量是Qgc：

(19)

热量Qh可表示为：

(20)

对流换热量Qd可表示为：

(21)

将式(19)，(20)，(21)代入(18)得：

(22)

根据傅里叶定律可知：

(23)

将(23)代入式(22)得：

(24)

式中，λg为采空区气体导热系数，W/(m·℃)；Ke为对流换热系数，J/(m2·s·℃)；Sn为单元体内冒落岩石与气体对流换热的表面积，m2；tg为气体温度，℃；ts为固体温度，℃；a为单元体的比表面积，a=Sn/ΔxΔyΔz，L/m；vx，vy为气流沿x轴、y轴方向的分量，m/s；ρg为气体的密度，kg/m3；cg为气体的比热，J/(kg·℃)。

式(24)是移动坐标下采空区内气体温度场能量方程。

2 软件开发

2.1 软件界面

图形显示模块的任务是把整个模拟过程处理存储的数据用图形显示出来，具有一个功能齐全的图形显示界面是有限体积软件中最基本的一个要求。对于大多数的用户而言，图形显示功能是至关重要的，因为图形显示能够有效地减小大量数据结构的堆积，进而减小了结构计算的错误。本设计软件图形显示模块分为3部分：登陆系统；参数输入界面；结果汇总显示。

显示界面是Tecplot10.0作为绘图以及数据分析情况下通用软件，该软件作为图形显示的优点是易学易用，界面友好，能够马上读入*.cas、*.dat文件等。因此基于Tecplot软件的这些优点，使其作为该软件的图形显示模块效果是非常可观的，而且该软件较成熟，解决了不能准确地处理数据的弊端。

2.2 软件功能

该软件的开发，不仅对采空区自然发火防治的研究起到了极大的作用，更是对注氮防灭火这一灭火技术有了更深入的研究。该软件的开发可应用到实际的生产中，通过各煤矿的实际情况，解算出注氮的最佳位置、最佳流量、最佳个数等实质性问题；还可计算采空区注氮条件下“三场”(氧浓度场、流场、温度场)的分布情况；避免了实际操作的困难，取而代之的是一次次的计算机上的操作，相对比人工操作精确而且可循环，试验性强，计算精度高，速度快，对煤矿解决火灾这个难题起到了极大的帮助。

3 数值模拟

3.1 网格划分

首先通过计算一个控制体，达到计算出构成采空区工作面所有单元的相关参数。为了清楚地表示出网格划分的方法，用图1所示MNPQ代表中心体积法中的矩形网格，四边形abcd表示选取的控制体。沿ab方向为采空区深度，沿ad方向为工作面方向，根据采空区自然发火物理模型ab方向为采空区的进风巷，dc为回风巷，在采空区中“两道两线”处为自燃易发部位，因此在对控制体abcd进行网格剖分时在自燃易发部位处的网格加密。采空区其余工作面部位采取疏密不一的网格线进行剖分，为了使得绘制的网格标准化和规范化，网格之间的距离引入数学理论中的等比数列且系数大于1，如图2所示。

图1 有限体积法矩形网格

图2 控制体网格划分

然后在采空区中注入氮气，选取一个注氮点，假设在该点注入氮气，氮气进入采空区后会在该点向四周扩散，划分了如图3所示的第二层网格。如图所示坐标为(e，f)的点的氮气含量系数为(e2+f2)1/2，用DI表示注入的氮气量，则该点G的氮气含量为DIG=DI×(e2+f2)1/2。同理，其余点的氮气量也可据此求出，这样就得出了采空区内所有点的氮气量。

图3 注氮含量计算示意

通过将两个网格叠加之后，从而进行模型的数值计算。由Tecplot作为图显得到图4的网格划分形式。

图4 网格划分

3.2 参数设置

在林南仓煤矿井下采集具有代表性的煤样进行工业分析、元素分析等基础实验，从而确定采空区注氮软件开发基础参数的设定，如表1所示。

表1 模型基础参数

由采空区自然发火原理中的煤氧复合理论得出，空气中的氧气与煤的表面结合后会发生大量的化学反应，并释放大量的热量，同时该过程更会加快煤与氧气的反应进程，使得火灾或爆炸迅速发生。因此，针对这一理论知识，通过升温氧化实验，确定煤氧化放热参数，如表2所示。

表2 煤氧化放热参数

3.3 模拟结果

3.3.1 氧浓度场

对比未注氮(如图5(a))和注氮(如图6(a))的氧浓度图可知，注氮使氧化带变窄，窒息带变宽，且氧化带向工作面方向推进。因为随着漏风气流流入采空区中，氧气持续向深部扩散，但氮气的注入不断稀释氧浓度，使氧浓度下降。且温度在注氮后逐渐降低，减轻了热量积聚放热量降低，使得耗氧速率减小，氧化带宽度则变窄。

图5 未注氮氧浓度模拟

对比未注氮和注氮后的氧浓度场三维图，发现两者的氧浓度均是逐渐降低的，但氧浓度变化趋势则由注氮前的凸形变为注氮后的凹形。观察图5(b)和图6(b)可知，三维图中反曲线的凸率逐渐减小，说明随着距回风巷距离的增大，氧浓度改变的速率逐渐减小。证明了未注氮时在进风巷周围由于漏风流将氧气带入采空区深部，而到达回风巷时，由于逆风的作用，使得氧气浓度的变化不再明显；而注氮后，靠近进风巷的位置处氧浓度随深度的改变呈凹形变化，且随着与进风巷距离的加大，凹率逐渐减小。

图6 注氮后氧浓度场模拟

3.3.2 温度场

对比未注氮(如图7(a))和注氮(如图8(a))二维图的温度场可知，工作面附近的温度较低，且注氮前后的温度差只有2℃，说明注氮前后进风巷中与漏风流的温度相接近。但随着风流流向采空区深部，低温气体与采空区深部的冒落岩石及煤体进行热量交换，温度急剧降低，因此采空区深部的温度差较大。注氮前的最高温度为70℃，注氮后的最高温度仅为40℃，温差达到了30℃左右。且注氮前后采空区温度与深度呈正比例关系，采空区越深部温度越高。注氮后进风巷附近的温度有明显的降低，而回风巷周围的温度较总体有明显的降低，但与注氮前相比较，高温区域由注氮前的采空区深部位置变为回风巷部位。由于注入氮气后，氮气被风流带入采空区深部，发生热传导降低原有遗煤温度，热量则积聚到回风巷；又因回风巷处的风量较小，热量不易被带出采空区，导致回风巷深部位置温度最高。

图7 采空区未注氮温度场模拟

图8 采空区注氮后温度场模拟

图7(b)(注氮前)与8(b)(注氮后)三维图对比发现，不仅注氮与未注氮最高点温度不同，且在进风侧附近，注氮后的温度随着采空区的深度的增加缓慢上升，而在回风侧附近，上升速率较大。因为氮气进入采空区后，由于氮气的低温降低了采空区内遗煤温度上升的速率，而在回风巷附近，由于逆风的作用，氮气的温度较低，所以大量的热量积聚在上隅角的位置处，使得该处的温度很高。

3.4 工程应用

根据林南仓煤矿的技术条件及生产能力确定相应的注氮参数，矿井下通常运用埋管方式进行防灭火，沿回采工作面进风巷巷帮铺设管路，伴随着工作面推进，管段注氮口会逐渐伸入深部，注氮口位置开始进入氧化自燃带时进行注氮，每隔一段步距重新铺设一条管路，在注氮的过程中采用封闭式注氮，即利用密闭墙上的注氮管向采空区中注入氮气。进行现场注氮后发现温度和氧气浓度明显降低，林南仓煤矿在正常生产的条件下是不易发生自燃的，因此注氮达到了相应的效果，由此说明，该软件的模拟结果切实可行，为日后采空区注氮技术研究起到了很大支撑作用。

4 结 论

(1)基于有限体积法，建立采空区注氮“氧浓度场”和“温度场”数学模型。并利用VB语言编制注氮防灭火模拟系统，通过Tecplot软件作为图显，开发了采空区注氮防灭火软件。可运用到煤矿实际生产中，并有效解决煤矿火灾问题。

(2)根据模拟得到的氧浓度场，发现注氮前后氧化带宽度变窄，并且向工作面方向推进，窒息带变宽。氧浓度变化曲线由凸形变为凹形。随着回风巷距离的增大，氧浓度变化速率越来越小。

(3)在温度场中，注氮前后工作面附近的温度差仅2℃，而采空区深部的温度差为30℃，且采空区温度与深度呈正比例关系。注氮后，进风巷附近温度明显降低，而回风巷周围温度较总体明显下降，但与注氮前比较高温区域由深部位置变为回风巷位置。