不同采空区开挖滑坡灾害与结构体动力响应研究

毕钰璋 孙新坡 何思明 单 雨 王安辉

（1.东南大学岩土工程研究所，江苏南京210096；2.四川理工学院土木工程学院，四川自贡643000；3.中国科学院山地灾害与地表过程重点实验室，四川成都610041；4.中国科学院水利部山地灾害与环境研究所，四川成都610041；5.中国科学院青藏高原地球科学卓越创新中心，北京100101；6.西南交通大学土木工程学院，四川成都610041）

滑坡灾害的发生通常伴随着大量势能转化为动 能的过程，并且以其强冲击力对周围的结构体造成极大的破坏［1-4］。滑坡灾害的发生在采矿工程中更为常见，人为的扰动会打破边坡的力学平衡，并且造成局部的应力集中，进而造成滑坡灾害的产生［5-8］。尤其是露天开采和地下开采的转换过程中，采空区会触发滑坡或者造成古滑坡的复活，对边坡的稳定性造成不可估量的影响，并对周围的交通路线、村庄、各种民用公用结构造成难以估量的威胁［9-10］。

关于采空区对边坡稳定性的影响，国内外学者都作了大量的工作。Fujisawa等人［11］建立了相关的数学模型并模拟了二维条件下的采空区对边坡稳定性的影响；Shao Yong等［6］基于实际工况，采用三维有限元的方法分析了某矿采空区对周围山体稳定性的影响，并给出了不同开挖条件下的位移云图；柴保红等［12］基于强度折减法采用FLAC3D对不同开采条件的边坡进行了详细的研究，分析了开挖地点在边坡下部、边坡中部、边坡上部等工况下对边坡稳定性的影响；宋伟东等［13］基于有限差分方法和极限平衡法研究了大冶铁矿东露天边坡的工况，分析了露采转地采应该满足的安全条件；丁桂伶［14］基于现场实地的调查并采用有限元方法对大村涧村采空区上方的边坡稳定性进行了分析，不仅得到了天然状态和饱和状态的安全系数，而且分析了房屋结构的沉降情况。虽然大量的研究都是关于采空区和滑坡的，然而前人的研究主要侧重于采空区对边坡稳定性的影响，对灾害产生后对周边结构体的潜在影响并没有分析，而这个恰恰是防灾减灾工程中重要的一个指标。

关于灾害体和结构体之间的动力机理研究属于灾害动力学的范畴，相关的研究工作主要集中在研究岩崩、泥石流等灾害对周边结构体的影响。Teufelsbauer等［15］基于二维离散元方法研究了颗粒流和结构体之间的相互作用，并得出了边坡形状等因素对于颗粒流冲击力的影响；Zanuttigh等［16］基于浅水方程研究了泥石流和灾害体的相互作用，揭示了流深等因素对于促进泥石流冲击力的关系；毕钰璋等研究了灾害和新型结构体的相互的作用［2］，以及约束条件下粗细混合颗粒和结构体之间的动力响应［17］，揭示了颗粒分选对于灾害破坏力的影响。可见揭示灾害体和结构体之间的动力响应规律对于指导工程应用有着十分重要的意义。

本研究采用二维离散元方法，通过和前人工况算例对比来确定相应的数值模拟参数；进而建立物理模型对不同开挖工况进行模拟，分析其崩塌规律；最后分析不同采空区条件下滑坡灾害与结构体之间的动力响应规律。旨在得到相关的规律性结论，并对实际的采矿工程有一定指导意义。

1 数值模型及参数选取

离散单元法（discrete element method）由Cundall和Strack在1979年提出［18-19］。离散单元法的基本组成元素包括二维圆盘单元、三维圆球单元、和墙单元。该方法通过圆盘或者圆球的粘结进而研究固体材料的力学特性。本次研究采用PFC2D（二维颗粒流方法）来搭建数值模型，进而进行一系列的数值分析。

1.1 颗粒接触本构模型

本次研究中对颗粒间的接触模型采用平行粘结模型。平行粘结模型是颗粒流中常用的一种模型，它是颗粒间相互粘结的模型，可以通过对颗粒体的胶结来模拟岩土体［19］。平行粘结模型的粘结只能发生在颗粒接触点很小的范围内，通常用于模拟粘性摩擦质地的材料。从物理的角度来看，通常粘结岩石和胶结土的颗粒间的距离是一定的，胶结体相当于“颗粒键”的作用，不仅可以传递力还可以传递弯矩［20］。在本次研究中，DEM主要作为数值工具来模拟实际工况并对宏观的滑坡灾害进行分析，而微观颗粒间的力学性能并不是本文分析的重点。基于此，本研究中“颗粒键”模型的选取尽量简化，从而减少微观参数的数量。

一般而言，平行粘结模型主要由颗粒密度、颗粒形状、颗粒尺寸、颗粒分布等微观特性共同决定。目前关于颗粒微观参数和实际材料的宏观参数之间没有确切的数学模型来描述。因此关于颗粒的微观参数的确定一般采用反演的方法，即用颗粒构成的宏观材料和实际材料相对比，使二者宏观力学特性相吻合，进而确定用于离散元数值模拟的微观力学参数。数值双轴试验和剪切试验通常被用于评价材料的微观参数［21-22］。在PFC2D中使用以下关系计算颗粒刚度以及接触键参数：

式中，kn是颗粒法向刚度；t是颗粒沿着平面的厚度；Ec是颗粒接触的杨氏模量；φn和φs分别表示接触键的法向强度和切向强度；σn和σs分别表示材料的法向强度和切向强度；R为2个颗粒的平均半径。

方程（1）展示了在二维情况下，单个接触键的强度与颗粒半径成比例，并且当其承载的任一应力等于或超过材料强度时，接触键将断裂。

1.2 数值实验的参数确定

离散元模拟中的对象的宏观力学行为取决于微观颗粒的力学参数的设定。目前关于离散元模拟的参数选取主要采用参数反演的方法。有些学者通过物理试验和数值模拟的堆积形态的拟合来确定最终的参数［23］；有些学者通过实际的工程数据和数值模拟的结果来反演参数［1-2］。本次研究的目的是采空区对边坡动力响应的影响，因此如果选取堆积形态或者力学参数的话就不具备很强的说服力。基于此，通过对实际的隧道开挖工程进行模拟，进而研究开挖后产生的边坡扰动，对比前人的计算结果来反演本次研究所需的参数。

Koizumi等［24］在研究隧道开挖对边坡失稳的过程中，采用了不同的计算模型来研究边坡的最大水平位移和“隧道——软弱岩层面”距离的关系。如图1（a）所示是Koizumi的二维有限元模型，整个边坡分成3个部分组成：沙土沉积层（1）、软弱岩层（2）、以及基岩层（3）。其中，隧道距离软弱岩层面的距离分别表示为：0.5D、1.0D、1.5D、2.0D、3.0D，其中D表示的是隧道半径。图1（b）为相同尺寸的离散元数值模拟的二维模型。

图2显示的是边坡表面最大水平位移和“隧道——软弱岩层面”距离的关系曲线图。Koizumi分别用双线性模型和应变软化模型分别模拟了不同“隧道—软弱岩层面”距离值对边坡表面最大水平位移的影响。中间曲线代表本次离散元数值模拟的结果。从图中可以很好地反映出本次数值模拟的结果和Koizumi对工况模拟的结果不仅在规律上、数值上都很接近，进而可以确定本次研究所用的数值模拟参数见表1。

?

1.3 本次研究的物理模型

本次研究采用二维圆盘粘结并模拟滑坡体。滑坡体模型的构建采用“雨滴法”建模，如图3（a）所示。相对于滑坡体400 m的长度和75 m的高度，颗粒的半径设置成0.6 m和0.8 m 2种，这样使得颗粒尺寸相对于滑坡体而言很小，有利于更好地去研究边坡的失稳和破坏。而2种以上不同尺寸的颗粒可以有效地防止颗粒有序性，而颗粒有序性会限制边坡破坏面的产生。如图3（b）所示，一共生成的初始颗粒数量是50 000个。边坡模型分条添加一种颜色，方便数值模拟直观地得到沉降、崩塌的结果对比。当“雨滴法”完成图3（b）的模型后，在对其一侧进行开挖，删除多余的颗粒，进而得到如图3（c）所示的最终的边坡模型。颗粒间被设置了初始的粘结力，进而保证边坡在没有外力条件下是稳定的。然后通过对模型中所有颗粒键施加较低的拉伸强度来降低粘结强度，进而在重力影响下颗粒键断裂，颗粒因此破碎并造成滑动，形成崩塌、滑坡。在本次研究中还设定了采空区的工况（图3（d）），在边坡的坡脚处设置了采空区，研究不同位置采空区造成滑坡的最大冲击力的结果。“采空区—自由面”距离为图中所示的S。相关的物理模型几何参数如表2所示。

2 “灾害—结构体”动力响应结果分析

?

研究“灾害—结构体”之间的相互作用对于实际的防灾减灾有着十分重要的意义。而灾害和结构体之间的动力响应又以力学的性能反应出来。因此本次研究着重研究不同采空区开挖条件下滑坡灾害对结构体的冲击力，对采空区高度、采空区—自由面之间的距离、结构体与坡脚的距离等一系列参数进行详细的分析。

图4显示的是边坡滑移面的生成过程。从图4可以看出，边坡的破坏从坡脚a点处开始，进而扩散到b点和c点，最终贯通整个坡面形成一个完整的边坡滑移面。由于滑坡内部的破坏机理十分复杂，本次研究不考虑开挖对边坡内部所造成的破坏，而是直接考虑滑移面产生后，滑坡灾害体运动过程中对结构体造成的破坏。评判指标主要为灾害体与结构体之间的冲击力。

2.1 “灾害体—结构体”堆积结果影响分析

图5显示了不同坡脚条件下的滑坡初始情况和崩塌结果，从图5可以看出图5（a2）（坡角90°）情况下的崩塌量要明显大于图5（b2）（坡角70°）情况下的崩塌量。图5（a2）条件下的平均下沉量也大于图5（b2）条件下的平均下沉量。不仅如此，图5（a2）情况下的滑动距离L同样也大于图5（b2）的情况。这是因为坡角的增大会使得潜在的滑坡量会随之增加，当颗粒键力逐渐减小时，边坡在重力的作用下会随之产生崩塌，而边坡面越陡，潜在滑动块的势能就越大，因此其滑动的距离L也大于坡面不陡的情况。

图5的（c1）显示的是坡角70°，挡墙距离坡脚5 m时的灾害堆积结果，α为最终堆积角度，Hd为灾害体在结构体接触面的堆积高度。通过对比图5（c1）和图5（c2）以及图5（d1）和图5（d2）可以得出，随着挡墙距离的增大，堆积角度α随之减小，灾害体在结构体接触面上的堆积高度Hd也随之减小。通过对比图5（c1）和图5（d1）以及图5（c2）和图5（d2）可以得出，随着边坡角度的增大，堆积角度α随之增大，灾害体在结构体接触面上的堆积高度Hd也随之增大。

2.2 采空区到边坡自由面距离对灾害冲击力的影响

如图6（a）所示，显示了在不同采空区离边坡自由面距离的情况下，灾害体对结构体的冲击力随时间步长的演化规律。从图中同样可以看出主滑体的冲击力要远远大于崩塌和滚石所造成的冲击力，此外S=20 m时的冲击力的值几乎在各个时段（除了1.7×106～2.1×106这段时步）都要大于其他各个工况的冲击力值。

图6（b）中给出了冲击力随不同采空区距边坡自由面距离的演化规律，分别给出了灾害体与结构体之间的最大冲击力以及平均冲击力的演化规律。从图中可以很明显地看出，最大冲击力的值先是随着S值的增加而增加，当S值增加到某个值之后，继续增加S的值会导致最大冲击力的略微下降。同样的规律可以反映到其平均冲击力上。可以看出设置不同的采空区距边坡自由面的距离会对灾害的动力性能造成不同的结果。而且，S值的变化会影响灾害的崩塌、滚石的冲击力的演化规律，使其呈现主滑体的冲击力的演化规律（最大冲击力）。

2.3 滑坡和结构体相互作用动力分析

图7显示的是70°坡角条件下的滑坡的最终堆积结果，其中采空区距离坡面的距离为10 m，采空区的采高为5 m。相对于没有采空区的情况（图5），图7条件下的灾害体更加碎屑化，破坏区域更大，堆积范围更广。从图中可以看出，堆积区B的面积基本和滑源区A的面积相吻合。除此之外，还存在孤立的堆积区C，与堆积区B之间存在一定的距离。

图8显示的是灾害体对结构体的冲击力（包括最大冲击力和平均冲击力）随采空区高度增加而变化的曲线图。从图8（a）可以看出，最大冲击力随着采空区高度的增加而增长，并在15 m时达到最大值，之后基本稳定在最大值附近；平均冲击力的变化规律类似于最大冲击力的规律，其峰值同样出现在采高15 m处，之后尽管增加采高，其平均冲击力值也是维持在该值附近。图8（b）的最大冲击力和平均冲击力的演化规律和图8（a）中的相似，同样是随着采高的增加而增加，并且维持在一个峰值附近不再变化。不同的是图8（b）中的最大冲击力和平均冲击力的峰值都要小于图8（a）中的相应的峰值，这是因为随着结构体离坡脚的距离变远，灾害体在滑动过程中会存在能量的损耗。

图8（c）中的变化规律和图8（a）和图8（b）中的略有不同，其最大冲击力维持在一个较大的值附近波动，与此同时，其平均冲击力的规律也是如此。这是因为结构体所处的距离决定的。如图7所示，当结构体位于堆积区A以内，可以拦截到所有的灾害体；而当结构体位于堆积区A之外，只有部分的崩塌体与结构体相撞击。该部分的崩塌体只占滑坡体灾害的极少比例，因而对结构体的冲击力也小得多。

3 结语

本次研究使用了二维离散元方法来对不同采高条件下的边坡失稳、滑动后的堆积、与结构体撞击等过程进行了分析。着重研究了不同采高条件下，滑坡灾害与结构体之间的相互作用。研究过程分为3个部分进行，首先，通过对前人的工作的总结、对比来确定本次研究将要使用的微观力学参数；其次，通过“雨滴法”在PFC2D中构建二维离散元模型；最后，基于所选的参数和搭建的物理模型对所提出的问题进行详细的研究。

通过二维离散元数值模拟可以得出，在无结构体拦截的条件下，灾害体的堆积结果和边坡的角度有关，角度越大，其滑程越远；当有结构体拦截时，灾害体的堆积结果不仅和边坡的角度有关，而且和其结构体距离坡脚的距离有关。开挖条件下的滑坡体的堆积区主要分成2个部分（图6），即主要的堆积区B和次要的堆积区C。当结构体位于主要的堆积区B内时，灾害体对结构体的最大冲击力、平均冲击力都会随着采空区高度的增加而增加，并且最终稳定在一个固定值不变；而当结构体位于主要的堆积区C内时，灾害体对结构体的冲击力则几乎维持一个范围内波动，不会出现之前的变化规律。

［1］ 毕钰璋，付跃升，何思明，等.牛眠沟地震滑坡碎屑化全过程离散元模拟［J］.中国地质灾害与防治学报，2015（3）:17-25.Bi Yuzhang，Fu Yuesheng，He Siming，et al.Simulation of the whole process of Niumiangou creek rock avalanche triggered by the earthquake using a distinct element method［J］.The Chinese Journal of Geological Hazard and Control，2015（3）:17-25.

［2］ 毕钰璋，何思明，付跃升，等.基于离散元方法的高速远程滑坡碎屑流新型防护结构［J］.山地学报，2015（5）:560-570.Bi Yuzhang，He Siming，Fu Yuesheng，et al.Simulation of the dynamic response of new type rock avalanche impact defense structure and the mechanism of energy dissipation base on DEM［J］.Journal of Mountain Science，2015（5）:560-570.

［3］Mcsaveney M J.Recent rockfalls and rock avalanches in Mount Cook national park，New Zealand［J］.Reviews in Engineering Geology，2002，15:35-70.

［4］Valentino R，Barla G，Montrasio L.Experimental analysis and micromechanical modelling of dry granular flow and impacts in laboratory flume tests［J］.Rock Mechanics and Rock Engineering，2008，41（1）:153-177.

［5］ 崔 杰，王兰生，王 卫，等.采空区边坡变形破裂演化机制研究［J］.采矿与安全工程学报，2008（4）:409-414.Cui Jie，Wang Lansheng，Wang Wei，et al.Deformation and fracturing mechanism of goaf slope［J］.Journal of Mining&Safety Engineering，2008（4）:409-414.

［6］Shao Yong，Yan Changhong，Xu Baotian.The stability analysis of mountain under intensive goaf［J］.Electronic Journal of Geotechnical Engineering，2016，20（4）:1181-1196.

［7］Shikai F A N.A discussion on the slope stability on the goaf［J］.Resources Environment and Engineering，2006，20:617-627.

［8］Dong L，Li X.An application of grey-general regression neural network for predicting landslide deformation of Dahu Mine in China［J］.Advanced Science Letters，2012，6（1）:577-581.

［9］Yang X J，Hou D G，Hao Z L，et al.Fuzzy comprehensive evaluation of landslide caused by underground mining subsidence and its monitoring［J］.International Journal of Environment and Pollution，2016，59（2-4）:284-302.

［10］Chan N W.Responding to landslide hazards in rapidly developing Malaysia:a case of economics versus environmental protection［J］.Disaster Prevention and Management:An International Journal，1998，7（1）:14-27.

［11］Fujiaswa K，Kasai M，Okuda S，et al.Construction of a numerical analysis model to evaluate interaction between a landslide and a tunnel［C］//The Japan Landslide Society 48th Annual Congress.Niigata，Japan［.s.n.］2009:71-72.

［12］柴红保，曹 平，柴国武，等.采空区对边坡稳定性的影响［J］.中南大学学报:自然科学版，2010，41（4）:1528-1534.Chai Hongbao，Cao Ping，Chai Guowu，et al.Influence of goaf on slope stability［J］.Journal of Central South University:Science and Technology ，2010，41（4）:1528-1534.

［13］宋伟东，杜建华，谢正平.大冶铁矿深凹露天开采最终边坡稳定性分析［J］.北京科技大学学报，2005，27（4）:385-389.Song Weidong，Du Jianhua，Xie Zhengpinig.Stability of the finished slope of Daye iron open-pit-mine［J］.Journal of University and Technology Beijing，2005，27（4）:385-389.

［14］丁桂伶.采空区上方边坡稳定性及滑坡危险性分析［J］.中国地质灾害与防治学报，2012，23（2）:38-43.Ding Guiling.Stability and landslide risk analysis of slope above goaf［J］.The Chinese Journal of Geological Hazard and Control，2012，23（2）:38-43.

［15］Teufelsbauer H，Wang Y，Pudasaini S P，et al.DEM simulation of impact force exerted by granular flow on rigid structures［J］.Acta Geotechnica，2011，6（3）:119.

［16］Zanuttigh B，Lamberti A.Experimental analysis of the impact of dry avalanches on structures and implication for debris flows［J］.Journal of Hydraulic research，2006，44（4）:522-534.

［17］毕钰璋，何思明，李新坡，等.约束条件下粗细混合颗粒动力机理分析［J］.岩土工程学报，2016，38（3）:529-536.Bi Yuzhang，He Siming，Li Xinpo，et al.Kinetic mechanism of mixed particles under constraint conditions［J］.Chinese Journal of Geotechnical Engineering，2016，38（3）:529-536.

［18］Cundall P A，Strack O D.A discrete numerical model for granular assemblies［J］.Geotechnique，1979，29（1）:47-65.

［19］Potyondy D O，Cundall P A.A bonded-particle model for rock［J］.International Journal of Rock Mechanics and Mining Sciences，2004，41（8）:1329-1364.

［20］Utili S，Nova R.DEM analysis of bonded granular geomaterials［J］.International Journal for Numerical and Analytical Methods in Geomechanics，2008，32（17）:1997-2031.

［21］Villard P，Chevalier B，Le Hello B，et al.Coupling between finite and discrete element methods for the modelling of earth structures reinforced by geosynthetic［J］.Computers and Geotechnics，2009，36（5）:709-717.

［22］Zhang L，Thornton C.A numerical examination of the direct shear tes［tJ］.Geotechnique，2007，57（4）:343-354.

［23］Bi Y，He S，Li X，et al.Geo-engineered buffer capacity of two-layered absorbing system under the impact of rock avalanches based on Discrete Element Method［J］.Journal of Mountain Science，2016，13（5）:917-929.

［24］Koizumi Y，Lee J，Yokota Y，et al.Numerical analysis of landslide behavior induced by tunnel excavation［C］//ISRM International Symposium-EUROCK 2010.Lausanne，Switzerland:International Society for Rock Mechanics，2010.