南芬露天礦滑坡機理及其動態過程

王長軍，曹春暉，陶志剛

(1.本鋼歪頭山鐵礦，遼寧 本溪 117005；2.深部巖土力學與地下工程國家重點實驗室，北京 100083；3.中國礦業大學(北京)力學與建筑工程學院，北京 100083)

隨著露天礦開采深度的增加，礦山邊坡面臨一系列地質災害問題[1-2]，比如滑坡、泥石流、崩塌、沉降、臺階失穩等。滑坡是露天礦邊坡常見的地質災害，往往會造成嚴重的經濟損失和人身財產損失。滑坡己經成為影響礦山安全生產和長久健康運營的重要災害問題，迫切需要對露天礦滑坡的機理和動力學過程做更深入的研究。

國內外學者對不同的露天礦滑坡機理進行了分析研究，開挖擾動、凍融循環、降雨、地震等是露天礦滑坡發生的主要原因[3-6]。多數研究僅考慮其中一種或少數幾種因素，滑坡機理缺少系統性綜合分析。另外，對露天礦滑坡的動力學全過程分析相對較少。受前人研究的啟發，本文以2019年1月4日發生在南芬露天礦下盤394～502 m臺階中部老滑坡體區域的滑坡為例，對其機理進行研究。結合現場勘查、地質資料，并采用光滑粒子流體動力學(SPH)方法結合Mohr-Coulomb彈塑性模型，建立三維數值計算模型，對露天礦滑坡進行動力學全過程分析，以了解其運動過程。SPH方法[7]最初用來研究天體物理的運行規律，并逐漸應用于流體力學、固體力學、沖擊爆炸等領域。作為一種無網格純拉格朗日方法，SPH方法在處理大變形問題如滑坡時具有獨特優勢，在巖土工程領域得到了廣泛應用，并取得了一系列進展[8-9]。相對于有限元等方法，SPH方法由于無網格，可以避免大變形過程中產生的網格畸形等問題。而一些離散元方法，在處理像露天礦滑坡這類工程尺度的大變形問題時，對計算機內存等要求較高，計算效率面臨挑戰。SPH方法以其高效率和可靠性，可以很好地再現露天礦滑坡的全過程，并根據數值結果為露天礦邊坡劃定危險區域，制定防災減災措施，為露天礦安全生產和運營提供指導性建議。

1 滑坡背景

1.1 滑坡介紹

本溪市南芬露天鐵礦是亞洲最大的單體露天礦山。近20年來，在特殊地形和長期礦山開采綜合影響下，采場下幫邊坡發生過多次滑坡，壓礦1 000多萬t，嚴重影響了礦山的正常開采。

2019年1月4日在南芬露天礦下盤394～502 m臺階中部老滑坡體區域發生滑坡，見圖1。滑體對礦山安全作業帶來了嚴重影響。由于及時預警，所幸未造成人身和設備損失，但滑落的碎石破壞了部分排水管路，阻擋了采礦車輛的正常出入，影響了正常的開采作業。本次滑坡體積大約為86 400 m3。露天礦邊坡整體較陡，滑坡體滑動方向約為264°。滑坡體厚度小于10 m，屬于淺層滑坡。滑坡源區上部寬度約為50 m，下部寬度約為140 m，上下垂直落差大約110 m。 滑坡傳播區垂直落差大約120 m，底部寬度約為240 m。露天鐵礦滑坡源區為老滑坡體堆積區，老滑坡發生后滑坡體整體呈松散狀態掛靠于邊坡上。老滑坡體堆積層主要為碎石土，整體穩定性較低，在爆破振動，開挖、降雨、凍融等因素下極易發生滑坡。

圖1 滑坡特征Fig.1 Feature of the landslide

1.2 地質和氣象水文資料

南芬礦區位于華北地臺遼東臺背斜營口-寬甸隆起的北緣太子河凹陷之中，新華夏系第二隆起與東西天山-陰山構造帶的交匯處。地貌為單斜構造，主要由變質巖地層組成，屬中高山侵蝕構造地貌。礦區內廣泛發育太古界鞍山群，其次為元古界遼河群和震旦系地層及新生界第四系地層。南芬露天鐵礦屬于大型鞍山式沉積變質鐵礦，含鐵巖段主要賦存于太古界鞍山群，表現為單斜構造。

該區域水文網絡不發育。冬天降雨不多，夏季降雨量增加，河水變得豐富。大氣降水是地表水和地下水的主要來源。年平均降雨量、最大降雨量和最小降雨量分別為848～856 mm、1 259.8 mm和518.5 mm。汛期降雨量占全年降雨量的70%。 該區地處北溫帶季風氣候區，四季分明。 年平均氣溫為7.7～8.2 ℃，最低氣溫為-24 ℃。 霜期(降雪期)一般為11月下旬至次年3月，最大凍結深度為1.49 m。

1.3 滑坡體土體成分

南芬露天礦下盤滑坡源主要為老滑坡堆積體，這些堆積的碎石土主要成分為礫石、塊石、礦石碎屑、細砂和黏土。老滑坡堆積體碎石土的主要巖性為底盤角閃巖、綠泥片巖、綠泥角閃片巖等。碎石土黏度較低，顆粒強度高，水穩定性好。經過現場勘查，結合室內和室外試驗，給出滑坡源松散堆積物的各級顆粒含量，見表1。

表1 不同級別顆粒組成情況Table 1 Composition of different grades of particles

2 滑坡機理分析

南芬露天礦年度剝巖量和鐵礦石開采量巨大，經過長年的開采，逐漸形成了超高邊坡。這些高邊坡在自重、降雨、爆破振動、機械開挖等因素作用下，其穩定性逐漸降低，并成為礦區安全生產的重大隱患[10]。露天礦上盤邊坡以逆傾向巖體為主，下盤邊坡則以順傾向巖體為主。在特殊地形和長期礦山開采綜合影響下，采場下幫邊坡形成了多處較大規模的滑坡體；其中，下盤394～502 m臺階中部老滑坡體區域，由于其為散體堆積體，穩定性較差，成為滑坡研究的重點區域。在2019年1月4日滑坡發生前，該區域曾發生明顯的變形；其中，在502臺階出現長度約為6 m的裂縫，最寬處20～30 cm(圖2(a))。主裂縫附近伴有數條寬度2～20 mm的細小裂縫，長度不等，這些裂縫使得降雨入滲變得容易。在降雨、振動及自重等作用下，這些裂紋逐漸擴展，使得降雨入滲進一步增強，邊坡穩定性逐漸降低。另外，在滑坡發生前，502 m臺階滑移量為7～8 m(圖2(b))。

“滑坡牛頓力空地監測預警系統”是深部巖土力學與地下工程國家重點實驗室何滿潮院士團隊自主研發的滑坡監測預警系統。該系統集滑坡監測、預警、支護加固等技術為一體。牛頓力[11]的變化可以作為判斷災變發生的條件來對滑坡進行預警，可以得到較為準確的滑坡臨滑信息。從本文滑坡發生前的牛頓力監測曲線(圖2(c))變化可以看出，在2019年1月3日，牛頓力出現突變，指示滑坡可能即將發生；1月4日滑坡發生，驗證了牛頓力監測系統的可靠性。

圖2 滑坡前兆與滑坡牛頓力監測曲線Fig.2 Landslide precursor and Newton force monitoring curve of landslide

露天礦邊坡的穩定性往往受到開挖巖石結構地質的顯著影響，如層理、節理和斷層，這些不連續面經過爆破振動、風化等活動逐漸破碎成更小的塊體顆粒。這些破碎的碎石土一部分堆積在邊坡臺階上，一部分成為老滑坡體的堆積體。老滑坡體堆積在露天礦臺階上，結構松散，穩定性差。

巖質邊坡的開挖通常涉及爆破。露天礦爆破的基本要求是產生一個碎渣堆，以適應挖掘作業以及運輸設備。在多層理、節理和高裂隙巖石中，爆炸氣體噴入、楔入并擴展原有的裂縫。總體破碎程度受構造地質控制。在爆破過程中，受施工技術等多方面制約，臺階完整性和穩定性會隨著爆破振動而降低。而老堆積體的松散巖體，受爆破振動影響，內部構造會受到擾動，并處在一個不斷變化調整的過程。在爆破振動作用下，老堆積體內的碎石土會逐漸松散，其穩定性亦隨之逐漸降低，當土體內部塊石之間的互鎖作用減弱，土體黏聚力和摩擦力降低到某一值時，滑坡就會發生。

坡腳開挖及臺階碎石土的清理，會降低露天礦邊坡穩定性。在老滑坡體下部開挖或清理臺階碎石土時(清理浮石，圖3)，上部土體失去一部分阻擋和支持作用，坡腳開挖處的部分臨空土體受到擾動而滑落，并牽拉上部土體，當上部土體的自重產生的下滑力大于抵擋下滑作用的摩擦阻力時，滑坡隨之發生。

圖3 浮石清理Fig.3 Cleaning of float stone

降雨通常被認為是滑坡最關鍵的誘發因素[16]。老滑坡體內含有細顆粒物，這些細顆粒物可以吸收水分，部分顆粒遇水膨脹，使得土體的持水能力增強。老滑坡堆積體內的水分可以降低顆粒間的摩擦力。在沖刷和滲透作用下，細顆粒會在骨架中遷移并被雨水帶走，粗顆粒物之間的間隙增大，從而改變原土體結構，而降雨入滲變得更加容易。當有持續降水時，土體內部的水壓力不斷增大，導致排土場土體抗剪強度和穩定性下降，并誘發滑坡的發生。此外，凍融也是導致滑坡發生一個原因。凍土層的膨脹會使土體局部隆起；當凍土層融化時，會引起土體收縮。這種膨脹和收縮導致排土場局部損傷、薄弱面和裂縫的增長和擴展。最后，坡面出現不均勻的崩塌和沉降。土體內部的某些凍土層在某些條件下成為軟弱面或者潛在滑動面，在自重或外部荷載擾動下，滑坡可能會發生。綜上所述，南芬露天礦排土場滑坡是爆破振動、坡腳開挖、降雨、凍融等綜合作用下發生的。綜合理解滑坡機理，對于制定露天礦防災減災措施具有重要意義。

3 滑坡動力學全過程分析

3.1 SPH理論和Mohr-Coulomb模型

SPH法[13]是將計算域離散成一定數量的粒子，每個粒子占據一定質量和體積，并攜帶密度、速度、加速度和壓力、應力等物理參數。對于域內任意的函數f(x)，其核近似形式可以表示為式(1)[17]。

(1)

式中：h為平滑長度；W為平滑函數。

滑坡過程中，粒子的運動和質量守恒的控制方程可以寫成式(2)和式(3)。

(2)

(3)

式中：ρ、vα、t、σαβ、xβ和fα分別為密度、速度、時間、應力張量、位置和體力；α和β為笛卡爾坐標方向。

老滑坡體為松散巖土顆粒，本研究采用彈塑性Mohr-Coulomb本構模型來研究其力學特性，其本構模型可以用三個不變量來表示，見式(4)。

f(I1,J2,J3)=Rmcq-ptanφ-c=0

(4)

其流動勢函數g可以表示為式(5)[18]。

(5)

式中：I1為應力第一不變量；J2和J3分別為偏應力第二不變量和第三不變量；φ為材料摩擦角；c為黏聚力；p為等效壓應力；q為Mises等效應力；為子午線偏心率，它定義了雙曲函數接近漸近線的速率；c0為初始黏聚力屈服應力；ψ為在高圍壓下在p-Rmwq平面上測得的膨脹角；Rmc和Rmw為分別為和φ以及ψ有關的函數。通過引入應力張量和應變率張量的加法分解以及塑性屈服理論，并經過一系列計算(推導過程參考[19])，SPH形式的應力應變關系可以表示為式(6)。

(6)

3.2 滑坡動力學過程數值模擬

露天礦邊坡因其規模較大，研究人員常常采用數值方法研究其穩定性。為了清楚地了解露天礦滑坡動力學過程，本文選取2019年1月4日發生在南芬露天礦下盤394～502 m臺階中部老滑坡體區域的滑坡為例，建立數值模型并進行全過程分析。邊坡巖體相對比較完整，其黏聚力為20 740 kPa，摩擦角為37°。而附著在邊坡表面的老滑坡堆積體為松散的碎石土，其黏聚力大大降低，而摩擦角也適當減小。南芬露天礦排土場滑坡土體形態、粒徑和該露天礦邊坡老滑坡堆積體亦較為接近。參考南芬露天礦排土場發生的滑坡土體的黏聚力和摩擦角，并根據中國科學院武漢巖土力學研究所進行的室內和現場試驗，結合露天礦工程地質資料，并根據試錯法獲得使數值模擬結果和實際滑坡相符的力學參數(表2)。

表2 SPH滑坡模型材料屬性Table 2 Material properties of the SPH model

露天礦滑坡SPH數值模型尺寸如圖4所示。滑坡體上部寬50 m，下部寬140 m，上下垂直落差為105 m。整個模型尺寸為：長412 m，寬300 m，高288 m。圖4中深色部分為滑坡源，在數值計算開始時自動轉換成粒子，粒子總數為15 839。采用ABAQUS/Explicit計算，計算時長為100 s。

圖4 SPH滑坡模型尺寸Fig.4 Model size of the landslide by the SPH method

模型邊界條件施加如下：限制模型左右兩側Z方向的位移；限制模型前后兩側X方向位移；限制模型下方三個方向的位移。整個模型施加重力荷載。滑坡土體和邊坡巖體直接設置通用接觸，并設置摩擦系數為0.75，以模擬實際接觸情況。

圖5為露天礦滑坡的真實形態和SPH數值模擬得到的形態。由圖5可知，數值模擬形態和真實形態比較吻合。松散巖土體在露天礦邊坡臺階上向下滑動，滑動形態呈倒V字形。由數值結果可知，SPH模型可以很好地模擬露天礦滑坡，并可以根據模擬結果滑動危險區域范圍，為露天礦邊坡制定防災減災策略。

圖5 真實滑坡與SPH模擬對比Fig.5 Comparison between real landslide and SPH simulation

為了更好地了解露天礦滑坡動力學全過程，選取t=0 s、10 s、20 s、30 s、50 s和100 s來分析不同時刻的滑坡形態，如圖6所示。由圖6可知，土體顆粒大概在t=20～30 s時滑動到最低處臺階上；從t=50 s起，滑坡形態不再有很大變化，表明滑坡逐漸趨于停止。另外，從圖6可以清晰看出，粒子在滑坡運動過程中逐漸向兩側擴散，這使得滑坡的最終形態呈現近似三角形，這與實際情況相符。

為了了解滑坡過程中土體顆粒的位移信息，露天礦滑坡在t=0 s、10 s、20 s、30 s、50 s和100 s的位移云圖如圖7所示。圖7顯示位移在中間線位置處為最大，并呈現出左右對稱的分布形式。

圖7 露天礦滑坡的總位移云圖Fig.7 Total displacement(U-Magnitude) contour of the landslides for the full high waste dump

圖8給出了露天礦滑坡從t=1 s到t=100 s的速度云圖。由圖8可知，粒子經歷了先加速后減速至零的過程。從圖8(d)和8(e)對比發現，粒子的速度云圖形態基本一致，但是圖8(e)上顯示帶有速度的粒子個數明顯比圖8(d)的數量少。其原因是，粒子在露天礦邊坡上向下滑動時，一部分粒子存留在了臺階上，這使得滑動層變薄，進而使得圖8(e)中的顯示帶有速度的粒子數目明顯減少。在實際露天礦滑坡過程中，土體一邊下滑一邊在各級臺階上堆積。數值模擬結果和實際滑坡表現出了很好的一致性。

圖8 露天礦滑坡的速度云圖Fig.8 Velocity (V-Magnitude) contour of the landslides for the full high waste dump

總的來說，SPH方法憑借其自身的高計算效率和占用較低計算機內存等優勢，對于描述露天礦滑坡等大變形問題是很適用的。在本文中，結合SPH方法與Mohr-Coulomb彈塑性理論，將南芬露天礦邊坡老滑坡松散體離散成一系列粒子，這些粒子攜帶了代表散體顆粒的位置、密度和速度等信息，可以完整地描述老滑坡體在自重、降雨、凍融及振動等作用下的動態滑動過程。該方法可以獲得任意時刻、任意位置粒子的速度、位移、能量耗散等信息，借此可以全面了解露天礦邊坡滑坡的動態全過程。數值結果也可以幫助更好地了解滑坡機理。SPH方法屬于無網格法，可以避免出現網格畸形而導致的計算不收斂問題。除了計算效率高之外，其可靠性亦能得到保障，如圖6所示，SPH計算得到的滑坡形態和范圍非常接近實際滑坡形態。另外，滑動距離和速度等信息同樣可以按要求在數值模擬結果中予以顯示(圖7和圖8)。數值計算結果可以為露天礦邊坡滑坡災害提供一些防治減災策略，并可以根據數值結果中滑坡形態和范圍為滑坡災害劃定范圍，提前制定防治減災措施，比如對邊坡進行加固支護等，以保障露天礦的安全運營。不過，在SPH方法中，露天礦邊坡老滑坡體被視為單一物質，未考慮老滑坡體內的弱面、空隙和水壓力等影響。在以后的研究中，如果考慮流固耦合等問題，可能會使得計算精度進一步提高，但計算效率可能會受到一定程度的影響。因次，在SPH方法中更加真實地描述露天礦邊坡物質組成，以及如何保障高計算效率，將是下一步的研究重點。

4 結 論

1) 本研究以露天礦下盤394～502 m臺階中部老滑坡體區域的滑坡為例，對其機理進行了研究。老堆積體的自身穩定性較低，開挖或爆破振動、降雨和凍融作用是導致滑坡的最主要原因。

2) 采用SPH方法結合彈塑性Mohr-Coulomb模型對露天礦滑坡進行動力學全過程模擬，數值結果在滑動距離、形態和范圍等方面與實際滑坡表現出較好的一致性。SPH方法可以用來描述露天礦滑坡動力學過程。

3) SPH數值結果可以為露天礦滑坡劃定危險范圍，為礦山安全運營提供防治減災的指導性建議。