基于MatDEM對高陡邊坡裂隙巖體失穩演化機制的研究

2023-02-18 13:14:04蔣佩伶王志松

中國礦業 2023年2期

蔣佩伶，王志松，蔣越，劉浩

(1.江西理工大學電氣工程與自動化學院，江西贛州 341000；2.江西理工大學應急管理與安全工程學院，江西贛州 341000；3.江西理工大學信息工程學院，江西贛州 341000)

通過對我國離子型重稀土儲量進行統計分析可知，贛州市重稀土儲量在國內同類型礦山儲量排行中位居第一。為了推動贛州地區經濟持續發展，大量企業對贛州市稀土礦產資源進行了大規模開采，直接導致礦區的地質環境越來越差。選用池浸、堆浸和原地浸礦的工藝開采方式使礦區土壤越來越松散，在降雨等惡劣天氣影響下，極其容易滑坡，使當地生產發展和居民生活的安全難以得到保障。近幾年，國家出臺相關法規對耕地進行保護，許多礦區為減少占地面積，不得不在原來的尾礦上進行二次堆積，長時間開采所形成的尾礦廢渣堆積也越來越高，據統計，贛州市超過200 m的尾礦達30多處[1]。

贛州市的尾礦主要位于山區，地勢復雜，礦區分布分散，導致監測成本高、難度大，不能集中統一監測。由于稀土利潤巨大導致難以監管，許多私企在偏僻地段的開采屢禁不止。贛州市稀土開采的規模持續擴大，直接導致礦山裸露地表區的土壤侵蝕越來越嚴重。與一般土質邊坡的性質相比，尾礦的廢渣堆積存在結構無序、弱膠結或無膠結的性質，規律和破壞模式差異較大[2-3]，給離子型稀土礦環境綜合治理工作埋下了重大的安全隱患。贛南地區雨水豐沛，特別是在每年的梅雨季節，受持續強降雨等極端天氣影響，雨水會深層次滲入邊坡，極易產生局部垮塌甚至大規模的滑坡災害[4-6]。如圖1所示為贛州市某尾礦滑坡實拍，如果類似的尾礦邊坡處于上游地區，滑坡災害所造成的損失還將成倍增加。

圖1 贛州市某尾礦堆積實拍Fig.1 A real shot of a tailings accumulation in Ganzhou City

邊坡的滑落破壞除了降雨以外還由很多其他因素導致，如邊坡的坡角大小、干灘比、活動荷載、內摩擦角、裂隙大小和深度等，對于這些因素前人進行了很多研究，也進行過相關數值模擬。但對于在降雨作用下單獨分析裂隙這一因素的數值模擬還比較少，加之國產數值模擬軟件近幾年處于推廣過程中，也缺少相關的研究分析。因此，本文以60°邊坡作為高陡邊坡中的一個特例，對其在降雨條件下的邊坡含水率和位移的演化規律進行相關研究。

1 數值模擬實驗的相關理論與過程

1.1 MatDEM簡介

離散元法(DEM)基于分子動力學納米技術，最早主要應用于物理學和流體動力學。該方法自從Cundall引入研究顆粒集合的力學行為以來開始受到關注。經過幾十年的發展證明，離散元已經成為一種廣泛應用的數值計算方法，特別是在巖土力學研究中，其整體力學性質的數值模擬可以與物理實驗進行對比分析。本文采用國產自主研發的離散元軟件(MatDEM)進行數值模擬。MatDEM基于矩陣離散元法，在單元數量和運行速度方面與同類型軟件相比有明顯的優勢，極大地方便了數值模擬。

1.2 降雨入滲與水分傳輸過程

土顆粒、孔隙和孔隙水構成MatDEM中一個離散元單元，因此每個單元可以設置和傳輸一定量水分，并通過把有限差分思想運用到數值模擬中，計算單元水分遷移量，MatDEM就可以實現對水分場的模擬[7-8]。同時，考慮水分場對土體抗拉強度等力學性質的影響，可以最終得到水分場和位移場在一定時間內的演化規律。

模型中每個單元顆粒的含水量代表顆粒大小區域內的平均含水量。雨水落到邊坡表面再向下傳播，每個時間步單元增長含水率計算見式(1)。

(1)

式中：ω為顆粒每次計算時的含水率；mw為本次迭代增加的含水量；mw0為本次迭代計算時的初始含水量；ms為固體顆粒的重量。

根據達西定律描述土體滲流場內水頭的分布規律，水分由高水頭向低水頭遷移，進而實現降雨入滲的過程，計算見式(2)。

(2)

式中：Q為滲流的流量；H1-H2為水頭損失；A為流過顆粒的斷面積；L為滲流通過的長度；k為滲透系數，表示孔隙在透水方面的物理性質。只有當總水頭差Δh>0時，水才會發生從總水頭高的點向低的點流動。降雨入滲過程是非飽和土體變為飽和的過程。因此，本文將顆粒間總水頭簡化為位置水頭。以圖2中白色基本單元為例，其會向周圍相對其水頭較低單元傳遞水分。根據達西定律，當白色顆粒會傳遞水分給下方三個黑色顆粒時，若兩者間水力梯度值大則水分傳遞量越大。因此假定在一個計算步內，相鄰兩單元水分傳遞量與兩者含水量差值以及水力梯度值有關，白色單元向傳遞給1號單元的含水量計算見式(3)。

(3)

式中：i01、i02以及i03分別為白色顆粒與三個黑色量之間的水力梯度值；ω0、ω1為白色顆粒與1號單元的含水量。1號單元在被傳遞dM01后，白色單元會減少相應的含水量，實現了水分運輸守恒，通過以上理論和方法繼而實現降雨水分場演化。

圖2 水分傳遞示意圖Fig.2 Schematic diagram of moisture transfer

1.3 離散單元的材料

在離散元建模中，建模重點和難點之一就是確定模型中每個單元力學參數。因此，MatDEM中存在一個轉換公式使堆積模型的宏觀力學性質和微觀力學參數相互轉換，使離散元堆積材料具有特定彈性模量和強度性質，并用于邊坡巖土不同地層之中，通過使用訓練后的材料就可以實現相同類型模型的建模自動化[9]。如果模擬模型是線彈性模型，微觀力學參數就可以通過五個宏觀力學參數性質計算得到，分別是楊氏模量、泊松比、抗壓強度、抗拉強度和內摩擦系數。該模型通過不斷的材料數值測試和自動調整，得到所需要的五個宏觀力學參數。為了減少模型的力學性質與設定值之間的差異，函數還會自動將設定值乘以相關系數再代入相關公式，本文對尾礦的不同材料設置詳見表1。

2 基于滑坡的離散單元模型建立

1) 隨機堆積模型。廢棄稀土尾礦堆積的邊坡土壤主要表現為一種各向異性、結構無序、弱膠結或無膠結的特殊地質體。因此，模型通過在堆積模型時設置分散系數來實現這一現象[10]，在模型箱內生成平均半徑為0.03 m的隨機單元并壓實這些單元，來替代實際中巖土顆粒的沉積過程。如圖3所示，建立一個高為8.5 m、寬為12 m的長方形初始堆積模型。

表1 尾礦的材料Table 1 Material of tailings

圖3 初始堆積模型Fig.3 Initial stacked model plot

2) MatDEM分層和賦予材料。MatDEM可以利用導入的坐標點文件生成不同劃分地層，只要確定地層邊界坐標點就可以對任何復雜二維模型進行建模，極大地降低了使用者的建模難度。編寫邊坡模型代碼導入MatDEM后，生成一個坡頂長為2 m、坡底長為4.73 m的尾礦堆積體，然后對完成分層的模型進行裂隙設置。由于裂隙生成比較復雜，所以模型中使用軟弱層來替代裂隙的影響效果。又因為單元之間孔隙的存在，邊坡裂隙含水率高于其他邊坡內部區域含水率。具體設置為邊坡裂隙含水率為常數0.3，非裂隙區域含水率為常數0.1。

3) 設置降雨。降雨為發生滑坡的重要原因之一[11-13]。在MatDEM代碼中，選取邊坡最外層為雨水層，并設置其含水率不變且為最高，具體設置雨水層的含水率為常數0.9，只要雨水層的含水率不發生變化，就可以表示降雨在持續進行。

4) 建立觀測區。為了便于分析，建立四個觀測區，每個觀測區中有300多個單元和兩組數據(圖4)。觀測區域1取z值為5.75～6.25 m，且x值為1.75～2.25 m的所有單元顆粒，觀測區域2取z值為3.75～4.25 m，且x值為1.75～2.25 m的所有單元顆粒，觀測區域3取z值為3.75～4.25 m，且x值為3.75～4.25 m的所有單元顆粒，觀測區域4取z值為1.75～2.25 m，且x值為3.75～4.25 m的所有單元顆粒。觀測區1和觀測區2的高差與觀測區3和觀測區4的高差相同，均為2 m。此外，觀測區2和觀測區3的水平距離也為2 m，這一設置除了方便進行數據定量分析外，還可以讓一個觀測區內的單元同時擁有裂隙區單元和非裂隙區單元。

3 數值實驗結果

本文完成了強降雨作用下邊坡從失穩、單元顆粒連接斷裂到整體滑移的全過程模擬，該過程共歷時約36 s，計算機模擬耗時約5 h，共生成12個MatDEM文件和8個供分析的數據表，主要數據為含水率、位移變化、配位數以及單元速度。通過對數據提取、分類、篩選、分析，得出位移場、水分場和裂隙演變規律。

3.1 水分場演變規律

為了對持續降雨環境進行模擬，設定最頂層顆粒單元的含水率為常數0.9。模擬結果顯示，在強降雨條件下，水分持續入滲，并導致邊坡滑移。圖5為坡體內的含水率時程分布圖。由圖5可知，在3 s、15 s、24 s、27 s時刻，分別入滲到坡體內約為1 m、2.5 m、3 m、4 m位置。坡體內含裂隙區域含水率的變化速率明顯較高，在這幾個時刻，分別入滲到距離破頂面約為1 m、3.5 m、4.5 m、6 m位置，其滲流速率約為非裂隙區域2倍左右。

圖4 邊坡模型示意圖Fig.4 Schematic diagram of slope model

圖5 含水率的時程分布圖Fig.5 Time-history distribution of moisture content

圖6展示了單元含水率變化情況。由圖6可知，A觀測點、B觀測點、C觀測點、D觀測點依次為觀測區1、觀測區2、觀測區3和觀測區4的非裂隙區單元觀測點；E觀測點、F觀測點、G觀測點、H觀測點依次為觀測區1、觀測區2、觀測區3和觀測區4的裂隙區單元觀測點。隨著時間的推移，四個觀測區所有單元的含水率均會逐漸增加。在圖6(a)中，非裂隙單元含水率從0 s時就開始變化。除了A觀測點與D觀測點的含水率變化曲線在6.1 s出現增速變化外，四個觀測區中的單元含水率曲線并未發生較大的增速變化。6.1 s之后，四個觀測點的含水率曲線并未發生較大的波動，從高到低依次為ωC、ωA、ωD、ωB。在圖6(b)中，裂隙區內單元含水率變化出現在6.1 s之后。在10.7 s之前，單元含水率變化由高到低排列分別為ωG、ωE、ωH、ωF，其中，ωH和ωF含水率相等。10.7 s之后，裂隙區單元含水率由高到低排列分別為ωE、ωG、ωF、ωH。其中，F觀測點的含水率變化從15 s開始逐漸增加，H觀測點的含水率變化在21.1 s之后開始增加。

為了研究三個觀測區的單元含水率變化，分別對觀測區中的八個觀測點的單元含水率差值進行分析。如圖7(a)所示，在非裂隙單元中，ωC-ωD差值曲線始終高于ωA-ωB差值曲線，且分別在15 s時達到含水率差值峰值，出現這類現象的原因主要是C觀測點處在坡底靠近邊坡面，而A觀測點與B觀測點的位置靠近裂隙，水分傳輸較快，含水率差值反而較小。與圖5所示15 s時含水率變化曲線進行對比分析可以發現，在15 s時，水分一直沿著裂隙帶遷移，水分的傳導層達到觀測區1附近，濕潤鋒達到觀測區2附近，而在觀測區4沒有水分尚未入滲。在圖7(b)中，又對裂隙單元的含水率差值進行了對比，可以發現，在11 s前，ωG-ωH差值曲線高于ωE-ωF差值曲線。出現這一現象的主要原因是G觀測點的位置要比E觀測點更靠近坡面，水分遷移更快。在10～18 s之間，兩條差值曲線高低發生改變，水分開始沿著裂隙帶向下加速遷至E觀測點，所以ωE-ωF差值曲線高于ωG-ωH差值曲線，ωE-ωF差值曲線在15 s出現最高點，最高點的含水率差值為常數0.375。在18 s時，出現兩條曲線的第二個交點。18 s之后，ωG-ωH差值曲線高于ωE-ωF差值曲線，并且ωG-ωH差值曲線在26 s達到最高點，其含水率差值約為常數0.31。出現第二個交點的原因主要是：在18 s時，裂隙帶上層的單元滑落，G觀測點與邊坡表面的水分傳輸距離縮短。整體來看，隨著時間增加，四個觀測區單元含水率都逐漸趨于飽和，增速也明顯變緩，在分別達到含水率差值峰值之后，差值曲線都出現了下降的趨勢，意味著四個觀測區單元含水率將逐漸趨于飽和。

圖6 單元含水率變化曲線Fig.6 Change curves of unit water content

圖7 含水率差值變化曲線Fig.7 Change curves of water content difference

單獨取觀測區1、觀測區2和觀測區3中的六個觀測點，進行單元水分橫向遷移與縱向遷移的快慢分析，并探究裂隙是否對水分遷移有一定的作用。因為前面已經設置三個觀測區的直線距離為2 m，可以對三個觀測區進行相同水分運距和時間、不同方向的含水率差值對比分析，如圖8所示。由圖8可知，15 s附近三個觀測區的含水率差值同時達到峰值。并且從圖8(a)中還可以看出，ωA-ωB差值曲線始終高于ωC-ωB差值曲線，說明縱向水力梯度要大于橫向水力梯度。相比于非裂隙區單元而言，在圖8(b)中，裂隙區單元含水率差值卻在10.1 s出現一個交點。在10.1 s之前，裂隙區單元與非裂隙區單元含水率差值變化一樣，但在10.1 s之后，裂隙對雨水遷移開始影響含水率的差值變化，ωG-ωF差值曲線高于ωE-ωH差值曲線，不同方向的水力梯度大小變為橫向水力梯度大于的縱向水力梯度。這一變化與非裂隙區有明顯差別，由此可以得出裂隙對水力梯度的大小有一定的影響。

圖8 不同方向含水率差值變化曲線Fig.8 Change curves of water content difference in different directions

3.2 位移演變規律

3.2.1 觀測區滑落位移分析

在高降雨強度條件下，邊坡的滑落破壞極其容易產生[14-16]。與整體破壞滑坡有所不同，因降雨入滲裂隙導致孔隙水壓增加的非整體破壞滑坡，往往是從裂隙引發的局部滑坡，并在雨水和自重作用下滑落至邊坡底部[17-18]。由圖9可知，當降雨入滲9 s后，分別在距離坡底1 m、2.5 m、6.5 m的位置出現了小規模的下沉破壞。當降雨入滲18 s后，邊坡沿著裂隙出現分層滑坡，水分使裂隙區域單元連接斷開，連接斷開的單元顆粒在重力和孔隙水壓的雙重作用下產生長約5 m的斷裂帶。隨著降雨繼續進行至27 s，破壞范圍逐漸擴大，滑體部分的顆粒含水率基本達到飽和，滑體最底端由2 m的高度沖擊到地面，滑體最上端的高度由7 m降至為4 m。當降雨入滲36 s后，滑落部分完全下降至底層地面，也代表著尾礦邊坡滑落完全結束。

3.2.2 移動滑體量分析

降雨入滲使排土場體內含水率飽和區域迅速增加，并不斷下移擴散，排土體自重增大[19]。在MatDEM中，篩選X軸方向移動速度大于0的單元，每0.1 s計算其面積之和作為移動滑體量，如圖10所示。由圖10可知，在整個邊坡失穩過程中，受降雨影響的移動滑體量一直圍繞在0 s后均線12.7 m2附近上下變化，并未發生太大的變化。在0 s時，邊坡未受到雨水侵蝕，單元也沒有任何移動，在1～6 s和12～22.5 s之間，受單元自身摩擦力的影響，產生兩個移動滑體量波動范圍較大的區域，特別是在16.5 s附近，移動滑體量出現一個最大振幅，由8 m2增加到17.5 m2。在21 s之后，移動滑體量的波動范圍趨于穩定，單元移動滑體量仍在均線12.7 m2附近。

圖9 位移變化圖Fig.9 Diagram of displacement change

圖10 移動滑體量Fig.10 Move slider amount

3.3 裂隙的演變規律

3.3.1 單元連接演變規律

裂隙作為優勢入滲通道，使得雨水能快速沿著裂隙滲入到坡體內部。同時，裂隙內水體入滲屬于有壓入滲[20]，與邊坡表面的入滲相比，裂隙入滲的速度更快，對邊坡的破壞作用也較大[21-23]。在MatDEM中，單元之間連接如圖11所示。其中，邊坡內部的非空白區域表示單元連接，增加的空白區域表示單元連接斷開。隨著水分開始入滲邊坡內部，裂隙單元連接逐漸斷開，并開始出現滑坡。與邊坡移動過程圖不同，降雨對裂隙單元間連接的破壞并不是36 s后才結束。由圖11(d)可知，單元之間連接在21 s前就完全斷開，并且單元連接沿著裂隙向上和向下都有斷開。以最內側裂隙帶為例進行分析：在6 s時，上方空白區域長約為1.5 m，下方空白區域長也約為1.5 m。但在18 s時，上方空白區域長約為5 m，下方空白區域長約為2 m，說明單元連接向下斷開的速度比向上斷開的速度更快。

圖11 單元連接變化Fig.11 Changes of unit connection

圖12 邊坡配位數過程圖Fig.12 Process map of slope coordination number

3.3.2 配位數演變規律

在MatDEM中，一個單元顆粒周圍連接了許多顆粒，而配位數就是某一單元顆粒周圍顆粒的數量之和，它反映單元密實程度[24]，如圖12所示為單元配位數變化過程。配位數主要受單元連接斷開影響而發生變化，如果單元連接斷開，則單元配位數會下降。在本文模型中，只要裂隙單元含水率超過常數0.5，則裂隙單元連接就會斷開，單元配位數也同時會減小。因此，配位數的變化還間接反映降雨對單元密實程度的影響。除此之外，非裂隙單元在受到外界力量或單元自身重力時，也會造成配位數發生改變。由圖12可知，在18 s之前，配位數主要受降雨的影響在裂隙處變化。通過配位數的變化，裂隙單元之間的連接力下降，裂隙和邊坡堆積巖層相互之間的作用力也同時減小，使邊坡更容易滑落。在18 s之后，配位數的變化主要受自身重力的沖擊影響而增加，滑體向下沖擊，使部分滑體前端部分的配位數下降至2個以下。

由于無法預先確定滑體的位置和滑體量，因此無法求得在滑落過程中的配位數，但是可以通過篩選配位數小于2個的單元，求得單元面積和來表示配位數的改變量，從而反映配位數量的變化趨勢。如圖13所示，0～36 s之間的配位數面積和一直在增加，反映出邊坡密實程度一直減小。0～9 s之間，受降雨影響，小于2個的配位數量的單元面積和增長比較快，說明單元的密實程度變化較快。在滑落中間階段，即21～27 s，受滑體沖擊地面的影響，配位數小于2個的單元面積和出現第二次快速增加，證實了在沖擊力的影響下，單元的密實程度也會進一步減小。在27 s之后，滑體大部分落至地面，滑落基本結束，配位數量變化趨于穩定，說明單元密實程度最終趨于穩定。