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

煤岩体瓦斯、水渗流耦合过程数值模型 及其在矿山工程中的应用 东北大学 杨天鸿 3 主要物理过程 气体压缩过程、气体吸附和解析过程、扩散过气体压缩过程、气体吸附和解析过程、扩散过 程、渗流过程、应力程、渗流过程、应力- -渗流耦合过程等渗流耦合过程等 4 物理数学方程 (1)(1)瓦斯渗流方程：瓦斯渗流方程： (2)(2)气体状态方程：气体状态方程： (3) Langmuir吸附 解析方程 4 物理数学方程 (4)(4)气体压力有效应力方程：气体压力有效应力方程： 4 物理数学方程 (5)K(5)Klinkenberglinkenberg 方程：方程： （滑脱效应）（滑脱效应） (6)渗透耦合方程： 瓦斯渗流-应力耦合方程： lFEMLAB为专门求解耦合偏微分方程组的有限元分析工具。 具有强大的处理功能： 它含有一些内嵌的经典物理模型， 包括单物理场和多物理场模型，可以直接用于分析。 l 功能最强大、最灵活的还是其偏微分方程组模式：系数形 式、通式与弱形式。这三个数学应用模式中：系数形式 (Coefficient form)，适宜求解线性问题；通式(General form) ，适宜求解非线性问题；弱形式(Weak form)，最为灵活，对 于边界条件、时间序列复杂模型尤为适宜，但应用也相对复 杂些。一般地，大多数物理问题均可采用通式模式进行求解 。 求解方法 FEMLAB简介 l 对于不同物理场中交叉耦合项的处理简单有效。一 方面，在各物理场的偏微分方程中考虑了不同场的影 响；另一方面，各物理场中的计算变量可以直接用于 耦合关系的定义。 l 该软件带有Script语言并兼容Matlab语言，具有强 大的二次开发功能，对于创新性理论研究尤为适合。 此外，FEMLAB还有强大的后处理功能。 求解方法 FEMLAB简介 一、基于数字图像处理技术的煤层瓦斯渗 流过程数值模拟 l在数字图像处理技术中，人们通常采用HSI颜 色空间来表述数字图像，因为该图像空间有利 于人肉眼的识别。HSI色彩空间中，颜色用色 度（Hue）、饱和度（Saturation）和亮度（ Intensity）来表示。其中H表示了肉眼看到的 颜色，S表示该颜色相对于白色的饱和度，I表 示的亮度。HSI的颜色空间的数值可以从RGB 的数据转换而来。其中I的数值是R、G和B的 算术平均值。 （1）煤的细观数字图像 (2) 图像I值的分布 煤样的细观扫描照片及基于数值图像技术获得的孔隙率分布 l图1（a）是文献18给出的具有突出倾向性煤 样的显微照片，从中可以看出煤样中的叶状的 碎斑结构。灰度低的部分可视为裂隙。图1（b ）给出了该图像I值的分布。由此可以看出I值 较好地反映了煤体中的结构特征。 煤样中各组分中的孔隙率、渗透率和初始瓦斯压力 像素的I值 组分孔隙率弹性模量 / GPa 渗透率 / m2 初始瓦斯压 力 / MPa 0.0 I 0.6裂隙带0.22.01.0e-163.0 0.6 I 1.0基质0.0520.01.0e-181.0 （a）初始孔隙率分布 （b）弹性模量分布 基于数值图像技术获得的孔隙率和渗透率分布 0x = 2 MPa n(p) = 0 0x = 2 MPa n(p) = 0 内部边界: pi = 0.1 MPa s0y = 1 MPa n(p) = 0 X Y n(p) = 0 初始条件: p0 A(14,10) B (6,10) E (7.5,12.5) C (12.5,10) D(12.5,7.5) 瓦斯运移的数值模型 瓦斯压力图 l下图给出的瓦斯压力图可以看出瓦斯从中间井孔不断 释放的整个过程。由于这里没有考虑煤层的补给，故 瓦斯压力不断降低，最终煤样中的瓦斯压力降低到到 内部孔边界的瓦斯压力，故瓦斯的运移过程停止。在 初始条件下，由于裂隙带中的瓦斯压力比煤基质中的 高，故瓦斯由裂隙带向外围的煤机制中不断扩散，使 得裂隙带和基质间的瓦斯压力剃度不断降低。直到时 间 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 瓦斯压力及渗流速度随时间的变化过程图 l下图给出了试样中5个点（其位置见瓦斯运移的数值 模型 ）的瓦斯压力时间曲线。其中，点C和D位 于裂隙带中，具有较大的初始瓦斯压力（为3.0 MPa ）。而点A、B和E位于基质中，其初始瓦斯压力为1.0 MPa。在t = 10 s前，随着瓦斯的抽放，裂隙带中的瓦 斯压力逐渐降低，而煤基质中的瓦斯压力不断增加， 这说明瓦斯从裂隙带向基质中不断渗流。当t = 10 s时 ，试样中的瓦斯开始集中往抽放孔处流动，最终在t = 1e5 s时，试样中的瓦斯压力达到抽放孔的压力值0.1 MPa。 试样中五个点（A、B、C、D和E）处的瓦斯压力时间曲线 l下图给出了5个特征点的孔隙率的变化曲线。 与瓦斯压力的分布曲线类似，在裂隙带的点C 和D处，受到外部边界应力的作用后，孔隙率 从初值0.2降低到了0.1855和0.188。随后，随 着瓦斯的不断释放，瓦斯压力降低，故有效应 力增加，所以，孔隙被压缩，所以孔隙率随着 瓦斯压力的降低而不断降低。 瓦斯渗流过程中的煤样的孔隙率变化（断裂带中的点C和D） l但是，如下图所示，在煤基质中，t = 1e-02s 时，首先是由于煤的变形导致了孔隙率下降， 但是，随着瓦斯从裂隙带向着煤基质中渗流， 故基质中的瓦斯压力会逐渐，所以导致了孔隙 率的增加。此后，由于瓦斯开始集中向抽放孔 中流动，整个煤样中的瓦斯压力不断降低会引 起有效应力的增加，故孔隙率会不断降低。 瓦斯渗流过程中的煤样的孔隙率变化(煤基质中的点A、B和E) 主要结论如下： l（1）由于裂隙带中的瓦斯压力比煤基质中的高，故 瓦斯由裂隙带向外围的煤机制中不断扩散，使得裂隙 带和基质间的瓦斯压力剃度不断降低。在某个时间后 ，瓦斯开始集中向抽放孔中流动。 l（2）在本文给定的模拟参数下，应力场引起的煤体 压缩是对其渗透性的影响要大于Klinkerberg效应和瓦 斯解吸效应所引起的渗透率变化。故总体上煤层中表 现出渗透率的降低，随着瓦斯的不断抽放，渗透率更 是不断降低。 二、冒落区瓦斯浓度扩散-对流及 风流场数值模拟 l瓦斯浓度扩散-对流场基本方程： l l DarcyDarcy方程 冒落区破碎岩体压实过程中的气体流场 l非DarcyDarcy方程(堆石体、土石坝流场) 冒落区瓦斯浓度扩散-对流 及风流场数值模拟 l通风流场基本方程： l Navier-Stokes equations: Free flow适合通风 巷道风流场 式中：v流体流速，m/s；p流体压力，Pa； 流体密度，kg/ m3；I单位矢量；F流体阻力 冒落区瓦斯浓度扩散-对流 及风流场数值模拟 l通风流场基本方程： lBrinkman equations: Fast flow in porous media，适 合冒落区风流场 模型建立 l计算模型和方案 l参照综采工作面的具体尺寸，建立如下图所示 的二维计算模型，不考虑势能，模型东西向即 工作面推进方向取400m,南北向宽240m。划 分成均质（5个分区，沿工作向采空区方向的 区域（1）至（5）宽度均取为80m）；即大约 1周的推进时间瞬时完成。假设这个期间内每 个区域的透气系数分布均不同，具体采用的计 算模型相关条件如下： 模型建立 不 透 气 边 界 不 透 气 边 界 浓度补给边界 采空冒落区 （5） （4） （3） （2） （1） 进风口 出风口 推进方向 进风口进风口进风口进风口 出风口出风口出风口出风口 浓度补给边界 模型建立 l边界条件： l（1）通风条件：左侧下20m为进风口边界， 左侧上20m为回风口边界，压力差100Pa,其他 边界为不透气边界。 l（2）扩散条件：右侧边界为绝缘对称边界， 其他边界有补给，推进区域1时，上下边界补 给量为1.2e-6mol/m2s，汇源项为3e-6 mol/m2s，以描述瓦斯通量随开采动态过程而 增加的瓦斯量。 模型建立 l初始条件： l域内具有一个大气压，瓦斯初始浓度3mol/m3 ；推进新的工作面区域时，旧区域瓦斯浓度模 拟结果C0再附加一初始补充浓度Cbc （取 3mol/m3）作为旧推进工作面的初始浓度 。 模型建立 l时间步长： l按照非均等的积数步长增大，初始值7s,终止 指7e5s，设定100个中间时间值。 l计算参数： l动粘系数=1.8e-5pas，流体密度=1200kg/m3 ，扩散系数D=2e-5m2 ，瞬态时间比例系数 =0.55 模型建立 12345 1 3.24e-8 2 2.88e-83.24e-8 3 2.52e-82.88e-83.24e-8 4 2.16e-82.52e-82.88e-83.24e-8 5 1.8e-82.16e-82.52e-82.88e-83.24e-8 已推区域 透气率 （m2） 新推区域 随推进进行透气率变化表 模拟结果分析 l推进区域1时的计算结果 l下图分别是当时间为：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时瓦斯浓度 模拟结果分析 Time=7e5s时瓦斯浓度和流线 A1A1 模拟结果分析 第一步推进时A1-A1切面浓度变化曲线 模拟结果分析 l推进区域2时的计算结果 l下图分别是当时间为：7s，7e3s，1.1e5s， 7e5s瞬态时间区域2采空区瓦斯浓度分布图 模拟结果分析 Time=7s时瓦斯浓度 模拟结果分析 Time=7e3s时瓦斯浓度 模拟结果分析 Time=1.1e5s时瓦斯浓度 模拟结果分析 第二步推进时Time=7e5s时瓦斯浓度和流线 A2A2 模拟结果分析 第二步推进时A2-A2切面浓度变化曲线 模拟结果分析 l推进区域3,4时的计算结果 模拟结果分析 第三步推进时Time=7e5s时瓦斯浓度和流线 A3A3 模拟结果分析 第三步推进时A3-A3切面浓度变化曲线 模拟结果分析 第四步推进时Time=7e5s时瓦斯浓度和流线 A4A4 模拟结果分析 第四步推进时A4-A4切面浓度变化曲线 讨论 l（1）在区域1右侧不透气边界附近，随着工作 面向左推进，风流作用逐渐减弱，该处瓦斯浓 度不断集聚增大，每个推进步瓦斯浓度分别为 7mol，12mol，17 mol，24 mol，表明该处距 离通风口越远，瓦斯集聚程度越大。 讨论 l（2）对于每个推进步，当时间达到7e5s（8d 左右）时，风流基本上可以把本区域瓦斯浓度 降低到安全范围内（接近0mol），表明风流 在本区域的流速明显大于前一个区域。而前一 个区域的下左半部分（占整个区域的1/4）的 瓦斯浓度也能得到有效降低，而对其他区域瓦 斯浓度驱散作用影响有限。 讨论 l（3）图9可见，瓦斯浓度沿工作面推进方向从 近到远瓦斯浓度成台阶状渐次增大，这表明通 风量只能在一定范围内降低瓦斯浓度，同时也 看出距离通风口越远，瓦斯补给的时间和补给 量越大。 讨论 l（4）虽然模型上下边界都是瓦斯补给边界， 但由于模型左下边界是进风口，瓦斯随风流对 流作用得到有效降低；而模型左上边界是出风 口，瓦斯容易在回风隅角大量聚集，所以该处 瓦斯浓度较高。 三、三维渗流耦合模型及瓦斯抽放 l根据实际的三维煤层瓦斯抽放过程，建立如图 1所示的长宽高为10m10m10m理想化的三 维计算模型，模型左下部边界假设为巷道，三 个瓦斯抽放孔K1、K 2、K 3之间的距离 1m，与水平面呈45倾角展布，K 1、K 3与K 2呈16夹角分布在孔2两侧，三个瓦斯抽放孔 长度7m，瓦斯抽放孔按照“以缝代孔”16原则 简化为定压力边界（25kPa）。 巷道所在位置 45 16 L1 L2 K1 K3 K2 外部载荷 l(1) 边界条件：四周和上下面均为零通量不透 气边界；四周和下部都约束法线方向的位移， 上部自由，作用有P分别为10MPa、20MPa、 2MPa、0.1MPa的外部载荷，同时模型具有自 重载荷。 l(2) 初始条件：内部有1MPa的初始瓦斯压力， 三个抽放孔的压力为0.25e5Pa。 l(3) 时间步长：按照非均等的积数步长增大， 初始值1s，终止指1e7s（10d左右），设定 100个中间时间步长。 l(4) 计算方案：模拟不同外部载荷条件下（ 2MPa、10MPa、20MPa、0.1MPa），瓦斯 抽放效果和渗透性变化规律。 l(5) 计算参数：相关参数列表1所列。 杨氏模量 E （Pa） 泊松比 煤层密度 s （kg/m3） 流体密度 1 （kg/m3） 流体粘滞性 （Pas） 饱和渗透率 （m2） 7.41e9 0.33 1250 1000 1e-3 5e-6 模拟结果分析 l推进区域2时的计算结果 l下图分别是当时间为：1s， 1.28e5s ， 1.08e6s，1e7s瞬态时间渗透性系数分布图 和压力等表面图 渗透性系数分布图 Time1s时渗透性系数分布图 渗透性系数分布图 Time1.28e5s时渗透性系数分布图 渗透性系数分布图 Time1.08e6s时渗透性系数分布图 渗透性系数分布图 Time1e7s时渗透性系数分布图 压力等表面图 Time1e7s时渗透性系数分布图 Time1s时压力等表面分布图 压力等表面图 Time1s时压力等表面分布图 Time1.28e5s时压力等表面分布图 压力等表面图 Time1.08e6s时压力等表面分布图 压力等表面图 Time1e7s时压力等表面分布图 Z=0m切面位置瓦斯压力随时间变化分布 外部载荷为0.1MPa时 Z=0m切面位置瓦斯压力随时间变化分布 外部载荷为2MPa时 Z=0m切面位置瓦斯压力随时间变化分布 外部载荷为10MPa时 Z=0m切面位置瓦斯压力随时间变化分布 外部载荷为20MPa时 L1位置瓦斯压力随时间变化曲线 1 2 3 4 5 6 1Time1s 2Time1e5s 3Time5e5s 4Time1e6s 5Time5e6s 6Time1e7s 外部载荷为0.1MPa时 L1位置瓦斯压力随时间变化曲线 外部载荷为2MPa时 1 2 3 4 5 6 1Time1s 2Time1e5s 3Time5e5s 4Time1e6s 5Time5e6s 6Time1e7s L1位置瓦斯压力随时间变化曲线 1 2 3 4 5 6 1Time1s 2Time1e5s 3Time5e5s 4Time1e6s 5Time5e6s 6Time1e7s 外部载荷为10MPa时 L1位置瓦斯压力随时间变化曲线 1Time1s 2Time1e5s 3Time5e5s 4Time1e6s 5Time5e6s 6Time1e7s 外部载荷为10MPa时 1 2 3 4 5 6 L2位