不同采空區開挖滑坡災害與結構體動力響應研究

畢鈺璋 孫新坡 何思明 單 雨 王安輝

（1.東南大學巖土工程研究所，江蘇南京210096；2.四川理工學院土木工程學院，四川自貢643000；3.中國科學院山地災害與地表過程重點實驗室，四川成都610041；4.中國科學院水利部山地災害與環境研究所，四川成都610041；5.中國科學院青藏高原地球科學卓越創新中心，北京100101；6.西南交通大學土木工程學院，四川成都610041）

滑坡災害的發生通常伴隨著大量勢能轉化為動 能的過程，并且以其強沖擊力對周圍的結構體造成極大的破壞［1-4］?；聻暮Φ陌l生在采礦工程中更為常見，人為的擾動會打破邊坡的力學平衡，并且造成局部的應力集中，進而造成滑坡災害的產生［5-8］。尤其是露天開采和地下開采的轉換過程中，采空區會觸發滑坡或者造成古滑坡的復活，對邊坡的穩定性造成不可估量的影響，并對周圍的交通路線、村莊、各種民用公用結構造成難以估量的威脅［9-10］。

關于采空區對邊坡穩定性的影響，國內外學者都作了大量的工作。Fujisawa等人［11］建立了相關的數學模型并模擬了二維條件下的采空區對邊坡穩定性的影響；Shao Yong等［6］基于實際工況，采用三維有限元的方法分析了某礦采空區對周圍山體穩定性的影響，并給出了不同開挖條件下的位移云圖；柴保紅等［12］基于強度折減法采用FLAC3D對不同開采條件的邊坡進行了詳細的研究，分析了開挖地點在邊坡下部、邊坡中部、邊坡上部等工況下對邊坡穩定性的影響；宋偉東等［13］基于有限差分方法和極限平衡法研究了大冶鐵礦東露天邊坡的工況，分析了露采轉地采應該滿足的安全條件；丁桂伶［14］基于現場實地的調查并采用有限元方法對大村澗村采空區上方的邊坡穩定性進行了分析，不僅得到了天然狀態和飽和狀態的安全系數，而且分析了房屋結構的沉降情況。雖然大量的研究都是關于采空區和滑坡的，然而前人的研究主要側重于采空區對邊坡穩定性的影響，對災害產生后對周邊結構體的潛在影響并沒有分析，而這個恰恰是防災減災工程中重要的一個指標。

關于災害體和結構體之間的動力機理研究屬于災害動力學的范疇，相關的研究工作主要集中在研究巖崩、泥石流等災害對周邊結構體的影響。Teufelsbauer等［15］基于二維離散元方法研究了顆粒流和結構體之間的相互作用，并得出了邊坡形狀等因素對于顆粒流沖擊力的影響；Zanuttigh等［16］基于淺水方程研究了泥石流和災害體的相互作用，揭示了流深等因素對于促進泥石流沖擊力的關系；畢鈺璋等研究了災害和新型結構體的相互的作用［2］，以及約束條件下粗細混合顆粒和結構體之間的動力響應［17］，揭示了顆粒分選對于災害破壞力的影響?？梢娊沂緸暮w和結構體之間的動力響應規律對于指導工程應用有著十分重要的意義。

本研究采用二維離散元方法，通過和前人工況算例對比來確定相應的數值模擬參數；進而建立物理模型對不同開挖工況進行模擬，分析其崩塌規律；最后分析不同采空區條件下滑坡災害與結構體之間的動力響應規律。旨在得到相關的規律性結論，并對實際的采礦工程有一定指導意義。

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種以上不同尺寸的顆?？梢杂行У胤乐诡w粒有序性，而顆粒有序性會限制邊坡破壞面的產生。如圖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.