带压开采底板破坏规律的三维数值模拟研究.pdf

第 35 卷 第 4 期 2004 年 7月 太原理工大学学报 JOURNAL OF TAIYUAN UNIVERSITY OF T ECHNOLOGY Vol. 35No. 4 July2004 文章编号 10079432 2004 04040004 带压开采底板破坏规律的三维数值模拟研究 刘志军, 胡耀青 太原理工大学 采矿工艺研究所, 山西 太原 030024 摘 要 从带压开采的角度建立了带压开采的三维固流耦合数学模型及其数值解法, 并针对太 原市东山煤矿带压开采, 进行了三维固流耦合数值模拟, 对带压开采底板的应力分布规律进行了详 细的探讨, 为进一步防治矿井突水提供理论依据。 关键词 带压开采; 固流耦合; 有限元; 数值模拟 中图分类号 TD82383 文献标识码 A 我国煤炭资源丰富, 是能源组成的重要部分, 然 而我国灰岩岩溶 Karst 分布广泛, 随着浅部煤层的 枯竭, 开采强度的增加, 矿井的深度不断延伸, 来自 煤层底部奥陶纪灰岩高承压水的危害日趋加剧。据 统计有 60 的煤矿不同程度地受到底板岩溶承压 水的威胁, 仅华北地区受其危害的矿井就有 230 多 个, 造成 40 左右的煤炭不能正常开采。因此带压 开采的研究多年来一直是我国乃至世界矿业工作者 关注的问题。关于带压开采的研究, 前人已作了大 量的研究, 从现场实测, 实验模拟到数值模拟都进行 了研究, 就数值模拟而言, 大多是二维模拟, 且没有 考虑耦合作用, 所以本文从固流耦合的观点出发, 建 立固流耦合三维数学模型及其相应的数值解法, 并 对开采过程中底板的应力分布规律进行探讨。 1 三维固流耦合数学模型及数值解法 1. 1 基本假设 在建立模型之前, 先引入如下的假设 1 固体围岩场的应力应变服从线性规律, 即简 化为弹性过程, 应力与变形遵守广义虎克定律, 其主 要原因是由于三维固流耦合分析计算工作量较大, 计算时间较长, 如采用非线性规律计算, 对微型机来 说难以承受, 考虑到工程的要求, 以分段线性的方法 来代替非线性, 其解足可以满足工程的要求。 ij ije 2 ∀ij. 式中 ij为总应力张量; ij为 kronecker 记号; , 为拉梅常数; ∀ij为总应变; e 为体积应变。 2 流场水压力对围岩的作用按内力处理, 即把 水压力分别作用到对应的结点上。 3 流场中的围岩被单相水饱和, 且水为不可压 缩体。 4 围岩体内微段压力梯度上 , 流体渗流服从 线性达西定律 , 即 q k , p p x . 式中 q 为流量; k , p 为渗透系数; 为体积应 力; p 为水压。 5 围岩固体骨架变形的有效应力规律遵守修 正的太沙基有效应力规律, 即 ij ij ∃ p ij. 式中 ∃为等效孔隙压力系数; ij为有效应力张量; ij为总应力张量; ij为 kronecker 记号。 6 裂隙渗流服从裂隙流的达西定律, 即 q kf p s . 式中 kf b2/ 12, b 为裂隙宽度, 是水压和应力的 函数; s 为裂隙切向坐标。 7 裂隙单元的变形服从 Goodman 节理单元变 形控制方程[ 1] n Kn∀n, s Ks∀s, n n- p. 8 裂隙单元有效应力规律遵循太沙基有效应 力规律, 即 - ∃ p, 收稿日期200309 02 基金项目 省基金 20021076, 20001012 和省青年基金资助项目 190- 12903306 作者简介 刘志军 1978- , 男, 内蒙集宁人, 在读硕士, 主要从事带压开采及水害治理研究。 式中 , 分别为裂隙法向总应力与有效应力。 1. 2 数学模型 在以上假设的基础上, 可以给出完整围岩的固 流耦合模型, 渗流方程 Kx 2p 2x Ky 2p 2y Kz 2p 2z S p t e t W. 1 变形方程 Uj , j i Ui,j j Fi ∃ p , i 0. 2 式中 s 为贮水系数; w 为源汇项; e 为体积变形; p 为水压; ∃为等效空隙压系数。 同样, 可得裂隙介质的固流耦合模型, 渗流方程 Kx 2p 2x Ky 2p 2y Sn p t e t W. 3 变形方程 n n- p b/ Dn, x sb/ Ds. 4 式中 e n s为裂隙的体积变形; Sn为裂隙的 贮水系数; Ds, Dn为无厚度节理裂隙单元单位长度 的切向刚度和法向刚度。 方程1 4组成带压开采三维固流耦合的数 学模型, 可以看到, 方程 1包含了体积变形, 只能从 方程 2中得到, 而方程 2 中包含了流场中的水压 力, 只能从方程1 中得到, 方程3, 4 也一样, 所 以为固流耦合方程。 1. 3 数值解法 关于固体变形方程的数值方法, 可参见文献 [ 2] , 不作具体分析, 以下结合文献[ 3] 的方法重点考 虑渗流方程的数值方法。 对方程 Kx 2p 2x Ky 2p 2y Kz 2p 2z S p t e t W, 采用剩余试函数法, 令其试探解为 p [ N ] { pn} Nipi, 其中 i 为结点数。 则上式可变为 Kx 2[ N] { pn} 2x Ky 2[ N] { pn} 2y Kz 2[ N] { pn} 2z - S [ N] { pn} t - e t - W Re 0. 5 式中, Re为方程的剩余值, 它是由于试探解引起的 误差, 现要求剩余值 Re的加权平均 Re在单元 e 内 域 上的任一点为零, 若令 W 为权函数, 则 Re WRed WRed 0.6 式6等于 0 的条件是分子为 0, 即要求 WRed 0.7 由于权函数和剩余函数是共轭的, 对第 i 个结 点而言, 权函数为 Wi, 选定 Wi的任意性很大, 与有 限元相结合, 可取 Wi Ni, 对典型单元 e 来说, 上 式可写为 e NiRed 0.8 式8为有限元法与剩余试函数法相结合的表达式, 把式5代入式8得 e Ni[ Kx 2[ N ] { pn} 2x Ky 2[ N ] { pn} 2y Kz 2[ N] { pn} 2z - S [ N] { pn} t - e t - W] d 0. 9 由分步积分和格林公式得 e NiKx x p x dx dydz Gij ∀ e SNiNjdxdydz ; Qi e q e t WNidx dy dz . 计算时, 可采用有限差分法, 即11式可写为 t[ ∋] { pt t} [ G] {pt t- pt} t{ Q}. 对于 e/ t 项也可以采用同样得方法, 其它计算方 法可参见文献[ 4] 。 2 有限元计算实例 采用如上的固流耦合数学模型及计算方法, 对 太原市东山煤矿一采区带压开采进行数值模拟, 其 计算模型简化为图 1, 材料力学特性参见文献[ 5] . 图 1 计算模型图 2. 1 边界条件的处理 根据计算区域及现场实际条件, 顶面取均布载 荷, q 100 kg/ cm2 相当于 400 m 的采深 , 底面为 固定边界条件, 前面与左面为侧向约束条件, 如在 XOZ 面, Y 方向的位移固定, 而其它方向可以自由 移动。另外两个面为给定应力边界条件, 如后面给 定其 X 向的载荷为 50 kg/ cm2, 即取侧压系数为 0. 5, 右面给定其 Y 方向的载荷也为 50 kg/ cm2. 2. 2 计算方法 根据实验模型, 把图 1 所示的计算模型进行网 格. 剖分, 共剖分 15 000 个单元, 15 350 个结点, 根 据实验情况, 开采 1 号, 2 号, 3 号工作面, 工作面的 长度分别为 100m, 80m, 60m, 工作面的走向长度都 为200 m, 每个工作面的开采步距都为 20 m, 即每开 采 20 m, 计算一次围岩应力, 整个区域开采完毕, 需 计算 30次, 而在此期间, 流场可按时间步长进行多 次循环, 如每开采 4 m, 那么 20 m 需 5 d 的时间, 按 10 小时的等步长计算, 需计算 12 次, 即围岩应力场 计算一次, 渗流场计算 12次作为一个大循环。 2. 3 底板围岩应力分布规律 篇幅所限, 只从一个角度分析顶底板的应力变 化规律, 即在 1 号工作面的中部取一纵剖面 AA 图 2 分析顶底板应力随工作面开采的变化规律。 图 2 AA 剖面示意图 图 3 6为 1 号工作面开采各阶段纵向中剖面 的应力分布云图。而图 7画出了 1 号工作面中部纵 剖面AA 随 2 号工作面开采的应力分布图。 图 3 1 号面开采20 m, 纵向中剖面 z应力分布 kg/ cm2 图 4 1 号面开采60 m, 纵向中剖面 z应力分布 kg/ cm2 402太 原 理 工 大 学 学 报 第 35卷 图 5 1 号面开采 140 m, 纵向中剖面 z应力分布 kg/ cm2 图 6 1 号面开采 180 m, 纵向中剖面 z应力分布 kg/ cm2 从图 3 7 中可得出如下几点结论 1 顶底板的应力分布规律基本上是一致的, 但 当工作面推进到 80 m 左右时, 顶底板的应力出现一 些细小的差别, 即 80 m 之前, 顶底板的拉应力区是 一致的, 底板拉应力的范围是从工作面后方 5 m 左 右开始一直到采空区, 最大拉应力区紧接煤层底面, 深度在 8 m 左右; 80 m 之后, 底板最大拉应力区下 移, 即在距煤层底面 10 m 以下。 2 从总体来看, 底板最大拉应力区出现在距 图 7 2 号面开采 120 m, 纵向中剖面 z应力分布 kg/ cm2 工作面后方 15 m 左右, 且随着工作的推进其应力不 可恢复, 即一直处于受拉状态, 这与前人计算的结果 有所区别[6]。 3 底板拉应力的深度在 80 m 之前是随着工作 面的推进而逐渐加深, 80 m 之后基本稳定在 20 m 左右。 4 按最大拉应力准则判断, 计算方案底板的破 坏深度在 10 15 m 之间。 5 在整个开采过程中, 顶底板的应力一直处于 动态变化, 也就是说, 本工作面的开采会影响其临近 区域 包括临近采空区 应力的分布, 所以突水不一 定发生在开采的工作面, 在其临近区域也可能诱发 突水事故。 3 结 论 l 承压水上采煤问题是一个应力场和渗流场 相互作用问题, 应该用固流耦合的方法加以解决。 2 通过分析底板应力分布、 破坏深度与规律, 揭示了带压开采中底板应力变化的运移规律。 3 通过计算, 提出了一些承压水上采煤的有益 建议, 并成功地应用于太原东山煤矿。 参考文献 [ 1] Cook N G W. natural joints in rock[ J] . Int J rock Mech Min Sci 2. The Third Geology Reconnaissance Department of Shanxi Province, Taiyuan 030600, China AbstractFrom the point of view of GIS, this paper analyzes the mineral resources programming diagram and basic elctronical data request of mineral resources, probes into the mineral resources pro gramming work based on the geography ination system GIS , discusses the technique that combining intellignet reasoning and figure digital ination, realizes the conjunction of digital infor mation and graphic ination, is good to intelligent communicating of the mineral resources program ming. Key wordsmineral resources; programming diagram; protracting 编辑 张爱绒 上接第 403 页 Study on 3D Numerical Simulation for the Failure Laws of Floor in Coal Mining with Water Pressure LIU Zhijun , HU Yao qing Mining Technology Institute, T aiyuan University of Technology , Taiyuan 030024, China Abstract In this paper, 3D solid fluid coupling mathematical model and numerical solution is built from the point of view of mining with water pressure, and numerical simulationson 3D solidfluid cou pling for mining with water pressure in Dongshan coal mine of T aiyuan is carried out . T he paper also discussed stress distributing laws of floor in the process of mining with water pressure, and provided a theoritical base for the prevention of water outburst from coal mine. Key wordsmining with water pressure; solidfluid coupling; finite element; numerical simulation 编辑 张爱绒 406太 原 理 工 大 学 学 报 第 35卷