FLAC和数值分析在矿山地表沉降预测中的应用.pdf

第5 7 卷第3 期 2005 年8 月 有色金属 N o n f e r r o t l M e t a l s V 0 1 ．5 7 ．N o ．3 A u g u s t 2 0 0 5 F L A C 和数值分析在矿山地表沉降预测中的应用 黄志安1 ，李示波1 ，赵永祥2 ，张英华1 ，倪文1 ，杨晓光1 1 ．北京科技大学，北京 10 0 0 8 3 ； 2 ．山东省乳山市人事局，山东乳山2 6 4 5 0 0 摘要使用F L A C 软件预测出采空区形成后引起的地表垂直和水平位移，然后使用数值分析方法，计算出地表各个沉降点 的斜率和曲率，预测结果对矿山建设有参考作用。 关键词采矿工程；地表沉降；预测；F L A C ；数值分析 中图分类号T D 3 2 5 ．2文献标识码A文章编号1 0 0 1 0 2 1 1 2 0 0 5 0 3 0 0 9 5 0 4 矿山地表沉降是一个不容忽视的问题，如果解 决得不好，轻则引起地表设施或建筑物的破坏，重则 有可能引起地表突然坍塌，带来灾难性的后果。地 下矿开采后形成采空区，在周围形成了一个移动带， 移动带内可能存在各种设施或建筑物，如果在已知 地下采空区大小的情况下，可以预测出地表沉降参 数，则可对地表设施的布置起到事先参考的作用，提 前采取相应措施，以确保采空区实际形成后，不会造 成大的安全事故。基于这一考虑，通过F L A C 软件 预测出采空区形成后引起的地表垂直和水平位移， 然后使用数值分析方法，计算出地表各个沉降点的 斜率和曲率。 1使用F L A C 和数值分析进行地表 预测的过程 F L A C F a s tL a g r a n g i a nA n a l y s i so fC o n t i n u a 程 序是一个有限差分程序，主要是为岩土工程应用而 开发的岩石力学计算程序，程序中包括了反映岩土 材料力学效应的特殊计算功能，可用于边坡、坝体、 隧道、硐室等岩土介质的应力和变形模拟与分析。 F L A C 程序建立在拉格朗日算法基础上，特别适合 模拟大变形和扭曲。F L A C 采用显式算法来获得模 型全部运动方程 包括内变量 的时间步长解，从而 可以追踪材料的渐进破坏和跨落。该软件是目前国 际岩土工程界十分推崇的计算机软件。使用F L A C 和数值分析进行地表沉降预测有如下步骤。 1 搜集矿山的地质综合柱状图，获取各岩层和 矿层的岩性及相关参数，并获取将开采的矿床整体 收稿日期2 0 0 4 0 2 2 6 作者简介黄志安 1 9 7 3 一 ，男，四川眉山市人，博士生。主要从事 采矿和安全技术工程方面的研究。 空间赋存情况。 2 在F L A C 程序中选择摩尔一库仑模型，将相 关数据输入到该模型中，通过该模型的计算处理后， 可提取所取截面各相关点的纵向位移和横向位移。 3 使用第 2 步所提取的地表各离散点的纵向 位移和横向位移值，通过数值分析方法中的数值微 分来求得各相应点的斜率和曲率，得出结论。 2 地表预测实例 某矿区属于地质条件简单型的矿区。所采煤层 产状近水平，倾角4 。，煤层距地表1 5 0 ～1 7 0 m ，上部 岩层夹若干薄煤层，厚度0 ．6 ～1 ．O m 。顶板为细及 中砂岩，灰白色，以长石、石英为主，含云母等暗色矿 物，并含植物化石，底板为粉及细砂岩。煤层顶、底 板岩层较为稳定。工作面处于矿井疏干区内，水文 条件简单，涌水量小。工作面走向长1 6 3 1 m ，沿倾向 长1 3 0 m ，开采的煤层厚度1 0 m 。所以该采空区上部 岩层移动为典型的平面应变问题。 1 由于柱状图的篇幅过大，现对其进行简化， 得出该矿区的煤层横断面图，如图1 所示。 2 选择F L A C 程序中的计算模型，然后划分网 格进行计算。 由于煤层长度方向远大于宽度方向，属于平面 应变问题，选取其中的一个平面进行计算。矿区涉 及的岩体有砂岩和煤层，风化性较强，属弹塑性材 料，因此选择摩尔．库仑 M o h r C o u l o m b 模型。摩尔 一库仑屈服准则为工 d 1 一d 3X 1 s i n 9 ／ 1 一 s i n | 9 0 2 c [ 1 s i n 9 ／ 1 一s i n 9 ] 1 ／2 ，式中仃1 ， d ，一最大和最小主应力；C 一岩土体的粘结力；9 一 岩土体的内摩擦角。 当f 0 时，岩土体将发生剪切破坏。岩土体在 万方数据 有色金属第5 7 卷 国一灰白色砂岩；圈一烁岩及砂烁岩；目一砂岩；●一煤 岩层厚度／m a 一5 0 ；b 一1 0 ；C 一1 0 ；d 一6 ；e 一1 8 ．6 ； f 一0 ．4 ；g 一2 6 ．5 ；h 一0 ．3 ；i 一3 0 ．7 ；j 一0 ．6 ； k 一4 0 ．8 ；l 一0 ．3 ；m 一1 5 ；n 一0 ．8 ；o 一1 0 。 图1 煤层横断面 F i g ．1 C r o s ss e c t i o no fc o a ls e a m 达到屈服极限后，在恒定的应力下将产生塑性变形， 在拉应力状态下，如果拉应力超过岩土体的抗拉强 度，岩土体将会产生拉破坏。各岩层的力学参数如 下。砂岩弹性模量E 1 9 ．3 1 0 3 M P a ；泊松比岸 0 ．3 8 ；粘结力c 0 ．7 M P a ；内摩擦角9 2 8 。；平均 容重’， 2 ．6 1 0 3 k g ／m s 。煤弹性模量E 4 ．1x 1 0 3 M P a ；泊松比肛 0 ．3 3 ；粘结力c 0 ．0 2 M P a ；内 摩擦角9 1 5 。；平均容重y 1 ．4 1 0 3 k g ／m 3 。砾 岩弹性模量E 2 0 t 0 3 M P a ；泊松比口 0 ．2 5 ； 粘结力C 0 ．5 M P a ；内摩擦角9 2 6 。；平均容重y 1 ．8 1 0 3k g ／m 3 。 划分F L A C 网格，并在地表设置历史记录点，记 录在运算后地表的水平和竖直位移值。网格的节点 间距以及这些记录点的间距都跟煤层的平均开采深 度有关。开采深度越大，选取的点的间距也越大。 相反，开采深度越小，开采煤层距离地表越近，选取 的点的间距就越小，点取得越密。在平均开采深度 1 0 0 ～2 0 0 m 时，选取的点间距应为1 5 m ，所以在地表 每1 5 m 设置一个点，记录其位移值。 将岩石和煤层相关参数输入到模型中，运行程 序，得出位移矢量图、Y 方向位移等值线图和z 方向 位移等值线图，如图 2 ～图 4 所示。 图2 位移矢量图 F i g ．2D i s p l a c e m e n tv e c t o rf i g u r e 图3Y 方向位移等值线图 F i g ．3D i s p l a c e m e n ti s o l i n ef i g u r eo fYo r i e n t a t i o n 图4x 方向位移等值线图 F i g ．4D i s p l a c e m e n ti s o l i n ef i g u r eo fzo r i e n t a t i o n 对网格图上的地表相应点从左到右依次进行编 号，以最左边的第一个点为基点，并将其取作坐标原 点，右上方向为两坐标轴的正方向。通过程序计算 后生成相应的z 、y 方向位移值，例如3 1 号点的z 和 Y 方向的位移值如图5 所示。最后得出的地表所有 点的结果如表1 所示。 表1 地表的水平和竖直位移值 T a b l e1H o r i z o n t a la n dv e r t i c a ld i s p l a c e m e n tv a l u eo fS u l { a c e 点号z ，点号 z ， 点号 z ， 1O ．O2 61 72 2 4 ．24 5 03 31 6 ．6一1 8 21 ．62 71 81 8 2 ．06 1 63 41 3 ．31 6 33 ．22 71 98 9 ．87 3 93 5一1 0 ．61 3 44 ．62 42 07 ．08 0 63 68 ．8一1 1 56 ．92 42 15 9 ．07 9 83 77 ．39 68 ．72 22 21 2 6 ．27 4 23 86 ．09 71 0 ．01 92 32 0 1 ．06 4 63 95 ．08 82 1 ．11 32 42 3 6 ．75 1 84 04 ．08 93 0 ．62 2 52 3 6 ．0 3 7 94 13 ．18 1 01 1 2 ．51 4 2 62 3 0 ．4 2 7 24 22 ．77 1 11 1 0 ．O一3 62 72 0 1 ．62 0 l4 31 ．77 1 21 1 6 ．45 92 8 1 7 5 ．2 1 5 84 4一1 ．18 1 31 1 9 ．99 0 2 91 0 2 ．5 1 0 44 50 ．78 1 4 1 8 4 ．31 4 63 05 8 ．46 1 4 60 ．47 1 52 1 2 ．72 0 4 3 13 7 ．5 3 54 70 ．07 1 62 1 8 ．12 9 7 3 2 2 1 ．82 4 万方数据 第3 期黄志安等F L A C 和数值分析在矿山地表沉降预测中的应用9 7 8 k j 雨 。 色 ≤ h 赫 魁 - 3 ．0 0 【 0 ’4 ．0 0 0 0 - 5 ． o 【 0 - 6 ．o 0 0 ～7 ．0 0 0 0 81 01 21 41 61 82 02 22 42 6 时间步数／1 0 2 81 01 21 41 61 82 02 22 42 6 时问步数／1 0 2 图53 1 点的水平和竖直位移图 F i g ．5 H o r i z o n t a la n dv e r t i c a ld i s p l a c e m e n t f i g u r eo fN o ．3 1p o i n t 3 由于有限差分是按照一定网度划分出网格， 即F L A C 得出的数值解是一些离散的点的数值，而 地表移动引起的对地面构筑物的破坏主要表现在地 表倾斜变形的影响、地表弯曲变形的影响、地表水平 拉伸或压缩变形的影响。这就需要求出地表相应点 位移变化的一阶导数或二阶导数。对离散点求导 数，使用数值微分方法比较方便。 若函数f z 由表格给出，求f z 在结点上的 微商，通常称为数值微分。最简单的数值微分公式是 用差商代替微商。例如由向前差商或向后差商，分别 得到厂 z o [ f x o h 一f z o ] ／ O h ， 厂 z o [ f z o 一f z o ～h ] ／ O h 。将以上 两个公式取算术平均值，即得到中心差分公式 中点 公式 ，见式 1 。 厂 z o [ “z o h 一f x o h ] ／2 h O h 2 1 由于高阶导数可用高阶差商近似。用二阶中心差 商近似厂 z ，便得到／， z 的数值微分公式 2 。 厂 z o ≈艿2 ／h 2 厂 z [ f z o h 一2 f x o f x o h ] ／h 2 2 显然，从理论上看，h 越小，微分近似越精确。但 从计算误差上看，h 越小，f z o 与f x o h 越接 近，根据误差理论，两个相近的数相减，将会大量丢 失有效数字。为了克服上述弊病，一个很自然的想法 是求式 1 或式 2 误差渐进展开式，再利用外推法 提高精度。根据T a y l o r 公式f x 八z 髓厂 z o h 2 ／2 厂 z 3 ／31 ，3 ’ z ⋯和f x h f x 一厂 z o h 2 ／2 ／， z 一 3 ／3 1 ，3 z ⋯，两式相减，除以2 ，移项得厂 z [ 厂 z 一f x h ] ／2 h 一乏，2 ‘ ” z h 2 ‘／ 2 i 十1 i 1 ～O O 。类似的，还可以得到 ／， z o [ 八z o h 一2 f z o f x o h ] ／h 2 ∑ 2 1 2 z h 2 i ／ 2 i 2 i 1 ～o O 。 令G h 为上面两式等号左端的第二项，在 R i c h a r d s o n 外推法中，取q 1 ／2 ，则对于一阶和二 阶导数可建立外推公式G l G h ，G 。 1 h [ G 。 h ／2 一4 - r a g 。 h ] ／ 1 4 一 m 1 ，2 ， ⋯ ，其中厂 z 一G 。 l h O h 2 m 1 ’ ， ／， z o 一G 。 l h O h 2 卅 1 ’ 。 计算时，步长分别取3 0 m 和1 5 m ，即可按照所得 的G 。 1 h 值计算出相应的一阶导数和二阶导数， 其中m 1 ，即外推一次即可。忽略高阶无穷小，相 应的G 。 1 h 即可取为一阶导数或二阶导数，计算 后所得的一阶导数和二阶导数值见表2 ，为节省篇 幅，表2 只选取了1 5 个点。 表2 相应点的一阶导数和二阶导数值 T a b l e 2F i r s td e r i v a t i v ea n ds e c o n dd e r i v a t i v eo fc o r r e s p o n d i n gp o i n t s 瑚瑚瑚瑚啪瑚瑚似|蒡喊弧呲言耄咖。加■o吨c；j矗 万方数据 9 8有色金属第5 7 卷 表2 中的对应值就是所需要的结果。比如移动 带内的建筑物和移动带内埋藏的管道，就需要考虑 建筑物和管道所在点的水平方向的斜率，看它是否 满足拉伸或压缩要求，还需要考虑竖直方向的曲率 是否满足建筑物和管道的弯曲要求。还比如移动带 内的高压线铁塔，则需要考虑所在点的竖直方向上 的斜率是否满足架线倾斜要求，移动带内的铁路则 需要考虑所在点的竖直方向上的斜率是否满足铁路 参考文献 的倾斜要求等，依此类推。 3结语 矿山实例经过简化处理，在矿山实际预测中，各 项工作都需要细化，才能得出理想的结果。使用三 维模型原理上与二维模型相同，即首先求出相关点 的水平位移和竖直位移，再求出相应点的一阶导数 和二阶导数。 [ 1 ] 刘钦圣，张晓丹，王兵团．数值计算方法教程[ M ] ．北京冶金工业出版社，1 9 9 8 2 2 9 2 3 0 ． [ 2 ] 李先炜．岩体力学性质[ M ] ．北京煤炭工业出版社，1 9 9 0 2 4 0 2 8 0 ． [ 3 ] 何国清，杨伦，凌赓娣，等．矿山开采沉陷学[ M ] ．北京中国矿业大学出版社，1 9 9 1 5 1 8 ，4 7 5 2 ． A p p l i c a t i o no fF L A Ca n dN u m e r i c a lV a l u eA n a l y s i st oM i n eS u r f a c eS e d i m e n t a t i o nP r e d i c t i o n H U A N GZ h i a n l ，L IS h i ．b 0 1 ，Z H A OY o n g - x i a n 9 2 ，Z H A N GY i n g - h u a l ，N IW e n l ，Y A N GX i a 伊g u a n 9 1 1 ．U n i v e r s i t yo fS c i e n c ea n dT e c h n o l o g yB e i j i n g ，B e i j i n g1 0 0 0 8 3 ，C h i n a ； 2 ．P e r s o n n e lB u r e a u ，R u s h a n2 6 4 5 0 0 ，S h a n d o n g ，C h i n a A b s t r a c t T h es u r f a c ev e r t i c a la n dh o r i z o n t a ld i s p l a c e m e n t sr e s u l t e df r o mt h em i n e d ．o u tp i ta r ep r e d i c t e db yu s eo ft h e F L A Cs o f t w a r ea n dt h es l o p ea n dc u r v a t u r eo fs u r f a c es e d i m e n t a t i o na te v e r yl o c u sa r ef i g u r e do u tw i t hn u m e r i c a l v a l u ea n a l y s i st h e o r y ．T h ep r e d i c t i o nr e s u l t sc a nb er e f e r r e df o rt h em i n ec o n s t r u c t i o n ． K e y w o r d s m i n i n ge n g i n e e r i n g ；g r o u n ds u b s i d i n gm i n e ；p r e d i c t i o n ；F L A C ；n u m e r i c a lv a l u ea n a l y s i s 万方数据