基于離散-連續耦合計算的露天礦邊坡失穩特征與影響研究

田 涯，涂立嘯，周 偉，王知樂，才慶祥，陸 翔

（1.中國礦業大學 礦業工程學院，江蘇 徐州 221116；2.中國礦業大學 露天礦山高新技術研究中心，江蘇 徐州 221116；3.中國礦業大學 煤炭資源與安全開采國家重點實驗室，江蘇 徐州 221116）

邊坡安全穩定是露天礦重要研究課題，滑坡災害發生將嚴重威脅人身安全與能源保障[1-2]。在實際生產過程中，企業往往會結合露天礦推進特點、空間形態，在規程要求范圍內犧牲邊坡部分穩定性以獲取經濟效益[3]、減少土地占用[4]。而降雨、凍融循環損傷、鹽堿蝕破壞等環境因素對巖體力學性質的弱化不可忽略[5-7]，邊坡失穩仍存在一定可能。邊坡災變過程的分析及滑體影響范圍研究對礦山應急管理與風險管控有著重要意義。

露天礦滑坡的研究中，楊天鴻等[8]通過構建失穩案例等數據庫結合巖體力學評價，從力學機理與案例匹配確定邊坡位移時間序列閾值等，實現對滑坡預測預警；李聰等[9-10]同樣以建立滑坡數據庫的方式，歸納滑坡不同階段變形特征、裂縫發育等信息來預測滑坡，并將滑坡演變過程分為初始變形、勻速變形、加速變形和急劇變形4 個階段；王立文[11]通過邊坡雷達監測獲取撫順西露天礦位移云圖、速度曲線等變化規律，對邊坡巖體變形階段進行劃分，通過速度閾值的設定對滑坡概率進行預測。在滑體形態研究中，王長軍等[12]以南芬露天礦為背景，采用光滑粒子流體動力學方法（SPH）的數值模擬再現了老滑坡區松散堆積體復活過程，滑動距離、形態等呈現較好一致性，但研究是建立在邊坡巖體和滑體力學強度差異較大且滑體形態已經給定的基礎上，普適性相對較弱；霍文等[13]則采用離散元方法對典型巖質臺階破壞過程、剝離臺階與內排臺階安全距離等進行研究，與現場吻合度較好。在算法方面，有限元法或有限差分法基于連續介質假定，網格節點變形量受限，能夠較好得對邊坡應力-應變關系進行表述，但無法較好反應邊坡位移、破壞的內容。離散元法對大變形問題有著較強處理能力，但計算過程要求顆粒數量繁多、顆粒間接觸冗余，計算精度與時間均受限制。因此離散元對邊坡研究往往集中在局部細節而非整體。

離散-連續耦合計算方法的出現一定程度彌補了離散元法與有限差分法的缺陷。研究表明耦合計算過程中邊坡應力、位移在離散域與連續域有著較高一致性，離散域細觀破裂機制可以有效表征邊坡宏觀塑性變形[14-15]。因此，以某巨厚煤層露天煤礦靠幫開采后端幫邊坡為研究背景，潛在失穩區域采用離散元方法分析滑坡特征、破壞過程、影響范圍等；基巖等變形量較小部分建立有限差分法網格，傳遞坡體內部應力、提高計算速度。

1 巖體宏細觀力學參數

所研究露天礦煤層巨厚緩傾斜，地層傾向南東-南，生產能力30 Mt/a。該礦在靠幫開采后端幫幫坡角由28°提高至31.84°，已接近極限平衡狀態。其中邊坡高度257.12 m，煤層厚度74.55 m，地層逆傾6.42°，單臺階高度15 m、臺階坡面角60°。

通過鉆探與現場踏勘將南幫地層主要分為第四系、灰色泥巖、泥巖、煤及灰色粉砂巖5 種，對煤巖試樣進行單軸抗壓強度試驗與直剪試驗，基于RMR評價法得到該礦巖體宏觀力學參數[16]。離散元中，巖體所呈現出的宏觀力學特性受離散單元間細觀力學參數影響，通過室內試驗所得到巖體力學參數無法直接應用于數值計算。常用方法是在離散元程序中復現室內試驗過程，以試錯法不斷調整接觸模型細觀參數使數值計算與室內試驗分別得到的應力-應變曲線相近，以此實現細觀參數標定[17]。采用直剪試驗獲得巖體細觀力學參數，剪切盒尺寸為30 m×30 m×30 m，顆粒數量107 202 個，顆粒半徑660～1 000 mm（均勻分布），加載速度0.15 m/s。顆粒間接觸默認為線性模型；顆粒膠結選用平行黏結模型，該模型黏結強度較高，能夠承受彎矩等載荷，一般多用此模型模擬巖石、混凝土等材料。巖體宏細觀力學參數見表1。

表1 巖體宏細觀力學參數Table 1 Macroscopic and mesoscopic parameters of rock mass

2 基于有限差分法的南幫邊坡分析

采用有限差分法對南幫端幫邊坡進行分析,通過有限差分法的計算，初步得到了邊坡失穩后的滑體范圍、應力分布與位移等信息。端幫邊坡屈服狀態如圖1，端幫邊坡最大剪應變增量云圖如圖2，端幫邊坡位移云圖如圖3。

圖2 端幫邊坡最大剪應變增量云圖Fig.2 Diagram of maximum shear strain increment of end slope

圖3 端幫邊坡位移云圖Fig.3 Diagram of end slope displacement

圖1 為邊坡在強度折減算法計算結束后單元體的屈服狀態，表示該區域巖體是否在應力作用下發生破壞，以及發生破壞的類型（紅色部分指區域巖體正受剪切作用發生破壞、藍色區域為單位正受張拉作用而破壞）。從單元區域狀態可以看出，邊坡潛滑面呈圓弧形，以剪切破壞為主，在滑體剪出口、滑體后緣處存在一定張拉破壞。其中，第四系地層受拉伸破壞最為明顯，第四系下方的灰色泥巖一定程度阻礙了剪切破壞區域與張拉破壞區域進一步貫通。

圖1 端幫邊坡屈服狀態Fig.1 Yield state of end slope

圖2 為邊坡模型計算結束后最大剪應變增量云圖，用來反映邊坡潛滑面[18]。通過最大剪應變增量同樣反映了端幫圓弧形滑動的趨勢，但滑體范圍相比圖1 有一定收束。

圖3 為邊坡位移云圖，滑體范圍與最大剪應變增量所反映出的相近，坡腳處局部位移達到29.4 m。而限于程序連續介質的假定，該位移僅為不平衡力作用在單元體后應變增量的累加，參考價值較小。

3 離散-連續耦合計算方法和模型

3.1 計算方法概述

離散-連續耦合計算基于FLAC3D6.0 程序，該版本支持PFC 離散顆粒模型與FLAC3D連續單元模型之間物理力學信息進行傳輸與交換。其中，物理力學信息的傳輸、交換是通過在離散域與連續域之間生成“耦合墻體”，耦合墻體頂點與連續域單元節點重合。單元節點在不平衡力作用下產生位移，耦合墻體頂點記錄位移過程并與單元節點保持同步，以此將連續域的信息傳輸至離散域中。離散元計算中，墻體運動對計算結果起到決定性作用。在連續域將單元節點的物理力學信息傳遞至耦合墻體頂點后，耦合墻體推動顆粒模型運動，運動后力學信息更新，墻體所受力與力矩通過等效力系統[19]分配至各個耦合墻體頂點，以此為媒介傳遞回連續域，此時1 個計算循環完成，程序最終通過顯式積分得到耦合計算結果。離散-連續耦合計算過程如圖4。

圖4 離散-連續耦合計算過程Fig.4 Discrete-continuous coupling calculation process

3.2 耦合計算模型

以圖1 的滑體范圍為基準，向邊坡內部偏移50 m，視為潛在的滑體范圍。將范圍內單元刪除，替換為離散元模型，進行離散-連續耦合計算。耦合計算模型如圖5，其中連續域內單元采用Mohr-Coulomb本構模型，離散域內顆粒接觸采用平行黏結模型，黏結失效后退化為線性模型，滑體形態由有限差分法計算結果給出。計算模型長930.25 m、高410.80 m、厚度5 m，邊坡參數與工程背景介紹相同。模型共生成單元14 924 個，顆粒302 298 個。

圖5 離散-連續耦合計算模型Fig.5 Discrete-continuous coupling calculation model

3.3 強度折減法

通過折減巖體強度反復計算，使邊坡達到臨界破壞狀態的方法稱為強度折減法[20]。邊坡在臨界破壞狀態時的強度折減系數視為其穩定性系數，是巖體原有力學強度和臨界破壞狀態強度的比值。該方法在有限差分法和離散元中均有研究[21-22]。通過有限差分法的計算已知邊坡穩定性系數為1.117，為了分析邊坡失穩過程的演化規律，且考慮到宏細觀力學參數間存有一定偏差，故在耦合計算過程中對巖體宏細觀力學參數進行1.150 的強度折減。

其中，有限差分法中強度折減法指邊坡某單元巖體抗剪強度τ 為[21]：

式中：C 為巖體黏聚力，Pa；φ 為巖體內摩擦角，（°）；σ 為法向應力，Pa。

對巖體進行強度折減，折減系數為F，此時抗剪強度τF為：

4 基于離散-連續耦合計算法的南幫邊坡分析

4.1 邊坡位移演化

邊坡水平位移演化過程如圖6，邊坡豎直方向位移云圖如圖7。

圖7 邊坡豎直方向位移云圖Fig.7 Diagram of vertical displacement of slope

圖6 可知，邊坡水平方向位移演化呈圓弧形滑動，與FLAC 計算結果相近。隨著計算步數的增大，端幫最大水平位移量不斷增大。邊坡受自重產生側向變形影響，坡面顆粒最先出現黏結破壞，向下滾落并堆積至坡腳，滾落顆粒最大位移值超過40 m。計算至4 000 步時（圖6（b）），受自重影響，滑體滑移，后緣出現張拉裂縫。邊坡水平方向位移在5～10 m，坡頂裂隙和坡腳水平位移明顯，呈現圓弧形滑動趨勢。8 000～24 000 步時（圖6（c）～圖6（g）），邊坡各區域位移隨著計算步數增大，其中水平位移大小一定程度與巖性有關，能夠體現邊坡滑體非均質特性。坡腳和520 m 水平的位移增幅較大，邊坡滑體清晰，坐落明顯，下部滑體向臨空面推擠，具有一定傾倒趨勢。計算結束時（圖6（h）），滑坡過程中煤臺階與灰色泥巖部分破壞嚴重，近似為2 個圓弧形滑面。中部灰色粉砂巖灰色泥巖互層在一定程度阻礙了上下2 部分滑面的貫通。與邊坡初始輪廓對比，邊坡整體破壞嚴重，原布置于端幫的各運輸通道均存在較大變形。

圖6 邊坡水平位移演化過程Fig.6 The evolution of horizontal displacement of slope

由圖7 可知，由最下部煤臺階裂隙，邊坡滑體可大致分為2 個部分。滑體在沿滑床向外剪出時，受坡腳顆粒堆積與逆傾構造影響而出現較大錯動，滑體大致分為上下2 個部分。其中上部分滑體后緣豎直位移較大，在25～45 m 之間，接近滑面處位移較大，上部分滑體豎直位移在向中部過程中逐漸減小直至滑體前緣；滑體前緣受下部煤層錯動影響，滑體前緣呈楔體向下滑移，位移在20～25 m 之間。下部分滑體位移較小，逆傾構造阻礙了下部分滑體進一步運動，對滑坡影響范圍起重要作用。

4.2 邊坡破壞過程

邊坡破壞過程如圖8。

由圖8 可知，隨著計算步數增大，邊坡破壞的巖體數量、破壞程度增大。在計算至4 000 步時（圖8（a）），邊坡開始顯現圓弧型滑動趨勢，巖體較明顯但相對完整，破碎巖體數量為2 453；12 000 步時（圖8（b）），下側巖體破壞嚴重，破碎巖體數量增幅較大，此時上側巖體相對完整，煤層及頂板有加劇破壞趨勢，破碎巖體數量增加至3 398；20 000 步時（圖8（c）），下側巖體破壞趨向穩定，上側巖體逐層破壞，破壞區域上移；計算結束時（圖8（d）），破碎巖體數量為4 412。

圖8 邊坡破壞過程Fig.8 Slope damage process

4.3 邊坡失穩影響范圍

邊坡失穩影響范圍如圖9。

圖9 邊坡失穩影響范圍Fig.9 Impact scope of slope slide

圖9 可知，在滑坡過程中，最上部580 m 水平出現拉裂縫，距坡肩61.37 m 內有明顯下沉，550～580 m 水平滑塌約30 m。邊坡巖層錯動明顯，下部煤臺階及頂板向臨空面處推擠，下部水平位移達42.39 m。通過耦合計算結果，初步給出了所述邊坡在滑坡災害發生后將影響的范圍，為應急疏散及邊坡安全距離確定提供一定參考。

5 結 語

1）通過室內試驗與離散元數值模擬得到所研究邊坡巖體的宏細觀力學參數，通過有限差分法軟件FLAC 計算得到連續介質下邊坡滑體的范圍。在此基礎上建立邊坡離散-連續耦合計算模型，以離散元方法建立潛在失穩區域模型，以有限差分網格建立基巖等模型，計算得到滑面位置與有限差分法近似。

2）離散-連續耦合計算得到邊坡呈圓弧形滑坡，滑坡過程最先為坡面顆粒黏結受側向變形影響發生少量破壞，黏結破壞后顆粒向下滾落并堆積在邊坡坡腳；隨后滑體后緣出現張拉裂縫，滑面逐漸清晰，滑體前緣煤層錯動，并在邊坡逆傾構造影響下阻礙滑坡進一步發展。

3）滑坡過程中滑體破壞嚴重，其中下側巖體優先破壞，集中在煤層及頂板處，隨后破壞區域上移，上側巖體逐漸破壞。滑坡結束后，邊坡最上部580 m水平出現較多張拉裂縫，距坡肩61.37 m 內有明顯下沉，約為30 m。邊坡巖層明顯錯動，下部煤臺階及頂板外臨空面推擠，下部水平位移達42.39 m。

4）離散-連續耦合計算能夠有效分析邊坡大變形問題，有著較高計算效率。計算初步給出了邊坡滑坡發生后將影響范圍，能夠為礦山應急疏散及邊坡安全距離確定提供一定參考。但限于時間與技術能力，僅在二維層面對邊坡滑坡特征及影響進行分析，未能體現露天礦工作幫與內排土場對端幫的支撐；采用強度折減法對邊坡巖體粘聚力與內摩擦角同時折減，未能體現特定環境對力學參數的不同影響。以上內容將在后續研究中陸續完善。