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

煤岩体瓦斯、水渗流耦合过程数值模型及其在矿山工程中的应用东北大学杨天鸿东北大学杨天鸿13主要物理过程气体压缩过程、气体吸附和解析过程、扩散过程、渗流过程、应力-渗流耦合过程等气体压缩过程、气体吸附和解析过程、扩散过程、渗流过程、应力-24物理数学方程(1)瓦斯渗流方程：(2)气体状态方程：(3)Langmuir吸附解析方程

(1)瓦斯渗流方程：(2)气体状态方程：(3)Langmu34物理数学方程(4)气体压力有效应力方程：(4)气体压力有效应力方程：44物理数学方程(5)Klinkenberg

方程：（滑脱效应）(6)渗透耦合方程：

(5)Klinkenberg方程：(6)渗透耦合方程：5瓦斯渗流-应力耦合方程：瓦斯渗流-应力耦合方程：6FEMLAB为专门求解耦合偏微分方程组的有限元分析工具。具有强大的处理功能：⑴它含有一些内嵌的经典物理模型，包括单物理场和多物理场模型，可以直接用于分析。⑵功能最强大、最灵活的还是其偏微分方程组模式：系数形式、通式与弱形式。这三个数学应用模式中：系数形式(Coefficientform)，适宜求解线性问题；通式(Generalform)，适宜求解非线性问题；弱形式(Weakform)，最为灵活，对于边界条件、时间序列复杂模型尤为适宜，但应用也相对复杂些。一般地，大多数物理问题均可采用通式模式进行求解。求解方法–FEMLAB简介FEMLAB为专门求解耦合偏微分方程组的有限元分析工具。具有7⑶对于不同物理场中交叉耦合项的处理简单有效。一方面，在各物理场的偏微分方程中考虑了不同场的影响；另一方面，各物理场中的计算变量可以直接用于耦合关系的定义。⑷该软件带有Script语言并兼容Matlab语言，具有强大的二次开发功能，对于创新性理论研究尤为适合。此外，FEMLAB还有强大的后处理功能。求解方法–FEMLAB简介⑶对于不同物理场中交叉耦合项的处理简单有效。一方面，在各物8一、基于数字图像处理技术的煤层瓦斯渗流过程数值模拟在数字图像处理技术中，人们通常采用HSI颜色空间来表述数字图像，因为该图像空间有利于人肉眼的识别。HSI色彩空间中，颜色用色度（Hue）、饱和度（Saturation）和亮度（Intensity）来表示。其中H表示了肉眼看到的颜色，S表示该颜色相对于白色的饱和度，I表示的亮度。HSI的颜色空间的数值可以从RGB的数据转换而来。其中I的数值是R、G和B的算术平均值。一、基于数字图像处理技术的煤层瓦斯渗流过程数值模拟在数字图9（1）煤的细观数字图像(2)图像I值的分布煤样的细观扫描照片及基于数值图像技术获得的孔隙率分布

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

煤样中各组分中的孔隙率、渗透率和初始瓦斯压力像素的I值组分孔隙率弹性模量/GPa渗透率/m2初始瓦斯压力/MPa0.0≤I<0.6裂隙带0.22.01.0e-163.00.6≤I≤1.0基质0.0520.01.0e-181.0煤样中各组分中的孔隙率、渗透率和初始瓦斯压力组分孔隙率弹性12

（a）初始孔隙率分布

（b）弹性模量分布基于数值图像技术获得的孔隙率和渗透率分布（a）初始孔隙率分布 13σ0x=2MPa

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

n(p)=0内部边界:pi=0.1MPas0y=1MPan(p)=0XYn(p)=0初始条件:p0A(14,10)B(6,10)E(7.5,12.5)C(12.5,10)D(12.5,7.5)瓦斯运移的数值模型

σ0x=2MPa

n(p)=0σ0x=214瓦斯压力图下图给出的瓦斯压力图可以看出瓦斯从中间井孔不断释放的整个过程。由于这里没有考虑煤层的补给，故瓦斯压力不断降低，最终煤样中的瓦斯压力降低到到内部孔边界的瓦斯压力，故瓦斯的运移过程停止。在初始条件下，由于裂隙带中的瓦斯压力比煤基质中的高，故瓦斯由裂隙带向外围的煤机制中不断扩散，使得裂隙带和基质间的瓦斯压力剃度不断降低。直到时间t=1e04s后，瓦斯开始集中向抽放孔中流动。最终在t=1e06s左右时，煤样内的瓦斯压力和抽放孔中的给定压力相同，瓦斯流动过程停止。瓦斯压力图下图给出的瓦斯压力图可以看出瓦斯从中间井孔不断释放15t=1e-01st=1e0st=1e01st=1e02st=1e03st=1e04st=1e05st=1e06st=1e07s瓦斯压力及渗流速度随时间的变化过程图t=1e-01s瓦斯压力及渗流速度随时间的变化过程图16下图给出了试样中5个点（其位置见瓦斯运移的数值模型）的瓦斯压力——时间曲线。其中，点C和D位于裂隙带中，具有较大的初始瓦斯压力（为3.0MPa）。而点A、B和E位于基质中，其初始瓦斯压力为1.0MPa。在t=10s前，随着瓦斯的抽放，裂隙带中的瓦斯压力逐渐降低，而煤基质中的瓦斯压力不断增加，这说明瓦斯从裂隙带向基质中不断渗流。当t=10s时，试样中的瓦斯开始集中往抽放孔处流动，最终在t=1e5s时，试样中的瓦斯压力达到抽放孔的压力值0.1MPa。下图给出了试样中5个点（其位置见瓦斯运移的数值模型）的瓦斯17试样中五个点（A、B、C、D和E）处的瓦斯压力—时间曲线试样中五个点（A、B、C、D和E）处的瓦斯压力—时间曲线18下图给出了5个特征点的孔隙率的变化曲线。与瓦斯压力的分布曲线类似，在裂隙带的点C和D处，受到外部边界应力的作用后，孔隙率从初值0.2降低到了0.1855和0.188。随后，随着瓦斯的不断释放，瓦斯压力降低，故有效应力增加，所以，孔隙被压缩，所以孔隙率随着瓦斯压力的降低而不断降低。下图给出了5个特征点的孔隙率的变化曲线。与瓦斯压力的分布曲线19瓦斯渗流过程中的煤样的孔隙率变化（断裂带中的点C和D）瓦斯渗流过程中的煤样的孔隙率变化（断裂带中的点C和D）20但是，如下图所示，在煤基质中，t=1e-02s时，首先是由于煤的变形导致了孔隙率下降，但是，随着瓦斯从裂隙带向着煤基质中渗流，故基质中的瓦斯压力会逐渐，所以导致了孔隙率的增加。此后，由于瓦斯开始集中向抽放孔中流动，整个煤样中的瓦斯压力不断降低会引起有效应力的增加，故孔隙率会不断降低。但是，如下图所示，在煤基质中，t=1e-02s时，首先是21瓦斯渗流过程中的煤样的孔隙率变化(煤基质中的点A、B和E)瓦斯渗流过程中的煤样的孔隙率变化(煤基质中的点A、B和E)22主要结论如下：（1）由于裂隙带中的瓦斯压力比煤基质中的高，故瓦斯由裂隙带向外围的煤机制中不断扩散，使得裂隙带和基质间的瓦斯压力剃度不断降低。在某个时间后，瓦斯开始集中向抽放孔中流动。（2）在本文给定的模拟参数下，应力场引起的煤体压缩是对其渗透性的影响要大于Klinkerberg效应和瓦斯解吸效应所引起的渗透率变化。故总体上煤层中表现出渗透率的降低，随着瓦斯的不断抽放，渗透率更是不断降低。主要结论如下：23二、冒落区瓦斯浓度扩散-对流及

风流场数值模拟瓦斯浓度扩散-对流场基本方程：Darcy方程二、冒落区瓦斯浓度扩散-对流及

风流场数值模拟瓦斯浓度扩散-24冒落区破碎岩体压实过程中的气体流场非Darcy方程(堆石体、土石坝流场)冒落区破碎岩体压实过程中的气体流场非Darcy方程(堆石体、25冒落区瓦斯浓度扩散-对流

及风流场数值模拟通风流场基本方程：

Navier-Stokesequations:Freeflow适合通风巷道风流场式中：v——流体流速，m/s；p——流体压力，Pa；ρ——流体密度，kg/m3；I——单位矢量；F——流体阻力冒落区瓦斯浓度扩散-对流

及风流场数值模拟通风流场基本方程：26冒落区瓦斯浓度扩散-对流

及风流场数值模拟通风流场基本方程：Brinkmanequations:Fastflowinporousmedia，适合冒落区风流场冒落区瓦斯浓度扩散-对流

及风流场数值模拟通风流场基本方程：27模型建立计算模型和方案参照综采工作面的具体尺寸，建立如下图所示的二维计算模型，不考虑势能，模型东西向即工作面推进方向取400m,南北向宽240m。划分成均质（5个分区，沿工作向采空区方向的区域（1）至（5）宽度均取为80m）；即大约1周的推进时间瞬时完成。假设这个期间内每个区域的透气系数分布均不同，具体采用的计算模型相关条件如下：模型建立计算模型和方案28模型建立不透气边界不透气边界浓度补给边界采空冒落区

（5）（4）（3）（2）（1）进风口出风口推进方向进风口进风口进风口进风口出风口出风口出风口出风口浓度补给边界模型建立不透气不透气浓度补给边界采空冒落区进风口出风口推进方29模型建立边界条件：（1）通风条件：左侧下20m为进风口边界，左侧上20m为回风口边界，压力差100Pa,其他边界为不透气边界。（2）扩散条件：右侧边界为绝缘对称边界，其他边界有补给，推进区域1时，上下边界补给量为1.2e-6mol/m2·s，汇源项为3e-6mol/m2·s，以描述瓦斯通量随开采动态过程而增加的瓦斯量。模型建立边界条件：30模型建立初始条件：域内具有一个大气压，瓦斯初始浓度3mol/m3；推进新的工作面区域时，旧区域瓦斯浓度模拟结果C0再附加一初始补充浓度Cbc（取3mol/m3）作为旧推进工作面的初始浓度。模型建立初始条件：31模型建立时间步长：按照非均等的积数步长增大，初始值7s,终止指7e5s，设定100个中间时间值。计算参数：动粘系数=1.8e-5pa·s，流体密度=1200kg/m3，扩散系数D=2e-5m2，瞬态时间比例系数=0.55模型建立时间步长：32模型建立1234513.24e-8

22.88e-83.24e-832.52e-82.88e-83.24e-842.16e-82.52e-82.88e-83.24e-851.8e-82.16e-82.52e-82.88e-83.24e-8已推区域透气率（m2）新推区域随推进进行透气率变化表模型建立1234513.24e-822.88e-83.2433模拟结果分析推进区域1时的计算结果下图分别是当时间为：7s，7e3s，7e4s，1.1e5s，2.5e5s，3.5e5s，7e5s瞬态时间区域1采空区瓦斯浓度分布图模拟结果分析推进区域1时的计算结果34模拟结果分析Time=7s时瓦斯浓度

模拟结果分析Time=7s时瓦斯浓度35模拟结果分析Time=7e3s时瓦斯浓度

模拟结果分析Time=7e3s时瓦斯浓度36模拟结果分析Time=7e4s时瓦斯浓度

模拟结果分析Time=7e4s时瓦斯浓度37模拟结果分析Time=1.1e5s时瓦斯浓度

模拟结果分析Time=1.1e5s时瓦斯浓度38模拟结果分析Time=2.5e5s时瓦斯浓度

模拟结果分析Time=2.5e5s时瓦斯浓度39模拟结果分析Time=3.5e5s时瓦斯浓度

模拟结果分析Time=3.5e5s时瓦斯浓度40模拟结果分析Time=7e5s时瓦斯浓度和流线

A1A1’模拟结果分析Time=7e5s时瓦斯浓度和流线A1A1’41模拟结果分析第一步推进时A1-A1’切面浓度变化曲线

模拟结果分析第一步推进时A1-A1’切面浓度变化曲线42模拟结果分析推进区域2时的计算结果下图分别是当时间为：7s，7e3s，1.1e5s，7e5s瞬态时间区域2采空区瓦斯浓度分布图模拟结果分析推进区域2时的计算结果43模拟结果分析Time=7s时瓦斯浓度模拟结果分析Time=7s时瓦斯浓度44模拟结果分析Time=7e3s时瓦斯浓度模拟结果分析Time=7e3s时瓦斯浓度45模拟结果分析Time=1.1e5s时瓦斯浓度模拟结果分析Time=1.1e5s时瓦斯浓度46模拟结果分析第二步推进时Time=7e5s时瓦斯浓度和流线

A2A2’模拟结果分析第二步推进时Time=7e5s时瓦斯浓度和流线47模拟结果分析第二步推进时A2-A2’切面浓度变化曲线

模拟结果分析第二步推进时A2-A2’切面浓度变化曲线48模拟结果分析推进区域3,4时的计算结果模拟结果分析推进区域3,4时的计算结果49模拟结果分析第三步推进时Time=7e5s时瓦斯浓度和流线

A3A3’模拟结果分析第三步推进时Time=7e5s时瓦斯浓度和流线50模拟结果分析第三步推进时A3-A3’切面浓度变化曲线

模拟结果分析第三步推进时A3-A3’切面浓度变化曲线51模拟结果分析第四步推进时Time=7e5s时瓦斯浓度和流线

A4A4’模拟结果分析第四步推进时Time=7e5s时瓦斯浓度和流线52模拟结果分析第四步推进时A4-A4’切面浓度变化曲线

模拟结果分析第四步推进时A4-A4’切面浓度变化曲线53讨论（1）在区域1右侧不透气边界附近，随着工作面向左推进，风流作用逐渐减弱，该处瓦斯浓度不断集聚增大，每个推进步瓦斯浓度分别为7mol，12mol，17mol，24mol，表明该处距离通风口越远，瓦斯集聚程度越大。讨论（1）在区域1右侧不透气边界附近，随着工作面向左推进，风54讨论（2）对于每个推进步，当时间达到7e5s（8d左右）时，风流基本上可以把本区域瓦斯浓度降低到安全范围内（接近0mol），表明风流在本区域的流速明显大于前一个区域。而前一个区域的下左半部分（占整个区域的1/4）的瓦斯浓度也能得到有效降低，而对其他区域瓦斯浓度驱散作用影响有限。讨论55讨论（3）图9可见，瓦斯浓度沿工作面推进方向从近到远瓦斯浓度成台阶状渐次增大，这表明通风量只能在一定范围内降低瓦斯浓度，同时也看出距离通风口越远，瓦斯补给的时间和补给量越大。讨论56讨论（4）虽然模型上下边界都是瓦斯补给边界，但由于模型左下边界是进风口，瓦斯随风流对流作用得到有效降低；而模型左上边界是出风口，瓦斯容易在回风隅角大量聚集，所以该处瓦斯浓度较高。讨论57三、三维渗流耦合模型及瓦斯抽放根据实际的三维煤层瓦斯抽放过程，建立如图1所示的长宽高为10m×10m×10m理想化的三维计算模型，模型左下部边界假设为巷道，三个瓦斯抽放孔——K1、K2、K3之间的距离1m，与水平面呈45°倾角展布，K1、K3与K2呈16°夹角分布在孔2两侧，三个瓦斯抽放孔长度7m，瓦斯抽放孔按照“以缝代孔”[16]原则简化为定压力边界（25kPa）。三、三维渗流耦合模型及瓦斯抽放根据实际的三维煤层瓦斯抽放过程58巷道所在位置45º16ºL1L2K1K3K2外部载荷巷道所在位置45º16ºL1L2K1K3K2外部载荷59(1)边界条件：四周和上下面均为零通量不透气边界；四周和下部都约束法线方向的位移，上部自由，作用有P分别为10MPa、20MPa、2MPa、0.1MPa的外部载荷，同时模型具有自重载荷。(2)初始条件：内部有1MPa的初始瓦斯压力，三个抽放孔的压力为0.25e5Pa。(1)边界条件：四周和上下面均为零通量不透气边界；四周和下60(3)时间步长：按照非均等的积数步长增大，初始值1s，终止指1e7s（10d左右），设定100个中间时间步长。(4)计算方案：模拟不同外部载荷条件下（2MPa、10MPa、20MPa、0.1MPa），瓦斯抽放效果和渗透性变化规律。(5)计算参数：相关参数列表1所列。(3)时间步长：按照非均等的积数步长增大，初始值1s，终止61杨氏模量E（Pa）泊松比ν煤层密度ρs（kg/m3）流体密度ρ1（kg/m3）流体粘滞性（Pa·s）

饱和渗透率（m2）7.41e90.33125010001e-35e-6杨氏模量泊松比煤层密度流体密度流体粘滞性饱和渗透率7.41e62模拟结果分析推进区域2时的计算结果下图分别是当时间为：1s，1.28e5s，1.08e6s，1e7s瞬态时间渗透性系数分布图和压力等表面图模拟结果分析推进区域2时的计算结果63渗透性系数分布图Time＝1s时渗透性系数分布图渗透性系数分布图Time＝1s时渗透性系数分布图64渗透性系数分布图Time＝1.28e5s时渗透性系数分布图渗透性系数分布图Time＝1.28e5s时渗透性系数分布图65渗透性系数分布图Time＝1.08e6s时渗透性系数分布图渗透性系数分布图Time＝1.08e6s时渗透性系数分布图66渗透性系数分布图Time＝1e7s时渗透性系数分布图渗透性系数分布图Time＝1e7s时渗透性系数分布图67压力等表面图Time＝1e7s时渗透性系数分布图Time＝1s时压力等表面分布图压力等表面图Time＝1e7s时渗透性系数分布图Time＝68压力等表面图Time＝1s时压力等表面分布图Time＝1.28e5s时压力等表面分布图压力等表面图Time＝1s时压力等表面分布图Time＝1.69压力等表面图Time＝1.08e6s时压力等表面分布图压力等表面图Time＝1.08e6s时压力等表面分布图70压力等表面图Time＝1e7s时压力等表面分布图压力等表面图Time＝1e7s时压力等表面分布图71Z=0m切面位置瓦斯压力随时间变化分布外部载荷为0.1MPa时Z=0m切面位置瓦斯压力随时间变化分布外部载荷为0.1MP72Z=0m切面位置瓦斯压力随时间变化分布外部载荷为2MPa时Z=0m切面位置瓦斯压力随时间变化分布外部载荷为2MPa时73Z=0m切面位置瓦斯压力随时间变化分布外部载荷为10MPa时Z=0m切面位置瓦斯压力随时间变化分布外部载荷为10MPa74Z=0m切面位置瓦斯压力随时间变化分布外部载荷为20MPa时Z=0m切面位置瓦斯压力随时间变化分布外部载荷为20MPa75L1位置瓦斯压力随时间变化曲线1234561－Time＝1s

2－Time＝1e5s

3－Time＝5e5s

4－Time＝1e6s

5－Time＝5e6s

6－Time＝1e7s外部载荷为0.1MPa时L1位置瓦斯压力随时间变化曲线1234561－Time＝176L1位置瓦斯压力随时间变化曲线外部载荷为2MPa时1234561－Time＝1s

2－Time＝1e5s

3－Time＝5e5s

4－Time＝1e6s

5－Time＝5e6s

6－Time＝1e7sL1位置瓦斯压力随时间变化曲线外部载荷为2MPa时123477L1位置瓦斯压力随时间变化曲线1234561－Time＝1s

2－Time＝1e5s

3－Time＝5e5s

4－Time＝1e6s

5－Time＝5e6s

6－Time＝1e7s外部载荷为10MPa时L1位置瓦斯压力随时间变化曲线1234561－Time＝178L1位置瓦斯压力随时间变化曲线1－Time＝1s

2－Time＝1e5s

3－Time＝5e5s

4－Time＝1e6s

5－Time＝5e6s

6－Time＝1e7s外部载荷为10MPa时123456L1位置瓦斯压力随时间变化曲线1－Time＝1s

2－Ti79L2位置瓦斯压力随时间变化曲线1－Time＝1s

2－Time＝1e5s

3－Time＝5e5s

4－Time＝1e6s

5－Time＝5e6s

6－Time＝1e7s外部载荷为0.1MPa时123456L2位置瓦斯压力随时间变化曲线1－Time＝1s

2－Ti80L2位置瓦斯压力随时间变化曲线1－Time＝1s

2－Time＝1e5s

3－Time＝5e5s

4－Time＝1e6s

5－Time＝5e6s

6－Time＝1e7s外部载荷为2MPa时123456L2位置瓦斯压力随时间变化曲线1－Time＝1s

2－Ti81L2位置瓦斯压力随时间变化曲线1－Time＝1s

2－Time＝1e5s

3－Time＝5e5s

4－Time＝1e6s

5－Time＝5e6s

6－Time＝1e7s外部载荷为10MPa时213456L2位置瓦斯压力随时间变化曲线1－Time＝1s

2－Ti82L2位置瓦斯压力随时间变化曲线1－Time＝1s

2－Time＝1e5s

3－Time＝5e5s

4－Time＝1e6s

5－Time＝5e6s

6－Time＝1e7s外部载荷为20MPa时134562L2位置瓦斯压力随时间变化曲线1－Time＝1s

2－Ti83不同围压条件下瓦斯压力分布曲线

（L1处；Time=1e7s）00.10.20.30.40.50.60.70.80246810Y方向/m围压为0.1MPa围压为2MPa围压为10MPa围压为20MPa瓦斯压力/MPa不同围压条件下瓦斯压力分布曲线

（L1处；Time=1e7s8400.10.20.30.40.50.60.70.80.90246810Y方向/m围压为0.1MPa围压为2MPa围压为10MPa围压为20MPa瓦斯压力/Mpa不同围压条件下瓦斯压力分布曲线

（L2处；Time=1e7s）00.10.20.30.40.50.60.70.80.902852007-05-018THANKS谢谢感谢国家自然科学基金委材料学部冶金与矿业学科2007-05-018THANKS谢谢感谢国家自然科学基金委86煤岩体瓦斯、水渗流耦合过程数值模型及其在矿山工程中的应用东北大学杨天鸿东北大学杨天鸿873主要物理过程气体压缩过程、气体吸附和解析过程、扩散过程、渗流过程、应力-渗流耦合过程等气体压缩过程、气体吸附和解析过程、扩散过程、渗流过程、应力-884物理数学方程(1)瓦斯渗流方程：(2)气体状态方程：(3)Langmuir吸附解析方程

(1)瓦斯渗流方程：(2)气体状态方程：(3)Langmu894物理数学方程(4)气体压力有效应力方程：(4)气体压力有效应力方程：904物理数学方程(5)Klinkenberg

方程：（滑脱效应）(6)渗透耦合方程：

(5)Klinkenberg方程：(6)渗透耦合方程：91瓦斯渗流-应力耦合方程：瓦斯渗流-应力耦合方程：92FEMLAB为专门求解耦合偏微分方程组的有限元分析工具。具有强大的处理功能：⑴它含有一些内嵌的经典物理模型，包括单物理场和多物理场模型，可以直接用于分析。⑵功能最强大、最灵活的还是其偏微分方程组模式：系数形式、通式与弱形式。这三个数学应用模式中：系数形式(Coefficientform)，适宜求解线性问题；通式(Generalform)，适宜求解非线性问题；弱形式(Weakform)，最为灵活，对于边界条件、时间序列复杂模型尤为适宜，但应用也相对复杂些。一般地，大多数物理问题均可采用通式模式进行求解。求解方法–FEMLAB简介FEMLAB为专门求解耦合偏微分方程组的有限元分析工具。具有93⑶对于不同物理场中交叉耦合项的处理简单有效。一方面，在各物理场的偏微分方程中考虑了不同场的影响；另一方面，各物理场中的计算变量可以直接用于耦合关系的定义。⑷该软件带有Script语言并兼容Matlab语言，具有强大的二次开发功能，对于创新性理论研究尤为适合。此外，FEMLAB还有强大的后处理功能。求解方法–FEMLAB简介⑶对于不同物理场中交叉耦合项的处理简单有效。一方面，在各物94一、基于数字图像处理技术的煤层瓦斯渗流过程数值模拟在数字图像处理技术中，人们通常采用HSI颜色空间来表述数字图像，因为该图像空间有利于人肉眼的识别。HSI色彩空间中，颜色用色度（Hue）、饱和度（Saturation）和亮度（Intensity）来表示。其中H表示了肉眼看到的颜色，S表示该颜色相对于白色的饱和度，I表示的亮度。HSI的颜色空间的数值可以从RGB的数据转换而来。其中I的数值是R、G和B的算术平均值。一、基于数字图像处理技术的煤层瓦斯渗流过程数值模拟在数字图95（1）煤的细观数字图像(2)图像I值的分布煤样的细观扫描照片及基于数值图像技术获得的孔隙率分布

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

煤样中各组分中的孔隙率、渗透率和初始瓦斯压力像素的I值组分孔隙率弹性模量/GPa渗透率/m2初始瓦斯压力/MPa0.0≤I<0.6裂隙带0.22.01.0e-163.00.6≤I≤1.0基质0.0520.01.0e-181.0煤样中各组分中的孔隙率、渗透率和初始瓦斯压力组分孔隙率弹性98

（a）初始孔隙率分布

（b）弹性模量分布基于数值图像技术获得的孔隙率和渗透率分布（a）初始孔隙率分布 99σ0x=2MPa

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

n(p)=0内部边界:pi=0.1MPas0y=1MPan(p)=0XYn(p)=0初始条件:p0A(14,10)B(6,10)E(7.5,12.5)C(12.5,10)D(12.5,7.5)瓦斯运移的数值模型

σ0x=2MPa

n(p)=0σ0x=2100瓦斯压力图下图给出的瓦斯压力图可以看出瓦斯从中间井孔不断释放的整个过程。由于这里没有考虑煤层的补给，故瓦斯压力不断降低，最终煤样中的瓦斯压力降低到到内部孔边界的瓦斯压力，故瓦斯的运移过程停止。在初始条件下，由于裂隙带中的瓦斯压力比煤基质中的高，故瓦斯由裂隙带向外围的煤机制中不断扩散，使得裂隙带和基质间的瓦斯压力剃度不断降低。直到时间t=1e04s后，瓦斯开始集中向抽放孔中流动。最终在t=1e06s左右时，煤样内的瓦斯压力和抽放孔中的给定压力相同，瓦斯流动过程停止。瓦斯压力图下图给出的瓦斯压力图可以看出瓦斯从中间井孔不断释放101t=1e-01st=1e0st=1e01st=1e02st=1e03st=1e04st=1e05st=1e06st=1e07s瓦斯压力及渗流速度随时间的变化过程图t=1e-01s瓦斯压力及渗流速度随时间的变化过程图102下图给出了试样中5个点（其位置见瓦斯运移的数值模型）的瓦斯压力——时间曲线。其中，点C和D位于裂隙带中，具有较大的初始瓦斯压力（为3.0MPa）。而点A、B和E位于基质中，其初始瓦斯压力为1.0MPa。在t=10s前，随着瓦斯的抽放，裂隙带中的瓦斯压力逐渐降低，而煤基质中的瓦斯压力不断增加，这说明瓦斯从裂隙带向基质中不断渗流。当t=10s时，试样中的瓦斯开始集中往抽放孔处流动，最终在t=1e5s时，试样中的瓦斯压力达到抽放孔的压力值0.1MPa。下图给出了试样中5个点（其位置见瓦斯运移的数值模型）的瓦斯103试样中五个点（A、B、C、D和E）处的瓦斯压力—时间曲线试样中五个点（A、B、C、D和E）处的瓦斯压力—时间曲线104下图给出了5个特征点的孔隙率的变化曲线。与瓦斯压力的分布曲线类似，在裂隙带的点C和D处，受到外部边界应力的作用后，孔隙率从初值0.2降低到了0.1855和0.188。随后，随着瓦斯的不断释放，瓦斯压力降低，故有效应力增加，所以，孔隙被压缩，所以孔隙率随着瓦斯压力的降低而不断降低。下图给出了5个特征点的孔隙率的变化曲线。与瓦斯压力的分布曲线105瓦斯渗流过程中的煤样的孔隙率变化（断裂带中的点C和D）瓦斯渗流过程中的煤样的孔隙率变化（断裂带中的点C和D）106但是，如下图所示，在煤基质中，t=1e-02s时，首先是由于煤的变形导致了孔隙率下降，但是，随着瓦斯从裂隙带向着煤基质中渗流，故基质中的瓦斯压力会逐渐，所以导致了孔隙率的增加。此后，由于瓦斯开始集中向抽放孔中流动，整个煤样中的瓦斯压力不断降低会引起有效应力的增加，故孔隙率会不断降低。但是，如下图所示，在煤基质中，t=1e-02s时，首先是107瓦斯渗流过程中的煤样的孔隙率变化(煤基质中的点A、B和E)瓦斯渗流过程中的煤样的孔隙率变化(煤基质中的点A、B和E)108主要结论如下：（1）由于裂隙带中的瓦斯压力比煤基质中的高，故瓦斯由裂隙带向外围的煤机制中不断扩散，使得裂隙带和基质间的瓦斯压力剃度不断降低。在某个时间后，瓦斯开始集中向抽放孔中流动。（2）在本文给定的模拟参数下，应力场引起的煤体压缩是对其渗透性的影响要大于Klinkerberg效应和瓦斯解吸效应所引起的渗透率变化。故总体上煤层中表现出渗透率的降低，随着瓦斯的不断抽放，渗透率更是不断降低。主要结论如下：109二、冒落区瓦斯浓度扩散-对流及

风流场数值模拟瓦斯浓度扩散-对流场基本方程：Darcy方程二、冒落区瓦斯浓度扩散-对流及

风流场数值模拟瓦斯浓度扩散-110冒落区破碎岩体压实过程中的气体流场非Darcy方程(堆石体、土石坝流场)冒落区破碎岩体压实过程中的气体流场非Darcy方程(堆石体、111冒落区瓦斯浓度扩散-对流

及风流场数值模拟通风流场基本方程：

Navier-Stokesequations:Freeflow适合通风巷道风流场式中：v——流体流速，m/s；p——流体压力，Pa；ρ——流体密度，kg/m3；I——单位矢量；F——流体阻力冒落区瓦斯浓度扩散-对流

及风流场数值模拟通风流场基本方程：112冒落区瓦斯浓度扩散-对流

及风流场数值模拟通风流场基本方程：Brinkmanequations:Fastflowinporousmedia，适合冒落区风流场冒落区瓦斯浓度扩散-对流

及风流场数值模拟通风流场基本方程：113模型建立计算模型和方案参照综采工作面的具体尺寸，建立如下图所示的二维计算模型，不考虑势能，模型东西向即工作面推进方向取400m,南北向宽240m。划分成均质（5个分区，沿工作向采空区方向的区域（1）至（5）宽度均取为80m）；即大约1周的推进时间瞬时完成。假设这个期间内每个区域的透气系数分布均不同，具体采用的计算模型相关条件如下：模型建立计算模型和方案114模型建立不透气边界不透气边界浓度补给边界采空冒落区

（5）（4）（3）（2）（1）进风口出风口推进方向进风口进风口进风口进风口出风口出风口出风口出风口浓度补给边界模型建立不透气不透气浓度补给边界采空冒落区进风口出风口推进方115模型建立边界条件：（1）通风条件：左侧下20m为进风口边界，左侧上20m为回风口边界，压力差100Pa,其他边界为不透气边界。（2）扩散条件：右侧边界为绝缘对称边界，其他边界有补给，推进区域1时，上下边界补给量为1.2e-6mol/m2·s，汇源项为3e-6mol/m2·s，以描述瓦斯通量随开采动态过程而增加的瓦斯量。模型建立边界条件：116模型建立初始条件：域内具有一个大气压，瓦斯初始浓度3mol/m3；推进新的工作面区域时，旧区域瓦斯浓度模拟结果C0再附加一初始补充浓度Cbc（取3mol/m3）作为旧推进工作面的初始浓度。模型建立初始条件：117模型建立时间步长：按照非均等的积数步长增大，初始值7s,终止指7e5s，设定100个中间时间值。计算参数：动粘系数=1.8e-5pa·s，流体密度=1200kg/m3，扩散系数D=2e-5m2，瞬态时间比例系数=0.55模型建立时间步长：118模型建立1234513.24e-8

22.88e-83.24e-832.52e-82.88e-83.24e-842.16e-82.52e-82.88e-83.24e-851.8e-82.16e-82.52e-82.88e-83.24e-8已推区域透气率（m2）新推区域随推进进行透气率变化表模型建立1234513.24e-822.88e-83.24119模拟结果分析推进区域1时的计算结果下图分别是当时间为：7s，7e3s，7e4s，1.1e5s，2.5e5s，3.5e5s，7e5s瞬态时间区域1采空区瓦斯浓度分布图模拟结果分析推进区域1时的计算结果120模拟结果分析Time=7s时瓦斯浓度

模拟结果分析Time=7s时瓦斯浓度121模拟结果分析Time=7e3s时瓦斯浓度

模拟结果分析Time=7e3s时瓦斯浓度122模拟结果分析Time=7e4s时瓦斯浓度

模拟结果分析Time=7e4s时瓦斯浓度123模拟结果分析Time=1.1e5s时瓦斯浓度

模拟结果分析Time=1.1e5s时瓦斯浓度124模拟结果分析Time=2.5e5s时瓦斯浓度

模拟结果分析Time=2.5e5s时瓦斯浓度125模拟结果分析Time=3.5e5s时瓦斯浓度

模拟结果分析Time=3.5e5s时瓦斯浓度126模拟结果分析Time=7e5s时瓦斯浓度和流线

A1A1’模拟结果分析Time=7e5s时瓦斯浓度和流线A1A1’127模拟结果分析第一步推进时A1-A1’切面浓度变化曲线

模拟结果分析第一步推进时A1-A1’切面浓度变化曲线128模拟结果分析推进区域2时的计算结果下图分别是当时间为：7s，7e3s，1.1e5s，7e5s瞬态时间区域2采空区瓦斯浓度分布图模拟结果分析推进区域2时的计算结果129模拟结果分析Time=7s时瓦斯浓度模拟结果分析Time=7s时瓦斯浓度130模拟结果分析Time=7e3s时瓦斯浓度模拟结果分析Time=7e3s时瓦斯浓度131模拟结果分析Time=1.1e5s时瓦斯浓度模拟结果分析Time=1.1e5s时瓦斯浓度132模拟结果分析第二步推进时Time=7e5s时瓦斯浓度和流线

A2A2’模拟结果分析第二步推进时Time=7e5s时瓦斯浓度和流线133模拟结果分析第二步推进时A2-A2’切面浓度变化曲线

模拟结果分析第二步推进时A2-A2’切面浓度变化曲线134模拟结果分析推进区域3,4时的计算结果模拟结果分析推进区域3,4时的计算结果135模拟结果分析第三步推进时Time=7e5s时瓦斯浓度和流线

A3A3’模拟结果分析第三步推进时Time=7e5s时瓦斯浓度和流线136模拟结果分析第三步推进时A3-A3’切面浓度变化曲线

模拟结果分析第三步推进时A3-A3’切面浓度变化曲线137模拟结果分析第四步推进时Time=7e5s时瓦斯浓度和流线

A4A4’模拟结果分析第四步推进时Time=7e5s时瓦斯浓度和流线138模拟结果分析第四步推进时A4-A4’切面浓度变化曲线

模拟结果分析第四步推进时A4-A4’切面浓度变化曲线139讨论（1）在区域1右侧不透气边界附近，随着工作面向左推进，风流作用逐渐减弱，该处瓦斯浓度不断集聚增大，每个推进步瓦斯浓度分别为7mol，12mol，17mol，24mol，表明该处距离通风口越远，瓦斯集聚程度越大。讨论（1）在区域1右侧不透气边界附近，随着工作面向左推进，风140讨论（2）对于每个推进步，当时间达到7e5s（8d左右）时，风流基本上可以把本区域瓦斯浓度降低到安全范围内（接近0mol），表明风流在本区域的流速明显大于前一个区域。而前一个区域的下左半部分（占整个区域的1/4）的瓦斯浓度也能得到有效降低，而对其他区域瓦斯浓度驱散作用影响有限。讨论141讨论（3）图9可见，瓦斯浓度沿工作面推进方向从近到远瓦斯浓度成台阶状渐次增大，这表明通风量只能在一定范围内降低瓦斯浓度，同时也看出距离通风口越远，瓦斯补给的时间和补给量越大。讨论142讨论（4）虽然模型上下边界都是瓦斯补给边界，但由于模型左下边界是进风口，瓦斯随风流对流作用得到有效降低；而模型左上边界是出风口，瓦斯容易在回风隅角大量聚集，所以该处瓦斯浓度较高。讨论143三、三维渗流耦合模型及瓦斯抽放根据实际的三维煤层瓦斯抽放过程，建立如图1所示的长宽高为10m×10m×10m理想化的三维计算模型，模型左下部边界假设为巷道，三个瓦斯抽放孔——K1、K2、K3之间的距离1m，与水平面呈45°倾角展布，K1、K3与K2呈16°夹角分布在孔2两侧，三个瓦斯抽放孔长度7m，瓦斯抽放孔按照“以缝代孔”[16]原则简化为定压力边界（25kPa）。三、三维渗流耦合模型及瓦斯抽放根据实际的三维煤层瓦斯抽放过程144巷道所在位置45º16ºL1L2K1K3K2外部载荷巷道所在位置45º16ºL1L2K1K3K2外部载荷145(1)边界条件：四周和上下面均为零通量不透气边界；四周和下部都约束法线方向的位移，上部自由，作用有P分别为10MPa、20MPa、2MPa、0.1MPa的外部载荷，同时模型具有自重载荷。(2)初始条件：内部有1MPa的初始瓦斯压力，三个抽放孔的压力为0.25e5Pa。(1)边界条件：四周和上下面均为零通量不透气边界；四周和下146(3)时间步长：按照非均等的积数步长增大，初始值1s，终止指1e7s（10d左右），设定100个中间时间步长。(4)计算方案：模拟不同外部载荷条件下（2MPa、10MPa、20MPa、0.1MPa），瓦斯抽放效果和渗透性变化规律。(5)计算参数：相关参数列表1所列。(3)时间步长：按照非均等的积数步长增大，初始值1s，终止147杨氏模量E（Pa）泊松比ν煤层密度ρs（kg/m3）流体密度ρ1（kg/m3）流体粘滞性（Pa·s）

饱和渗透率（m2）7.41e90.33125010001e-35e-6杨氏模量泊松比煤层密度流体密度流体粘滞性饱和渗透率7.41e148模拟结果分析推进区域2时的计算结果下图分别是当时间为：1s，1.28e5s，1.08e6s，1e7s瞬态时间渗透性系数分布图和压力等表面图模拟结果分析推进区域2时的计算结果149渗透性系数分布图Time＝1s时渗透性系数分布图渗透性系数分布图Time＝1s时渗透性系数分布图150渗透性系数分布图Time＝1.28e5s时渗透性系数分布图渗透性系数分布图Time＝1.28e5s时渗透性系数分布图151渗透性系数分布图Time＝1.08e6s时渗透性系数分布图渗透性系数分布图Time＝1.08e6s时渗透性系数分布图152渗透性系数分布图Time＝1e7s时渗透性系数分布图渗透性系数分布图Time＝1e7s时渗透性系数分布图153压力等表面图Time＝1e7s时渗透性系数分布图Time＝1s时压力等表面分布图压力等表面图Time＝1e7s时渗透性系数分布图Time＝154压力等表面图Time＝1s时压力等表面分布图Time＝1.28e5s时压力等表面分布图压力等表面图Time＝1s时压力等表面分布图Time＝