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

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

【技术实现步骤摘要】

本专利技术涉及一种井巷瞬态风流稳定性计算方法，尤其涉及一种计算含高浓度瓦斯井巷瞬态风流稳定性的方法。
技术介绍
高浓度瓦斯对井巷风流的稳定性产生一定的影响；以煤与瓦斯突出为例，当突出发生以后，高浓度瓦斯积聚在井巷，引起井巷风流逆转，甚至诱导瓦斯爆炸等二次灾害的发生。前人对火灾时期井巷瞬态风流稳定性进行了较为深入的研究，主要考虑在火风压作用下，井巷风流的灾变规律。周延，王省身分别采用静态分析法（布德雷克法）和动态分析法对上行风流火灾旁侧支路风流逆转的临界条件进行了对比分析；张兴凯，王振财等根据矿井火灾时期的燃烧规律以及通风网路风流流动特征，分析了矿井火灾时期风流非稳态过程的通风原理，建立了矿井火灾燃烧数学模型及风流非稳定流动方程组，并给出了解算方法，同时通过编制计算机软件实现了对矿井火灾燃烧过程及其风流变化过程的模拟。瓦斯气流在井巷中的运移是瞬态的，其对井巷分支风流的影响也是动态的。目前国内外学者对含高浓度瓦斯井巷瞬态风流稳定性研究较少，尚需提出一种可靠的计算方法。研究结果对矿井抗灾救护、灾变通风以及有效防止二次灾害事故的发生均具有重要的理论和现实意义。具体内容鉴于现有计算方法的不足，本专利技术提供了一种含高浓度瓦斯井巷瞬态风流稳定性计算方法，该方法考虑高浓度瓦斯对井巷分支瞬态风流的影响，利用该计算方法可以实现求解通风网络任意井巷分支在不同时刻的风量变化规律。为解决上述技术问题，本专利技术专利采用如下技术方案：其计算模型由瓦斯弥散方程和井巷通风网络方程构成，其中式（1）中ui为井巷分支i的风速，xi为沿井巷分支i的空间距离，t为时间，ci为井巷分支i的瓦斯体积分数，Ex为瓦斯弥散系数，式（2）中qk(k＝1，2，…，b)为分支风量，Ri为分支i风阻，为分支i的通风机风压，为分支i的瓦斯风压，cji为回路矩阵，cik为回路矩阵cji的转置矩阵，其中瓦斯弥散方程（1）求解方法是首先采用Crank-Nicholson方法进行有限差分，然后利用追赶法得到井巷分瓦斯浓度，同时井巷通风网络方程（2）求解方法是采用龙格-库塔法求解井巷分支的风量（风速），求解过程中需对方程（1）和方程（2）进行耦合求解。所述的瓦斯弥散系数Ex，其特征在于瓦斯弥散系数其中r为圆管半径，α为摩擦阻力系数，u为井巷分支风流速度。所述的分支风量，其特征在于，分支风量qk(k＝1,2,…,b)为待求参数，也是时间的函数，即所求的值是b个未知函数，b为通风网络中独立回路的个数；qk是评判井巷风流稳定性的指标，qk＝0表示风流停滞，qk＜0表示风流逆转，qk＞0表示风流方向保持不变。所述的通风机风压，其特征在于，风机风压a0，a1，…，an为曲线拟合系数；给定三个或三个以上风机的风压、风量数据，通过最小二乘法求解拟合系数。所述的瓦斯风压其特征在于，为分支i的平均密度，为tn时刻分支第j个空间步长的气体密度，为tn时刻分支第j个空间步长的瓦斯体积分数，ρm和ρa分别为甲烷和空气的密度，h为空间步长，Li为井巷分支i的长度，N(i)为分支i的空间步长数，zi(0)和zi(Li)分别为分支i的始末点标高。所述的瓦斯弥散方程（1）的求解方法，其特征在于，采用Crank-Nicholson方法进行有限差分，得到三对角方程组，通风系统中风流的流动为一维不可压缩流，通风网络中各巷道最后一个空间步长在某时刻的瓦斯浓度与上一时刻前一步长的瓦斯浓度相等，瓦斯释放后立即与风流混合，且巷道交汇处各风流中的瓦斯立即相互混合，利用追赶法求解三对角方程组，得到井巷分瓦斯浓度。井巷通风网络方程(2)的求解方法，其特征在于，令将各变量表示成矩阵形式为：A＝[ajk]，C＝[cji]，K＝[kii]，CT＝[cik]，D＝(D1,D2,…Db)T，q＝[q1,q2,…,qk]T(k＝1,2,…,b)，其中kii=Ki，且A＝CKCT，式（2）表示成矩阵形式为选用改进的欧拉公式（龙格-库塔法）求解各井巷分支的风量，在此基础上，根据井巷分支的断面积，得出井巷分支的风速。所述的求解过程中需对方程（1）和方程（2）进行耦合求解，其特征在于，对所述的式（1）和式（2）进行耦合求解，t0时刻各参数的值已知，先求解式（1）t1时刻瓦斯浓度，然后将解得的值代入式（2）求解t1时刻风量，得到的风量代入式（1）求解t2时刻的瓦斯浓度，依次类推即可得到各个时间步长巷道的浓度分布和风量大小。附图说明图1是本专利技术的含高浓度瓦斯井巷瞬态风流稳定性计算方法整体示意图；具体实施方式下面将对本专利技术实施例作进一步地详细描述；由于对于任意井巷分支，瓦斯浓度弥散方程为：由于式（1）中u和Ex并非常数，而是时间的函数，该模型为非线性偏微分方程，其解析解难于计算，所以本专利技术采用有限差分的方法计算上式的数值解；即在微分方程中用差商代替偏导数，得到相应的差分方程，通过解差分方程得到微分方程解的近似值。在进行差分之前，先将各巷道按时间和空间位置划分差分网格。假设巷道的长度为L，沿巷道风流方向划分N段区间，则每段区间长度h=L/N，xj为第j个空间步长距始节点的距离，表示为jh；同时将时间按同一时间步长τ划分区间，则tn时刻对应的时间为nτ，采用Crank-Nicholson差分格式。对于井巷中瓦斯浓度弥散方程按Crank-Nicholson格式差分为化简可得：其中：上面的瓦斯弥散模型建立在单一巷道基础上，而任何一个矿井通风系统都存在着许多条巷道，风流在各相关联的巷道运移过程中，前后分支存在着质的连续传递。因此需要建立各分支的瓦斯浓度的网络联合模型。为此，假设：1）通风系统中风流的流动为一维不可压缩流；2）通风网络中各巷道最后一个空间步长在某时刻的瓦斯浓度与上一时刻前一步长的瓦斯浓度相等；3）瓦斯释放后立即与风流混合，且巷道交汇处各风流中的瓦斯立即相互混合。根据上述假设，若通风网络总分支数为m，第i个分支划分的空间步长数为N（i），则通风网络中任意分支i的瓦斯浓度可按下式计算，即式中，i＝1,2,…,m；j＝1,2,…,N(i)-1。任意分支i的瓦斯浓度边界条件为式中：表示tn+1时刻分支i最后一个空间步长的瓦斯浓度；表示tn时刻分支i倒数第二个空间步长的瓦斯浓度；和Qi'分别表示tn+1时刻，由相邻分支流入分支i的单位体积瓦斯量和单位体积风量。显然要得到边界条件（5），就需要先求解出和Qi'。为此，定义矩阵表示通风网络中分支i与分支k的风流关系，且时，记δ=1，表示分支k的风流流入分支i；时，记δ=0，表示分支k的风流流出分支i或者分支k和分支i不相关联。则有式中，表示tn+1时刻分支k的风量。将式（6）和式（7）代入式（5），得由于本专利技术不考虑同一巷道的截面变化，所以对于同一条巷道其风速u和弥散系数Ex不随空间位置发生变化，仅是时间的函数。所以可以令A[xj]＝a，B[xj]＝b，C[xj]＝c，D[xj]＝dj，且式中分别表示tn时刻分支i的风流速度和弥散系数。对于分支i，若tn时刻的风速和弥散系数已知，结合空间步长和时间步长则可计算出a、b、c，同时根据该分支中各空间步长的瓦斯浓度，可以得出dj(2≤j≤N-1)，又根据分支i的瓦斯浓度边界条件，可以得到tn+1时刻巷道初末两个步长的瓦斯浓度x1和xN，将巷道中2≤j本文档来自技高网...

【技术保护点】
一种含高浓度瓦斯井巷瞬态风流稳定性计算方法，其计算模型由瓦斯弥散方程∂ci∂t+ui∂ci∂xi=Ex∂2ci∂xi2---(1)]]>和井巷通风网络方程Σi=1m[cjiKi(Σk=1bcikdqkdt)]=Σi=1mcji(hFi+hMi-Ri|Σk=1bcikqk|Σk=1bcikqk)---(2)]]>构成，其中式(1)中ui为井巷分支i的风速，xi为沿井巷分支i的空间距离，t为时间，ci为井巷分支i的瓦斯体积分数，Ex为瓦斯弥散系数，式(2)中为分支风量，Ri为分支i风阻，为分支i的通风机风压，为分支i的瓦斯风压，cji为回路矩阵，cik为回路矩阵cji的转置矩阵，其中瓦斯弥散方程(1)求解方法是首先采用Crank‑Nicholson方法进行有限差分，然后利用追赶法得到井巷分瓦斯浓度，同时井巷通风网络方程(2)求解方法是采用龙格‑库塔法求解井巷分支的风量(风速)，求解过程中需对方程(1)和方程(2)进行耦合求解。

【技术特征摘要】
1.一种含高浓度瓦斯井巷瞬态风流稳定性计算方法，其计算模型由瓦斯弥散方程∂ci∂t+ui∂ci∂xi=Ex∂2ci∂xi2---(1)]]>和井巷通风网络方程Σi=1m[cjiKi(Σk=1bcikdqkdt)]=Σi=1mcji(hFi+hMi-Ri|Σk=1bcikqk|Σk=1bcikqk)---(2)]]>构成，其中式(1)中ui为井巷分支i的风速，xi为沿井巷分支i的空间距离，t为时间，ci为井巷分支i的瓦斯体积分数，Ex为瓦斯弥散系数，式(2)中为分支风量，Ri为分支i风阻，为分支i的通风机风压，为分支i的瓦斯风压，cji为回路矩阵，cik为回路矩阵cji的转置矩阵，其中瓦斯弥散方程(1)求解方法是首先采用Crank-Nicholson方法进行有限差分，然后利用追赶法得到井巷分瓦斯浓度，同时井巷通风网络方程(2)求解方法是采用龙格-库塔法求解井巷分支的风量(风速)，求解过程中需对方程(1)和方程(2)进行耦合求解。2.根据权利要求1所述的一种含高浓度瓦斯井巷瞬态风流稳定性计算方法，其特征在于：所述的瓦斯弥散系数其中r为圆管半径，α为摩擦阻力系数，u为井巷分支的风速。3.根据权利要求1所述的一种含高浓度瓦斯井巷瞬态风流稳定性计算方法，其特征在于：所述的分支风量qk(k＝1,2,…,b)为待求参数，也是时间的函数，即所求的值是b个未知函数，b为通风网络中独立回路的个数；qk是评判井巷风流稳定性的指标，qk＝0表示风流停滞，qk＜0表示风流逆转，qk＞0表示风流方向保持不变。4.根据权利要求1所述的一种含高浓度瓦斯井巷瞬态风流稳定性计算方法，其特征在于：所述的通风机风压a0，a1，…，an为曲线拟合系数；给定三个或三个以上风机的风压、风量数据，...

【专利技术属性】
技术研发人员：周爱桃，王凯，范灵鹏，刘昂，
申请(专利权)人：中国矿业大学北京，
类型：发明
国别省市：北京;11