一种含高浓度瓦斯井巷瞬态风流稳定性计算方法.pdf

本发明专利公开了一种含高浓度瓦斯井巷瞬态风流稳定性计算方法，其计算模型由瓦斯弥散方程和井巷通风网络方程构成，瓦斯弥散方程采用Crank‑Nicholson方法进行有限差分，利用追赶法得到井巷分瓦斯浓度；井巷通风网络方程采用龙格‑库塔法求解井巷分支的风量(风速)；在此基础上对两方程进行耦合求解；本发明计算模型考虑了高浓度瓦斯对井巷风流稳定性的影响，实现了含高浓度瓦斯井巷分支瞬态风流的定量计算。

1.一种含高浓度瓦斯井巷瞬态风流稳定性计算方法，其计算模型由瓦斯弥散方程
$\frac{\∂c_i}{\∂t}+u_i\frac{\∂c_i}{\∂x_i}=E_x\frac{{\∂}^2{c}_i}{\∂x_i^2}---\left(1\right)$
和井巷通风网络方程
$\underset{i=1}{\overset{m}{\Σ}}\[c_{ji}{K}_i\left(\underset{k=1}{\overset{b}{\Σ}}{c}_{ik}\frac{{\mathrm{dq}}_k}{dt}\right)\]=\underset{i=1}{\overset{m}{\Σ}}{c}_{ji}\left(h_{F_i}+h_{M_i}-R_i\left|\underset{k=1}{\overset{b}{\Σ}}{c}_{ik}{q}_k\right|\underset{k=1}{\overset{b}{\Σ}}{c}_{ik}{q}_k\right)---\left(2\right)$
构成，其中式(1)中u_i为井巷分支i的风速，x_i为沿井巷分支i的空间距离，t为时间，c_i为
井巷分支i的瓦斯体积分数，E_x为瓦斯弥散系数，式(2)中为分支风量，R_i为分支
i风阻，为分支i的通风机风压，为分支i的瓦斯风压，c_ji为回路矩阵，c_ik为回路矩阵c_ji
的转置矩阵，其中瓦斯弥散方程(1)求解方法是首先采用Crank-Nicholson方法进行有限差
分，然后利用追赶法得到井巷分瓦斯浓度，同时井巷通风网络方程(2)求解方法是采用龙
格-库塔法求解井巷分支的风量(风速)，求解过程中需对方程(1)和方程(2)进行耦合求解。
2.根据权利要求1所述的一种含高浓度瓦斯井巷瞬态风流稳定性计算方法，其特征在
于：所述的瓦斯弥散系数其中r为圆管半径，α为摩擦阻力系数，u为井巷分支
的风速。
3.根据权利要求1所述的一种含高浓度瓦斯井巷瞬态风流稳定性计算方法，其特征在
于：所述的分支风量q_k(k＝1,2,…,b)为待求参数，也是时间的函数，即所求的值是b个未知
函数，b为通风网络中独立回路的个数；q_k是评判井巷风流稳定性的指标，q_k＝0表示风流停
滞，q_k＜0表示风流逆转，q_k＞0表示风流方向保持不变。
4.根据权利要求1所述的一种含高浓度瓦斯井巷瞬态风流稳定性计算方法，其特征在
于：所述的通风机风压a₀，a₁，…，a_n为曲线拟合系数；给定三个或三个以
上风机的风压、风量数据，通过最小二乘法求解拟合系数。
5.根据权利要求1所述的一种含高浓度瓦斯井巷瞬态风流稳定性计算方法，其特征在
于：所述的瓦斯风压为分支i的平均密度，为t_n时刻分
支第j个空间步长的气体密度，为t_n时刻分支第j个空
间步长的瓦斯体积分数，ρ_m和ρ_a分别为甲烷和空气的密度，h为空间步长，L_i为井巷分支i的
长度，N(i)为分支i的空间步长数，z_i(0)和z_i(L_i)分别为分支i的始末节点标高。
6.根据权利要求1所述的一种含高浓度瓦斯井巷瞬态风流稳定性计算方法，其特征在
于：所述瓦斯弥散方程(1)的求解方法采用Crank-Nicholson方法进行有限差分，得到三对
角方程组，通风系统中风流的流动为一维不可压缩流，通风网络中各井巷分支最后一个空
间步长在某时刻的瓦斯浓度与上一时刻前一步长的瓦斯浓度相等，瓦斯释放后立即与风流
混合，且巷道交汇处各风流中的瓦斯立即相互混合，利用追赶法求解三对角方程组，得到井
巷分支瓦斯浓度。
7.根据权利要求1所述的一种含高浓度瓦斯井巷瞬态风流稳定性计算方法，其特征在
于：所述的井巷通风网络方程(2)的求解方法式(2)，令
将各变量表示成矩阵形式为：A＝[a_jk]，C＝
[c_ji]，K＝[k_ii]，C^T＝[c_ik]，D＝(D₁,D₂,…D_b)^T，q＝[q₁,q₂,…,q_k]^T(k＝1,2,…,b)，其中k_ii＝
K_i，且A＝CKC^T，式(2)表示成矩阵形式为选用改进的欧拉公式(龙格-库塔法)求解
各井巷分支的风量，在此基础上，根据井巷分支的断面积，得出井巷分支的风速。
8.根据权利要求1所述的一种含高浓度瓦斯井巷瞬态风流稳定性计算方法，其特征在
于：所述的求解过程中需对方程(1)和方程(2)进行耦合求解，t₀时刻各参数的值已知，先求
解式(1)t₁时刻瓦斯浓度，然后将解得的值代入式(2)求解t₁时刻风量，得到的风量代入式
(1)求解t₂时刻的瓦斯浓度，依次类推即可得到各个时间步长巷道的浓度分布和风量大小。

一种含高浓度瓦斯井巷瞬态风流稳定性计算方法

技术领域

本发明涉及一种井巷瞬态风流稳定性计算方法，尤其涉及一种计算含高浓度瓦斯
井巷瞬态风流稳定性的方法。

背景技术

高浓度瓦斯对井巷风流的稳定性产生一定的影响；以煤与瓦斯突出为例，当突出
发生以后，高浓度瓦斯积聚在井巷，引起井巷风流逆转，甚至诱导瓦斯爆炸等二次灾害的发
生。前人对火灾时期井巷瞬态风流稳定性进行了较为深入的研究，主要考虑在火风压作用
下，井巷风流的灾变规律。周延，王省身分别采用静态分析法（布德雷克法）和动态分析法对
上行风流火灾旁侧支路风流逆转的临界条件进行了对比分析；张兴凯，王振财等根据矿井
火灾时期的燃烧规律以及通风网路风流流动特征，分析了矿井火灾时期风流非稳态过程的
通风原理，建立了矿井火灾燃烧数学模型及风流非稳定流动方程组，并给出了解算方法，同
时通过编制计算机软件实现了对矿井火灾燃烧过程及其风流变化过程的模拟。瓦斯气流在
井巷中的运移是瞬态的，其对井巷分支风流的影响也是动态的。目前国内外学者对含高浓
度瓦斯井巷瞬态风流稳定性研究较少，尚需提出一种可靠的计算方法。研究结果对矿井抗
灾救护、灾变通风以及有效防止二次灾害事故的发生均具有重要的理论和现实意义。

具体内容

鉴于现有计算方法的不足，本发明提供了一种含高浓度瓦斯井巷瞬态风流稳定性
计算方法，该方法考虑高浓度瓦斯对井巷分支瞬态风流的影响，利用该计算方法可以实现
求解通风网络任意井巷分支在不同时刻的风量变化规律。

为解决上述技术问题，本发明专利采用如下技术方案：其计算模型由瓦斯弥散方
程

和井巷通风网络方程

构成，其中式（1）中u_i为井巷分支i的风速，x_i为沿井巷分支i的空间距离，t为时间，c_i为
井巷分支i的瓦斯体积分数，E_x为瓦斯弥散系数，式（2）中q_k(k＝1，2，…，b)为分支风量，R_i为
分支i风阻，为分支i的通风机风压，为分支i的瓦斯风压，c_ji为回路矩阵，c_ik为回路矩
阵c_ji的转置矩阵，其中瓦斯弥散方程（1）求解方法是首先采用Crank-Nicholson方法进行有
限差分，然后利用追赶法得到井巷分瓦斯浓度，同时井巷通风网络方程（2）求解方法是采用
龙格-库塔法求解井巷分支的风量（风速），求解过程中需对方程（1）和方程（2）进行耦合求
解。

所述的瓦斯弥散系数E_x，其特征在于瓦斯弥散系数其中r为圆管半
径，α为摩擦阻力系数，u为井巷分支风流速度。

所述的分支风量，其特征在于，分支风量q_k(k＝1,2,…,b)为待求参数，也是时间
的函数，即所求的值是b个未知函数，b为通风网络中独立回路的个数；q_k是评判井巷风流稳
定性的指标，q_k＝0表示风流停滞，q_k＜0表示风流逆转，q_k＞0表示风流方向保持不变。

所述的通风机风压，其特征在于，风机风压a₀，a₁，…，a_n为曲线
拟合系数；给定三个或三个以上风机的风压、风量数据，通过最小二乘法求解拟合系数。

所述的瓦斯风压其特征在于，为分支i的平均密度，
为t_n时刻分支第j个空间步长的气体密度，
为t_n时刻分支第j个空间步长的瓦斯体积分数，ρ_m和ρ_a分别为甲烷和空气的密度，h为空间
步长，L_i为井巷分支i的长度，N(i)为分支i的空间步长数，z_i(0)和z_i(L_i)分别为分支i的始
末点标高。

所述的瓦斯弥散方程（1）的求解方法，其特征在于，采用Crank-Nicholson方法进
行有限差分，得到三对角方程组，通风系统中风流的流动为一维不可压缩流，通风网络中各
巷道最后一个空间步长在某时刻的瓦斯浓度与上一时刻前一步长的瓦斯浓度相等，瓦斯释
放后立即与风流混合，且巷道交汇处各风流中的瓦斯立即相互混合，利用追赶法求解三对
角方程组，得到井巷分瓦斯浓度。

井巷通风网络方程(2)的求解方法，其特征在于，令
将各变量表示成矩阵形式为：A＝[a_jk]，C＝
[c_ji]，K＝[k_ii]，C^T＝[c_ik]，D＝(D₁,D₂,…D_b)^T，q＝[q₁,q₂,…,q_k]^T(k＝1,2,…,b)，其中k_ii=
K_i，且A＝CKC^T，式（2）表示成矩阵形式为选用改进的欧拉公式（龙格-库塔法）求解
各井巷分支的风量，在此基础上，根据井巷分支的断面积，得出井巷分支的风速。

所述的求解过程中需对方程（1）和方程（2）进行耦合求解，其特征在于，对所述的
式（1）和式（2）进行耦合求解，t₀时刻各参数的值已知，先求解式（1）t₁时刻瓦斯浓度，然后将
解得的值代入式（2）求解t₁时刻风量，得到的风量代入式（1）求解t₂时刻的瓦斯浓度，依次类
推即可得到各个时间步长巷道的浓度分布和风量大小。

附图说明

图1是本发明的含高浓度瓦斯井巷瞬态风流稳定性计算方法整体示意图；

具体实施方式

下面将对本发明实施例作进一步地详细描述；由于对于任意井巷分支，瓦斯浓度
弥散方程为：由于式（1）中u和E_x并非常数，而是时间的函数，该模型为
非线性偏微分方程，其解析解难于计算，所以本发明采用有限差分的方法计算上式的数值
解；即在微分方程中用差商代替偏导数，得到相应的差分方程，通过解差分方程得到微分方
程解的近似值。在进行差分之前，先将各巷道按时间和空间位置划分差分网格。假设巷道的
长度为L，沿巷道风流方向划分N段区间，则每段区间长度h=L/N，x_j为第j个空间步长距始节
点的距离，表示为jh；同时将时间按同一时间步长τ划分区间，则t_n时刻对应的时间为nτ，采
用Crank-Nicholson差分格式。对于井巷中瓦斯浓度弥散方程按Crank-Nicholson格式差分
为

化简可得：
其中：

上面的瓦斯弥散模型建立在
单一巷道基础上，而任何一个矿井通风系统都存在着许多条巷道，风流在各相关联的巷道
运移过程中，前后分支存在着质的连续传递。因此需要建立各分支的瓦斯浓度的网络联合
模型。为此，假设：1）通风系统中风流的流动为一维不可压缩流；2）通风网络中各巷道最后
一个空间步长在某时刻的瓦斯浓度与上一时刻前一步长的瓦斯浓度相等；3）瓦斯释放后立
即与风流混合，且巷道交汇处各风流中的瓦斯立即相互混合。根据上述假设，若通风网络总
分支数为m，第i个分支划分的空间步长数为N（i），则通风网络中任意分支i的瓦斯浓度可按
下式计算，即

式中，i＝1,2,…,m；j＝1,2,…,N(i)-1。任意分支i的瓦斯浓度边界条件为

式中：表示t_n+1时刻分支i最后一个空间步长的瓦斯浓度；表示t_n时刻分支i倒
数第二个空间步长的瓦斯浓度；和Q_i'分别表示t_n+1时刻，由相邻分支流入分支i的单位体
积瓦斯量和单位体积风量。显然要得到边界条件（5），就需要先求解出和Q_i'。为此，定义
矩阵表示通风网络中分支i与分支k的风流关系，且

时，记δ=1，表示分支k的风流流入分支i；时，记δ=0，表示分支k的
风流流出分支i或者分支k和分支i不相关联。则有

式中，表示t_n+1时刻分支k的风量。将式（6）和式（7）代入式（5），得

由于本发明不考虑同一巷道的截面变化，所以对于同一条巷道其风速u和弥散系数E_x
不随空间位置发生变化，仅是时间的函数。所以可以令A[x_j]＝a，B[x_j]＝b，C[x_j]＝c，D[x_j]
＝d_j，且

式中分别表示t_n时刻分支i的风流速度和弥散系数。对于分支i，若t_n时刻的风
速和弥散系数已知，结合空间步长和时间步长则可计算出a、b、c，同时根据该分支中各
空间步长的瓦斯浓度，可以得出d_j(2≤j≤N-1)，又根据分支i的瓦斯浓度边界条件，可以得
到t_n+1时刻巷道初末两个步长的瓦斯浓度x₁和x_N，将巷道中2≤j≤N-1空间步长的数据代入
式（3）并展开得

该方程组为三对角方程组，可以用追赶法进行求解，从而得到t_n+1时刻分支i各空间步
长的瓦斯浓度。井巷通风网络方程表达式为：

式（2）中R_i为分支i风阻，为分支i的通风机风压，可以通过风机特性曲线方程计算得
到；为分支i的瓦斯风压，与井巷分支中的瓦斯浓度、巷道高差相关。是时间和风
量的函数，风量q_k(k＝1,2,…,b)为待求参数，也是时间的函数，即所求的值是b个未知函
数。式（2）较为复杂，难于直接对其进行求解，因此需要根据相关数值解法对其进行转化，本
发明采用的是改进的欧拉公式（龙格-库塔法）方法。令

为方便求解，将各变量表示成矩阵形式为：A＝[a_jk]，C＝[c_ji]，K＝[k_ii]，C^T＝[c_ik]，D＝
(D₁,D₂,…D_b)^T，q＝[q₁,q₂,…,q_k]^T(k＝1,2,…,b)，其中k_ii=K_i，且A＝CKC^T。所以式（2）表示成
矩阵形式为

即

由式（11）可以列出b个独立的微分方程，方程可解。采用有限差分法，选用改进的欧拉
公式（龙格-库塔法）求解上式能达到足够的精度。根据有限差分理论，先将时间按同一时间
步长τ划分区间，则t_n时刻对应的时间为nτ。由上述分析可知，式（11）实际上就是风量的偏
导关于时间和基本分支风量的函数，所以可表示成q'_k＝f_k(t,q₁,q₂,…,q_b)(k＝1,2,…,b)，
则根据改进的欧拉公式可得

f_k(t,q₁,q₂,…,q_b)(k＝1,2,…,b)的表达式应当是关于时间和各基本分支风量的函
数，所以要正确求解上式，需要把其他变量转换成关于时间和各基本分支风量的函数。根据
前面的推导，该表达式中还包含各分支风机风压和瓦斯风压两个变量。风机的风压可
以通过风机特性曲线方程计算得到，风压特性曲线反映的是通风机的风压与风量的关系。
通常可用下面多项式拟合：h_F＝a₀+a₁q+a₂q²+a₃q³+...+a_nqⁿ，式中a₀，a₁，…，a_n为曲线拟合系
数，q为风机所在分支风量。曲线的多项式次数根据计算精度要求确定，本文取到2即可。即
h_F＝a₀+a₁q+a₂q²。为了拟合风机性能曲线的工作段，给定其工作段上三个点(q₁,h₁),(q₂,
h₂),(q₃,h₃)，代入上式，即可得到含三个未知数的线性方程。解此联立方程组就可以得到风
压特性曲线方程中的各项拟合系数。

所以对于任意分支i，若已知风机所在分支风量Q_i，则该分支风机风压可以表示成：

此外，欲求解瓦斯风压应先计算各分支瓦斯的平均密度即

由于该气体为甲烷和空气的混合气体，所以

式中，M_i为分支i的空气总质量，kg；V_i为分支i的空气总体量，m³；A_i为分支i的截面积，
m²；h为空间步长，m；为t_n时刻分支第j个空间步长的气体密度，kg/m³；为t_n时刻分支第j
个空间步长的瓦斯浓度（体积分数）；ρ_m为甲烷密度，kg/m³；ρ_a为空气密度，kg/m³。式(1)中的
速度u由风量q计算得到，而在采用改进的欧拉方法求解式（2）风量q时又要用到式（1）中的
浓度c，所以本发明将式(1)和式（2）进行耦合求解。t₀时刻各参数的值已知，先求解式（1）t₁
时刻瓦斯浓度，然后将解得的值代入式（2）求解t₁时刻风量，得到的风量代入式（1）求解t₂时
刻的瓦斯浓度，依次类推即可得到各个时间步长巷道的浓度分布和风量大小。