基于最小流动单元的采空区氮气充注压力损失计算方法

技术领域

本发明属于矿井通风和采空区防灾减灾技术领域，具体是涉及一种基于最小流动单元的采空区氮气充注压力损失计算方法。

背景技术

70％以上的煤矿火灾事故发生在毗邻采区的采空区。采空区遗煤自燃发火及其灭火救灾一直是行业内的热点问题，采用注氮治理自燃发火是最常见和行之有效的灭火措施。但是，现有计算方法和实施细则，无法准确计算灭火用所注氮气在采空区内流动的距离，从而无法准确预测所注氮气能影响的范围，导致液氮浪费和灭火困难。

在确定采空区注氮参数方法及其措施的研究进展方面，李东等人提出了一种采空区U型通风工作面全断面帷幕注氮的防灭火方法，提高了氮气惰化效率。朱红青发明了自动控制旋转牵引式注氮防灭火装置，实现了采空区注氮点在空间上的连续，提升了采空区注氮的惰化效果。秦波涛和鲁义提出了一种高效治理浅埋藏煤层大面积采空区遗煤自燃的方法，集堵漏控风与快速惰化降温为一体，高效防治了遗煤自燃。郭君柳探讨了采空区“三带”的观测以及注氮防灭火的实用方法，提高了煤炭开采的安全性和经济效益。根据质量、动量和组分守恒方程，朱红青建立采空区气体体积分数变化的理论流体力学模型，通过初始化、赋边界值和迭代计算，数值模拟确定最佳注氮位置和注氮量。李宗翔结合段王煤矿进行了“两进一回”复杂采场采空区注氮防灭火模拟，确定了最佳的注氮量和注氮位置。何俊忠通过SF₆示踪气体测定宏岩煤矿4405采空区的漏风量，并确定主要漏风方向，在此基础上优化了采空区注氮工艺。张琪结合大同煤矿集团公司煤峪矿现场实际情况，通过对综采放顶煤工作面火灾隐患治理技术的研究，提出采用“旁路式”注氮法对综放工作面采空区内的火灾隐患进行治理。针对高瓦斯易自燃煤层，罗新荣以某矿综采面为原型，二次开发了带有抽采孔的计算流体力学代码。刘星魁利用Fluent软件模拟了注氮前后采空区自燃带范围的改变情况，并分析了注氮位置和注氮量对采空区氧化带位置分布的影响，从中拟合出了最佳的注氮参数。董军为确定5-2S-2综放工作面采空区连续注氮防灭火的最佳工艺参数，依据计算流体力学基本理论，构建了采空区渗流场数学物理模型，对不同注氮条件下的采空区氧浓度分布及“三带”划分进行数值模拟研究，得到最佳注氮参数。上述研究是防治采空区自燃发火的重要基础，但是，若将这些方法和装置应用于实际工程中还存在问题。首先，由于不同煤矿的采空区物性参数和生产工艺参数都不一样，采用数值模拟确定注氮参数的方法，耗费时间多，技术难度大，难以普及；其次，上述研究都不能明确氮气流动过程中压损与注氮参数之间的量化关系，因此，无法快速确定注氮的具体参数。为了克服上述不足，本发明提出了一种基于最小流动单元的采空区氮气充注参数优化方法。

发明内容

为了解决上述技术问题，本发明提供一种基于最小流动单元的采空区氮气充注压力损失计算方法。

本发明采用的技术方案是：一种基于最小流动单元的采空区氮气充注压力损失计算方法，包括如下步骤：

1)建立采空区最小流动模型，并确定其最小流动单元；

2)以最小流动单元为研究对象，根据断面面积随流动距离变化关系式的不同，将流动过程划分为四个不同的阶段，并确定每个阶段的过流断面积表达式；

3)根据过流断面积表达式和连续性方程，确定氮气的流速，初始速度为u₀，最小流动单元入口处的过流断面积为A₀，则流过氮气的总体积为A₀u₀，下游断面的流速即用氮气体积流量除以过流断面面积获得；

4)根据Navier-Stocks方程，确定氮气在最小流动单元流动的压力损失值，具体操作如下：

a)基于注入氮气时氮气的流动为稳定流，且最小流动单元沿z轴方向所受的力只有压力和粘滞力，于是Navier-Stocks方程变为：

式中：方程左边三项分别为X、Y、Z方向上的惯性力分量，方程右边两项分别为静压力梯度项和粘性力项，u_x、u_y、u_z是流体质点单位时间内在X、Y、Z方向上的位置变化，x、y和z分别为X轴、Y轴和Z轴的长度变量，和分别为X轴、Y轴和Z轴的偏导数，ρ为氮气密度，P为静压强，为静压强在Z轴上的偏导数，υ为氮气运动粘度，为拉普拉斯算子符号，为粘性力散度；

b)Z轴方向为氮气进气断面的法线方向，以Z轴初始坐标轴，依据笛卡尔坐标系准则，分别确定出X轴方向和Y轴方向；由于X轴、Y轴方向以Z轴为对称轴对称，故X、Y方向的流速大小相等；有dx＝u_xdt、dy＝u_ydt、dz＝u_zdt、u_x＝tgθ₁u_z和u_y＝tgθ₂u_z，式中：dx、dy和dz分别为X轴、Y轴和Z轴的长度微分量，dt为时间微分量，θ₁、θ₂分别为速度分量u_x与速度分量u_z之间的夹角、速度分量u_y与速度分量u_z之间的夹角，将其代入欧拉流体力学的中，能得出tgθ₁≈tgθ₂<<1，继而，得到X轴、Y轴方向上的压损远远小于Z轴方向上的压损，故只考虑Z轴方向上的流动过程；

c)利用平均速度的概念，计算出各个阶段的最小流动单元的平均速度，并将其代入Navier-Stocks方程，得到单位长度上的压力损失计算公式，进而求出流动过程中的压力损失；

式中，为u_z在有限容积V_z内的平均值，dV为容积的微分量，A_z为坐标值z位置的过流断面面积。

上述的基于最小流动单元的采空区氮气充注压力损失计算方法中，步骤1)中的最小流动模型包括8个直径为d的球体，1个直径为小径球体和6个直径为的小径半球；8个直径为d的球体分别位于边长为d的立方体八个顶点处，每个直径为d的球体与相邻的直径为d的球体相切；所述的直径为小径球体位于8个直径为d的球体的中心处，与8个直径为d的球体分别相切；6个直径为的小径半球分别位于立方体的六个面的中心位置，直径为的小径半球与周围球面相切。

上述的基于最小流动单元的采空区氮气充注压力损失计算方法中，步骤1)中最小流动单元确定包括如下步骤：以边长为d的立方体的各条棱边的中心点处为切分点，将最小流动模型切分呈八个等分的小单元，每个小单元中由一个直径为d的球体和4个1/8小径球体组成，将切割出来的小单元视为最小流动模型的最小流动单元。

上述的基于最小流动单元的采空区氮气充注压力损失计算方法中，步骤2)具体操作如下：

对于最小流动单元而言，在Z轴方向，坐标原点为最小流动模型中心处小径球体的球心，z＝0点处过流断面面积为其中：d为球体直径，π为圆周率；

当满足约束条件时，z为过流断面在Z轴上的坐标，随着z的增加，过流断面的面积等于边长为d的方形断面减去1个随z变化的小径半圆面积和2个随z变化的大径半圆面积，即：

同理：当时，

当时，

当时，

在一个周期内，过流断面变化，由上述四个阶段构成；当分别满足公式(1)(2)、(3)或(4)的断面关系式时，分别称之为第一流动阶段、第二流动阶段、第三流动阶段或第四流动阶段。

上述的基于最小流动单元的采空区氮气充注压力损失计算方法中，步骤4)中步骤c)所得的个流动阶段的压力损失计算公式如下：

第一流动阶段的最小流动单元在单位长度上的压力损失计算公式为：

J＝Au₀+Bu₀²，其中：

式中：J为压力梯度，A和B为组合变量，u₀为最小流动单元(z＝0)的形心处沿Z轴方向的断面平均流速，α'＝α₁+α₂，α₁、α₂分别为X轴与Z轴的粘性力比值、Y轴与Z轴的粘性力比值，β'＝β₁+β₂，β₁、β₂分别为X轴与Z轴的惯性力比值、Y轴与Z轴的惯性力比值，μ为氮气动力粘度；

第二流动阶段的最小流动单元在单位长度上的压力损失计算公式为：

第三流动阶段的最小流动单元在单位长度上的压力损失计算公式为：

第四流动阶段的最小流动单元在单位长度上的压力损失计算公式为：

与现有技术相比，本发明的有益效果是：

本发明能快速计算出氮气流动过程中的压力损失，为注入氮气时氮气所需的初始压力提供理论支持，能实现对充注范围的准确计算，又能使计算过程变得简单方便。

附图说明

图1是最小流动模型正视图。

图2是最小流动模型左视图。

图3是最小流动模型俯视图。

图4是最小流动单元的结构图。

图5是氮气流经最小流动单元第一流动阶段的压损与温度的关系图。

图6是氮气流经最小流动单元第二流动阶段的压损与温度的关系图。

图7是氮气流经最小流动单元第三流动阶段的压损与温度的关系图。

图8是氮气流经最小流动单元第四流动阶段的压损与温度的关系图。

图9是氮气流经最小流动单元一个完整阶段压损与温度的关系图。

具体实施方式

下面结合附图对本发明作进一步详细的描述。

本发明包括如下步骤：

1)建立最小流动模型，确定最小流动模型的最小流动单元。对于一个边长为2d的立方体，恰好能容纳下8个直径为d的球体，其最松散的排列方式为球体之间及球体与立方体的侧面均相切，则此包含8个球体的立方体孔隙率为1-π/6，约为0.4764；继而，由8个直径为d的大径球体和4个直径为的小径球体组成的最小流动模型孔隙率为0.3737。

因此，建立如图1-3所示的最小流动模型，最小流动模型包括8个直径为d的球体，1个直径为小径球体和6个直径为的小径半球；8个直径为d的球体分别位于边长为d的立方体八个顶点处，每个直径为d的球体与相邻的直径为d的球体相切；所述的直径为小径球体位于8个直径为d的球体的中心处，与8个直径为d的球体分别相切；6个直径为的小径半球分别位于立方体的六个面的中心位置，直径为的小径半球与周围球面相切。

为了简化理论推导过程，在最小流动模型的基础上，切割出最小流动单元。所述最小流动单元，具体如下：

以边长为d的立方体的各条棱边的中心点处为切分点，将最小流动模型切分呈八个等分的小单元，每个小单元中由一个直径为d的球体和4个1/8小径球体组成，将切割出来的小单元视为最小流动模型的最小流动单元，如图4所示。每个最小流动单元的孔隙率为0.3737，与八个大球和四个小球组成的模型孔隙率相等。在最小流动单元上的流体流动，有一个完整的流动周期。

2)Z轴方向为氮气进气断面的法线方向，以Z轴初始坐标轴，依据笛卡尔坐标系准则，分别确定出X轴方向和Y轴方向；以最小流动模型的中心处的小径球体的球心为原点建立坐标系。

对于最小流动单元而言，在Z轴方向，从z＝0点处过流断面面积为：(A₀为初始过流断面面积，d球体直径，π为圆周率)。

当满足约束条件时，z为断面在Z轴方向上的坐标，随着z的增加，过流断面的面积等于边长为d的方形断面减去1个随z变化的小径半圆面积和2个随z变化的大径半圆面积，即：

式中，A_z为坐标值z位置的过流断面面积。

同理，当时，

当时，

当时，

根据过流断面面积随z的关系式不同，划分出不同流动阶段,当分别满足公式(1)、(2)、(3)或(4)的过流断面关系式时，称之为第一流动阶段、第二流动阶段、第三流动阶段或第四流动阶段。在一个周期内，过流断面变化，由上述四个阶段构成；在此基础上，形成了一个周期内最小流动单元的微流动过程。

3)根据Navier-Stocks方程，确定氮气在最小流动单元流动的压力损失值，具体操作如下：

由于不可压缩流体的运动规律符合Navier-Stocks方程，有：

式中：方程左边三项分别为X、Y、Z方向上的惯性力分量，方程右边两项分别为静压力梯度项和粘性力项，u_x、u_y、u_z是流体质点单位时间内在X、Y、Z方向上的位置变化，t为时间，为时间偏导数，为u_z的偏导数，和分别为X轴、Y轴和Z轴的偏导数，ρ为氮气密度，f_z为Z轴的体积力分量，P为静压强，为静压强在Z轴上的偏导数，μ为氮气动力粘度，为拉普拉斯算子符号，为粘性力散度。

基于注入氮气时氮气的流动为稳定流，且最小流动单元沿z轴方向所受的力只有压力和粘滞力，因此f_z＝0；且，根据物性参数之间的物理学基本关系，有μ＝ρυ；则氮气的Navier-Stocks方程，即公式(5)简化为：

式中，υ为氮气运动粘度。

在公式(6)中，方程的左端为惯性力作用项，方程右端第一项为压力损失项、第二项为粘滞力作用项。只要分别将方程中的各项对最小流动单元作体积分就可以得出流体在典型最小流动单元中的运动方程。

建立四个阶段的微流动模型：

a)第一流动阶段的模型建立

在最小流动单元内，当z＝0时，设形心处z轴方向点速度等于过流断面的平均流速，且为u₀；其中，该过流断面的面积为A₀；沿Z轴方向，坐标为z处氮气的流速为u_z，过流断面面积为A_z，根据连续性方程则有：

A₀u₀＝A_zu_z>

当z＝0时，过流断面面积坐标为z处的过流断面面积

那么，坐标为z处的过流断面平均流速为：

因此，沿Z轴方向的当地加速度为：

由式(8)和式(9)，得迁移加速度为：

得出了坐标为z的处的过流断面所受的惯性力项和粘滞力项在Z轴方向的分量。由于实际流动只能发生在孔隙部分，要得到整个立方体单元在Z轴方向所受的合力就要将孔隙部分所受的力平均到整个最小流动单元的过流断面上。设立方体的体积为V₀，孔隙部分的体积为V_z，根据质量守恒、动量守恒可得：

式中：代表氮气在不同流动阶段的速度u_z的平均值，为不同流动阶段同一时刻由于在Z轴方向上位置不同引起的单位长度上速度的变化的平均值，为不同流动阶段上Z轴方向上粘性力作用的平均值。

把式(8)代入式(12)，得：

其中，第一流动阶段的空隙部分体积为则由公式(15)得第一流动阶段的平均断面过流速度为：

把式(10)代入(13)式，得第一流动阶段的平均惯性力为：

把式(11)代入式(14)，得第一流动阶段的平均粘性力为：

对于惯性力和粘滞力在X、Y方向的分量，由对称性，得：

式中：分别为的平均值。

式中，分别为的平均值。

由于氮气从孔腹流至孔喉或者从孔喉流至孔腹的过程中，流线也在不断的收缩或者放大，在球面附近则会发生边界层分离并形成漩涡。因此，渗透流速沿X坐标轴和Y坐标轴方向的速度梯度分布函数难以获取，无法通过理论直接推导出惯性力和粘滞力在X、Y方向的分量。但由对称性可知，Z轴方向的惯性力和粘滞力在X、Y轴方向的分量具有一定的相关性，Irmay就曾用数学估计的方法对此完成了证明。沿Y轴方向Δy处取一截面，设u_x和u_z分别为X和Z方向上的分速度，且：

式中，θ₁为速度分量u_x与速度分量u_z之间的夹角；

进一步，得到：

根据量纲相同的实验流体力学原理，可得：

设：

式中，β₁、β₂分别为X轴与Z轴的惯性力比值、Y轴与Z轴的惯性力比值。

同理，可得：

设：

式中，α₁、α₂分别为X轴与Z轴的粘性力比值、Y轴与Z轴的粘性力比值。

令β'＝β₁+β₂，则过流断面上的平均惯性力可以表示为：

式中：表示不同阶段同一时刻由于在X、Y、Z方向上位置不同引起的单位长度上速度的变化的平均值。

再令α'＝α₁+α₂，则过流断面上的平均粘滞力项可以表示为：

把公式(6)应用于过流断面上，有：

式中，代表粘性力散度，为过流断面上的平均压力梯度，P为过流断面中心点的氮气压力。

移项整理公式(29)，得：

把μ＝ρ·υ，代入公式(30)，有：

令并将公式(27)和(28)，代入公式(31)中，整理可得：

式中：J为单位长度上的压力损失。

此外，据前述推导过程有第一流动阶段的即：第一流动阶段的氮气的平均流速约为孔腹流速u₀的1.8069倍。

令：

再令：

把公式(33)和(34)代入公式(32)，得：

J＝Au₀+Bu₀²>

b)同理，在第二流动阶段中，单位长度上的压力损失模型为：

c)同理，在第三流动阶段中，单位长度上的压力损失模型为：

d)同理，在第四流动阶段中，单位长度上的压力损失模型为：

e)根据欧拉流体力学，有：

根据牛顿经典力学，

根据公式(21)，且类似的定义有：

式中，θ₂为速度分量u_y与速度分量u_z之间的夹角。

将公式(49)和(50)带入公式(48)，得：

由于初速度u₀的方向与Z轴一致，故由氮气流动引起的压强主要作用于Z轴方向并影响Z轴方向上的流动，且tgθ₁≈tgθ₂<<1，故可忽略在X，Y轴方向上的压力损失，只考虑主流方向Z轴上的压力损失，因此步骤e)所得的压力损失即为Z轴上的压力损失。

由前面的推导结果可知，在模型确定的情况下，球直径d，初速度u₀均为常量，动力粘度μ和氮气密度ρ为影响压力损失J的两个变量。以球直径d为5mm，初始速度u₀为1m/s为例，由于采空区温度变化范围集中在20℃到60℃之间，且在不同的温度下对应的动力粘度μ和氮气密度ρ均不同，所以应在不同温度条件下研究J与μ、ρ的关系。分别取20℃、30℃、40℃、50℃、60℃下的动力粘度μ和密度ρ带入四个流动阶段的式子中，算出不同粘度系数下每米的压力损失J值和模型四个流动阶段的压力损失，并将运动粘度υ和压力损失J的关系用图表表示出来，四个不同阶段的单位长度下的压损J与温度T之间的关系图，如图5-8所示。

在图5-8中，水平轴为υ×10^-5，单位为1m²/s；纵轴为单位长度上的压力损失，单位为Pa/m。通过分析具体实施方案，所归纳出的规律如下：①氮气在模型内流动的过程中，压损与初始速度满足：的二次关系式。②随着温度变化，氮气动力粘度和密度也将发生变化，两者共同影响单位长度上的压损。③随着运动粘度的增加，整个模型上的压损将减小。由图9可知，氮气流动过程中，在粘性系数不变的情况下，第一流动阶段的压损最小，第二流动阶段的压损最大。第三第四流动阶段的压损很接近。而四个流动阶段的压损也随着运动粘度系数的增加而增加。本发明与现有方法相比，能快速计算出氮气流动过程中的压力损失，能实现对充注范围的准确计算，又能使计算过程变得简单方便。