CN101806224B - 井下煤层气抽采产能预测方法 - Google Patents

Abstract

本发明提出了一种井下煤层气抽采产能预测方法，考虑了煤层渗透率的实际变化，并能对抽采钻孔组(群)的抽采产能进行预测，包括如下步骤：1)测取待测煤层原始参数；2)通过求解本发明提出的模型，获得抽采后的平均渗透率、抽采后的煤层气压力、抽采后的煤层气平均含量以及煤层气抽采量；本发明提出了的井下煤层气抽采产能预测方法，以达西渗流定律为基础，建立了基于抽采过程中渗透率变化的井下煤层气产能预测的数学模型，由于数学模型中考虑了抽采过程中渗透率的实际变化，因此能正确反映煤层气的渗流状况。在进一步的求解步骤中，采用差分离散的方法对数学模型进行求解，收敛速度快，精度更高，稳定性好，提高计算的精度和速度。

Description

井下煤层气抽采产能预测方法

技术领域

本发明涉及煤气开采技术领域，具体涉及一种井下煤层气抽采产能的预测方法。

背景技术

煤层气，即煤矿瓦斯以吸附状态或游离状态赋存于煤层中，通常吸附状态的瓦斯量占90％以上，游离状态的瓦斯量占10％左右。当煤层开采时，存在于煤层中的瓦斯随煤层开采的扰动不断释放，从而造成瓦斯爆炸、煤与瓦斯突出、瓦斯窒息等危害。对井下煤层气进行抽采，不仅是治理煤矿瓦斯灾害的有效方法，同时也是煤层气资源开发的一种不可或缺的方式。

而对井下煤层气产能进行预测，能引导对特定区块煤层气资源进行有序开发，并可为井下煤层气抽采钻孔的合理布置提供依据，提高抽放效果、缩短抽放时间，进而实现矿井的安全生产。

目前井下煤层气产能预测方法有许多缺陷，如以下几个方面：仅针对单孔煤层气抽采量，而不能对抽采钻孔组(群)的产能进行预测；煤层气抽采模型是建立在抽采过程中渗透率不变的情况下，但根据现场生产实际表明，在抽采过程中煤层的渗透率是不断变化的，现有模型均不能够满足生产实际。

发明内容

有鉴于此，为了解决上述问题，本发明提出了一种井下煤层气抽采产能预测方法，考虑了煤层渗透率的实际变化，并能对抽采钻孔组(群)的抽采产能进行预测。

本发明的目的是这样实现的：井下煤层气抽采产能预测方法，包括如下步骤：

井下煤层气抽采产能预测方法，其特征在于：包括如下步骤：

1)测取待测煤层原始参数；

2)通过求解如下模型，获得抽采后的平均渗透率、抽采后的煤层气压力、抽采后的煤层气平均含量以及煤层气抽采量：

$\frac{\∂}{\∂x}\left(K_x\frac{\∂P}{\∂x}\right)+\frac{\∂}{\∂y}\left(K_y\frac{\∂P}{\∂y}\right)+\frac{\∂}{\∂z}\left(K_z\frac{\∂P}{\∂z}\right)=A\left(P\right)\frac{\∂P}{\∂t};$

其中：

$K=K_0(1+{(1-\mathrm{Rlg}\frac{\mathrm{bp}+1}{{\mathrm{bp}}_0+1})}^2-e^{3C_{\mathrm{\φ}}(p-p_0)});$

其中：p₀为初始煤层气压力，单位MPa；

p为煤层气压力，单位MPa；

b为煤的吸附常数，单位1/MPa；

R为煤基质收缩系数；

K₀为初始煤层气渗透率；

$A\left(P\right)=\frac{2\mathrm{\φ\μ}}{p_N}+\frac{100-W-A}{100*(1+0.31W)}\·\mathrm{\μ}\·{\mathrm{\γ}}_m[\frac{\mathrm{ab}}{(1+b\sqrt{P})}+\frac{\mathrm{ab}}{{(1+b\sqrt{P})}^2}]$

K为煤层渗透率，m²；

μ为煤层气的绝对粘度，Pa·s；

a为煤的最大煤层吸附量，单位m³/t

b为煤的吸附常数，单位1/MPa；

W为煤层中单位体积煤中煤层气总量；

A为单位体积煤中游离煤层气含量，m³/m³；

P为煤层气压力，单位MPa；

γ_m为煤的容重，单位t/m³；

φ为煤层孔隙率；

p_N为1标准大气压，0.1013MPa。

进一步，所述步骤2)的求解过程具体包括如下步骤：

21)将待测煤层原始参数作为初始值计算抽采范围，对抽采范围进行网格划分；

22)根据步骤21)的网格划分，计算系数矩阵并建立求解矩阵：

系数矩阵一：

$A_{i,j}=\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}$ $B_{i,j}=-[\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}+\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}+A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}],$ $C_{i,j}=\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2},$

$D_{i,j}=-\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j-1}^k+[\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}+\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}-A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}]P_{i,j}^k-\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j+1}^k$

$\left(\begin{array}{cccccc}{B}_{1,j}& C_{1,j}& ...& & & 0\\ A_{2,j}& B_{2,j}& C_{2,j}& ...& & 0\\ 0& A_{3,j}& B_{3,j}& C_{3,j}& ...& 0\\ ......& & & & & \\ 0& 0& 0& ...& A_{m-1,j}& B_{m-1,j}\end{array}\right)$

系数矩阵二：

$A_{i,j}^{\′}=\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2},$ $B_{i,j}^{\′}=-[\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}+\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}+A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}],$

$C_{i,j}^{\′}=\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2},$ $D_{i,j}^{\′}=-\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i-1,j}^{k+\frac{1}{2}}+[\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}+\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}-A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}]P_{i,j}^{k+\frac{1}{2}}-\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i+1,j}^{k+\frac{1}{2}}$

$\left(\begin{array}{cccccc}{B}_{i,1}^{\′}& C_{i,1}^{\′}& ...& & & 0\\ A_{i,2}^{\′}& B_{i,2}^{\′}& C_{i,2}^{\′}& ...& & 0\\ 0& A_{i,3}^{\′}& B_{i,3}^{\′}& C_{i,3}^{\′}& ...& 0\\ ......& & & & & \\ 0& 0& 0& ...& A_{i,n-1}^{\′}& B_{i,n-1}^{\′}\end{array}\right)$

求解矩阵：

$\left(\begin{array}{cccccc}{B}_{1,j}& C_{1,j}& ...& & & 0\\ A_{2,j}& B_{2,j}& C_{2,j}& ...& & 0\\ 0& A_{3,j}& B_{3,j}& C_{3,j}& ...& 0\\ ......& & & & & \\ 0& 0& 0& ...& A_{m-1,j}& B_{m-1,j}\end{array}\right)\left(\begin{array}{c}{P}_{1,j}\\ P_2{,}_j\\ P_{3,j}\\ .\\ .\\ .\\ P_{m-1,j}\end{array}\right)=\left(\begin{array}{c}{D}_{1,j}\\ D_{2,j}\\ D_{3,j}\\ .\\ .\\ .\\ D_{m-1,j}\end{array}\right)$

$\left(\begin{array}{cccccc}{B}_{i,1}^{\′}& C_{i,1}^{\′}& ...& & & 0\\ A_{i,2}^{\′}& B_{i,2}^{\′}& C_{i,2}^{\′}& ...& & 0\\ 0& A_{i,3}^{\′}& B_{i,3}^{\′}& C_{i,3}^{\′}& ...& 0\\ ......& & & & & \\ 0& 0& 0& ...& A_{i,n-1}^{\′}& B_{i,n-1}^{\′}\end{array}\right)\left(\begin{array}{c}{P}_{i,1}\\ P_{i,2}\\ P_{i,3}\\ .\\ .\\ .\\ P_{i,n-1}\end{array}\right)=\left(\begin{array}{c}{D}_{i,1}^{\′}\\ D_{i,2}^{\′}\\ D_{i,3}^{\′}\\ .\\ .\\ .\\ D_{i,n-1}^{\′}\end{array}\right))$

23)计算第一个时间步长的各结点煤层气压力和渗透率；

24)将计算所得的各结点煤层气压力和渗透率赋为初值；

25)利用步骤24)所得的初值重新计算抽采范围并进行网格划分；

26)根据步骤25)的网格划分，计算系数矩阵并建立求解矩阵；

27)计算下一时间步长的各结点煤层气压力和渗透率，并返回步骤24)，直到计算完成，计算完成的条件为

本发明提出了的井下煤层气抽采产能预测方法，以达西渗流定律为基础，建立了基于抽采过程中渗透率变化的井下煤层气产能预测的数学模型，由于数学模型中考虑了抽采过程中渗透率的实际变化，因此能正确反映煤层气的渗流状况。在进一步的求解步骤中，采用差分离散的方法对数学模型进行求解，收敛速度快，精度更高，稳定性好，提高计算的精度和速度。

本发明的其他优点、目标，和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书，权利要求书，以及附图中所特别指出的结构来实现和获得。

附图说明

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步的详细描述：

图1示出了本发明井下煤层气抽采产能预测方法的流程示意图；

图2示出了网格划分示意图；

图3示出了煤层透气性系数分别为4m2/MPa2.d的顺层单孔抽采量随时间变化图；

图4示出了煤层透气性系数分别为60m2/MPa2.d的顺层单孔抽采量随时间变化图；

图5示出了煤层透气性系数分别为40m2/MPa2.d的顺层平行钻孔抽采量随时间变化图；

图6示出了煤层透气性系数分别为60m2/MPa2.d的顺层平行钻孔抽采量随时间变化图；

图7示出了煤层透气性系数分别为4m2/MPa2.d的穿层单孔抽采量随时间变化图；

图8示出了煤层透气性系数分别为10m2/MPa2.d的穿层单孔抽采量随时间变化图；

图9示出了煤层透气性系数分别为4m2/MPa2.d的穿层平行钻孔抽采量随时间变化图；

图10示出了煤层透气性系数分别为10m2/MPa2.d的穿层平行钻孔抽采量随时间变化图。

具体实施方式

以下将对本发明的优选实施例进行详细的描述。

影响煤层气抽采过程中渗透率变化的主要因素是有效应力和煤基质收缩效应。但两者对煤层渗透率的影响正好相反，一正一负。为此，如何衡量两种因素对煤层渗透率影响的综合效应就成为一个必须解决的问题。

在煤层气抽采过程中，随着煤层气压力的下降，有效应力增大，使得煤层渗透率呈减小趋势，而此时煤基质收缩，煤层渗透率呈增大趋势。两种因素对煤层渗透率的影响程度可分别用以下两式来表示：

${\mathrm{\ΔK}}_1=K_0{e}^{3C_{\mathrm{\φ}}(p-p_0)}-K_0...(1-1)$

${\mathrm{\ΔK}}_2=K_0{(1-\mathrm{\θlg}\frac{\mathrm{bp}+1}{{\mathrm{bp}}_0+1})}^2-K_0...(1-2)$

式中：

ΔK₁为有效应力作用下的渗透率变化量；

ΔK₂为煤基质收缩作用下的渗透率变化量；

P₀为初始煤层气压力，单位MPa；

P为煤层气压力，单位MPa；

b为吸附常数，为1/MPa；

θ为煤基质收缩系数；

C_φ为孔隙压缩系数；

K₀为初始煤层气渗透率。

则煤层渗透率最终变化量为两种效应之差，即：

${\mathrm{\ΔK}}_2-{\mathrm{\ΔK}}_1=K_0{(1-\mathrm{\θlg}\frac{\mathrm{bp}+1}{{\mathrm{bp}}_0+1})}^2-K_0{e}^{3C_{\mathrm{\φ}}(p-p_0)}$

因此，当煤层气压力从p₀→p时，煤层的渗透率从K₀→K。而根据上式有：

$K=K_0+K_0[{(1-\mathrm{\θlg}\frac{\mathrm{bp}+1}{{\mathrm{bp}}_0+1})}^2-e^{3C_{\mathrm{\φ}}(p-p_0)}]$

即：

$K=K_0(1+{(1-\mathrm{\θlg}\frac{\mathrm{bp}+1}{{\mathrm{bp}}_0+1})}^2-e^{3C_{\mathrm{\φ}}(p-p_0)})...(1-3)$

式(1-3)即为在煤层气抽采过程中渗透率的变化规律。

煤在形成过程中受到很多因素的影响，使得煤层具有非均质性，但从宏观上看，在一个较大的区域内(除断层等地质构造以外)可以看作是均质的。煤层气的流固耦合渗流规律是一个复杂的问题，涉及流体力学、岩石力学等诸多学科，为此为了使问题简化，突出重点，按如下假设来推导煤层气流动方程：

1)煤层中的煤层气含量满足朗格缪尔方程：

煤层中煤层气含量由两部分组成，一部分是处在吸附状态下，而另一部分是以游离形式赋存的。吸附煤层气量由朗格缪尔公式给出：

${\mathrm{\Γ}}_1=\frac{100-W-A}{100(1+0.31W)}\·\left(\frac{\mathrm{abp}}{1+\mathrm{bp}}\right){\mathrm{\γ}}_m...(2-1)$

式中：

Γ₁-单位体积煤中吸附煤层气含量，m³/m³；

a-煤的最大煤层气吸附量(m³/t)；

b-煤的吸附常数(m³/t或1/MPa)；

p-煤层气(瓦斯)压力，MPa；

γ_m-煤的容重，(t/m³)。

煤中游离的煤层气由下列公式给出：

${\mathrm{\Γ}}_2=\mathrm{\φ}\·\frac{p}{p_N}...(2-2)$

式中：

Γ₂-单位体积煤中游离煤层气含量，m³/m³；

φ-煤层孔隙率；

p_N-1标准大气压，0.1013MPa。

煤层中单位体积煤中煤层气总量(体积含量)W为：

$\mathrm{\Γ}={\mathrm{\Γ}}_1+{\mathrm{\Γ}}_2=\mathrm{\φ}\frac{p}{p_N}+\frac{100-W-A}{100*(1+0.31W)}\·\frac{\mathrm{abp}}{1+\mathrm{bp}}{\mathrm{\γ}}_m...(2-3)$

通常吸附煤层气占80％～90％，而游离煤层气仅占10％～20％。

2)煤层气为理想气体，流动过程为等温过程，有如下煤层气状态方程：

$\mathrm{\ρ}=\frac{p}{\mathrm{RT}}...(2-4)$

式中：

T-温度；

R-煤层气气体常数。

3)一定范围内的煤层具有相同的原始煤层气压力；

4)煤层顶、底板与煤层相比含煤层气量少，因此可以认为顶底板不含煤层气；

5)煤层气在煤层中的流动符合达西定律，即：

$\left.\begin{array}{c}{q}_x=-\frac{K_x}{\mathrm{\μ}}\frac{\∂p}{\∂x}\\ q_y=-\frac{K_x}{\mathrm{\μ}}\frac{\∂p}{\∂y}\\ q_z=-\frac{K_x}{\mathrm{\μ}}\frac{\∂p}{\∂z}\end{array}\right\}...(2-5)$

式中：

K-煤层渗透率，m²；

μ-煤层气的绝对粘度，Pa·s；

对于煤体单元来说，煤层气含量的流入流出的质量流量等于煤体单元煤层气质量含量的变化量，即：

$\mathrm{div}\left(\mathrm{\ρ}\stackrel{\→}{q}\right)=-\frac{\∂M}{\∂t}...(2-6)$

式中：

ρ-在温度为293K，煤层气压力为p时的煤层气密度，(g/cm³)；

M-单位煤体所含的煤层气含量，(g/cm³)，M＝ρΓ；

t-时间，s。

在dt时间内煤体单元在相互垂直X、Y、Z方向流入流出的净质量流量为：

$-\frac{\∂\left(\mathrm{\ρ}{q}_x\right)}{\∂x}\mathrm{dxdydzdt},-\frac{\∂\left(\mathrm{\ρ}{q}_y\right)}{\∂y}\mathrm{dxdydzdt},-\frac{\∂\left(\mathrm{\ρ}{q}_z\right)}{\∂z}\mathrm{dxdydzdt}$

在dt时间内煤体单元流入流出的质量流量为：

$-[\frac{\∂\left(\mathrm{\ρ}{q}_x\right)}{\∂x}+\frac{\∂\left(\mathrm{\ρ}{q}_y\right)}{\∂y}+\frac{\∂\left(\mathrm{\ρ}{q}_z\right)}{\∂z}]\mathrm{dxdydzdt}$

在dt时间内煤体单元内煤层气质量变化量为：

$-\frac{\∂M}{\∂t}\mathrm{dxdydzdt}$

$-[\frac{\∂\left(\mathrm{\ρ}{q}_x\right)}{\∂x}+\frac{\∂\left(\mathrm{\ρ}{q}_y\right)}{\∂y}+\frac{\∂\left(\mathrm{\ρ}{q}_z\right)}{\∂z}]\mathrm{dxdydzdt}=\frac{\∂M}{\∂t}\mathrm{dxdydzdt}$

式中：

M-单元体ΔxΔyΔz内含有的煤层气质量，即：

$M=(\mathrm{\φ}\frac{p}{p_N}\mathrm{\ρ}+\frac{100-W-A}{100*(1+0.31W)}\·\frac{\mathrm{abp}}{1+\mathrm{bp}}{\mathrm{\γ}}_m\mathrm{\ρ})\mathrm{\Δx\Δy\Δz}$

将式(2-4)代入上式有：

$M=(\mathrm{\φ}\frac{p}{p_N}+\frac{100-W-A}{100*(1+0.31W)}\·\frac{\mathrm{abp}}{1+\mathrm{bp}}{\mathrm{\γ}}_m)\mathrm{\Δx\Δy\Δz}\frac{p}{\mathrm{RT}}$

$=(\frac{\mathrm{\φ}}{p_N}+\frac{100-W-A}{100*(1+0.31W)}\·\frac{\mathrm{ab}}{1+\mathrm{bp}}{\mathrm{\γ}}_m)\mathrm{\Δx\Δy\Δz}\frac{p^2}{\mathrm{RT}}$

因此有：

$\frac{\∂M}{\∂t}=\mathrm{\Δx\Δy\Δz}\frac{1}{\mathrm{RT}}[\frac{\mathrm{\φ}}{p_N}+\frac{100-W-A}{100*(1+0.31W)}\·{\mathrm{\γ}}_m(\frac{\mathrm{ab}}{2(1+\mathrm{bp})}+\frac{\mathrm{ab}}{2{(1+\mathrm{bp})}^2})]\frac{{\∂p}^2}{\∂t}...(2-7)$

而：

$\frac{\∂}{\∂x}\left(\mathrm{\ρ}{q}_x\right)=-\frac{\∂}{\∂x}\left(\frac{p}{\mathrm{RT}}\frac{K_x}{\mathrm{\μ}}\frac{\∂p}{\∂x}\right)=-\frac{1}{2\mathrm{\μRT}}\frac{\∂}{\∂x}\left(K_x\frac{{\∂p}^2}{\∂x}\right)...(2-8)$

同理有：

$\frac{\∂}{\∂x}\left(\mathrm{\ρ}{q}_y\right)=-\frac{\∂}{\∂y}\left(\frac{p}{\mathrm{RT}}\frac{K_y}{\mathrm{\μ}}\frac{\∂p}{\∂x}\right)=-\frac{1}{2\mathrm{\μRT}}\frac{\∂}{\∂y}\left(K_y\frac{{\∂p}^2}{\∂y}\right)...(2-9)$

$\frac{\∂}{\∂z}\left(\mathrm{\ρ}{q}_z\right)=-\frac{\∂}{\∂z}\left(\frac{p}{\mathrm{RT}}\frac{K_z}{\mathrm{\μ}}\frac{\∂p}{\∂z}\right)=-\frac{1}{2\mathrm{\μRT}}\frac{\∂}{\∂z}\left(K_z\frac{{\∂p}^2}{\∂z}\right)...(2-10)$

将式(2-7)、(2-8)、(2-9)、(2-10)代入方程(2-6)可得：

$\frac{\∂}{\∂x}\left(K_x\frac{{\∂p}^2}{\∂x}\right)+\frac{\∂}{\∂y}\left(K_y\frac{{\∂p}^2}{\∂y}\right)+\frac{\∂}{\∂z}\left(K_z\frac{{\∂p}^2}{\∂z}\right)$

$=[\frac{2\mathrm{\φ\μ}}{p_N}+\frac{100-W-A}{100*(1+0.31W)}\·\mathrm{\μ}\·{\mathrm{\γ}}_m[\frac{\mathrm{ab}}{(1+\mathrm{bp})}+\frac{\mathrm{ab}}{{(1+\mathrm{bp})}^2}\left]\right]\frac{{\∂p}^2}{\∂t}$

令：p²＝P，且令 $A\left(P\right)=\frac{2\mathrm{\φ\μ}}{p_N}+\frac{100-W-A}{100*(1+0.31W)}\·\mathrm{\μ}\·{\mathrm{\γ}}_m[\frac{\mathrm{ab}}{(1+b\sqrt{P})}+\frac{\mathrm{ab}}{{(1+b\sqrt{P})}^2}]$ 则有：

$\frac{\∂}{\∂x}\left(K_x\frac{\∂P}{\∂x}\right)+\frac{\∂}{\∂y}\left(K_y\frac{\∂P}{\∂y}\right)+\frac{\∂}{\∂z}\left(K_z\frac{\∂P}{\∂z}\right)=A\left(P\right)\frac{\∂P}{\∂t}...(2-11)$

式(2-11)即为抽采钻孔周围煤层气的渗流模型。

对于二维渗流微分方程：

$A\left(P\right)=\frac{\∂P}{\∂t}=[\frac{\∂}{\∂x}\left(K\frac{\∂P}{\∂x}\right)+\frac{\∂}{\∂y}\left(K\frac{\∂P}{\∂y}\right)]$

首先对二维空间D进行网格剖分，用Δx_i＝x_i-x_i-1表示沿x方向第i单元的步长；Δy_j＝y_j-y_j-1表示沿y方向第j个区间的步长；Δt为t方向的步长。根据上述差分代替微分的方法有：

$\frac{\∂p}{\∂t}\≈\frac{P_{i,j}^{(k+1)}-P_{i,j}^{\left(k\right)}}{\mathrm{\Δt}}$

$\frac{\∂}{\∂x}{\left(K\frac{\∂P}{\∂x}\right)}_{(i,j)}\≈[{\left(K\frac{\∂P}{\∂x}\right)}_{i+\frac{1}{2},j}-{\left(K\frac{\∂P}{\∂y}\right)}_{i-\frac{1}{2},j}]/({\mathrm{\Δx}}_i+{\mathrm{\Δx}}_{i+1})/2$

$\≈\frac{2}{{\mathrm{\Δx}}_i+{\mathrm{\Δx}}_{i+1}}(K_{i+\frac{1}{2},j}\frac{P_{i+1,j}-P_{i,j}}{\mathrm{\Δ}{x}_{i+1}}+K_{i-\frac{1}{2},j}\frac{P_{i-1,j}-P_{i,j}}{{\mathrm{\Δx}}_i})$

同理可得：

$\frac{\∂}{\∂y}{\left(K\frac{\∂P}{\∂y}\right)}_{(i,j)}\≈\frac{2}{{\mathrm{\Δy}}_i+{\mathrm{\Δy}}_{i+1}}(K_{i,j+\frac{1}{2}}\frac{P_{i,j+1}-P_{i,j}}{{\mathrm{\Δy}}_{i+1}}+K_{i,j-\frac{1}{2}}\frac{P_{i,j-1}-P_{i,j}}{\mathrm{\Δ}{y}_i})$

将上述几式代入：

$A\left(P_{i,j}\right)\frac{P_{i,j}^{(k+1)}-P_{i,j}^{\left(k\right)}}{\mathrm{\Δt}}=\frac{2}{{\mathrm{\Δx}}_i+{\mathrm{\Δx}}_{i+1}}(K_{i+\frac{1}{2},j}\frac{P_{i+1,j}-P_{i,j}}{{\mathrm{\Δx}}_{i+1}}+K_{i-\frac{1}{2},j}\frac{P_{i-1,j}-P_{i,j}}{{\mathrm{\Δx}}_i})+$

$\frac{2}{{\mathrm{\Δy}}_i+{\mathrm{\Δy}}_{i+1}}(K_{i,j+\frac{1}{2}}\frac{P_{i,j+1}-P_{i,j}}{{\mathrm{\Δy}}_{i+1}}+K_{i,j-\frac{1}{2}}\frac{P_{i,j-1}-P_{i,j}}{{\mathrm{\Δy}}_i})$

整理后可以得到：

$A\left(P_{i,j}\right)\frac{P_{i,j}^{(k+1)}-P_{i,j}^{\left(k\right)}}{\mathrm{\Δt}}=\frac{{2K}_{i-\frac{1}{2},j}}{({\mathrm{\Δx}}_i+{\mathrm{\Δx}}_{i+1}){\mathrm{\Δx}}_i}(P_{i-1,j}-P_{i,j})+\frac{{2K}_{i+\frac{1}{2},j}}{({\mathrm{\Δx}}_i+{\mathrm{\Δx}}_{i+1}){\mathrm{\Δx}}_{i+1}}(P_{i+1,j}-P_{i,j})+$

$\frac{2K_{i,j-\frac{1}{2}}}{({\mathrm{\Δy}}_i+{\mathrm{\Δy}}_{i+1}){\mathrm{\Δy}}_i}(P_{i,j-1}-P_{i,j})+\frac{{2K}_{i,j+\frac{1}{2}}}{({\mathrm{\Δy}}_i+{\mathrm{\Δy}}_{i+1}){\mathrm{\Δy}}_{i+1}}(P_{i,j+1}-P_{i,j})$

引入记号：

$A_{i,j}=\frac{{2K}_{i-\frac{1}{2},j}}{({\mathrm{\Δx}}_i+{\mathrm{\Δx}}_{i+1}){\mathrm{\Δx}}_i};$ $B_{i,j}=\frac{{2K}_{i+\frac{1}{2},j}}{({\mathrm{\Δx}}_i+{\mathrm{\Δx}}_{i+1}){\mathrm{\Δx}}_{i+1}}$

$C_{i,j}=\frac{{2K}_{i,j-\frac{1}{2}}}{({\mathrm{\Δy}}_i+{\mathrm{\Δy}}_{i+1}){\mathrm{\Δy}}_i};$ $D_{i,j}=\frac{{2K}_{i,j+\frac{1}{2}}}{({\mathrm{\Δy}}_i+{\mathrm{\Δy}}_{i+1}){\mathrm{\Δy}}_{i+1}}$

则原式可以写为：

$A\left(P_{i,j}\right)\frac{P_{i,j}^{(k+1)}-P_{i,j}^{\left(k\right)}}{\mathrm{\Δt}}=A_{i,j}(P_{i-1,j}-P_{i,j})+B_{i,j}(P_{i+1,j}-P_{i,j})+C_{i,j}(P_{i,j-1}-P_{i,j})+D_{i,j}(P_{i,j+1}-P_{i,j})$

将上式中由第K层分成两步计算，在第K层到第K+1层中间引进一个过渡层，即第K+1/2层，对

使用隐式，对

使用显式，得到：

$A\left(P_{i,j}\right)\frac{P_{i,j}^{(k+\frac{1}{2})}-P_{i,j}^{\left(k\right)}}{\mathrm{\Δt}/2}=A_{i,j}(P_{i-1,j}^{k+\frac{1}{2}}-P_{i,j}^{k+\frac{1}{2}})+B_{i,j}(P_{i+1,j}^{k+\frac{1}{2}}-P_{i,j}^{k+\frac{1}{2}})+C_{i,j}(P_{i,j-1}^k-P_{i,j}^k)+D_{i,j}(P_{i,j+1}^k-P_{i,j}^k)$

式中P_i，j ^k为K层已知值，

是过渡层第(K+1/2)层的待求值。求出第(K+1/2)层之后，从(K+1/2)层到K+1层，对

使用显式，对

使用隐式，得到：

$A\left(P_{i,j}\right)\frac{P_{i,j}^{(k+1)}-P_{i,j}^{(k+\frac{1}{2})}}{\mathrm{\Δt}/2}=A_{i,j}(P_{i-1,j}^{k+\frac{1}{2}}-P_{i,j}^{k+\frac{1}{2}})+B_{i,j}(P_{i+1,j}^{k+\frac{1}{2}}-P_{i,j}^{k+\frac{1}{2}})+C_{i,j}(P_{i,j-1}^{k+1}-P_{i,j}^{k+1})+D_{i,j}(P_{i,j+1}^{k+1}-P_{i,j}^{k+1})$

根据上式，即可用求得的

求出第K+1层的P_i，j ^k+1。

将上述两式整理可得交替方向隐式格式，即ADI格式：

$A_{i,j}{P}_{i-1,j}^{k+\frac{1}{2}}-(A_{i,j}+B_{i,j}+A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}})P_{i,j}^{k+\frac{1}{2}}+B_{i,j}{P}_{i+1,j}^{k+\frac{1}{2}}=-C_{i,j}{P}_{i,j-1}^k+(C_{i,j}+D_{i,j}-A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}})P_{i,j}^k-D_{i,j}{P}_{i,j+1}^k$

$-C_{i,j}{P}_{i,j-1}^{k+1}+(C_{i,j}+D_{i,j}+A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}})P_{i,j}^{k+1}-D_{i,j}{P}_{i,j+1}^{k+1}=A_{i,j}{P}_{i-1,j}^{k+\frac{1}{2}}-(A_{i,j}+B_{i,j}-A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}})P_{i,j}^{k+\frac{1}{2}}+B_{i,j}{P}_{i+1,j}^{k+\frac{1}{2}}$

以上两式均为三对角线性方程组，计算量和存储量均较小。

如果Δx_i＝x_i+1，Δy_i＝y_i+1，x、y方向的步长分别相等，则上两式的变为：

$\frac{K_{i-\frac{1}{2},j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i-1,j}^{k+\frac{1}{2}}-[\frac{K_{i-\frac{1}{2},j}}{{\left(\mathrm{\Δx}\right)}^2}+\frac{K_{i+\frac{1}{2},j}}{{\left(\mathrm{\Δx}\right)}^2}+A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δ}}]P_{i,j}^{k+\frac{1}{2}}+\frac{K_{i+\frac{1}{2},j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i+1,j}^{k+\frac{1}{2}}=$

$\frac{K_{i,j-\frac{1}{2}}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j-1}^k+[\frac{K_{i,j-\frac{1}{2}}}{{\left(\mathrm{\Δy}\right)}^2}+\frac{K_{i,j+\frac{1}{2}}}{{\left(\mathrm{\Δy}\right)}^2}-A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}]P_{i,j}^k-\frac{K_{i,j+\frac{1}{2}}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j+1}^k$

$\frac{K_{i,j-\frac{1}{2}}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j-1}^{k+1}-[\frac{K_{i,j-\frac{1}{2}}}{{\left(\mathrm{\Δy}\right)}^2}+\frac{K_{i,j+\frac{1}{2}}}{{\left(\mathrm{\Δy}\right)}^2}+A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}]P_{i,j}^{k+1}+\frac{K_{i,j+\frac{1}{2}}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j+1}^{k+1}=$

$-\frac{K_{i-\frac{1}{2},j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i-1,j}^{k+\frac{1}{2}}+[\frac{K_{i-\frac{1}{2},j}}{{\left(\mathrm{\Δx}\right)}^2}+\frac{K_{i+\frac{1}{2},j}}{{\left(\mathrm{\Δx}\right)}^2}-A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}]P_{i,j}^{k+\frac{1}{2}}-\frac{K_{i+\frac{1}{2},j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i+1,j}^{k+\frac{1}{2}}$

令： $K_{i,j}=K_{i-\frac{1}{2},j},$ $K_{i+1,j}=K_{i+\frac{1}{2},j},$ $K_{i,j}=K_{i,j-\frac{1}{2}},$ $K_{i,j+1}=K_{i,j+\frac{1}{2}}$

则上述两式变为：

$\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i-1,j}^{k+\frac{1}{2}}-[\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}+\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}+A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}]P_{i,j}^{k+\frac{1}{2}}+\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i+1,j}^{k+\frac{1}{2}}=$

$-\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j-1}^k+[\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}+\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}-A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}]P_{i,j}^k-\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j+1}^k...(3-1)$

$\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j-1}^{k+1}-[\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}+\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}+A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}]P_{i,j}^{k+1}+\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j+1}^{k+1}=$

$-\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i-1,j}^{k+\frac{1}{2}}+[\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}+\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}-A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}]P_{i,j}^{k+\frac{1}{2}}-\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i+1,j}^{k+\frac{1}{2}}...(3-2)$

令：

$A_{i,j}=\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}$ $B_{i,j}=-[\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}+\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}+A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}],$ $C_{i,j}=\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2},$

$D_{i,j}=-\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j-1}^k+[\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}+\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}-A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}]P_{i,j}^k-\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j+1}^k$

则式(4-1)变为：

$A_{i,j}{P}_{i-1,j}^{k+\frac{1}{2}}+B_{i,j}{P}_{i,j}^{k+\frac{1}{2}}+C_{i,j}{P}_{i+1,j}^{k+\frac{1}{2}}=D_{i,j},j=(\mathrm{1,2,3,4},,....m-1)...(3-3)$

同理令：

$A_{i,j}^{\′}=\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2},B_{i,j}^{\′}=-[\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}+\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}+A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}],$

$C_{i,j}^{\′}=\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2},$ $D_{i,j}^{\′}=-\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i-1,j}^{k+\frac{1}{2}}+[\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}+\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}-A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}]P_{i,j}^{k+\frac{1}{2}}-\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i+1,j}^{k+\frac{1}{2}}$

则式(4-2)变为：

$A_{i,j}^{\′}{P}_{i,j-1}^{k+1}+B_{i,j}^{\′}{P}_{i,j}^{k+1}+C_{i,j}^{\′}{P}_{i,j+1}^{k+1}=D_{i,j}^{\′},i=(\mathrm{1,2,3,4},....n-1)...(3-4)$

上述两式分别写成矩阵形式有：

$\left(\begin{array}{cccccc}{B}_{1,j}& C_{1,j}& ...& & & 0\\ A_{2,j}& B_{2,j}& C_{2,j}& ...& & 0\\ 0& A_{3,j}& B_{3,j}& C_{3,j}& ...& 0\\ ......& & & & & \\ 0& 0& 0& ...& A_{m-1,j}& B_{m-1,j}\end{array}\right)\left(\begin{array}{c}{P}_{1,j}\\ P_2{,}_j\\ P_{3,j}\\ .\\ .\\ .\\ P_{m-1,j}\end{array}\right)=\left(\begin{array}{c}{D}_{1,j}\\ D_{2,j}\\ D_{3,j}\\ .\\ .\\ .\\ D_{m-1,j}\end{array}\right)$

$\left(\begin{array}{cccccc}{B}_{i,1}^{\′}& C_{i,1}^{\′}& ...& & & 0\\ A_{i,2}^{\′}& B_{i,2}^{\′}& C_{i,2}^{\′}& ...& & 0\\ 0& A_{i,3}^{\′}& B_{i,3}^{\′}& C_{i,3}^{\′}& ...& 0\\ ......& & & & & \\ 0& 0& 0& ...& A_{i,n-1}^{\′}& B_{i,n-1}^{\′}\end{array}\right)\left(\begin{array}{c}{P}_{i,1}\\ P_{i,2}\\ P_{i,3}\\ .\\ .\\ .\\ P_{i,n-1}\end{array}\right)=\left(\begin{array}{c}{D}_{i,1}^{\′}\\ D_{i,2}^{\′}\\ D_{i,3}^{\′}\\ .\\ .\\ .\\ D_{i,n-1}^{\′}\end{array}\right)$

方程的初始条件和边界条件变为：

$P_{i,j}^{\left(0\right)}=P_0(i=\mathrm{1,2},...m-1;j=\mathrm{1,2},...,n-1)$

$P_{0,j}^{\left(k\right)}=P_{1,j}^{\left(k\right)},P_{m,j}^{\left(k\right)}=P_{m-1,j}^{\left(k\right)}(k=\mathrm{1,2},...,l)$

$P_{i,0}^{\left(k\right)}=0,P_{i,n}^{\left(k\right)}=0(k=\mathrm{1,2},...,l)$

上述两个矩阵均可用数值分析方法进行求解。

同理可以对三维渗流方程进行离散求解，不同之处在于三维方程最终离散为三个矩阵，且对于顺煤层钻孔和穿层钻孔应根据情况确定不同的边界条件，但同样是可以用数值方法进行求解的。

在矿井中不同情况下的煤层气流动形式的边界条件和初始条件见下表。但对于井下煤层气抽采钻孔来说，由于钻孔布置方式不同，因此其对应的定解条件也不同。

根据所建立的不同布置形式抽采钻孔的物理模型，结合达西定律，确定抽采钻孔的初始条件为：

t＝0，P＝P₀＝p₀ ²

顺层钻孔的边界条件：

x、y≤R₀，P＝P₁＝p₁ ²

x、y＝R，P＝P₀＝p₀ ²

$\frac{\∂P}{\∂x}=0,$ x＝0，x＝R

$\frac{\∂P}{\∂y}=0,$ $y=\±\frac{m}{2}$

$\frac{\∂P}{\∂z}=0,$ z＝0，z＝L

穿层钻孔的边界条件：

x、y≤R₀，P＝P₁＝p₁ ²

x、y＝R，P＝P₀＝p₀ ²

$\frac{\∂P}{\∂x}=0,$ x＝0，x＝R

$\frac{\∂P}{\∂y}=0,$ y＝0，y＝R

$\frac{\∂P}{\∂z}=0,$ z＝0，z＝L

式中：

P₀-煤层气原始压力(p₀)的平方，MPa2；

P₁-大气压或钻孔抽放负压的平方，MPa2；

R₀-钻孔半径，m；

R-未扰动煤层气边界的距离，即抽采半径，m；

m-煤层厚度，m；

L-钻孔长度，m。

以如上的表述为依据，本实施例的井下煤层气抽采产能预测方法参见图1，包括如下步骤：

1)测取待测煤层原始参数，包括如下表所示的原始参数：

2)通过求解如下模型，获得抽采后的平均渗透率、抽采后的煤层气压力、抽采后的煤层气平均含量以及煤层气抽采量：

$\frac{\∂}{\∂x}\left(K_x\frac{\∂P}{\∂x}\right)+\frac{\∂}{\∂y}\left(K_y\frac{\∂P}{\∂y}\right)+\frac{\∂}{\∂z}\left(K_z\frac{\∂P}{\∂z}\right)=A\left(P\right)\frac{\∂P}{\∂t};$

其中：

$K=K_0(1+{(1-\mathrm{\θlg}\frac{\mathrm{bp}+1}{{\mathrm{bp}}_0+1})}^2{e}^{3C_{\mathrm{\φ}}(p-p_0)});$

其中：P₀为初始煤层气压力，单位MPa；

P为煤层气压力，单位MPa；

b为煤的吸附常数，单位1/MPa；

θ为煤基质收缩系数；

K₀为初始煤层气渗透率；

$A\left(P\right)=\frac{2\mathrm{\φ\μ}}{p_N}+\frac{100-W-A}{100*(1+0.31W)}\·\mathrm{\μ}\·{\mathrm{\γ}}_m[\frac{\mathrm{ab}}{(1+b\sqrt{P})}+\frac{\mathrm{ab}}{{(1+b\sqrt{P})}^2}]$

K-煤层渗透率，m²；

μ-煤层气的绝对粘度，Pa·s；

a为煤的最大煤层吸附量，单位m³/t

b为煤的吸附常数，单位1/MPa；

W为煤的水分，单位％；

A为煤的灰分，单位％；

P为煤层气压力，单位MPa；

γ_m为煤的容重，单位t/m³；

φ为煤层孔隙率；

p_N为1标准大气压，0.1013MPa。

所述步骤2)的求解过程具体包括如下步骤：

21)将待测煤层原始参数作为初始值计算抽采范围，对抽采范围进行网格划分；参见图2，网格划分是用一组有限个离散的点来代替原有的连续空间，即是从原点沿x方向用Δx_i＝x_i-x_i-1(表示沿x方向第i单元的步长)划分；从原点沿y方向用Δy_j＝y_j-y_j-1(表示第j个区间的步长)划分；Δt为t方向的步长进行时间的划分。

22)根据步骤21)的网格划分，计算系数矩阵并建立求解矩阵：

系数矩阵一：

$A_{i,j}=\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}$ $B_{i,j}=-[\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}+\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}+A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}],$ $C_{i,j}=\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2},$

$D_{i,j}=-\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j-1}^k+[\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}+\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}-A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}]P_{i,j}^k-\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j+1}^k$

$\left(\begin{array}{cccccc}{B}_{1,j}& C_{1,j}& ...& & & 0\\ A_{2,j}& B_{2,j}& C_{2,j}& ...& & 0\\ 0& A_{3,j}& B_{3,j}& C_{3,j}& ...& 0\\ ......& & & & & \\ 0& 0& 0& ...& A_{m-1,j}& B_{m-1,j}\end{array}\right)$

系数矩阵二：

$A_{i,j}^{\′}=\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2},$ $B_{i,j}^{\′}=-[\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}+\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}+A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}],$

$C_{i,j}^{\′}=\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2},$ $D_{i,j}^{\′}=-\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i-1,j}^{k+\frac{1}{2}}+[\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}+\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}-A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}]P_{i,j}^{k+\frac{1}{2}}-\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i+1,j}^{k+\frac{1}{2}}$

$\left(\begin{array}{cccccc}{B}_{i,1}^{\′}& C_{i,1}^{\′}& ...& & & 0\\ A_{i,2}^{\′}& B_{i,2}^{\′}& C_{i,2}^{\′}& ...& & 0\\ 0& A_{i,3}^{\′}& B_{i,3}^{\′}& C_{i,3}^{\′}& ...& 0\\ ......& & & & & \\ 0& 0& 0& ...& A_{i,n-1}^{\′}& B_{i,n-1}^{\′}\end{array}\right)$

求解矩阵：

$\left(\begin{array}{cccccc}{B}_{1,j}& C_{1,j}& ...& & & 0\\ A_{2,j}& B_{2,j}& C_{2,j}& ...& & 0\\ 0& A_{3,j}& B_{3,j}& C_{3,j}& ...& 0\\ ......& & & & & \\ 0& 0& 0& ...& A_{m-1,j}& B_{m-1,j}\end{array}\right)\left(\begin{array}{c}{P}_{1,j}\\ P_2{,}_j\\ P_{3,j}\\ .\\ .\\ .\\ P_{m-1,j}\end{array}\right)=\left(\begin{array}{c}{D}_{1,j}\\ D_{2,j}\\ D_{3,j}\\ .\\ .\\ .\\ D_{m-1,j}\end{array}\right)$

$\left(\begin{array}{cccccc}{B}_{i,1}^{\′}& C_{i,1}^{\′}& ...& & & 0\\ A_{i,2}^{\′}& B_{i,2}^{\′}& C_{i,2}^{\′}& ...& & 0\\ 0& A_{i,3}^{\′}& B_{i,3}^{\′}& C_{i,3}^{\′}& ...& 0\\ ......& & & & & \\ 0& 0& 0& ...& A_{i,n-1}^{\′}& B_{i,n-1}^{\′}\end{array}\right)\left(\begin{array}{c}{P}_{i,1}\\ P_{i,2}\\ P_{i,3}\\ .\\ .\\ .\\ P_{i,n-1}\end{array}\right)=\left(\begin{array}{c}{D}_{i,1}^{\′}\\ D_{i,2}^{\′}\\ D_{i,3}^{\′}\\ .\\ .\\ .\\ D_{i,n-1}^{\′}\end{array}\right)$

23)计算第一个时间步长的各结点煤层气压力和渗透率；

24)将计算所得的各结点煤层气压力和渗透率赋为初值；

25)利用步骤24)所得的初值重新计算抽采范围并进行网格划分；

26)根据步骤25)的网格划分，计算系数矩阵并建立求解矩阵；

27)计算下一时间步长的各结点煤层气压力和渗透率，并返回步骤24)，直到计算完成，计算完成的条件为

以某矿为例，抽采原始参数如下表所示：

采用上述方法编制模拟软件，对煤层厚度是6m、煤层透气性系数分别为40m2/MPa2.d、60m2/MPa2.d(对应的渗透率为1md、1.5md)的顺层单孔和顺层平行钻孔的抽采产能进行了模拟分析，结果如图2、3所示，从图2、3可以看出，煤层气单孔抽采量随时间的增长，流量在不断衰减，而单孔抽采总量在不断增加，且呈线性增长；透气性系数(渗透率)大的煤层，其单孔初始流量较大。

顺层平行钻孔抽采产能数值模拟结果见图4、5所示。由图4、5可以看出，随着抽采时间的不断增加，单孔平均百米抽采总量在不断增加，而单孔百米平均抽采量在不断衰减。

将图4、5与图2、3相比发现，顺层平行钻孔单孔抽采量的变化规律与顺层单孔抽采量变化规律不同，这两种不同透气性系数(渗透率)煤层的抽采钻孔其抽采量前后衰减规律不同，且均有“拐点”出现，将该点处的时间、抽采半径对比分析发现，该点即为钻孔间距的中点，可见模拟结果与前述分析是一致的。

通过对比还可以发现，在“拐点”前单孔衰减规律和平行钻孔衰减规律相同，“拐点”后平行钻孔的衰减速度较单孔衰减的快，分析认为由于平行钻孔两两“截抽”，对于平行钻孔中的某一个钻孔来说，其抽采区域是一定的，同时根据数值模拟结果发现，同样的抽采时间内，单孔抽采范围的透气性系数(渗透率)和平行钻孔单孔抽采影响范围内的平均透气性系数(渗透率)均有所增加，但平行钻孔单孔抽采影响范围内煤层的平均透气性系数(渗透率)增长的幅度较单个钻孔抽采影响范围内煤层平均透气性系数(渗透率)增长幅度大，导致煤层气抽采量下降的速度较快。

采用上述方法编制模拟软件，对煤层厚度是6m、煤层透气性系数分别为4m2/MPa2.d、10m2/MPa2.d(对应的渗透率为0.1md、0.25md)的穿层单孔和穿层平行钻孔的抽采产能进行了模拟分析。

穿层单孔抽采产能数值模拟结果见图6、7。由图6、7可知，穿层单孔百米抽采量随抽采时间的增加，抽采量也在不断增加，这与顺层单孔的预测结果不同，这是由于煤层气的“补给”范围不同造成的，穿层钻孔的煤层气“补给”范围随抽采时间的增加在X、Y方向不受限制，而顺层钻孔煤层气的“补给”范围随时间的增加在Y方向受到煤层顶、底板的限制。

穿层平行钻孔抽采产能数值模拟结果见图8、9。由图8、9可知，穿层平行钻孔与顺层平行钻孔的产能预测曲线类似，均有“拐点”，且拐点前后穿层平行钻孔的产能曲线与顺层平行钻孔的产能曲线不同，这是由于两者的煤层气“补给”方式不同，但拐点出现的原因是相同的，均是由于周围钻孔的“截抽”造成的。

以上所述仅为本发明的优选实施例，并不用于限制本发明，显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (2)

1.井下煤层气抽采产能预测方法，其特征在于：包括如下步骤：

1)测取待测煤层原始参数，包括初始煤层气压力、煤的吸附常数、初始煤层气渗透率、煤基质收缩系数、煤的容重、孔隙率、孔隙压缩系数、灰分和水分；

2)将测得的煤层原始参数代入模型中进行求解，获得抽采后的平均渗透率、抽采后的煤层气压力、抽采后的煤层气平均含量以及煤层气抽采量：

$\frac{\∂}{\∂x}\left(K_x\frac{\∂P}{\∂x}\right)+\frac{\∂}{\∂y}\left(K_y\frac{\∂P}{\∂y}\right)+\frac{\∂}{\∂z}\left(K_z\frac{\∂P}{\∂z}\right)=A\left(P\right)\frac{\∂P}{\∂t};$

其中：

$K=K_0(1+{(1-\mathrm{Rlg}\frac{\mathrm{bp}+1}{{\mathrm{bp}}_0+1})}^2-e^{3C_{\mathrm{\φ}}(p-p_0)});$

其中：p₀为初始煤层气压力，单位MPa；

p为煤层气压力，单位MPa；

b为煤的吸附常数，单位1/MPa；

R为煤基质收缩系数；

K₀为初始煤层气渗透率；

$A\left(P\right)=\frac{2\mathrm{\φ\μ}}{p_N}+\frac{100-W-A}{100*(1+0.31W)}\·\mathrm{\μ}\·{\mathrm{\γ}}_m[\frac{\mathrm{ab}}{(1+b\sqrt{P})}+\frac{\mathrm{ab}}{{(1+b\sqrt{P})}^2}]$

K为煤层渗透率，m²；

μ为煤层气的绝对粘度，Pa·s；

a为煤的最大煤层吸附量，单位m³/t

b为煤的吸附常数，单位1/MPa；

W为煤的水分，单位％；

A为煤的灰分，单位％；

P为煤层气压力，单位MPa；

γ_m为煤的容重，单位t/m³；

φ为煤层孔隙率；

c_φ为孔隙压缩系数；

p_N为1标准大气压，0.1013MPa；

所述，抽采后的平均渗透率为抽采一定时间后，抽采影响范围内的煤体的渗透率；所述抽采后的煤层气压力为抽采一定时间后，抽采影响范围内的煤体内残余的煤层气的压力；所述，抽采后煤层气含量为抽采一定时间后，抽采影响范围内的煤体内残余的煤层气含量；所述煤层气抽采量为一定抽采时间段内，抽采钻孔所抽出的煤层气量。

2.如权利要求1所述的井下煤层气抽采产能预测方法，其特征在于：所述步骤2)的求解过程具体包括如下步骤：

21)将待测煤层原始参数作为初始值计算抽采范围，对抽采范围进行网格划分；

22)根据步骤21)的网格划分，计算系数矩阵并建立求解矩阵：

系数矩阵一：

$A_{i,j}=\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}$ $B_{i,j}=-[\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}+\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}+A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}],$ $C_{i,j}=\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2},$

$D_{i,j}=-\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j-1}^k+[\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}+\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}-A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}]P_{i,j}^k-\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}{P}_{i,j+1}^k$

$\left(\begin{array}{cccccc}{B}_{1,j}& C_{1,j}& ...& & & 0\\ A_{2,j}& B_{2,j}& C_{2,j}& ...& & 0\\ 0& A_{3,j}& B_{3,j}& C_{3,j}& ...& 0\\ ......& & & & & \\ 0& 0& 0& ...& A_{m-1,j}& B_{m-1,j}\end{array}\right)$

系数矩阵二：

$A_{i,j}^{\′}=\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2},$ $B_{i,j}^{\′}=-[\frac{K_{i,j-1}}{{\left(\mathrm{\Δy}\right)}^2}+\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2}+A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}],$

$C_{i,j}^{\′}=\frac{K_{i,j+1}}{{\left(\mathrm{\Δy}\right)}^2},$ $D_{i,j}^{\′}=-\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i-1,j}^{k+\frac{1}{2}}+[\frac{K_{i-1,j}}{{\left(\mathrm{\Δx}\right)}^2}+\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}-A\left(P_{i,j}\right)\frac{2}{\mathrm{\Δt}}]P_{i,j}^{k+\frac{1}{2}}-\frac{K_{i+1,j}}{{\left(\mathrm{\Δx}\right)}^2}{P}_{i+1,j}^{k+\frac{1}{2}}$

$\left(\begin{array}{cccccc}{B}_{i,1}^{\′}& C_{i,1}^{\′}& ...& & & 0\\ A_{i,2}^{\′}& B_{i,2}^{\′}& C_{i,2}^{\′}& ...& & 0\\ 0& A_{i,3}^{\′}& B_{i,3}^{\′}& C_{i,3}^{\′}& ...& 0\\ ......& & & & & \\ 0& 0& 0& ...& A_{i,n-1}^{\′}& B_{i,n-1}^{\′}\end{array}\right)$

求解矩阵：

$\left(\begin{array}{cccccc}{B}_{1,j}& C_{1,j}& ...& & & 0\\ A_{2,j}& B_{2,j}& C_{2,j}& ...& & 0\\ 0& A_{3,j}& B_{3,j}& C_{3,j}& ...& 0\\ ......& & & & & \\ 0& 0& 0& ...& A_{m-1,j}& B_{m-1,j}\end{array}\right)\left(\begin{array}{c}{P}_{1,j}\\ P_{2,j}\\ P_{3,j}\\ .\\ .\\ .\\ P_{m-1,j}\end{array}\right)=\left(\begin{array}{c}{D}_{1,j}\\ D_{2,j}\\ D_{3,j}\\ .\\ .\\ .\\ D_{m-1,j}\end{array}\right)$

$\left(\begin{array}{cccccc}{B}_{i,1}^{\′}& C_{i,1}^{\′}& ...& & & 0\\ A_{i,2}^{\′}& B_{i,2}^{\′}& C_{i,2}^{\′}& ...& & 0\\ 0& A_{i,3}^{\′}& B_{i,3}^{\′}& C_{i,3}^{\′}& ...& 0\\ ......& & & & & \\ 0& 0& 0& ...& A_{i,n-1}^{\′}& B_{i,n-1}^{\′}\end{array}\right)\left(\begin{array}{c}{P}_{i,1}\\ P_{i,2}\\ P_{i,3}\\ .\\ .\\ .\\ P_{i,n-1}\end{array}\right)=\left(\begin{array}{c}{D}_{i,1}^{\′}\\ D_{i,2}^{\′}\\ D_{i,3}^{\′}\\ .\\ .\\ .\\ D_{i,n-1}^{\′}\end{array}\right))$

A，B，C，D均代指系数矩阵；

Δx表示沿x方向的步长；

Δt表示t方向的步长；

23)计算第一个时间步长的各结点煤层气压力和渗透率；

24)将计算所得的各结点煤层气压力和渗透率赋为初值；

25)利用步骤24)所得的初值重新计算抽采范围并进行网格划分；

26)根据步骤25)的网格划分，计算系数矩阵并建立求解矩阵；

27)计算下一时间步长的各结点煤层气压力和渗透率，并返回步骤24)，直到计算完成，计算完成的条件为

Images