当前位置:首页 >> 电力/水利 >>

煤岩体瓦斯、水渗流耦合过程数值模型及其在矿山工程中的应用


煤岩体瓦斯、水渗流耦合过程数值模型 及其在矿山工程中的应用

东北大学

杨天鸿

3 主要物理过程
气体压缩过程、气体吸附和解析过程、扩散过 程、渗流过程、应力-渗流耦合过程等

4 物理数学方程
?m (1)瓦斯渗流方程: + ? ? (ρ g q g ) = Q p ?t

kg qg = ? ? (?p + ρ g gD ) μg

(2)气体状态方程:

ρ g = βp
(3) Langmuir吸附 解析方程

β = M g (RT )
?φ a1 a 2 ρ s ? 2 ?p m = β? + ? p 1+ a p ? 2 ? 0 ?

4 物理数学方程
(4)气体压力有效应力方程:

2Gν σ ij = 2Gε ij + ε kk δ ij ? αpδ ij 1 ? 2v

4 物理数学方程
(5)Klinkenberg 方程: (滑脱效应) (6)渗透耦合方程:

qg = ?

kg

μg

? (?p + ρ g gD )

? ? α k k ∞0.36 ? ? k g = k ∞ ?1 + ? p ? ? ?

k ∞ = k ∞ 0 ? exp[22.2(φ φ 0 ? 1)]

φ = (φ 0 ? φ r ) exp(α φ ? σ v ) + φ r

瓦斯渗流-应力耦合方程:
2 ? φ a1a2 ρ s ? ?p 2 a1a2 ρ s ? ? p? β? + 2 ? p0 1 + a2 p 2(1 + a2 p ) ? ?t ? ? k ∞ ? α k k ∞ 0.36 ? 2 ? ?1 + ??p ? = Q p ? ? ? ?β ? ? p ? ? μg ? ? ? ?

Gu i , jj

G + u j , ji ? αp,i + Fi = 0 1 ? 2ν

求解方法 – FEMLAB简介
FEMLAB为专门求解耦合偏微分方程组的有限元分析工具。 具有强大的处理功能:⑴ 它含有一些内嵌的经典物理模型, 包括单物理场和多物理场模型,可以直接用于分析。 ⑵ 功能最强大、最灵活的还是其偏微分方程组模式:系数形 式、通式与弱形式。这三个数学应用模式中:系数形式 (Coefficient form),适宜求解线性问题;通式(General form),适宜求解非线性问题;弱形式(Weak form),最为灵 活,对于边界条件、时间序列复杂模型尤为适宜,但应用也 相对复杂些。一般地,大多数物理问题均可采用通式模式进 行求解。

求解方法 – FEMLAB简介

⑶ 对于不同物理场中交叉耦合项的处理简单有效。一 方面,在各物理场的偏微分方程中考虑了不同场的影 响;另一方面,各物理场中的计算变量可以直接用于 耦合关系的定义。 ⑷ 该软件带有Script语言并兼容Matlab语言,具有强 大的二次开发功能,对于创新性理论研究尤为适合。 此外,FEMLAB还有强大的后处理功能。

一、基于数字图像处理技术的煤层瓦斯渗 流过程数值模拟
在数字图像处理技术中,人们通常采用HSI颜 色空间来表述数字图像,因为该图像空间有利 于人肉眼的识别。HSI色彩空间中,颜色用色 度(Hue)、饱和度(Saturation)和亮度 (Intensity)来表示。其中H表示了肉眼看到 的颜色,S表示该颜色相对于白色的饱和度,I 表示的亮度。HSI的颜色空间的数值可以从 RGB的数据转换而来。其中I的数值是R、G和 B的算术平均值。

(1)煤的细观数字图像

(2) 图像I值的分布

煤样的细观扫描照片及基于数值图像技术获得的孔隙率分布

图1(a)是文献[18]给出的具有突出倾向性煤 样的显微照片,从中可以看出煤样中的叶状的 碎斑结构。灰度低的部分可视为裂隙。图1(b) 给出了该图像I值的分布。由此可以看出I值较 好地反映了煤体中的结构特征。

煤样中各组分中的孔隙率、渗透率和初始瓦斯压力

组分 像素的I值 0.0 ≤ I < 0.6 0.6 ≤ I ≤ 1.0 裂隙带 基质

孔隙率 0.2 0.05

弹性模量 / GPa 2.0 20.0

渗透率 初始瓦斯压力 / MPa / m2 1.0e-16 1.0e-18 3.0 1.0

(a)初始孔隙率分布 (b)弹性模量分布 基于数值图像技术获得的孔隙率和渗透率分布

σ0y = 1 MPa
n?(?p) = 0 内部边界: pi = 0.1 MPa

E (7.5,12.5)

Y C (12.5,10) A(14,10) X D(12.5,7.5)

σ0x = 2 MPa
n?(?p) = 0
B (6,10)

σ0x = 2 MPa
n?(?p) = 0

初始条件: p0

n?(?p) = 0

瓦斯运移的数值模型

瓦斯压力图
下图给出的瓦斯压力图可以看出瓦斯从中间井孔不断 释放的整个过程。由于这里没有考虑煤层的补给,故 瓦斯压力不断降低,最终煤样中的瓦斯压力降低到到 内部孔边界的瓦斯压力,故瓦斯的运移过程停止。在 初始条件下,由于裂隙带中的瓦斯压力比煤基质中的 高,故瓦斯由裂隙带向外围的煤机制中不断扩散,使 得裂隙带和基质间的瓦斯压力剃度不断降低。直到时 间 t = 1e04 s后,瓦斯开始集中向抽放孔中流动。最 终在t = 1e06 s左右时,煤样内的瓦斯压力和抽放孔中 的给定压力相同,瓦斯流动过程停止。

t = 1e-01 s t = 1e0 s t = 1e01 s t = 1e02 s t = 1e03 s t = 1e04 s t = 1e05 s t = 1e 06 s t = 1e 07 s

瓦斯压力及渗流速度随时间的变化过程图

下图给出了试样中5个点(其位置见瓦斯运移的数值 模型 )的瓦斯压力——时间曲线。其中,点C和D位 于裂隙带中,具有较大的初始瓦斯压力(为3.0 MPa)。而点A、B和E位于基质中,其初始瓦斯压力 为1.0 MPa。在t = 10 s前,随着瓦斯的抽放,裂隙带 中的瓦斯压力逐渐降低,而煤基质中的瓦斯压力不断 增加,这说明瓦斯从裂隙带向基质中不断渗流。当t = 10 s时,试样中的瓦斯开始集中往抽放孔处流动,最 终在t = 1e5 s时,试样中的瓦斯压力达到抽放孔的压 力值0.1 MPa。

3

x 10

6

2.5

2 p (Pa)

A (0.014,0.01) B (0.006,0.01) C (0.0125,0.01) D (0.0125,0.0075) E (0.0075,0.0125)

1.5

1

0.5

0 -2 10

10

0

10

2

10 Time (s)

4

10

6

10

8

试样中五个点(A、B、C、D和E)处的瓦斯压力—时间曲线

下图给出了5个特征点的孔隙率的变化曲线。 与瓦斯压力的分布曲线类似,在裂隙带的点C 和D处,受到外部边界应力的作用后,孔隙率 从初值0.2降低到了0.1855和0.188。随后,随 着瓦斯的不断释放,瓦斯压力降低,故有效应 力增加,所以,孔隙被压缩,所以孔隙率随着 瓦斯压力的降低而不断降低。

0.19

0.185

C (0.0125,0.01) D (0.0125,0.0075)

0.18 phi

0.175

0.17

0.165

0.16 -2 10

10

0

10

2

10 Time (s)

4

10

6

10

8

瓦斯渗流过程中的煤样的孔隙率变化(断裂带中的点C和D)

但是,如下图所示,在煤基质中,t = 1e-02s 时,首先是由于煤的变形导致了孔隙率下降, 但是,随着瓦斯从裂隙带向着煤基质中渗流, 故基质中的瓦斯压力会逐渐,所以导致了孔隙 率的增加。此后,由于瓦斯开始集中向抽放孔 中流动,整个煤样中的瓦斯压力不断降低会引 起有效应力的增加,故孔隙率会不断降低。

0.038 0.0375 0.037 0.0365 0.036 phi 0.0355 0.035 0.0345 0.034 0.0335 -2 10
0

A (0.014,0.01) B (0.006,0.01) E (0.0075,0.0125) 10 10
2

10 Time (s)

4

10

6

10

8

瓦斯渗流过程中的煤样的孔隙率变化(煤基质中的点A、B和E)

主要结论如下:
(1)由于裂隙带中的瓦斯压力比煤基质中的高,故 瓦斯由裂隙带向外围的煤机制中不断扩散,使得裂隙 带和基质间的瓦斯压力剃度不断降低。在某个时间 后,瓦斯开始集中向抽放孔中流动。 (2)在本文给定的模拟参数下,应力场引起的煤体 压缩是对其渗透性的影响要大于Klinkerberg效应和瓦 斯解吸效应所引起的渗透率变化。故总体上煤层中表 现出渗透率的降低,随着瓦斯的不断抽放,渗透率更 是不断降低。

二、冒落区瓦斯浓度扩散-对流及 风流场数值模拟

瓦斯浓度扩散-对流场基本方程:

? ?C ?C ( Dij + vi C ) = ?xi ?x j ?t
Darcy方程 Darcy

qg = ?

kg

μg

? (?p + ρ g gD )

冒落区破碎岩体压实过程中的气体流场
非Darcy方程(堆石体、土石坝流场) Darcy

冒落区瓦斯浓度扩散-对流 及风流场数值模拟
通风流场基本方程: Navier-Stokes equations: Free flow适合通风 巷道风流场

ρv ? ?v = ? ? (? pI + η (?v + (?v ) ) ) + F
T

式中:v——流体流速,m/s;p——流体压力,Pa;ρ—— 流体密度,kg/ m3;I——单位矢量;F——流体阻力

冒落区瓦斯浓度扩散-对流 及风流场数值模拟
通风流场基本方程: Brinkman equations: Fast flow in porous media,适 合冒落区风流场

(? ? ? η(?u + (?u ) ))-? η u + ?p ? F ?=0 ? ? k
T

?

?

? ?u = 0

模型建立
计算模型和方案 参照综采工作面的具体尺寸,建立如下图所示 的二维计算模型,不考虑势能,模型东西向即 工作面推进方向取400m,南北向宽240m。划 分成均质(5个分区,沿工作向采空区方向的 区域(1)至(5)宽度均取为80m);即大约 1周的推进时间瞬时完成。假设这个期间内每 个区域的透气系数分布均不同,具体采用的计 算模型相关条件如下:

模型建立
浓度补给边界

出风口

出风口

出风口

出风口

出风口

不 透 气 边 界

采空冒落区
(5) (4) (3) (2) (1)

不 透 气 边 界

进风口

进风口

进风口

进风口

进风口

浓度补给边界 推进方向

模型建立
边界条件: (1)通风条件:左侧下20m为进风口边界, 左侧上20m为回风口边界,压力差100Pa,其他 边界为不透气边界。 (2)扩散条件:右侧边界为绝缘对称边界, 其他边界有补给,推进区域1时,上下边界补 给量为1.2e-6mol/m2·s,汇源项为3e-6 mol/m2·s,以描述瓦斯通量随开采动态过程而 增加的瓦斯量。

模型建立
初始条件:

域内具有一个大气压,瓦斯初始浓度 3mol/m3;推进新的工作面区域时,旧区域瓦 斯浓度模拟结果C0再附加一初始补充浓度Cbc (取3mol/m3)作为旧推进工作面的初始浓度 。

δt

模型建立
时间步长: 按照非均等的积数步长增大,初始值7s,终止 指7e5s,设定100个中间时间值。 计算参数: 动粘系数=1.8e-5pa·s,流体密度 =1200kg/m3,扩散系数D=2e-5m2 ,瞬态时间 比例系数=0.55

模型建立
随推进进行透气率变化表
透气率 已推区域

1
3.24e-8

2

3

4

5

(m2)

新推区域

1 2 3 4 5

2.88e-8 3.24e-8 2.52e-8 2.88e-8 3.24e-8 2.16e-8 2.52e-8 2.88e-8 3.24e-8 1.8e-8 2.16e-8 2.52e-8 2.88e-8 3.24e-8

模拟结果分析
推进区域1时的计算结果

下图分别是当时间为:7s,7e3s,7e4s, 1.1e5s,2.5e5s,3.5e5s,7e5s瞬态时间区 域1采空区瓦斯浓度分布图

模拟结果分析

Time=7s时瓦斯浓度

模拟结果分析

Time=7e3s时瓦斯浓度

模拟结果分析

Time=7e4s时瓦斯浓度

模拟结果分析

Time=1.1e5s时瓦斯浓度

模拟结果分析

Time=2.5e5s时瓦斯浓度

模拟结果分析

Time=3.5e5s时瓦斯浓度

模拟结果分析

A1

A1 ’

Time=7e5s时瓦斯浓度和流线

模拟结果分析
6 5 瓦斯浓度/mol/m3 4 3 2 1 0 0 20 40 推进方向/m 60 80
Time=7s Time=7e3s Time=7e4s Time=1.1e5s Time=2.5e5s Time=3.5e5s Time=7e5s

第一步推进时A1-A1’切面浓度变化曲线

模拟结果分析
推进区域2时的计算结果

下图分别是当时间为:7s,7e3s,1.1e5s, 7e5s瞬态时间区域2采空区瓦斯浓度分布图

模拟结果分析

Time=7s时瓦斯浓度

模拟结果分析

Time=7e3s时瓦斯浓度

模拟结果分析

Time=1.1e5s时瓦斯浓度

模拟结果分析

A2

A2 ’

第二步推进时Time=7e5s时瓦斯浓度和流线

模拟结果分析
9 8 7 6 5 4 3 2 1 0 0

瓦斯浓度/mol/m3

Time=7s Time=7e3s Time=7e4s Time=1.1e4s Time=2.5e4s Time=3.5e4s Time7e5s

20

40

60

80 100 推进方向/m

120

140

160

第二步推进时A2-A2’切面浓度变化曲线

模拟结果分析
推进区域3,4时的计算结果

模拟结果分析

A3

A3 ’

第三步推进时Time=7e5s时瓦斯浓度和流线

模拟结果分析
14 12 瓦斯浓度/mol/m3 10 8 6 4 2 0 0 40 80 120 推进方向/m 160 200 240
Time=7s Time=7e3s Time=7e4s Time=1.1e5s Time=2.5e5s Time=3.5e5s Time=7e5s

第三步推进时A3-A3’切面浓度变化曲线

模拟结果分析

A4

A4 ’

第四步推进时Time=7e5s时瓦斯浓度和流线

模拟结果分析
20 18 16 14 12 10 8 6 4 2 0 0

瓦斯浓度/mol/m3

Time=7s Time=7e3s Time=7e4s Time=1.1e5s Time=2.5e5s Time=3.5e5s Time=7e5s

50

100

150 200 推进方向/m

250

300

第四步推进时A4-A4’切面浓度变化曲线

讨论

(1)在区域1右侧不透气边界附近,随着工作 面向左推进,风流作用逐渐减弱,该处瓦斯浓 度不断集聚增大,每个推进步瓦斯浓度分别为 7mol,12mol,17 mol,24 mol,表明该处距 离通风口越远,瓦斯集聚程度越大。

讨论

(2)对于每个推进步,当时间达到7e5s(8d 左右)时,风流基本上可以把本区域瓦斯浓度 降低到安全范围内(接近0mol),表明风流 在本区域的流速明显大于前一个区域。而前一 个区域的下左半部分(占整个区域的1/4)的 瓦斯浓度也能得到有效降低,而对其他区域瓦 斯浓度驱散作用影响有限。

讨论

(3)图9可见,瓦斯浓度沿工作面推进方向从 近到远瓦斯浓度成台阶状渐次增大,这表明通 风量只能在一定范围内降低瓦斯浓度,同时也 看出距离通风口越远,瓦斯补给的时间和补给 量越大。

讨论

(4)虽然模型上下边界都是瓦斯补给边界, 但由于模型左下边界是进风口,瓦斯随风流对 流作用得到有效降低;而模型左上边界是出风 口,瓦斯容易在回风隅角大量聚集,所以该处 瓦斯浓度较高。

三、三维渗流耦合模型及瓦斯抽放
根据实际的三维煤层瓦斯抽放过程,建立如图 1所示的长宽高为10m×10m×10m理想化的 三维计算模型,模型左下部边界假设为巷道, 三个瓦斯抽放孔——K1、K 2、K 3之间的距 离1m,与水平面呈45°倾角展布,K 1、K 3 与K 2呈16°夹角分布在孔2两侧,三个瓦斯 抽放孔长度7m,瓦斯抽放孔按照“以缝代 孔”[16]原则简化为定压力边界(25kPa)。

外部载荷

K1 K2 K3
16?

L2

45?

L1

巷道所在位置

(1) 边界条件:四周和上下面均为零通量不透 气边界;四周和下部都约束法线方向的位移, 上部自由,作用有P分别为10MPa、20MPa、 2MPa、0.1MPa的外部载荷,同时模型具有自 重载荷。 (2) 初始条件:内部有1MPa的初始瓦斯压力, 三个抽放孔的压力为0.25e5Pa。

(3) 时间步长:按照非均等的积数步长增大, 初始值1s,终止指1e7s(10d左右),设定 100个中间时间步长。 (4) 计算方案:模拟不同外部载荷条件下 (2MPa、10MPa、20MPa、0.1MPa),瓦 斯抽放效果和渗透性变化规律。 (5) 计算参数:相关参数列表1所列。

杨氏模量 E (Pa)

泊松比 ν

煤层密度 ρs (kg/m3)

流体密度 ρ1 (kg/m3)

流体粘滞性 饱和渗透率 (Pa·s) (m2)

7.41e9

0.33

1250

1000

1e-3

5e-6

模拟结果分析
推进区域2时的计算结果

下图分别是当时间为:1s, 1.28e5s , 1.08e6s,1e7s瞬态时间渗透性系数分布图 和压力等表面图

渗透性系数分布图

Time=1s时渗透性系数分布图

渗透性系数分布图

Time=1.28e5s时渗透性系数分布图

渗透性系数分布图

Time=1.08e6s时渗透性系数分布图

渗透性系数分布图

Time=1e7s时渗透性系数分布图

压力等表面图

Time=1e7s时渗透性系数分布图 Time=1s时压力等表面分布图

压力等表面图

Time=1.28e5s时压力等表面分布图 Time=1s时压力等表面分布图

压力等表面图

Time=1.08e6s时压力等表面分布图

压力等表面图

Time=1e7s时压力等表面分布图

Z=0m切面位置瓦斯压力随时间变化分布
外 部 载 荷 为 时

0.1MPa

Z=0m切面位置瓦斯压力随时间变化分布
外 部 载 荷 为 时

2MPa

Z=0m切面位置瓦斯压力随时间变化分布
外 部 载 荷 为 时

10MPa

Z=0m切面位置瓦斯压力随时间变化分布
外 部 载 荷 为 时

20MPa

L1位置瓦斯压力随时间变化曲线
外部载荷为0.1MPa时
1
1-Time=1s 2-Time=1e5s 3-Time=5e5s 4-Time=1e6s 5-Time=5e6s 6-Time=1e7s

2 3 4 5 6

L1位置瓦斯压力随时间变化曲线
外部载荷为2MPa时
1
1-Time=1s 2-Time=1e5s 3-Time=5e5s 4-Time=1e6s 5-Time=5e6s 6-Time=1e7s

2 3 4 5 6

L1位置瓦斯压力随时间变化曲线
外部载荷为10MPa时
1
1-Time=1s 2-Time=1e5s 3-Time=5e5s 4-Time=1e6s 5-Time=5e6s 6-Time=1e7s

2 3 4 5 6

L1位置瓦斯压力随时间变化曲线
外部载荷为10MPa时
1
1-Time=1s 2-Time=1e5s 3-Time=5e5s 4-Time=1e6s 5-Time=5e6s 6-Time=1e7s

2 3 4 5 6

L2位置瓦斯压力随时间变化曲线
外部载荷为0.1MPa时
1 3
1-Time=1s 2-Time=1e5s 3-Time=5e5s 4-Time=1e6s 5-Time=5e6s 6-Time=1e7s

2

4

5 6

L2位置瓦斯压力随时间变化曲线
外部载荷为2MPa时
1 2
1-Time=1s 2-Time=1e5s 3-Time=5e5s 4-Time=1e6s 5-Time=5e6s 6-Time=1e7s

3 4 5 6

L2位置瓦斯压力随时间变化曲线
外部载荷为10MPa时
1
1-Time=1s 2-Time=1e5s 3-Time=5e5s 4-Time=1e6s 5-Time=5e6s 6-Time=1e7s

2 3 4 5 6

L2位置瓦斯压力随时间变化曲线
外部载荷为20MPa时
1
1-Time=1s 2-Time=1e5s 3-Time=5e5s 4-Time=1e6s 5-Time=5e6s 6-Time=1e7s

2 3 4 5 6

不同围压条件下瓦斯压力分布曲线 (L1处;Time=1e7s)
0.8 0.7 0.6 瓦斯压力/MPa 0.5 0.4 0.3
围压为0.1MPa

0.2 0.1 0

围压为2MPa 围压为10MPa 围压为20MPa

0

2

4

Y方向/m

6

8

10

不同围压条件下瓦斯压力分布曲线 (L2处;Time=1e7s)
0.9 0.8 瓦斯压力/Mpa 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 2 4 Y方向/m 6 8 10
围压为0.1MPa 围压为2MPa 围压为10MPa 围压为20MPa

感谢国家自然科学基金委材料学部冶金与矿业学科

2007-05-018


相关文章:
COMSOL Multiphysics在岩土工程中的应用
普通的偏微分方程形式,也可以使用 COMSOLMultiphysics 提供的特定的 物理应用模型...比如煤岩体瓦斯水渗流耦合过程数值模 型及其在矿山工程中的应用、COMSOL 在...
基于应力-渗流耦合的保护层开采瓦斯渗流特征数值研究
瓦斯渗流特征数值研究_冶金/矿山/地质_工程科技_专业...2.1 模型建立及边界条件 根据平煤集团某矿煤层群...煤岩体瓦斯水渗流耦合... 86页 免费 ©...
瓦斯渗流文献综述
瓦斯的形 成、存储释放一直伴随着整个成煤过程。...一些数值 模拟软件的引入来对瓦斯在煤层中流动进行...体积应变方程应 作为瓦斯渗流耦合分析中主要的控制...
煤层气储层的数值模型及初步应用
煤岩体瓦斯水渗流耦合... 86页 免费 测井在煤层...模型,并推导了其有限 差分的数值模型,编写了数值...煤层气的数值模拟也在气井的实践中得到了广泛的应用...
博士前沿讲座论文
因此, 采动岩体力学是矿山工程力学的主干学科,甚至...的突水、瓦斯突出、突水溃沙等自然灾害是渗流失稳...变形与渗流耦合系统动力学行为;变质量岩体 渗流;非 ...
深部煤炭开发中煤与瓦斯共采理论
建立破断煤岩体中瓦斯非稳态流动数学模型; 应用密度...应变过程、 应力场-裂隙场-瓦斯耦合以及卸压条件...综合运用矿山工程力学、工程地质学、构造地质学、岩石...
煤矿采空区瓦斯渗流规律及其数值模拟研究
煤矿采空区瓦斯渗流规律及其数值模拟研究 煤矿采空区上覆岩层结构移动规律分析 综放工艺在开采高含量瓦斯厚煤层的推广应用中之所以遇到困难, 往往是由 于综放面...
饱和—非饱和土地基力学渗流耦合问题的分析与计算
问题的饱和-非饱和多孔多相介质的有限元数值模型。...工程建筑中,因此研究此类地基的力学和渗流耦合过程,...并将本文发展的有限元数值方法应用于挡水结构双层...
饱和―非饱和土地基力学渗流耦合问题的分析与计算
问题的饱和-非饱和多孔多相介质的有限元数值模型。 ...工程建筑中,因此研究此类地基的力学和渗流耦合过程,...并将本文发展的有限元数值方法应用于挡水结构双层...
受载含瓦斯气固耦合渗流规律实验及模型研究-秦恒洁
受载含瓦斯气固耦合渗流规律实验及模型研究-秦恒洁_...同时感谢姚邦华老师在论文数值模拟过程中给予我的指导...[64]进行了不同孔隙压力时,煤岩体在全应力-应变...
更多相关标签:
渗流应力耦合 | abaqus 渗流 流固耦合 | 瓦斯渗流 | 矿山岩体力学 | 矿山岩体力学 pdf | 矿山岩体力学重点 | 矿山岩体力学801 | 工程岩体分级标准 |