二維離散單元法在某舊金屬礦改造邊坡穩(wěn)定性研究中的應用

董夢龍，李運勝，李雨坤

（1.河海大學 地球科學與工程學院，江蘇 南京 210098；2.浙江省地礦勘察院，浙江 杭州 310013）

二維離散單元法在某舊金屬礦改造邊坡穩(wěn)定性研究中的應用

董夢龍1，李運勝2，李雨坤1

（1.河海大學 地球科學與工程學院，江蘇 南京 210098；2.浙江省地礦勘察院，浙江 杭州 310013）

金屬礦山露天開采后將形成規(guī)模較大的邊坡，如何處理利用廢棄礦山的邊坡資源是提高國土資源利用的重要途徑。本文以某金屬礦山廢棄采坑邊坡為例，利用二維離散單元法研究其穩(wěn)定性。該礦山邊坡為巖質邊坡與松散堆積體（尾礦）的結合體，局部為陡崖和陡坡，工程地質形態(tài)復雜。研究選取了2個具有代表性的剖面，分析計算了天然工況、暴雨工況和地震工況下的邊坡穩(wěn)定性，得出相應的穩(wěn)定系數(shù)，并預測了邊坡可能發(fā)生的失穩(wěn)模式，并提出了相應的治理意見。研究表明，二維離散單元法能適用于復雜形態(tài)下的邊坡穩(wěn)定性分析，并可作為條分法計算的重要補充。

采坑邊坡；穩(wěn)定分析；二維離散單元法

0 引言

邊坡在自然界以及人類工程中及其常見，是工程建設中最重要的地質環(huán)境之一，也是工程建設中最常見的工程形式。由于人類工程建設活動的不斷增加，在礦山、水利、交通以及城建等部門都涉及到大量的邊坡問題。我國是一個礦業(yè)大國，現(xiàn)有大量的金屬礦山采用露天開采，從而形成規(guī)模巨大的邊坡工程，有的甚至高達300 m以上。而這些金屬礦邊坡的穩(wěn)定性問題十分嚴重，一旦失穩(wěn)將對生產造成嚴重影響，并帶來巨大的經濟損失。近年來金屬礦邊坡失穩(wěn)滑坡事故時有發(fā)生，例如，本鋼南芬露天礦在1999年5月至2002年7月的短短3年時間內，就發(fā)生4次近百萬方的大型滑坡災害，造成了數(shù)千萬元的經濟損失；酒泉鋼鐵公司的黑溝鐵礦在其基建期間發(fā)生了一次極為嚴重滑坡和泥石流災害，該地質災害直接堵塞了酒泉市、嘉峪關市兩市唯一的水源北大河，造成直接經濟損失高達4 000余萬元；2004年四川省雅安市寶興縣隴東鎮(zhèn)宇通礦山的“10·18”特大巖體垮塌事故造成巨大的經濟損失和人員傷亡[1-2]。這些事故無一不說明了礦山邊坡工程穩(wěn)定性的重要性。另外，隨著一些礦山生產周期結束，廢棄露采礦山的綜合整治也是提高國土資源利用的重要工程途徑。

邊坡穩(wěn)定性分析的方法很多，但總的說來可分為兩大類，即以極限平衡理論為基礎的條分法和彈塑性理論為基礎的數(shù)值計算方法，其中離散單元法為現(xiàn)階段數(shù)值分析方法中發(fā)展較快的一種。離散單元法是分析離散介質變形和運動趨勢的有力工具，自1971年Cundall提出離散單元模型以來，就廣泛的應用于各類工程的數(shù)值模擬和分析研究中[3]。李作良等采用離散單元法對鎢礦采空區(qū)進行數(shù)值模擬，分析其穩(wěn)定狀況[4]。張國慶等利用UDEC二維離散元軟件對傾斜工作面開采導致的覆巖移動及裂隙發(fā)育過程進行數(shù)值模擬，模擬了開采過程中覆巖移動及裂隙發(fā)育規(guī)律[5]。

本文以某舊鐵礦采坑邊坡為例，研究了離散單元法對該類邊坡的穩(wěn)定評價的適用性。

1 工程概況

1.1 工程簡介

某景區(qū)擬利用某舊鐵礦的采坑進行改造修建景區(qū)建筑物（圖1），達到廢棄礦山綜合改造，提高土地資源利用的目的。擬利用的采坑深約30 m，采坑頂部至周邊山體頂部高差30～100 m，采坑底部至周邊山體的頂部高差60～130 m不等，坡度在20°～45°之間，巖層局部直立后反傾，該邊坡的穩(wěn)定性與項目主體建筑及附屬設施的安全直接相關，所以對該邊坡工程進行穩(wěn)定性分析和評價是十分必要，而且具有重要的工程意義。

圖1 舊鐵礦采坑全景圖Fig.1 Panoramic view of the old iron ore pit

1.2 區(qū)域地質構造

根據(jù)研究區(qū)域的區(qū)域地質資料，所研究的舊鐵礦邊坡位于揚子準地臺下?lián)P子臺褶帶，區(qū)域構造帶屬淮陽山字形反射弧之西段。區(qū)域內的構造運動主要受燕山-喜馬拉雅運動的影響。燕山構造期內使得區(qū)域內的褶皺活動和斷層的升降活動劇烈；而喜馬拉雅期內的構造運動相對較弱，對區(qū)域內的構造影響相對較小。所以研究區(qū)域內主要的構造運動為褶皺與斷層。

1.3 地層巖性

研究區(qū)域內地層屬寧蕪地層小區(qū)北部，以中生代侏羅系地層分布為主。出露地層的巖性主要為侏羅系中風化至全強風化凝灰?guī)r和第四系松散堆積層，上覆第四系主要以人工填土層（Q4ml）及坡積層（Q3dl）。

1.4 邊坡現(xiàn)狀

研究區(qū)域東西兩側的最高海拔標高約247.5 m與201.6 m，中部為連接兩峰的鞍形山脊，后鞍形山脊處開山采石，形成了南北長約462 m，東西寬約294 m的采礦坑，坑底高程在110 m，南北兩側高程在160～170 m之間，采坑四周均為坡度較陡的山坡。由此形成高約50～130 m不等的高邊坡，坡度平均在35°左右，局部直立或反傾。

礦坑西側、北側和礦坑南部高平臺西南側和東南側共有五個尾礦渣堆積體。由尾礦渣堆積而成，厚度在5～35 m，土質不均，顆粒成分從黏性土到漂石均有，結構從松散到致密不均，滲透性好。礦坑南側高平臺高程在157～163 m之間，呈南北向長條形狀，面積9 000 m2。局部尾礦渣堆積體有變形跡象。

所以該研究邊坡為巖質邊坡和礦渣堆積體相結合的邊坡，局部為開礦形成的陡崖和陡坡，離散狀態(tài)的坡體體尤其適合離散單元法分析穩(wěn)定狀況。

工程地質平面圖見圖2。

圖2 工程地質平面圖Fig.2 Engineering geological plan

2 計算原理與軟件介紹

離散元方法針對由不同的小塊體所構成的整體進行應力、應變的分析計算，而對于各不同塊體之間的連接是通過接觸點的耦合的模擬而互相連接在一起。對于巖體這一特殊的材料體分析，其主要的組成部分為巖石和結構面兩大部分。而巖體中的結構面的剛度和強度相比于巖石本身的強度要小得多?？紤]到巖體這一特殊性質，通常將巖體中的巖石假定為剛性體以減少不確定性（自由度）的數(shù)量，所以巖體產生的變形主要為巖體中結構面的總位移所產生，而結構面的位移是由各接觸點（面）的變形所引起[6-8]。

對于巖體的研究是將其假定為各種離散塊體的堆砌，而這些塊體之間的相互作用力是根據(jù)位移和力的關系式來求解，其中每一個塊體的運動都遵循牛頓運動的基本定律—即力和力矩的平衡定律。三大方程（平衡方程、變形協(xié)調方程和本構方程）是數(shù)值分析模型的建立的基礎，模型同時還必須滿足相應的邊界條件。而離散元單個塊體之間是不存在相互的變形約束，每個塊體僅需要滿足的是以下的物理方程和運動方程。

2.1 物理方程

離散單元法的物理方程是表述力與位移之間的相互關系。塊體之間的力的法向分力為Fn與它們之間的法向位移量Sn為正比關系，即：

式中：系數(shù)Kn為法向剛度系數(shù)。

塊體之間的力切向的分量Ft與切向的相對位移量St的關系：

式中：系數(shù)Kt為接觸面的切向剛度系數(shù)[5-6]。

2.2 運動方程

對應塊體所受的一組力F1，…，F(xiàn)n，產生合力F合和合力距M，由牛頓第二運動定律確定該塊體的加速度和角速度，然后可以確定時間步長Δt內的速度和角速度以及位移和轉動量。

式中：a為加速度，m為巖塊的質量。

上式應該理解為分別對X方向和Y方向求加速度。然后分別進行向前差分的數(shù)值積分，可以得到巖塊X方向和Y方向的速度和位移：

式中：V（t1）為t1時間的速度；U（t1）為t1時間的位移，t0為起始時間，Δt為時間步長。

計算時依據(jù)時步迭代并遍歷每個塊體，直到每個塊體不再出現(xiàn)不平衡力和不平衡力矩為止。

離散單元法區(qū)別于其他數(shù)值方法的主要特征是：（1）它能反映塊體之間接觸面的滑移、分離和翻轉等不連續(xù)變形，還能計算塊體內部的變形和應力分布狀態(tài)。（2）可以用于求解非線性大位移和動力穩(wěn)定問題，這一應用是根據(jù)顯示時間差分解法來求解動力平衡方程[9]。

2.3 軟件介紹

UDEC（Universal Distinct Element Code）是由美國Itasca軟件公司開發(fā)的，是在離散單元法基礎上編制的通用的離散元程序。UDEC將不連續(xù)介質視為離散塊體的集合，不連續(xù)性則看作塊體之間的邊界條件。本文主要采用針對非連續(xù)介質模擬的離散元數(shù)值計算程序UDEC4.0。

3 數(shù)值建模

3.1 計算工況

結合坡體特征與坡體所處地質環(huán)境，該邊坡穩(wěn)定性計算工況見表1。根據(jù)工程區(qū)氣候資料，區(qū)內雨量充沛，年降水1 200 mm，年平均降水量1 106 mm，降水最多季節(jié)為7月份，494.5 mm，最少為2月份，僅7.7 mm。所以暴雨工況下降水量計算值取160 mm/d。根據(jù)工程區(qū)地震資料，邊坡設計地震加速度值取0.10 g。

表1 邊坡穩(wěn)定性計算工況Tab.1 Slope stability calculation

3.2 計算參數(shù)

強度參數(shù)的選取合理與否，對邊坡穩(wěn)定性計算起關鍵性作用。邊坡失穩(wěn)形成后，在新的環(huán)境條件下，結構遭到破壞、抗剪強度降低至殘余強度的滑帶土，在上覆滑體自重應力作用下發(fā)生壓密，其殘余強度有一定的恢復特性；而在暴雨工況下，滑體土通常并不能達到飽和，對滑帶土的軟化作用也往往不會完全達到飽和軟化的程度。根據(jù)滑帶土室內試驗結果，綜合考慮邊坡失穩(wěn)時的土體狀態(tài)，以天然和飽和狀態(tài)下的峰值和殘余強度試驗值為基礎綜合選取。

（1）自然工況下，滑動面抗剪強度參數(shù)以滑帶土天然狀態(tài)抗剪峰值強度與殘余強度之間取值，接近于峰值強度；

（2）暴雨工況下，考慮到滑動面巖土一般并不會達到飽和，抗剪強度參數(shù)以滑帶土飽和狀態(tài)抗剪峰值強度與殘余強度之間、接近于殘余強度取值；

（3）當計算滑動面位于滑體土內時，自然、暴雨工況下，抗剪強度參數(shù)分別以滑體土天然和飽和抗剪強度計算；當計算滑動面沿滑帶土和部分穿越滑體土時，分段計算。

根據(jù)該邊坡的工程地質詳細勘察報告，結合巖體力學參數(shù)實驗，計算模型所選取的巖土體物理力學參數(shù)和結構面力學參數(shù)見表2和表3。

表2 不同巖層物理力學參數(shù)取值Tab.2 Physical and mechanical parameters of different strata

表3 結構面物理力學參數(shù)取值Tab.3 Physical and mechanical parameters of structural plane

3.3 計算模型

對該工程實例的二維離散元穩(wěn)定性分析主要是針對該礦坑的南部廢石堆積體和北部廢石堆積體，見圖2，所以計算選取了該邊坡具有代表性的剖面3-3'、6-6'。

圖3為DP3-3'剖面。由圖3可見，坡表主要由殘坡積土和全、強風化凝灰?guī)r組成，下覆基巖為凝灰?guī)r。巖石受通過場區(qū)的兩條構造影響，節(jié)理裂隙發(fā)育，主要方向為NNE向、NE向、NWW向，上述節(jié)理裂隙延伸長度較長，除NNW傾角較小外，其余各節(jié)理裂隙傾角均較陡。邊坡的淺表層，巖石的表面張性裂隙發(fā)育，局部直立或反傾。

根據(jù)DP3-3'剖面資料，建立相應數(shù)值模型。根據(jù)巖體風化程度，模型所選取巖體材料主要有全強風化凝灰?guī)r、中風化凝灰?guī)r和新鮮凝灰?guī)r三類。具體巖土體材料參數(shù)見表2。根據(jù)所模擬巖土體特性，模型所采用本構模型為摩爾-庫倫破壞模型。模型中設置結構面主要考慮有層面和非層面兩類結構面。結構面參數(shù)見表3。結構面同樣采用摩爾-庫倫破壞模型。模型尺寸高100 m，寬150 m；共劃分block 287個，zone 2 976個。模型底部、正面與背面邊界條件設置為固定位移邊界（即zero-velocity boundary）以控制模型底部與正面、背面的變形。頂部設置為自由邊界（即free boundary）。垂直方向初始地應力由巖土體自身重力決定；設置兩方向水平初始地應力為垂直向的一半。具體計算模型見圖4。計算中迭代至不平衡力趨近于0，該模型中設置迭代至不平衡力的數(shù)量級達到10-5。

圖4 DP3-3'剖面離散元計算模型Fig.4 DP3-3'discrete element model

圖5為DP6-6'剖面，坡表主要由殘坡積土和全、強風化凝灰?guī)r組成，下覆基巖為凝灰?guī)r。巖石受通過場區(qū)的兩條構造影響，節(jié)理裂隙發(fā)育，主要方向為NNE向、NE向、NWW向，上述節(jié)理裂隙延伸長度較長，除NNW傾角較小外，其余各節(jié)理裂隙傾角均較陡。邊坡的淺表層，巖石的表面張性裂隙發(fā)育，局部直立或反傾。

圖5 邊坡計算剖面DP6-6'Fig.5 Slope calculation section DP6-6'

根據(jù)DP6-6'剖面資料，建立相應數(shù)值模型。根據(jù)巖體風化程度，模型所選取巖體材料主要有全強風化凝灰?guī)r和新鮮凝灰?guī)r兩類。具體巖土體材料參數(shù)見表2。根據(jù)所模擬巖土體特性，模型所采用本構模型為摩爾-庫倫破壞模型。模型中設置結構面主要考慮有層面和非層面的結構面兩類。結構面參數(shù)見表3。結構面同樣采用摩爾-庫倫破壞模型。模型尺寸高135 m，寬195 m；共劃分block 372個，zone 3 615個。其他設置與DP3-3'剖面相同。具體計算模型見圖6。

圖6 DP6-6'剖面離散元計算模型Fig.6 DP6-6'discrete element model

4 結果分析

4.1 DP3-3'剖面計算結果分析

DP3-3'剖面在天然工況、暴雨工況（160 mm/d）和地震工況（水平加速度系數(shù)為0.1）下，邊坡穩(wěn)定性計算結果見圖7、圖8、圖9。

計算結果表明DP3-3'剖面潛在滑動面位于邊坡中部，全、強風化凝灰?guī)r與基巖分界面，水平埋深約9 m，豎直方向埋深約6 m。天然工況下，邊坡的穩(wěn)定系數(shù)為1.10，邊坡巖土體的最大位移約8 cm，邊坡整體穩(wěn)定性較好，不會發(fā)生滑動破壞。暴雨工況下，邊坡的穩(wěn)定系數(shù)為0.93，邊坡穩(wěn)定性差，極易發(fā)生滑坡。地震工況下，邊坡的穩(wěn)定系數(shù)為1.01，邊坡處于極限平衡狀態(tài)，考慮到邊坡表層巖體大量發(fā)育有風化裂隙、爆破裂隙等，邊坡易發(fā)生局部崩塌、塊體掉落等災害。

圖7 天然工況下局部失穩(wěn)矢量圖Fig.7 Local instability vector diagram under natural condition

圖8 暴雨工況局部失穩(wěn)矢量圖Fig.8 Local instability vector diagram under rainstorm condition

圖9 地震工況下局部失穩(wěn)矢量圖Fig.9 Local instability vector under seismic condition

4.2 DP6-6'剖面計算結果分析

DP6-6'剖面在天然工況、暴雨工況（160 mm/d）和地震工況（水平加速度系數(shù)為0.1）下，邊坡穩(wěn)定性計算結果見圖10、圖11、圖12。

圖10 天然工況下局部失穩(wěn)矢量圖Fig.10 Local instability vector diagram under natural condition

圖11 暴雨工況下局部失穩(wěn)矢量圖Fig.11 Local instability vector diagram under rainstorm condition

圖12 地震工況下局部失穩(wěn)矢量圖Fig.12 Local instability vector under seismic condition

計算結果表明，DP6-6'剖面的潛在失穩(wěn)區(qū)域集中在邊坡頂部，潛在滑動面的埋深約15 m。天然狀況下，邊坡的穩(wěn)定系數(shù)為1.08，邊坡整體穩(wěn)定，不會發(fā)生破壞。在暴雨條件下，由于雨水的入滲產生的向坡外的滲流力以及降水導致的巖土體力學參數(shù)的降低，邊坡穩(wěn)定性變差，穩(wěn)定系數(shù)為0.92，可能發(fā)生滑動破壞。地震工況下，邊坡的穩(wěn)定系數(shù)為0.97，邊坡處于極限平衡狀態(tài)，易發(fā)生傾倒、崩塌等地質災害。

4.3與條分法計算結果比較

根據(jù)有關邊坡工程技術規(guī)范[10-11]，研究同時采用摩根斯坦-普賴斯法（Mogenstem-price）和畢肖普法（Bishop）對兩個剖面進行穩(wěn)定性分析，用來驗證離散單元法對該類邊坡分析的合理性。計算結果見表4和表5。

表4 DP3-3'剖面不同計算方法穩(wěn)定系數(shù)比較Tab.4 Comparison of stability coefficient at DP3-3'profile

表5 DP6-6'剖面不同計算方法穩(wěn)定系數(shù)比較Tab.5 Comparison of stability coefficient at DP6-6'profile

通過計算結果比較表明，二維離散單元法得出的邊坡穩(wěn)定系數(shù)與二維極限平衡計算結果相近，說明離散單元法對該類邊坡的穩(wěn)定分析是合理的、適用的、也是更嚴格的。

根據(jù)計算結果和不同的失穩(wěn)模式為邊坡治理方案提供了依據(jù)。

5 結語

本文簡要介紹了二維離散單元法，并將該方法運用到舊金屬礦廢礦邊坡穩(wěn)定分析這一工程實際問題上。通過2個典型剖面的離散元法穩(wěn)定性計算分析，得出穩(wěn)定系數(shù)，并預測了邊坡可能發(fā)生的失穩(wěn)模式，解決了實際工程問題。可見二維離散單元法對于舊金屬礦邊坡的穩(wěn)定評價是適用的、可行的，這對目前很多舊金屬礦的二次開發(fā)利用具有一定的工程指導意義。

[1] 金屬非金屬礦山露天礦邊坡調查報告[R].馬鞍山礦山研究院，2006.12. Investigation report of open pit slope in metal and nonmetal mines [R].Sinosteel Maanshan Institute of Mining Research Co.，2006.12.

[2] 程東幸，劉大安.中國典型重大邊坡工程穩(wěn)定性與安全評價現(xiàn)狀研究[J].西部探工，2008（12）：33-37. CHEN Dongxing，LIU Daan.Studyon stability and safety evaluation of typical slope engineering in China [J].Western Exploration，2008（12）：33-37.

[3] Cundall PA（1988）Formulation of a three-dimensional distinct element model-Part 1.A scheme to detect and represent contacts in a system composed of many polyhedral blocks[J].Int J Rock Mech Min Sci&Geomech Abstr.，25：107-116.

[4] 李作良，李一帆.某鎢礦采空區(qū)穩(wěn)定性的離散元數(shù)值模擬[J].中國鎢業(yè)，2008，23（6）：1-4. LI Zuoliang，LI Yifan.Discrete element numerical simulation of stability in tungsten mine stope [J].China Tungsten Industry，2008，23（6）：1-4.

[5] 張國慶，黃 霆.某礦UDEC數(shù)值模擬研究[J].煤，2010，19（6）：21-23. ZHANG Guoqing，HUANG Ting.Study of numerical simulation on UDEC in one mine[J].Coal，2010，19（6）：21-23.

[6] FALLS S D，YOUNG R P.Acoustic emission and ultrasonic-velocity methods used to characterize the excavation disturbance associated with deep tunnels in hard rock [J].Tectonophysics，1998：178-181.

[7] SHENG Q，YUE Z Q，LEE C F.Estimating the excavation disturbed zone in the permanent shiplock slopes of the Three Gorges Project China[J].International Journal of Rock Mechanics &Mining Sciences，2002：86-89.

[8] GRIFFITHS DV，LANE PA.Slope stability analysis by finite elements[J].Geotechnique，1999：209-215.

[9] BROWN E T，HOEK E.Trends in relationships between measured in-situ stress and depth[J].Int.J.Rock Mech.Min.Sci.&Geomech. Abstr，1978（15）：211-215.

[10] 重慶市設計院.建筑邊坡工程設計規(guī)范[S].GB50330—2013.北京：中國建筑工業(yè)出版社，2014.

[11]水利部水利水電規(guī)劃設計總院.水利水電工程邊坡設計規(guī)范[S].北京：中國水利水電出版社，2009.

Application of Two Dimensional Discrete Element Method to the Slope Stability Analysis of an Exhausted Metallic Mine

DONG Menglong1,LI Yunsheng2,LI Yukun1
(1.School of Earth Sciences and Engineering,Hohai University,Nanjing 210098,Jiangsu,China;2.Zhejiang Geological Survey,Hangzhou 310013,Zhejiang, China)

The recycle of slope resources in the exhausted metallic mines is an important way to improve the use of land resources.This paper studies the stability of the abandoned slope in metallic mine by applying two-dimensional discrete element method.The mine slope is a combination of rock slope and loose accumulation body(tailings)with complex engineering geological form.After choosing two representative profiles,we calculate the corresponding slope stability coefficients under the conditions of norm,rainstorm and earthquake.The potential slope instability modes and corresponding treatment technologies are predicted.The results show that the two-dimensional discrete element method,as an important supplement to the calculation of strip method,can be applied to the stability analysis of slope under complex form.

metallic mine slope;stability analysis;two dimensional discrete element method

TD22；P642

A

10.3969/j.issn.1009－0622.2017.01.006

2017-01-11

董夢龍（1988-），男，安徽滁州人，博士研究生，主要從事巖體結構與工程穩(wěn)定研究。