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

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


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

东北大学

杨天鸿

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

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


qg ? ? kg

?g

? ??p ? ? g gD?

(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 ? ? ? ?? ?g ? p ? ? ? ? ? ? ?

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

渗透率 初始瓦斯压力 / m2 / MPa 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

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)

p (Pa)

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
?

K Darcy方程 ?h vi ? ?

kn ?xi g qg ? ? ? ??p ? ? g gD?

?g

冒落区破碎岩体压实过程中的气体流场
?

非Darcy方程(堆石体、土石坝流场)

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

通风流场基本方程: Navier-Stokes equations: Free flow适合通风 巷道风流场

?v ? ?v ? ? ? ?? pI ? ? (?v ? (?v)T )? ? F
式中: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

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

20

40 推进方向/m

60

80

第一步推进时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

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

40

80

120 推进方向/m

160

200

240

第三步推进时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

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

0

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


相关文章:
更多相关标签: