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

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

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

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

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

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

1 數(shù)值模擬實(shí)驗(yàn)的相關(guān)理論與過程

1.1 MatDEM簡(jiǎn)介

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

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

土顆粒、孔隙和孔隙水構(gòu)成MatDEM中一個(gè)離散元單元，因此每個(gè)單元可以設(shè)置和傳輸一定量水分，并通過把有限差分思想運(yùn)用到數(shù)值模擬中，計(jì)算單元水分遷移量，MatDEM就可以實(shí)現(xiàn)對(duì)水分場(chǎng)的模擬[7-8]。同時(shí)，考慮水分場(chǎng)對(duì)土體抗拉強(qiáng)度等力學(xué)性質(zhì)的影響，可以最終得到水分場(chǎng)和位移場(chǎng)在一定時(shí)間內(nèi)的演化規(guī)律。

模型中每個(gè)單元顆粒的含水量代表顆粒大小區(qū)域內(nèi)的平均含水量。雨水落到邊坡表面再向下傳播，每個(gè)時(shí)間步單元增長(zhǎng)含水率計(jì)算見式(1)。

(1)

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

根據(jù)達(dá)西定律描述土體滲流場(chǎng)內(nèi)水頭的分布規(guī)律，水分由高水頭向低水頭遷移，進(jìn)而實(shí)現(xiàn)降雨入滲的過程，計(jì)算見式(2)。

(2)

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

(3)

式中：i01、i02以及i03分別為白色顆粒與三個(gè)黑色量之間的水力梯度值；ω0、ω1為白色顆粒與1號(hào)單元的含水量。1號(hào)單元在被傳遞dM01后，白色單元會(huì)減少相應(yīng)的含水量，實(shí)現(xiàn)了水分運(yùn)輸守恒，通過以上理論和方法繼而實(shí)現(xiàn)降雨水分場(chǎng)演化。

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

1.3 離散單元的材料

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

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

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

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

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

2) MatDEM分層和賦予材料。MatDEM可以利用導(dǎo)入的坐標(biāo)點(diǎn)文件生成不同劃分地層，只要確定地層邊界坐標(biāo)點(diǎn)就可以對(duì)任何復(fù)雜二維模型進(jìn)行建模，極大地降低了使用者的建模難度。 編寫邊坡模型代碼導(dǎo)入MatDEM后，生成一個(gè)坡頂長(zhǎng)為2 m、坡底長(zhǎng)為4.73 m的尾礦堆積體，然后對(duì)完成分層的模型進(jìn)行裂隙設(shè)置。由于裂隙生成比較復(fù)雜，所以模型中使用軟弱層來替代裂隙的影響效果。又因?yàn)閱卧g孔隙的存在，邊坡裂隙含水率高于其他邊坡內(nèi)部區(qū)域含水率。具體設(shè)置為邊坡裂隙含水率為常數(shù)0.3，非裂隙區(qū)域含水率為常數(shù)0.1。

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

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

3 數(shù)值實(shí)驗(yàn)結(jié)果

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

3.1 水分場(chǎng)演變規(guī)律

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

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

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

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

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

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

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

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

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

3.2 位移演變規(guī)律

3.2.1 觀測(cè)區(qū)滑落位移分析

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

3.2.2 移動(dòng)滑體量分析

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

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

圖10 移動(dòng)滑體量Fig.10 Move slider amount

3.3 裂隙的演變規(guī)律

3.3.1 單元連接演變規(guī)律

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

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

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

3.3.2 配位數(shù)演變規(guī)律

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

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

圖13 配位數(shù)變化趨勢(shì)圖Fig.13 Variation trend graph of coordination number

4 結(jié) 論

1) 降雨入滲到裂隙區(qū)域并導(dǎo)致裂隙的單元連接斷開，邊坡穩(wěn)定性發(fā)生改變，裂隙處的配位數(shù)即密實(shí)程度下降，在持續(xù)降雨作用下，最終導(dǎo)致了滑坡的發(fā)生。

2) 0 s之后，受降雨影響的滑體移動(dòng)總量并沒發(fā)生太大改變，這說明降雨對(duì)滑體移動(dòng)的量并沒有太大影響。

3) 將模擬得到的水分遷移以及位移場(chǎng)變化規(guī)律與室內(nèi)試驗(yàn)及前人研究結(jié)論相對(duì)比，可以比較清楚地觀察降雨誘發(fā)滑坡的全過程，能夠較為有效地進(jìn)行降雨誘發(fā)滑坡模擬，為利用離散元法模擬降雨誘發(fā)滑坡過程提供了一個(gè)嶄新的思路。

4) 贛州市尾礦邊坡的預(yù)防和管理工作主要是以防止開采、標(biāo)準(zhǔn)化和風(fēng)險(xiǎn)發(fā)生時(shí)的早期警告為目的。特別是在贛州市的梅雨季節(jié)，降雨會(huì)持續(xù)相當(dāng)一段時(shí)間。了解邊坡不穩(wěn)定性機(jī)制，特別是對(duì)于高陡邊坡的數(shù)值模擬研究分析，有助于治理贛州市尾礦邊坡的治理，減小滑坡發(fā)生的可能性。