煤层瓦斯渗流与煤体变形的耦合数学模型及数值解法

第15卷　第 2期 岩石力学与工程学报 V ol . 15 N o . 2 1996 1996年 6月 Chinese J ournal of Rock Mechanics and Engineering 15( 1996) , 135—142 煤层瓦斯渗流与煤体变形的 耦合数学模型及数值解法* 　　　　　　　　梁　冰　章梦涛　　　　　　　　　　　王泳嘉　　　　　 　　　　　　　　　(阜新矿业学院　阜新　123000)　　　　　　　　(东北大学　沈阳　110006) 提要　根据固流耦合作用的基本理论, 考虑瓦斯渗流对煤体本构关系的影响 , 提出了煤和瓦 斯耦合作用的数学模型, 给出了用有限元求其数值解的方法及程序。在此基础上, 对采动影 响下考虑煤层变形影响时瓦斯在采空区的流动规律进行了数值模拟分析, 为采空区瓦斯抽放 提供了理论基础和科学依据。 关键词　瓦斯渗流, 固流耦合作用, 数学模型 1　引言 在煤炭开采过程中, 煤岩体变形与其中孔隙流体的流动相互影响、相互作用即固流耦 合作用的问题很多。实际上, 在含瓦斯煤的开采过程中煤层变形和瓦斯在煤层中的流动都 是在固流耦合作用下的煤层变形和瓦斯流动[ 1]。但是过去研究开采过程中的瓦斯流动和含 瓦斯煤的变形破坏等问题时, 只是根据传统的渗流力学或岩石力学的理论进行研究, 因此 研究结果与实际情况存在很大的误差。 近年来, 学者们开始应用固流耦合作用理论研究煤体变形及煤层中的瓦斯流动[ 2]。瓦 斯对煤层变形不仅有孔隙压力的力学作用, 而且还有煤对瓦斯的吸附解吸对煤的力学性质 所产生的影响 [ 3]。本文将在考虑瓦斯吸附变化对煤本构关系影响的基础上建立煤层瓦斯吸 附变化对煤体变形耦合作用的数学模型, 并给出其有限元的数值解法和程序。最后利用该 程序, 对四川打通二矿 7#煤层瓦斯在采空区的流动规律进行数值模拟分析, 其结果可为采 空区瓦斯抽放提供理论基础和科学依据。 2　瓦斯渗流与煤体变形耦合作用的数学模型 2. 1　瓦斯在煤层中流动的渗流场方程 煤层开采时, 一部分瓦斯以游离状态赋存于煤层的裂隙和孔隙中, 一部分以吸附状态 附着于煤微孔隙的 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 面。在一定的孔隙压力和温度下两者处于动平衡状态。目前研究瓦斯 在煤层中的运移有两种方式: 一是在裂隙和孔隙中的瓦斯渗流流动; 二是吸附瓦斯从微孔 1994年 11月 24日收到初稿, 1995年 4月 21日收到修改稿。 * 本项目得到国家自然科学基金的资助(批准号 19272019)。 面的解吸, 并由微孔中向裂隙和孔隙中的扩散。但从宏观角度看, 在开采过程中瓦斯在煤 层中的运动可以视为其在煤层中的渗流流动 [ 4]。瓦斯在煤层孔隙和裂隙中的流动极为复 杂, 可采用渗流力学的研究方法对其进行研究。根据大多数研究, 可将煤层简化为由裂隙 和孔隙共同组成的渗流通道的多重介质模型。由于裂隙渗透率远比孔隙渗透率高, 在采动 影响下这种情况更加严重, 因此瓦斯在可变形煤层中的流动可以视为在裂隙网络中的流 动。研究固流耦合作用下的瓦斯渗流, 显然其渗流方程是非线性的, 但是在孔隙压力变化 的微段内是可以视为线性的, 服从达西定律, 即有 qi = - ( K ijU) ′i ( 1) 式中 qi 为瓦斯渗流速度分量, K ij 为渗透系数张量, U为压力势。由于瓦斯密度很小, 重力 产生的势可忽略不计, 则有 U= p / (Qg) ( 2) 式中 p 为瓦斯压力, Q为瓦斯密度。 将式( 2)代入式( 1)中得 qi = - ( K ij Qg p ) ′i ( 3) 对于裂隙介质的渗透系数 K ij 尽管有许多理论计算公式, 但是在实际应用中特别在考虑变 形影响时, 只能通过宏观试验的方法确定。目前对此研究还只限于峰值强度前, 本文采用 的是文[ 2]中所提供的关系式: K ij = a0 exp( - a1 (′) ( 4) 式中 ( ′为有效体积应力, a0, a1为试验常数。 煤层中瓦斯含量为吸附瓦斯量 W c和游离瓦斯量W CB之和。吸附量W C由拉格谬尔公式 给出, 即 W c = abp 1 + bp r mQ 游离瓦斯量 W CH由下列公式给出 W CB = G p p 0 Q 则总含量W 为 W = W c + W CB = abp 1 + bp r mQ+ G p p 0 Q ( 5) 式中: W c为单位体积煤体中的吸附瓦斯量( g/ cm3 ) ; a为煤体的最大吸附瓦斯量( cm3 / g) ; b为煤体的吸附常数( cm2 / g ) ; r m为煤体的重率( g / cm3 ) ; W CB 为 单位体积中游离瓦斯的含 量( g/ cm3 ) ; p 0为 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 气压, p 0 = 0. 1013M Pa; G为孔隙率。 由于瓦斯为可压缩气体, 需要补充一个气体状态方程。设瓦斯为理想气体, 流动过程 为等温过程, 则状态方程为 Q= p / ( RT ) ( 6) 式中 T 为瓦斯气体的温度, R 为瓦斯气体常数。采用连续介质力学方法, 可以得到瓦斯质 量守恒方程为 (Qqi ) ′i = ŒWŒt ( 7) ·136· 岩石力学与工程学报 1996 年 　　将坐标系选在渗透主轴上, 并将上式综合可得可压缩瓦斯气体在煤层中流动的渗流方 程 Œ Œx ( K x Œp 2 Œx ) + Œ Œy ( K y Œp 2 Œy ) + Œ Œx ( K z Œp 2 Œz ) = [ 2G p 0 + abrm ( 1 + bp ) 2 + abr m ( 1 + bp ) ] Œp 2 Œt ( 8) 2. 2　变形场方程 变形场方程由以下三部分组成。 2. 2. 1　平衡方程 由表征单元体处于应力平衡状态, 可以得出按总应力 Rij 表示的平衡方程 Rij′j + f i = 0 ( 9) 式中 f i 为体积力。 设煤体骨架的有效应力由太沙基公式给出 Rij′= Rij - pDij ( 10) 　　将式( 10)代入式( 9)中可以得到按有效应力 Rij′表示的平衡方程为 ( R′ij )′j + ( p Dij )′j + f i = 0 ( 11) 2. 2. 2　几何变形方程 设岩体骨架所发生的变形是小变形, 则几何方程为 Eij = 12 ( ui, j + uj , i ) ( 12) 式中 Eij , ui 分别为应变分量和位移分量。 2. 2. 3　本构方程 考虑瓦斯对煤体的作用、煤体在塑性变形中显著的扩容现象以及含瓦斯煤体的应变软 化特性, 采用内时理论建立煤和瓦斯耦合作用下煤的内时本构方程[ 5]。 其中体积响应为 Rkk′= r 1Ehkk + r2 [ 1 - exp( - B1Ehkk ) ] + a1( Ekk a1 ) a2　( Ekk ≤ a1) ( 13a) Rkk′= r 1Ehkk + r2 [ 1 - exp( - B1Ehkk ) ] + a5( a1 - a4 ) a4 - a5 ( Ekk - a2 ) a4 + 　　　　　　 + R1　　　( a2 ≤ Ekk ≤ a1) ( 13b) 　　Rkk′= r 1Ehkk + r2 [ 1 - exp( - B1EhK K ) ] + a7 a6 [ ea6( Ekk- a2) - 1] + R 2　( EKK ≤ a2) ( 13c) 式中 Ehkk 为静水压力作用下的体积应变, r1、r 2、B1为材料常数, 由静水压缩三轴试验确定。 a1、a2 分别为剪胀发生及峰值强度处的体积应变值, R1、R 2为此两处的体积应力值。a3、 a4、a5、a6、a7为材料常数, 根据体积应力—体积应变曲线确定。 偏斜响应: S ij ′ =∫z soQ( z s - z′s) de pijdz′sdz′s ( 14) 式中, depij = deij - dS ij / ( 2G0 ) Q( Z) = ∑3 r = 1 C re - A r Z ·137·第 15 卷　第 2 期 梁冰等: 煤层瓦斯渗流与煤体变形的耦合数学模型及数值解法 　　取　dz s = ‖deij - dS ij / ( 2G0 )‖/ f s (Fs) 　　 f s(Fz ) = 1 + BFz 式中的所有常数 C1、C2、C3、A1、A2、A3、B表示含瓦斯煤在变形过程中瓦斯对煤的力学性质 的影响, 由单轴压缩试验确定。 以上数学方程, 加上边界条件和初始条件构成了煤体变形和瓦斯渗流流动的数学模 型。 3　数学模型的数值解法 煤体变形和瓦斯渗流的数学模型是相当复杂的, 对于具体工程实际问题其解析解也是 十分困难的, 甚至是不可能的, 因此只能借助于计算机求其数值解。 3. 1　瓦斯渗流场方程数值解的有限元方法 3. 1. 1　瓦斯渗流场方程的非线性处理 设　　　　　 　　　S( p ) = [ 2G p 0 + abr m ( 1 + bp ) 2 + abrm ( 1 + bp ) ] P ( 15) 则可将式( 8)改成如下形式 Œ Œx ( K xp Œp Œx ) + Œ Œy ( K yp Œp Œy ) + Œ Œz ( K z p Œp Œz ) = S( p ) Œp Œt ( 16) 　　采用迭代法对上式进行非线性处理。假设压力 { p } 的初值为{ p 0} , 将式( 16)展开得 ( K xp Œ2p Œx 2 + K yp Œ2p Œy 2 + K zp Œ2p Œz 2 ) ûp o + + [ ŒŒx ( K xp ) + ŒŒy ( K yp ) + ŒŒz ( K zp ) ] ûp o = S ( p 0) ŒpŒt ( 17) 　　忽略二阶微量 ŒŒx ( K xp ) ûp 0等, 且假设煤介质的渗透性为各向同性, 即 K x = K y = K z = K 代入式( 17)中得 K p 0 ( Œ2p Œx 2 + Œ2p Œy 2 + Œ2p Œz 2 ) = S ( p 0 ) Œp Œt ( 18) 　　对方程( 18)求解, 可得各节点的瓦斯压力 { p 1} , 再以{ p 1} 代替{ p 0} , 方程( 18) 的右端 项以{ p 0} 代入 S ( p ) 中, 重新求解方程( 18) , 反复多次, 直到满足要求的精度为止。 若瓦斯的渗流为非稳定的, 则需要对时间进行离散, 即把从初始时间 t0 到最终时间 tn 的这段时间分成 n段, 每段时间的步长为 $t , 当 $ t较小时, 可以认为Œp /Œt在该时间步长 内保持不变, 且等于( p i - p i- 1 ) / $ t, p i- 1是上一个时间步长终了时的瓦斯压力值, p i 是本 次时间步长终了时的瓦斯压力值, 则在 $ t时间内可将非稳定的瓦斯渗流变为稳定的渗流 问题。 3. 1. 2　瓦斯渗流场离散方程的建立 建立瓦斯渗流场的离散方程可利用泛函求极值的方法。设在边界 S1 上瓦斯的压力给 定, 在边界 S 2上瓦斯的流量给定, 即 p = p g　　在 S 1上 K x Œp Œx l x + K y Œp Œy l y = - q0　　在 S2上 ·138· 岩石力学与工程学报 1996 年 则其相应的泛函方程为 I ( p ) = 1 2∫v {Kp 0[ (ŒpŒx ) 2 + (ŒpŒy ) 2 + (ŒpŒz ) 2] + 2S 0ŒpŒt p } dv +∫sq0p dS ( 19) 　　对于空间问题, 一般采用空间 8节点等参单元。对 p 求一阶导数得 DI ( p ) = ∑M e (∫ve[ B c] T[ K c] [ B c] ûJûdNdGdF{ p } + + 2∫veS ( p 0) ŒpŒt [ N ] TûJ ûdNdGdF+∫S 2 q0 [ N ] T dS ) ( 20) 式中的求和是对单元进行的, M 为所划分区域的单元总数, [ K c] 为渗透系数矩阵。单元的 传导矩阵为 [ T ] e =∫∫∫ve[ Bc] T [ K c] [ Bc] ûJûdNdGdF ( 21) 单元的贮量矩阵为 [ S ] e =∫∫∫veS ( p 0) [ N ] TûJ ûdNdGdF ( 22) 单元的列矢量由边界条件形成。对于第一类边界条件, 即边界上的压力给定, 则 {H } e = [ N ] T p g = {N 1p g　N 2p g　⋯　⋯　N 8p g} T ( 23) 对于第二类边界条件, 即边界上流量给定, 设此边界为 N= ± 1的 GF面上, 则 {H } e =∫- 1- 1∫1- 1[ N ] TN= ±1q 0dA GF =∫- 1- 1∫1- 1q0[ N ] TN= ±1 ( EG%F- E2GF) N= ±1dGdF ( 24) 式中E G = (ŒxŒG) 2 + ( Œy ŒG) 2 + ( Œz ŒG) 2 E F= ( ŒxŒF) 2 + (ŒyŒF) 2 + (ŒzŒF) 2 E GF= ŒxŒGŒxŒF+ ŒyŒGŒyŒF+ ŒzŒGŒzŒF 离散后的代数方程组为 { [ T ] $ t + [ S] } { p } t j - 1 + $t = [ S ] { p } tj - 1 + {H }$ t ( 25) 当不考虑瓦斯压力随时间变化即瓦斯在煤层中的流动为稳定渗流时, 方程变为 [ T ] { p } = {H } ( 26) 3. 2　煤岩体变形场离散方程的建立 由于煤岩体中孔隙瓦斯的存在, 可将瓦斯对煤体的作用分为两部分处理。纯力学作用 以孔隙压力梯度分量的值作为体积力作用于煤体中, 非力学作用在煤岩体的内时本构方程 和瓦斯渗流方程的源汇项中已经予以考虑。体积力 {Q} 应包含有孔隙压力所产生的部分, 设 {Q} = {Q 1} + {Q 2} ( 27) 式中 {Q1} 为真正的体积力, {Q2} 为孔隙压力梯度产生的体积力。单元体内任一点的孔隙 压力为 p , 且 p = p ( x , y , z ) , 则由孔隙压力梯度ŒpŒx、 Œp Œy、 Œp Œz 所产生的体积力的等效节点 ·139·第 15 卷　第 2 期 梁冰等: 煤层瓦斯渗流与煤体变形的耦合数学模型及数值解法 载荷为 {R} e =∫∫∫ve[ N ] T{ŒpŒx ŒpŒy ŒpŒz } T dV ( 28) 　　煤岩体变形场的离散方程为 [ K ] { D} = {F} ( 29) 3. 3　耦合数值计算的程序框图 根据方程( 26)和方程( 29)可以得到煤体变形和瓦斯渗流耦合作用的离散方程为 [ T ] { p } = {H } [ K ] {D} = {F} ( 30) 对其进行耦合求解的程序框图如图 1所示。 图1 耦合程序的计算框图 Fig. 1 Flow chart of numerical calculation of the coupling program 开 始 读入计算所需的数据, 给变量赋初值 调用变形场子程序def, 计算煤岩体的变形和应力 调用渗流场子程序ry1, 计算瓦斯渗流场的压力 送代循环 改变初值 调用变形场子程度def 根据应力计算各个单元的渗透系数 调用渗流场子程度ry1 判断是否满足收敛要求 计算结果写入数据文件 ·140· 岩石力学与工程学报 1996 年 4　计算实例 根据所编制的三维空间等参有限元的耦合计算程序, 结合四川打通二矿在 7# 煤层 N1709工作面开采时对上、下邻近煤层卸压后瓦斯向开采层采空区流动状况进行了数值模 拟。该工作面已推进了 310m, 工作面总长 240m, 中间为进风巷, 在东西两侧为回风巷。模 拟空间示意图如图 2所示, 取 7# 上方 26m 的煤线为上边界, 取 12# 煤层底板为下边界。 图 2　采动空间瓦斯流动模拟空间示意图 Fig . 2　Diagr am of g as mig rat ion in mining influenced space 图 3　8# 煤层瓦斯压力变化曲线 F ig . 3　Curve of g as pressur e in No . 8 coa l seam 4　采动影响下邻近煤层瓦斯压力的分布及流向形态图 am o f gas mig r ation and ga s pr essure distr ibution near coal seam by mining influence 边界条件为 ( 1) 四周固体边界, 瓦斯无流入流出即 ( K Œp Œx i ) Giûs = 0 　　( 2) 下部边界取 12# 煤层的原始瓦斯压力为 0. 8M Pa, 到工作面煤壁开始卸压, 经 40m 后瓦斯压力降至 0. 3MPa; ·141·第 15 卷　第 2 期 梁冰等: 煤层瓦斯渗流与煤体变形的耦合数学模型及数值解法 ( 3) 上部边界取顶部煤线原始瓦斯压力为 0. 8M Pa; 取从工作面煤壁位置开始下降, 经 30m 后瓦斯压力下降至 0. 2MPa。 在上述条件下进行计算, 并使 8#煤层实测瓦斯压力变化曲线作为检验结果如图 3所 示。计算所得到的采动影响下邻近煤层瓦斯压力分布及流向形态图如图 4所示。 由采动影响下邻近层瓦斯压力的分布可见: ( 1) 在距工作面 30m 以内部位, 各个邻近层瓦斯压力从原始瓦斯压力值迅速下降, 且 距开采层越近, 压力变化及变化区域越小, 各层瓦斯压力梯度负方向指向距工作面 10m～ 50m 距离的采空区处。由压力分布与流线正交关系, 各邻近层中瓦斯渗透沿等压力线方向 流动。由于卸压初始阶段邻近煤层部位的瓦斯压力变化大, 所释放出来的瓦斯量多。这些 部位的瓦斯沿等压力线方向流向采空区, 在距工作面 10m～50m 处最为强烈, 构成瓦斯向 采空区流动的高峰区域。 ( 2) 距工作面 50m 以外的区域, 各邻近层瓦斯压力均下降至最小值, 且保持稳定。压 力等值线与煤层平行, 压力梯度变小为卸压区域压力值的 1/ 5左右, 这一区域的瓦斯流量 较小。 参 考 文 献 1　周世宁: 瓦斯在煤层流动的机理. 煤炭学报, 1990; 16( 1) : 17- 23 2　章梦涛, 梁冰, 梁栋: 采动影响下煤层内瓦斯流动状况的数学模型及数学分析. 见: 第二届全国岩石力学数值计算与 模型试验学术研讨会论文集, 上海: 同济大学出版社, 1990; 423—428 3　梁冰, 章梦涛, 王泳嘉: 瓦斯对煤的力学性质及力学响应影响的试验研究. 岩土工程学报, 1995; 17( 5) : 12- 18 4　Saghaf i A . : 煤层瓦斯流动的计算机模拟及其在预测瓦斯涌出和抽放瓦斯中的应用. 见: 第二十二届国际采矿安全 会议论文集, 北京: 煤炭工业出版社, 1987; 98- 103 5　梁冰, 章梦涛, 王泳嘉: 含瓦斯煤的内时本构模型. 岩土力学, 1995; 16( 3) : 22- 28 MATHEMATICAL MODEL AND NUMERICAL METHOD FOR COUPLED GAS FLOW IN COAL SEAMS AND COAL DEFORMATION 　　　　Liang Bing　Zhang M engtao　　　　　　　Wang Yongjia 　　　(F ux in M ining Institute,　Fux in　123000)　(N or theast Univ er sity , 　Shenyang　110006) Abstract　A coupling mathematical model of so lid and f luid for g as f low in coal seams and coal defo rmat ion is pr esented; and the numerical method of the mathematical model and FEM progr am are g iven. T he r egularity o f g as f low in the mined-out area inf luenced by the deformat ion of co al seam is obtained by using the numer ical method and the pro gram . The g iven results provide a theoret ical basis and scient ific data for the gas ext raction f rom the mined-out ar ea. Key words　gas f low , solid-fluid coupling effect , mathemat ical model ·142· 岩石力学与工程学报 1996 年