地下开采引起地表沉陷的数值模拟.pdf

总第 1 2 8期 2 0 0 6 年第 1 2 期 西部探矿工程 WES I 、 一CHI NA EXP L0RAT1 0N ENGI NEE RI NG s e r i e s No ．1 2 8 I e e ． 2 O 0 6 文章编号 1 0 0 4 --5 7 1 6 2 0 0 6 1 2 0 2 9 9 0 3 中图分类号 T U 4 3 3 文献标识码 B 地 下开 采引起地表沉 陷的数值模 拟 廉海， 魏秀泉， 甘德清 河北理工大学资源与环境学院， 河北 唐山0 6 3 0 0 9 摘要 以石人沟铁矿为研究对象， 运用 F I A C软件， 采用数值模拟的方法， 对石人沟铁矿地表沉陷数学模型的开采过 程进行分步开挖模拟， 并总结出了地下开采引起的上覆岩层和地表的移动规律 ， 给 出了不同开采阶段的地表动态移动 曲线 。 关键词 地下开采； 岩层移动； 地表沉陷； 数值模拟 1 国内外开采沉陷理论的研究发展 地下开采引起的岩层移动及地表下沉是一个复杂的物理、 力学变化过程， 它涉及到采矿工程的各个领域。在国外 ， 早在 1 5 世纪。 地表开采沉陷问题就已被人们所注意。但这个时期开采沉 陷理论发展较慢， 直到 1 9世纪， 相继出现了德国人依琴斯凯的 “ 二等分线理论” ， 耳西哈_ 3 0 的“ 自然斜面理论” ， 法国人法约尔 。 J F a y o 1 的开采沉陷“ 拱形理论” 等， 人们对开采影响问题才有了 一 些初步的认识和研究。 1 9世纪末， 沉陷理论研究逐步深入， 在研究各种沉陷分布与 各类角度的同时， 廿 土 开始了对开采引起的地表下沉与变形规律 的研究。1 9 5 3年波兰学者萨武斯托维奇 提出了将岩体视为均 质线弹性体， 提出计算岩体下沉的方法。随后李特维尼申等学者 应用非连续介质力学中的颗粒体介质力学研究岩层及地表移动 问题， 用概率论的方法建立 由地下单元开采所引起的岩层及地 表单元下沉盆地表达式、 单元水平移动表达式， 经迭加建立 ． F 沉 的剖面方程及其移动与变形分布表达式。随后通过实践应用到 2 o世纪 7 0年代末期， 形成了地表移动破坏的理论体系。 从 2 0世纪 7 0年代至今， 随着科学技术的发展和研究手段的 提高， 开采沉陷理论中渗透了其它相关学科的理论， 得到了进一 步的发展。在这一阶段， 陆续问世的有限元、 离散元、 边界元及各 种算法和程序充分发挥了各 自的长处， 在很大程度上促进 了岩 石软科学研究以及岩体地下工程的发展。 在我国， 地下开采引起的岩层移动及地表沉陷理论是建国 以后发展起来的。一些学者从不同的侧面用不同的方法对开采 沉陷理论进行了探讨研究 。1 9 8 1 年刘天泉与仲惟林等学者合 作， 提出了覆岩破坏的基本规律和冒落带高度的计算方法。李增 琪 1 9 8 5 采用 F o u r i e r 变换推出了岩层与地表移动表达式。杨 硕 1 9 9 0 提 出了开采沉陷 的力学预测模式。李永树 1 9 9 5 ～ 1 9 9 7 等建立了褶曲构造煤层三维空间开采时地表单元下沉盆地 和单元水平移动表达式。崔希民 1 9 9 6 对主断面的地表移动与 变形进行了实时位移的分析， 应用流变模型进行开采沉陷研究。 另外， 一些现场技术人员对采动滑坡、 采动对断层的影响、 巨厚松 散层下采矿对地表移动规律、 条带开采地表移动规律等方面也 进行了大量的研究。 总之， 经过几十年的研究， 开采沉陷理论已取得了很大发展 和完善。 2 石人沟铁矿简介 石人沟铁矿位于河北省遵化市境内。矿区所处 区域地质构 造位置为燕山沉降带马兰峪复背斜轴部。矿区为一单斜构造。 工程地质条件属简单型。 石人沟铁矿矿体走向长 2 ． 6 k in， 自南向北分布在 3 o号勘探 线至 5号勘探线之间。矿区浅部矿体采用露天开采， 露天开采境 界最低标高为0 m， 一般在 1 6 2 5 m之间。露天开采境界南北长 2 ． 8 k in。由南向北分为三个采区， 2 8 ～1 8线为南区， 1 8 ～8线为 中区， 8线以北为北区。 石人沟铁矿矿体赋存于太古界迁西群马兰峪组片麻岩中。 磁铁石英岩矿体与角闪斜长片麻岩岩层相互平行产出， 有 M 0 ， M1 ～M5计 6个矿体， 属多层矿。据统计， 各矿体之间夹层平均 厚度为 2 1 ． 3 m。矿体走向近南北， 向西倾斜， 倾角 5 0 。 ～7 O 。 ， 平均 倾角 6 o 。 ， 属急倾斜矿体。矿体顶底板围岩为角闪斜长片麻岩、 黑云角闪斜长 片麻 岩， 岩石 普 氏硬 度系 数为 6～ 1 0 ， 体重 2 ． 8, ／ m。 ， 松散系数 1 ． 5 ； 矿石为磁铁石英岩， 普氏硬度系数 1 0 ～ 1 2 ， 体重 3 ． 4 t ／ m。 ， 松散系数 1 ． 5 ， 矿体与围岩均比较稳固。矿区 示意 图如 图 1 所 示。 3 地表沉陷数学模型 的建立 本次研究分别以南区、 北区为例 ， 应用 F I AC软件进行数值 模拟， 模拟石人沟铁矿露天转地下后， 地下开采引起的岩层破坏 和地表移动的过程， 分析不同开采中段、 不同岩层的移动角， 并最 终得出地表下沉曲线模型。 3 ． 1 基本假设 在数值模拟过程中， 为了使计算结果比较接近实际情况， 对 岩体介质性质、 计算模型、 矿山地质条件、 受力条件、 采矿工艺及 采矿方法等都作了必要的假设。 1 对矿岩性质的假设。假设矿岩为各向同性、 均质且符合 摩尔库仑弹塑性模型的介质。 2 对计算模型的假设。对地下工程开挖来说， 地下矿山开 采是一个空间问题 ， 应采用三维空间计算模型更为合理。但一般 来说， 在同等条件下， 二维数值模拟结果与三维数值模拟的计算 维普资讯 3 0 0 西部探矿工程 De c ． 2 0 0 6 NO ． 1 2 结果比较接近。因此， 计算模型简化为二维平面模型。 3 矿房结构的简化。为模拟方便 ， 对巷道工程、 每一矿房的 开挖步数等不予考虑， 模拟时简化为实体。矿房结构参数 矿块 长度 6 0 m； 境界顶柱 1 6 m； 中段高度 4 4 m； 间柱宽度 1 0 m。 图 1 石人沟铁矿矿区现状示意图 4 计算不考虑时间有关的物理量。 3 ． 2 模型 建立 进行模拟计算的模型采用 F L AC 2 ． 2 5建立 ， 坐标原点设在 左下角 ， X轴水平， 由左向右方向为正方向； Y轴竖直， 由下向上 为正。在模型的左侧， 固定 X、 Y方向位移， 在其底部只固定 Y 方向位移， 在其右侧只固定 X方向位移。 3 ． 3 网格划分 由于此次模拟是二维模拟， 并考虑了计算机的内存 和软件 的版本问题， 所以网格数不能超过 1 5 0 0个。网格密度的划分应 视研究的内容和重点进行。开采范围内的矿体 以及围岩处网格 划分的比较密集， 其它部分的网格划分的比较大。 3 ． 4 矿岩物理力学参数选取 在本次模拟过程中， 力学模型选为线弹性模型， 依据分析需 要将整个模型划分为三种材料， 即 M1 矿体、 M2矿体、 黑云母角 闪斜长片麻岩， 其参数数值如表 1 所示。 表 1 矿岩物理 力学 参数 块体 抗拉抗压内 聚内摩 弹 性 泊 密 度 强 度 强 度 力 擦 角 模 量 松 依 据 称 1 0 4 M Pg ／ c m3 X MP a MP a MP a 0 M P a 比 ⋯ ’ tp 。 比 3 ． 5 边界条件的确定 边界条件按剖面所处位置的应力条件以及按照开采沉陷的 原理确定。边界条件为 水平方向上左右两侧取水平位移约束边 界 即 u 一O ， 下边界取垂直位移约束边界 垂直位移 u 0 ， 上 边界 地表 为自由边界。 由于研究区域范围不大， 可以认为研究区域处于均匀分布 的应力场中， 并且认为岩石中不存在构造应力， 则岩石中的应力 主要由上覆岩石的重量引起， 垂直应力和水平应力可以分别 由 下式求得 d 一 H 一 H 式中 上覆岩石的重力密度； Hl l一单元立方体所在的深度； 岩石的泊松比。 在模拟过程中， 模型的边界约束条件如图 2 、 图 3 所示。 ／ r ／ 二 ＼ ／ ＼ ／ 二 ＼ ， ’ ＼ 图 2 原岩应力模拟的边界约束条件 4 数值模拟结果分析 4 ． 1 矿房、 矿柱动态模拟过程及结果分析 从矿房分步开挖的模拟结果看 ， 最大、 最小主应力随开采的 进行， 在空区两侧应力值逐渐增加， 均高于原岩应力值。开挖区 域越大． 向空区两侧重分布的应力越大。北区开挖 8个矿房后最 大拉应力和最大压应 力值都没有达到相应的极限强度， 空区相 对来说比较稳定。但结合实际开采情况， 北区开挖 8 个矿房后就 回收矿柱， 处理了采空区。南 区在开挖 5 个矿房后， 拉应力和压 应力都达到了极限强度， 矿房顶、 底板受到了不lq程度的破坏。 所以当采完第 5 个矿房后就处理了空区， 即崩落岩层。在后面的 矿房开采中， 每开挖 4 ～5个矿房就处理一次采空区。在南区整 个模拟开采过程中， 进行了 3次矿柱回收。回收矿柱后， 应力再 次重新分布， 产生应力降低区。随着地下矿房的开挖 ， 应力又逐 渐增大。直到达到抗压强度和抗拉强度， 再次破坏顶板岩层。 4 ． 2 地表移动动态曲线的建立 由于此次模拟采用 F I AC软件模拟， 运用 F I AC程序可以 维普资讯 2 O 0 6定 第 1 2 期 廉海， 魏秀泉， 甘德清 地下开采引起地表沉陷的数值模拟 3 O1 对一些关键点的历史位移和最大主应力、 最小主应力进行跟踪 记录。在此次研究过程中， 每个模型都选取了采空区上方地表的 许多关键点进行记录。根据历史位移记录曲线得出竖直方向的 位移值， 画出每个模型在开采过程中地表沉陷的移动曲线。 4 ． 2 ． 1 北区地表沉陷动态移动过程曲线 根据前面的数值模拟结果和对地表关键点的历史位移的记 录， 得出了地表沉陷过程 图。由于北 区在 开挖 1个和 2个矿房 时， 移动量非常微小， 无法在图形上表示。所 以北区只画出了开 挖 8 个矿房且开挖矿柱后的移动曲线图。其动态移动曲线图见 图 4 。 4 ． 2 ． 2 南区地表沉陷动态移动过程曲线 根据数值模拟结果和历史位移记录可知， 南区由于上部为 废石松散层， 垂直位移量较大。在开采过程中， 当开挖 5个矿房 后， 就要处理矿柱和空区； 当开采到第 9 个矿房时， 第二次回收矿 柱和处理空区； 当开采到 1 1 个矿房时， 又进行了矿柱回收。因 此， 在此次模拟过程中， 共回收了三次矿柱。第一次回收矿柱后 地表最大下沉值为 1 8 m， 第二次回收矿柱后地表最大下沉值为 2 6 m， 第三次回收矿柱后地表最大下沉值为 2 7 m。 F LAC 2 ．25 S 七 日 P B pp 】 ied l 、 o ce s Hdx U c L u _ -一 6 ． 6 09E 06 O 、6 0 n - 1 2 0 n 地 表 图 3 南区边界约束条件 图 4 北区地表 沉陷移动过程 5 结论 应力状态分析 从矿房分步开挖的模拟结果看， 最大、 最小主 应力随开采的进行 ， 在空区两侧应力值逐渐增加， 均高于原岩应 力值。开挖区域越大， 向空区两侧重分布的应力越大。 垂直位移分析 在地下开采过程中， 顶 、 底板产生竖直向下的 位移， 这是由于顶、 底板处的拉应力引起的。随着矿房不断开采， 最大下沉值越来越大。最大下沉值发生在空区中部的上方， 矿房 的上部和矿房与矿柱接触的角点附近， 且偏向于开采前进 的一 侧． 呈不对称的拱形分布， 向两侧逐渐减小。 通过对模型结果的分析 ， 结合地表沉陷随开挖过程的动 态移动曲线 ， 得出了岩层 及地表移动破坏的规律 ， 并对确定 移动角和圈定移动范围提 出了特殊的方法 。在开采前 ， 先根 据工程类 比法和经验公式 ， 参照类似矿山的移动角确定该矿 山的移动角， 并根据此移动角从开采矿体的最深部划出该矿 的大致移动范围。在模拟过程中， 由于是分步开挖 ， 随着矿 房开采数 目、 开采时 间和开采空间的不 同， 随时分期 、 分布 划出地表移动曲线 ， 确定出不同岩层 的移动角 ， 并对参照的 移动角进行改进。 在地表沉陷理论中， 确定移动角和圈定移动范围是一个非 常重要的内容。本文没有研究出具体的移动角修正值， 还有待作 进一步的研究。 参考文献 E 1 ] 吴侃． 开采引起 的地表裂缝 深度 和宽度 预 汁[ J ] ．内卜新矿业 学院学 报 ， 1 9 9 7 ， 1 2 5 4 9 5 5 2 ． E 2 3 王树兀， 赵德勤． 矿{ h 开采破坏学[ M] ． 冶金_ i 业出版社， 1 9 9 3 ， 1 2 ． E 3 ] 贺跃光． 工程开挖 Ij l 起的地表移动与变形模型及监测技术研究 [ I ] ． 中南大学博 J __- 学位论文， 2 0 0 3 ， 4 l ～2 ． [ 4 ] 洪靖文． 采动覆岩动态移动破坏规律及开采沉陷预计系统研究 [ 1 j ． 中 矿业大学博士论文， 1 9 9 9 ． 维普资讯