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

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

Info

Publication number
CN108733923B
CN108733923B CN201810486161.4A CN201810486161A CN108733923B CN 108733923 B CN108733923 B CN 108733923B CN 201810486161 A CN201810486161 A CN 201810486161A CN 108733923 B CN108733923 B CN 108733923B
Authority
CN
China
Prior art keywords
flow
axis
diameter
nitrogen
pressure loss
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810486161.4A
Other languages
English (en)
Other versions
CN108733923A (zh
Inventor
陈世强
姜文
王海桥
鲁义
王鹏飞
田峰
于琦
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ping An Electric Co Ltd
Original Assignee
Hunan University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hunan University of Science and Technology filed Critical Hunan University of Science and Technology
Priority to CN201810486161.4A priority Critical patent/CN108733923B/zh
Publication of CN108733923A publication Critical patent/CN108733923A/zh
Priority to PCT/CN2019/087503 priority patent/WO2019223628A1/zh
Priority to KR1020207004816A priority patent/KR102307223B1/ko
Application granted granted Critical
Publication of CN108733923B publication Critical patent/CN108733923B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/28Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • EFIXED CONSTRUCTIONS
    • E21EARTH DRILLING; MINING
    • E21FSAFETY DEVICES, TRANSPORT, FILLING-UP, RESCUE, VENTILATION, OR DRAINING IN OR OF MINES OR TUNNELS
    • E21F5/00Means or methods for preventing, binding, depositing, or removing dust; Preventing explosions or fires
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/11Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00Details relating to the application field
    • G06F2113/08Fluids

Abstract

本发明公开了一种基于最小流动单元的采空区氮气充注压力损失计算方法,包括如下步骤:(1)建立采空区最小流动模型,并确定其最小流动单元;(2)以最小流动单元为研究对象,根据断面面积随流动距离变化关系式的不同,将流动过程划分为四个不同的阶段,确定每个阶段的过流断面积表达式;(3)根据Navier‑Stocks方程,确定氮气在最小流动单元流动的压力损失值。本发明能快速计算出氮气流动过程中的压力损失,为注入氮气时氮气所需的初始压力提供理论支持,能实现对充注范围的准确计算,又能使计算过程变得简单方便。

Description

基于最小流动单元的采空区氮气充注压力损失计算方法
技术领域
本发明属于矿井通风和采空区防灾减灾技术领域,具体是涉及一种基于最小流动单元的采空区氮气充注压力损失计算方法。
背景技术
70%以上的煤矿火灾事故发生在毗邻采区的采空区。采空区遗煤自燃发火及其灭火救灾一直是行业内的热点问题,采用注氮治理自燃发火是最常见和行之有效的灭火措施。但是,现有计算方法和实施细则,无法准确计算灭火用所注氮气在采空区内流动的距离,从而无法准确预测所注氮气能影响的范围,导致液氮浪费和灭火困难。
在确定采空区注氮参数方法及其措施的研究进展方面,李东等人提出了一种采空区U型通风工作面全断面帷幕注氮的防灭火方法,提高了氮气惰化效率。朱红青发明了自动控制旋转牵引式注氮防灭火装置,实现了采空区注氮点在空间上的连续,提升了采空区注氮的惰化效果。秦波涛和鲁义提出了一种高效治理浅埋藏煤层大面积采空区遗煤自燃的方法,集堵漏控风与快速惰化降温为一体,高效防治了遗煤自燃。郭君柳探讨了采空区“三带”的观测以及注氮防灭火的实用方法,提高了煤炭开采的安全性和经济效益。根据质量、动量和组分守恒方程,朱红青建立采空区气体体积分数变化的理论流体力学模型,通过初始化、赋边界值和迭代计算,数值模拟确定最佳注氮位置和注氮量。李宗翔结合段王煤矿进行了“两进一回”复杂采场采空区注氮防灭火模拟,确定了最佳的注氮量和注氮位置。何俊忠通过SF6示踪气体测定宏岩煤矿4405采空区的漏风量,并确定主要漏风方向,在此基础上优化了采空区注氮工艺。张琪结合大同煤矿集团公司煤峪矿现场实际情况,通过对综采放顶煤工作面火灾隐患治理技术的研究,提出采用“旁路式”注氮法对综放工作面采空区内的火灾隐患进行治理。针对高瓦斯易自燃煤层,罗新荣以某矿综采面为原型,二次开发了带有抽采孔的计算流体力学代码。刘星魁利用Fluent软件模拟了注氮前后采空区自燃带范围的改变情况,并分析了注氮位置和注氮量对采空区氧化带位置分布的影响,从中拟合出了最佳的注氮参数。董军为确定5-2S-2综放工作面采空区连续注氮防灭火的最佳工艺参数,依据计算流体力学基本理论,构建了采空区渗流场数学物理模型,对不同注氮条件下的采空区氧浓度分布及“三带”划分进行数值模拟研究,得到最佳注氮参数。上述研究是防治采空区自燃发火的重要基础,但是,若将这些方法和装置应用于实际工程中还存在问题。首先,由于不同煤矿的采空区物性参数和生产工艺参数都不一样,采用数值模拟确定注氮参数的方法,耗费时间多,技术难度大,难以普及;其次,上述研究都不能明确氮气流动过程中压损与注氮参数之间的量化关系,因此,无法快速确定注氮的具体参数。为了克服上述不足,本发明提出了一种基于最小流动单元的采空区氮气充注参数优化方法。
发明内容
为了解决上述技术问题,本发明提供一种基于最小流动单元的采空区氮气充注压力损失计算方法。
本发明采用的技术方案是:一种基于最小流动单元的采空区氮气充注压力损失计算方法,包括如下步骤:
1)建立采空区最小流动模型,并确定其最小流动单元;
2)以最小流动单元为研究对象,根据断面面积随流动距离变化关系式的不同,将流动过程划分为四个不同的阶段,并确定每个阶段的过流断面积表达式;
3)根据过流断面积表达式和连续性方程,确定氮气的流速,初始速度为u0,最小流动单元入口处的过流断面积为A0,则流过氮气的总体积为A0u0,下游断面的流速即用氮气体积流量除以过流断面面积获得;
4)根据Navier-Stocks方程,确定氮气在最小流动单元流动的压力损失值,具体操作如下:
a)基于注入氮气时氮气的流动为稳定流,且最小流动单元沿z轴方向所受的力只有压力和粘滞力,于是Navier-Stocks方程变为:
Figure GDA0002353356280000031
式中:方程左边三项分别为X、Y、Z方向上的惯性力分量,方程右边两项分别为静压力梯度项和粘性力项,ux、uy、uz是流体质点单位时间内在X、Y、Z方向上的位置变化,x、y和z分别为X轴、Y轴和Z轴的长度变量,
Figure GDA0002353356280000033
Figure GDA0002353356280000034
Figure GDA0002353356280000035
分别为X轴、Y轴和Z轴的偏导数,ρ为氮气密度,P为静压强,
Figure GDA0002353356280000036
为静压强在Z轴上的偏导数,υ为氮气运动粘度,
Figure GDA0002353356280000037
为拉普拉斯算子符号,
Figure GDA0002353356280000038
为粘性力散度;
b)Z轴方向为氮气进气断面的法线方向,以Z轴初始坐标轴,依据笛卡尔坐标系准则,分别确定出X轴方向和Y轴方向;由于X轴、Y轴方向以Z轴为对称轴对称,故X、Y方向的流速大小相等;有dx=uxdt、dy=uydt、dz=uzdt、ux=tgθ1uz和uy=tgθ2uz,式中:dx、dy和dz分别为X轴、Y轴和Z轴的长度微分量,dt为时间微分量,θ1、θ2分别为速度分量ux与速度分量uz之间的夹角、速度分量uy与速度分量uz之间的夹角,将其代入欧拉流体力学的
Figure GDA0002353356280000032
中,能得出tgθ1≈tgθ2<<1,继而,得到X轴、Y轴方向上的压损远远小于Z轴方向上的压损,故只考虑Z轴方向上的流动过程;
c)利用平均速度的概念,
Figure GDA0002353356280000041
计算出各个阶段的最小流动单元的平均速度,并将其代入Navier-Stocks方程,得到单位长度上的压力损失计算公式,进而求出流动过程中的压力损失;
式中,
Figure GDA0002353356280000042
为uz在有限容积Vz内的平均值,dV为容积的微分量,Az为坐标值z位置的过流断面面积。
上述的基于最小流动单元的采空区氮气充注压力损失计算方法中,步骤1)中的最小流动模型包括8个直径为d的球体,1个直径为
Figure GDA0002353356280000043
小径球体和6个直径为
Figure GDA0002353356280000044
的小径半球;8个直径为d的球体分别位于边长为d的立方体八个顶点处,每个直径为d的球体与相邻的直径为d的球体相切;所述的直径为
Figure GDA0002353356280000045
小径球体位于8个直径为d的球体的中心处,与8个直径为d的球体分别相切;6个直径为
Figure GDA0002353356280000046
的小径半球分别位于立方体的六个面的中心位置,直径为
Figure GDA0002353356280000047
的小径半球与周围球面相切。
上述的基于最小流动单元的采空区氮气充注压力损失计算方法中,步骤1)中最小流动单元确定包括如下步骤:以边长为d的立方体的各条棱边的中心点处为切分点,将最小流动模型切分呈八个等分的小单元,每个小单元中由一个直径为d的球体和4个1/8小径球体组成,将切割出来的小单元视为最小流动模型的最小流动单元。
上述的基于最小流动单元的采空区氮气充注压力损失计算方法中,步骤2)具体操作如下:
对于最小流动单元而言,在Z轴方向,坐标原点为最小流动模型中心处小径球体的球心,z=0点处过流断面面积为
Figure GDA0002353356280000048
其中:d为球体直径,π为圆周率;
当满足约束条件
Figure GDA0002353356280000051
时,z为过流断面在Z轴上的坐标,随着z的增加,过流断面的面积等于边长为d的方形断面减去1个随z变化的小径半圆面积和2个随z变化的大径半圆面积,即:
Figure GDA0002353356280000052
同理:当
Figure GDA0002353356280000053
时,
Figure GDA0002353356280000054
Figure GDA0002353356280000055
时,
Figure GDA0002353356280000056
Figure GDA0002353356280000057
时,
Figure GDA0002353356280000058
在一个周期内,过流断面变化,由上述四个阶段构成;当分别满足公式(1)(2)、(3)或(4)的断面关系式时,分别称之为第一流动阶段、第二流动阶段、第三流动阶段或第四流动阶段。
上述的基于最小流动单元的采空区氮气充注压力损失计算方法中,步骤4)中步骤c)所得的个流动阶段的压力损失计算公式如下:
第一流动阶段的最小流动单元在单位长度上的压力损失计算公式为:
J=Au0+Bu0 2,其中:
Figure GDA0002353356280000059
式中:J为压力梯度,A和B为组合变量,u0为最小流动单元(z=0)的形心处沿Z轴方向的断面平均流速,α'=α12,α1、α2分别为X轴与Z轴的粘性力比值、Y轴与Z轴的粘性力比值,β'=β12,β1、β2分别为X轴与Z轴的惯性力比值、Y轴与Z轴的惯性力比值,μ为氮气动力粘度;
第二流动阶段的最小流动单元在单位长度上的压力损失计算公式为:
Figure GDA0002353356280000061
第三流动阶段的最小流动单元在单位长度上的压力损失计算公式为:
Figure GDA0002353356280000062
第四流动阶段的最小流动单元在单位长度上的压力损失计算公式为:
Figure GDA0002353356280000063
与现有技术相比,本发明的有益效果是:
本发明能快速计算出氮气流动过程中的压力损失,为注入氮气时氮气所需的初始压力提供理论支持,能实现对充注范围的准确计算,又能使计算过程变得简单方便。
附图说明
图1是最小流动模型正视图。
图2是最小流动模型左视图。
图3是最小流动模型俯视图。
图4是最小流动单元的结构图。
图5是氮气流经最小流动单元第一流动阶段的压损与温度的关系图。
图6是氮气流经最小流动单元第二流动阶段的压损与温度的关系图。
图7是氮气流经最小流动单元第三流动阶段的压损与温度的关系图。
图8是氮气流经最小流动单元第四流动阶段的压损与温度的关系图。
图9是氮气流经最小流动单元一个完整阶段压损与温度的关系图。
具体实施方式
下面结合附图对本发明作进一步详细的描述。
本发明包括如下步骤:
1)建立最小流动模型,确定最小流动模型的最小流动单元。对于一个边长为2d的立方体,恰好能容纳下8个直径为d的球体,其最松散的排列方式为球体之间及球体与立方体的侧面均相切,则此包含8个球体的立方体孔隙率为1-π/6,约为0.4764;继而,由8个直径为d的大径球体和4个直径为
Figure GDA0002353356280000071
的小径球体组成的最小流动模型孔隙率为0.3737。
因此,建立如图1-3所示的最小流动模型,最小流动模型包括8个直径为d的球体,1个直径为
Figure GDA0002353356280000072
小径球体和6个直径为
Figure GDA0002353356280000073
的小径半球;8个直径为d的球体分别位于边长为d的立方体八个顶点处,每个直径为d的球体与相邻的直径为d的球体相切;所述的直径为
Figure GDA0002353356280000074
小径球体位于8个直径为d的球体的中心处,与8个直径为d的球体分别相切;6个直径为
Figure GDA0002353356280000075
的小径半球分别位于立方体的六个面的中心位置,直径为
Figure GDA0002353356280000076
的小径半球与周围球面相切。
为了简化理论推导过程,在最小流动模型的基础上,切割出最小流动单元。所述最小流动单元,具体如下:
以边长为d的立方体的各条棱边的中心点处为切分点,将最小流动模型切分呈八个等分的小单元,每个小单元中由一个直径为d的球体和4个1/8小径球体组成,将切割出来的小单元视为最小流动模型的最小流动单元,如图4所示。每个最小流动单元的孔隙率为0.3737,与八个大球和四个小球组成的模型孔隙率相等。在最小流动单元上的流体流动,有一个完整的流动周期。
2)Z轴方向为氮气进气断面的法线方向,以Z轴初始坐标轴,依据笛卡尔坐标系准则,分别确定出X轴方向和Y轴方向;以最小流动模型的中心处的小径球体的球心为原点建立坐标系。
对于最小流动单元而言,在Z轴方向,从z=0点处过流断面面积为:
Figure GDA0002353356280000081
(A0为初始过流断面面积,d球体直径,π为圆周率)。
当满足约束条件
Figure GDA0002353356280000082
时,z为断面在Z轴方向上的坐标,随着z的增加,过流断面的面积等于边长为d的方形断面减去1个随z变化的小径半圆面积和2个随z变化的大径半圆面积,即:
Figure GDA0002353356280000083
式中,Az为坐标值z位置的过流断面面积。
同理,当
Figure GDA0002353356280000084
时,
Figure GDA0002353356280000085
Figure GDA0002353356280000086
时,
Figure GDA0002353356280000087
Figure GDA0002353356280000088
时,
Figure GDA0002353356280000089
根据过流断面面积随z的关系式不同,划分出不同流动阶段,当分别满足公式(1)、(2)、(3)或(4)的过流断面关系式时,称之为第一流动阶段、第二流动阶段、第三流动阶段或第四流动阶段。在一个周期内,过流断面变化,由上述四个阶段构成;在此基础上,形成了一个周期内最小流动单元的微流动过程。
3)根据Navier-Stocks方程,确定氮气在最小流动单元流动的压力损失值,具体操作如下:
由于不可压缩流体的运动规律符合Navier-Stocks方程,有:
Figure GDA0002353356280000091
式中:方程左边三项分别为X、Y、Z方向上的惯性力分量,方程右边两项分别为静压力梯度项和粘性力项,ux、uy、uz是流体质点单位时间内在X、Y、Z方向上的位置变化,t为时间,
Figure GDA0002353356280000094
为时间偏导数,
Figure GDA0002353356280000095
为uz的偏导数,
Figure GDA0002353356280000096
Figure GDA0002353356280000097
Figure GDA0002353356280000098
分别为X轴、Y轴和Z轴的偏导数,ρ为氮气密度,fz为Z轴的体积力分量,P为静压强,
Figure GDA0002353356280000099
为静压强在Z轴上的偏导数,μ为氮气动力粘度,
Figure GDA00023533562800000910
为拉普拉斯算子符号,
Figure GDA00023533562800000911
为粘性力散度。
基于注入氮气时氮气的流动为稳定流,且最小流动单元沿z轴方向所受的力只有压力和粘滞力,因此
Figure GDA0002353356280000092
fz=0;且,根据物性参数之间的物理学基本关系,有μ=ρυ;则氮气的Navier-Stocks方程,即公式(5)简化为:
Figure GDA0002353356280000093
式中,υ为氮气运动粘度。
在公式(6)中,方程的左端为惯性力作用项,方程右端第一项为压力损失项、第二项为粘滞力作用项。只要分别将方程中的各项对最小流动单元作体积分就可以得出流体在典型最小流动单元中的运动方程。
建立四个阶段的微流动模型:
a)第一流动阶段的模型建立
在最小流动单元内,当z=0时,设形心处z轴方向点速度等于过流断面的平均流速,且为u0;其中,该过流断面的面积为A0;沿Z轴方向,坐标为z处氮气的流速为uz,过流断面面积为Az,根据连续性方程则有:
A0u0=Azuz (7);
当z=0时,过流断面面积
Figure GDA0002353356280000101
坐标为z处的过流断面面积为
Figure GDA0002353356280000102
那么,坐标为z处的过流断面平均流速为:
Figure GDA0002353356280000103
因此,沿Z轴方向的当地加速度为:
Figure GDA0002353356280000104
由式(8)和式(9),得迁移加速度为:
Figure GDA0002353356280000105
Figure GDA0002353356280000106
得出了坐标为z的处的过流断面所受的惯性力项和粘滞力项在Z轴方向的分量。由于实际流动只能发生在孔隙部分,要得到整个立方体单元在Z轴方向所受的合力就要将孔隙部分所受的力平均到整个最小流动单元的过流断面上。设立方体的体积为V0,孔隙部分的体积为Vz,根据质量守恒、动量守恒可得:
Figure GDA0002353356280000111
Figure GDA0002353356280000112
Figure GDA0002353356280000113
式中:
Figure GDA0002353356280000114
代表氮气在不同流动阶段的速度uz的平均值,
Figure GDA0002353356280000115
为不同流动阶段同一时刻由于在Z轴方向上位置不同引起的单位长度上速度的变化
Figure GDA0002353356280000116
的平均值,
Figure GDA0002353356280000117
为不同流动阶段上Z轴方向上粘性力作用
Figure GDA0002353356280000118
的平均值。
把式(8)代入式(12),得:
Figure GDA0002353356280000119
其中,第一流动阶段的空隙部分体积为
Figure GDA00023533562800001110
则由公式(15)得第一流动阶段的平均断面过流速度为:
Figure GDA00023533562800001111
把式(10)代入(13)式,得第一流动阶段的平均惯性力为:
Figure GDA00023533562800001112
把式(11)代入式(14),得第一流动阶段的平均粘性力为:
Figure GDA00023533562800001113
对于惯性力和粘滞力在X、Y方向的分量,由对称性,得:
Figure GDA0002353356280000121
式中:
Figure GDA0002353356280000122
分别为
Figure GDA0002353356280000123
的平均值。
Figure GDA0002353356280000124
式中,
Figure GDA0002353356280000125
分别为
Figure GDA0002353356280000126
的平均值。
由于氮气从孔腹流至孔喉或者从孔喉流至孔腹的过程中,流线也在不断的收缩或者放大,在球面附近则会发生边界层分离并形成漩涡。因此,渗透流速沿X坐标轴和Y坐标轴方向的速度梯度分布函数难以获取,无法通过理论直接推导出惯性力和粘滞力在X、Y方向的分量。但由对称性可知,Z轴方向的惯性力和粘滞力在X、Y轴方向的分量具有一定的相关性,Irmay就曾用数学估计的方法对此完成了证明。沿Y轴方向Δy处取一截面,设ux和uz分别为X和Z方向上的分速度,且:
Figure GDA0002353356280000127
式中,θ1为速度分量ux与速度分量uz之间的夹角;
进一步,得到:
Figure GDA0002353356280000128
根据量纲相同的实验流体力学原理,可得:
Figure GDA0002353356280000129
设:
Figure GDA0002353356280000131
式中,β1、β2分别为X轴与Z轴的惯性力比值、Y轴与Z轴的惯性力比值。
同理,可得:
Figure GDA0002353356280000132
设:
Figure GDA0002353356280000133
式中,α1、α2分别为X轴与Z轴的粘性力比值、Y轴与Z轴的粘性力比值。
令β'=β12,则过流断面上的平均惯性力可以表示为:
Figure GDA0002353356280000134
式中:
Figure GDA0002353356280000135
表示不同阶段同一时刻由于在X、Y、Z方向上位置不同引起的单位长度上速度的变化
Figure GDA0002353356280000136
的平均值。
再令α'=α12,则过流断面上的平均粘滞力项可以表示为:
Figure GDA0002353356280000137
把公式(6)应用于过流断面上,有:
Figure GDA0002353356280000138
式中,
Figure GDA0002353356280000141
代表粘性力散度,
Figure GDA0002353356280000142
为过流断面上的平均压力梯度,P为过流断面中心点的氮气压力。
移项整理公式(29),得:
Figure GDA0002353356280000143
把μ=ρ·υ,代入公式(30),有:
Figure GDA0002353356280000144
Figure GDA0002353356280000145
并将公式(27)和(28),代入公式(31)中,整理可得:
Figure GDA0002353356280000146
式中:J为单位长度上的压力损失。
此外,据前述推导过程有第一流动阶段的
Figure GDA0002353356280000147
即:第一流动阶段的氮气的平均流速
Figure GDA0002353356280000148
约为孔腹流速u0的1.8069倍。
令:
Figure GDA0002353356280000149
再令:
Figure GDA00023533562800001410
把公式(33)和(34)代入公式(32),得:
J=Au0+Bu0 2 (35);
b)同理,在第二流动阶段中,单位长度上的压力损失模型为:
Figure GDA00023533562800001411
Figure GDA00023533562800001412
Figure GDA0002353356280000151
Figure GDA0002353356280000152
c)同理,在第三流动阶段中,单位长度上的压力损失模型为:
Figure GDA0002353356280000153
Figure GDA0002353356280000154
Figure GDA0002353356280000155
Figure GDA0002353356280000156
d)同理,在第四流动阶段中,单位长度上的压力损失模型为:
Figure GDA0002353356280000157
Figure GDA0002353356280000158
Figure GDA0002353356280000159
Figure GDA00023533562800001510
e)根据欧拉流体力学,有:
Figure GDA00023533562800001511
根据牛顿经典力学,
Figure GDA00023533562800001513
根据公式(21),且类似的定义有:
Figure GDA00023533562800001514
式中,θ2为速度分量uy与速度分量uz之间的夹角。
将公式(49)和(50)带入公式(48),得:
Figure GDA00023533562800001512
由于初速度u0的方向与Z轴一致,故由氮气流动引起的压强主要作用于Z轴方向并影响Z轴方向上的流动,且tgθ1≈tgθ2<<1,故可忽略在X,Y轴方向上的压力损失,只考虑主流方向Z轴上的压力损失,因此步骤e)所得的压力损失即为Z轴上的压力损失。
由前面的推导结果可知,在模型确定的情况下,球直径d,初速度u0均为常量,动力粘度μ和氮气密度ρ为影响压力损失J的两个变量。以球直径d为5mm,初始速度u0为1m/s为例,由于采空区温度变化范围集中在20℃到60℃之间,且在不同的温度下对应的动力粘度μ和氮气密度ρ均不同,所以应在不同温度条件下研究J与μ、ρ的关系。分别取20℃、30℃、40℃、50℃、60℃下的动力粘度μ和密度ρ带入四个流动阶段的式子中,算出不同粘度系数下每米的压力损失J值和模型四个流动阶段的压力损失,并将运动粘度υ和压力损失J的关系用图表表示出来,四个不同阶段的单位长度下的压损J与温度T之间的关系图,如图5-8所示。
在图5-8中,水平轴为υ×10-5,单位为1m2/s;纵轴为单位长度上的压力损失,单位为Pa/m。通过分析具体实施方案,所归纳出的规律如下:①氮气在模型内流动的过程中,压损与初始速度满足:
Figure GDA0002353356280000161
的二次关系式。②随着温度变化,氮气动力粘度和密度也将发生变化,两者共同影响单位长度上的压损。③随着运动粘度的增加,整个模型上的压损将减小。由图9可知,氮气流动过程中,在粘性系数不变的情况下,第一流动阶段的压损最小,第二流动阶段的压损最大。第三第四流动阶段的压损很接近。而四个流动阶段的压损也随着运动粘度系数的增加而增加。本发明与现有方法相比,能快速计算出氮气流动过程中的压力损失,能实现对充注范围的准确计算,又能使计算过程变得简单方便。

Claims (5)

1.一种基于最小流动单元的采空区氮气充注压力损失计算方法,包括如下步骤:1)建立采空区最小流动模型,并确定其最小流动单元;
2)以最小流动单元为研究对象,根据断面面积随流动距离变化关系式的不同,将流动过程划分为四个不同的阶段,并确定每个阶段的过流断面积表达式;
3)根据过流断面积表达式和连续性方程,确定氮气的流速,初始速度为u0,最小流动单元入口处的过流断面积为A0,则流过氮气的总体积为A0u0,下游断面的流速即用氮气体积流量除以过流断面面积获得;
4)根据Navier-Stocks方程,确定氮气在最小流动单元流动的压力损失值,具体操作如下:
a)基于注入氮气时氮气的流动为稳定流,且最小流动单元沿z轴方向所受的力只有压力和粘滞力,于是Navier-Stocks方程变为:
Figure FDA0002353356270000011
式中:方程左边三项分别为X、Y、Z方向上的惯性力分量,方程右边两项分别为静压力梯度项和粘性力项,ux、uy、uz是流体质点单位时间内在X、Y、Z方向上的位置变化,x、y和z分别为X轴、Y轴和Z轴的长度变量,
Figure FDA0002353356270000012
Figure FDA0002353356270000013
Figure FDA0002353356270000014
分别为X轴、Y轴和Z轴的偏导数,ρ为氮气密度,P为静压强,
Figure FDA0002353356270000015
为静压强在Z轴上的偏导数,υ为氮气运动粘度,
Figure FDA0002353356270000016
为拉普拉斯算子符号,
Figure FDA0002353356270000017
为粘性力散度;
b)Z轴方向为氮气进气断面的法线方向,以Z轴初始坐标轴,依据笛卡尔坐标系准则,分别确定出X轴方向和Y轴方向;由于X轴、Y轴方向以Z轴为对称轴对称,故X、Y方向的流速大小相等;有dx=uxdt、dy=uydt、dz=uzdt、ux=tgθ1uz和uy=tgθ2uz式中:dx、dy和dz分别为X轴、Y轴和Z轴的长度微分量,dt为时间微分量,θ1、θ2分别为速度分量ux与速度分量uz之间的夹角、速度分量uy与速度分量uz之间的夹角,将其代入欧拉流体力学的
Figure FDA0002353356270000021
中,能得出tgθ1≈tgθ2<<1,继而得到X轴、Y轴方向上的压损远远小于Z轴方向上的压损,故只考虑Z轴方向上的流动过程;
c)利用平均速度的概念,
Figure FDA0002353356270000022
计算出各个阶段的最小流动单元的平均速度,并将其代入Navier-Stocks方程,得到单位长度上的压力损失计算公式,进而求出流动过程中的压力损失;
式中,
Figure FDA0002353356270000023
为uz在有限容积Vz内的平均值,dV为容积的微分量,Az为坐标值z位置的过流断面面积。
2.根据权利要求1所述的基于最小流动单元的采空区氮气充注压力损失计算方法,步骤1)中的最小流动模型包括8个直径为d的球体,1个直径为
Figure FDA0002353356270000024
小径球体和6个直径为
Figure FDA0002353356270000025
的小径半球;8个直径为d的球体分别位于边长为d的立方体八个顶点处,每个直径为d的球体与相邻的直径为d的球体相切;所述的直径为
Figure FDA0002353356270000026
小径球体位于8个直径为d的球体的中心处,与8个直径为d的球体分别相切;6个直径为
Figure FDA0002353356270000027
的小径半球分别位于立方体的六个面的中心位置,直径为
Figure FDA0002353356270000028
的小径半球与周围球面相切。
3.根据权利要求2所述的基于最小流动单元的采空区氮气充注压力损失计算方法,步骤1)中最小流动单元确定包括如下步骤:以边长为d的立方体的各条棱边的中心点处为切分点,将最小流动模型切分呈八个等分的小单元,每个小单元中由一个直径为d的球体和4个1/8小径球体组成,将切割出来的小单元视为最小流动模型的最小流动单元。
4.根据权利要求3所述的基于最小流动单元的采空区氮气充注压力损失计算方法,步骤2)中具体操作如下:
对于最小流动单元而言,在Z轴方向,坐标原点为最小流动模型中心处小径球体的球心,z=0点处过流断面面积为
Figure FDA0002353356270000031
其中:d为球体直径,π为圆周率;
当满足约束条件
Figure FDA0002353356270000032
时,z为过流断面在Z轴上的坐标,随着z的增加,过流断面的面积等于边长为d的方形断面减去1个随z变化的小径半圆面积和2个随z变化的大径半圆面积,即:
Figure FDA0002353356270000033
同理:当
Figure FDA0002353356270000034
时,
Figure FDA0002353356270000035
Figure FDA0002353356270000036
时,
Figure FDA0002353356270000037
Figure FDA0002353356270000038
时,
Figure FDA0002353356270000039
在一个周期内,过流断面变化,由上述四个阶段构成;当分别满足公式(1)(2)、(3)或(4)的断面关系式时,分别称之为第一流动阶段、第二流动阶段、第三流动阶段或第四流动阶段。
5.根据权利要求4所述的基于最小流动单元的采空区氮气充注压力损失计算方法,步骤4)中步骤c)所得的个流动阶段的压力损失计算公式如下:
第一流动阶段的最小流动单元在单位长度上的压力损失计算公式为:J=Au0+Bu0 2,其中:
Figure FDA0002353356270000041
J为单位长度上的压力损失,A和B为组合变量,u0为最小流动单元(z=0)的形心处沿Z轴方向的断面平均流速,α'=α12,α1、α2分别为X轴与Z轴的粘性力比值、Y轴与Z轴的粘性力比值,β'=β12,β1、β2分别为X轴与Z轴的惯性力比值、Y轴与Z轴的惯性力比值,μ为氮气动力粘度;
第二流动阶段的最小流动单元在单位长度上的压力损失计算公式为:
Figure FDA0002353356270000042
第三流动阶段的最小流动单元在单位长度上的压力损失计算公式为:
Figure FDA0002353356270000043
第四流动阶段的最小流动单元在单位长度上的压力损失计算公式为:
Figure FDA0002353356270000044
CN201810486161.4A 2018-05-21 2018-05-21 基于最小流动单元的采空区氮气充注压力损失计算方法 Active CN108733923B (zh)

Priority Applications (3)

Application Number Priority Date Filing Date Title
CN201810486161.4A CN108733923B (zh) 2018-05-21 2018-05-21 基于最小流动单元的采空区氮气充注压力损失计算方法
PCT/CN2019/087503 WO2019223628A1 (zh) 2018-05-21 2019-05-18 基于最小流动单元的采空区氮气充注压力损失计算方法
KR1020207004816A KR102307223B1 (ko) 2018-05-21 2019-05-18 최소 유동 유닛에 기초한 채굴된 공동 구역의 질소 가스 충전 압력 손실 계산 방법

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810486161.4A CN108733923B (zh) 2018-05-21 2018-05-21 基于最小流动单元的采空区氮气充注压力损失计算方法

Publications (2)

Publication Number Publication Date
CN108733923A CN108733923A (zh) 2018-11-02
CN108733923B true CN108733923B (zh) 2020-04-14

Family

ID=63938553

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810486161.4A Active CN108733923B (zh) 2018-05-21 2018-05-21 基于最小流动单元的采空区氮气充注压力损失计算方法

Country Status (3)

Country Link
KR (1) KR102307223B1 (zh)
CN (1) CN108733923B (zh)
WO (1) WO2019223628A1 (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108733923B (zh) * 2018-05-21 2020-04-14 湖南科技大学 基于最小流动单元的采空区氮气充注压力损失计算方法
CN109505642B (zh) * 2018-08-21 2020-01-17 湖南科技大学 特长公路隧道开式可控循环通风的节能量计算方法
CN112214943B (zh) * 2020-10-23 2022-04-22 湖南科技大学 一种拟三维煤矿采空区氮气流动压力损失计算方法
CN113673180B (zh) * 2021-08-10 2022-04-12 北京科技大学 膏体充填长距离水平输送管道泄漏位置的快速计算方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102828767A (zh) * 2012-07-19 2012-12-19 大同煤矿集团有限责任公司 采空区自然发火防治方法
CN106499431A (zh) * 2016-11-04 2017-03-15 大同煤矿集团有限责任公司 沿空巷道掘进期间相邻采空区注氮防火方法
CN106761893A (zh) * 2017-03-04 2017-05-31 辽宁工程技术大学 一种煤矿巷道防灭火系统

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
AU2012201821B2 (en) * 2011-05-20 2015-01-29 Sandvik Intellectual Property Ab Fire suppression valve improvements
JP5498523B2 (ja) * 2012-03-07 2014-05-21 住友ゴム工業株式会社 可塑性材料の押出シミュレーション方法及び装置
KR101472706B1 (ko) 2014-03-07 2014-12-15 서울대학교산학협력단 관내 다상 환상 유동 해석 방법 및 장치
CN108733923B (zh) * 2018-05-21 2020-04-14 湖南科技大学 基于最小流动单元的采空区氮气充注压力损失计算方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102828767A (zh) * 2012-07-19 2012-12-19 大同煤矿集团有限责任公司 采空区自然发火防治方法
CN106499431A (zh) * 2016-11-04 2017-03-15 大同煤矿集团有限责任公司 沿空巷道掘进期间相邻采空区注氮防火方法
CN106761893A (zh) * 2017-03-04 2017-05-31 辽宁工程技术大学 一种煤矿巷道防灭火系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
不同煤自燃特性参数下采空区"三带"分布规律的研究;朱红青 等;《矿业研究与开发》;20140630;第34卷(第03期);54-57 *
注氮控氧防火系统综述;姜文源 等;《给水排水》;20060110(第01期);63-66 *
采空区"三带"观测及注氮防灭火研究;郭柳君;《中国新技术新产品》;20150810(第15期);155 *

Also Published As

Publication number Publication date
KR102307223B1 (ko) 2021-09-29
KR20200032155A (ko) 2020-03-25
CN108733923A (zh) 2018-11-02
WO2019223628A1 (zh) 2019-11-28

Similar Documents

Publication Publication Date Title
CN108733923B (zh) 基于最小流动单元的采空区氮气充注压力损失计算方法
Brodny et al. Analysis of methane hazard conditions in mine headings
CN106528994B (zh) 一种基于气液交界面耦合的调压室通气洞风速模拟方法
Jiang et al. Ventilation control of tunnel drilling dust based on numerical simulation
Shi et al. Simulation of the wake vortex and trajectory characteristics of successively launched multiple projectiles
Wang et al. Modeling erosion process in elbows of petroleum pipelines using large eddy simulation
Fernandez-Feria Dam-break flow for arbitrary slopes of the bottom
Wang et al. Study of jet flow and dust motion in flat chambers based on theory of gas-soild two phase flow
Liu et al. Theoretical analysis and numerical simulation of mechanical energy loss and wall resistance of steady open channel flow
Wang et al. Large Eddy Simulation of Slurry Erosion in Submerged Impinging Jets
Yang et al. Simulation and numerical calculation on pipeline leakage process
Rassaei et al. Numerical flow model stepped spillways in order to maximize energy dissipation using FLUENT software
Erwanto et al. The effect of river flow velocity distribution on indications of the occurrence of degradation of the Tambong River basin
Liu et al. 3D simulation for dynamic characteristics of airflow and dust control in a laneway of coal mine
Uchida et al. Quasi 3D numerical simulation for flow and bed variation with various sand waves
Tinney et al. Free streamline theory for segmental jet deflectors
Li et al. Numerical Investigation on Gas Accumulation and Gas Migration in Enlarged Boreholes of Horizontal Gas Wells
Xin et al. Numerical analysis of the influence range of a tunnel ventilation resistance grid
Zhang et al. Study on influencing factors of tunnel ventilation
Chaochao et al. Hydraulic model and numerical simulation for bend stilling basin
Li et al. A case study for escape routes optimization after water inrush in a backward excavated karst tunnel
Pakgar et al. Numerical simulation of flow on a siphon spillway and investigation of the effect of a bottom/outlet angle on hydraulic parameters
Poza Sánchez Understanding flows over channel floor with cross-sectional joints
MOHAMMED et al. Side weir flow investigation in a circular channel using computational fluid dynamics (CFD) for subcritical flow condition
Dinçer et al. Effect of downstream channel slope on numerical modelling of dam break induced flows

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20220415

Address after: 411100, 12 Ping'an Road, Yuhu Bay, Yuhu District, Hunan, Xiangtan

Patentee after: PING'AN ELECTRIC Co.,Ltd.

Address before: 411201 Taoyuan Road, Yuhu District, Hunan, Xiangtan

Patentee before: HUNAN University OF SCIENCE AND TECHNOLOGY

PE01 Entry into force of the registration of the contract for pledge of patent right
PE01 Entry into force of the registration of the contract for pledge of patent right

Denomination of invention: Calculation method of nitrogen filling pressure loss in goaf based on minimum flow unit

Effective date of registration: 20221130

Granted publication date: 20200414

Pledgee: Hunan Tancheng Financing Guarantee Group Co.,Ltd.

Pledgor: PING'AN ELECTRIC Co.,Ltd.

Registration number: Y2022980024116