朱永東，王雪冬，孫延峰,王 翠,劉光偉

(1.遼寧工程技術(shù)大學(xué) 礦業(yè)學(xué)院，遼寧 阜新 123000；2.中冶沈勘工程技術(shù)有限公司，遼寧 沈陽 110169)

內(nèi)排土場(chǎng)是人工堆積的松散土石體，在降雨條件下容易引發(fā)滑坡，嚴(yán)重威脅人類的生命財(cái)產(chǎn)安全，因此，降雨型滑坡啟動(dòng)與運(yùn)動(dòng)特征研究對(duì)內(nèi)排土場(chǎng)防災(zāi)減災(zāi)工作意義重大。

降雨型滑坡方式和預(yù)測(cè)預(yù)報(bào)的研究成果較多，許強(qiáng)、曾裕平等[1-2]結(jié)合已發(fā)生的滑坡，證實(shí)了斜坡變形一般要經(jīng)歷初始、等速和加速變形階段，提出了預(yù)測(cè)預(yù)報(bào)方法。沈佳等[3]通過試驗(yàn)和模擬發(fā)現(xiàn)強(qiáng)降雨滑坡具有突發(fā)性和破壞迅速的特征。李聰?shù)萚4]采用統(tǒng)計(jì)分析法得出，當(dāng)滑面為硬性結(jié)構(gòu)時(shí)，滑坡加速與急劇變形階段歷時(shí)會(huì)更短。COLLINS等[5]提出了一種程序化的降雨入滲誘發(fā)滑坡預(yù)測(cè)公式。馮杭建等[6]利用GIS技術(shù)與確定性系數(shù)方法，分析了降雨型滑坡影響因子的敏感性。徐毅、吳火珍、SENTHILKUMAR等[7-9]以非飽和土力學(xué)理論為基礎(chǔ)，研究發(fā)現(xiàn)降雨作用下，基質(zhì)吸力降低導(dǎo)致土體抗剪強(qiáng)度降低，進(jìn)而誘發(fā)滑坡，為降雨型滑坡機(jī)理研究提供了思路。

以自然傾倒方式為主的內(nèi)排土場(chǎng)，坡度大且土料松散，在降雨入滲條件下坡體表面常發(fā)生溯源侵蝕破壞，多呈平行于坡表的淺層破壞特征[10]。張群等[11]認(rèn)為在降雨條件下，界面效應(yīng)是導(dǎo)致沿基覆界面淺層滑坡的原因。陶志剛等[12]發(fā)現(xiàn)內(nèi)排土場(chǎng)在降雨作用下泥質(zhì)膠結(jié)材料沖刷流失，排棄物間隙變大形成潛在滑動(dòng)面，進(jìn)而產(chǎn)生滑坡。WEI等[13]結(jié)合野外監(jiān)測(cè)和數(shù)值方法，建立了淺層滑坡觸發(fā)的水文條件模型。黃剛海等[14]借助ANSYS與FLAC3D給出了內(nèi)排土場(chǎng)產(chǎn)生滑坡危險(xiǎn)所需要的降雨量指標(biāo)。目前，通過試驗(yàn)、數(shù)值方法等綜合方法研究淺層滑坡過程較多，從動(dòng)力學(xué)角度探究滑動(dòng)過程的研究成果相對(duì)較少[15]。

MANENTI等[16]用流體動(dòng)力學(xué)模擬降雨型淺層滑坡來評(píng)估滑坡前緣速度與對(duì)下游墻體的沖擊力。MARTELLONI等[17]用離散元(DEM)的分子動(dòng)力學(xué)(MD)方法模擬并識(shí)別出淺層滑坡的動(dòng)力學(xué)特性，從而對(duì)邊坡失穩(wěn)時(shí)間進(jìn)行預(yù)測(cè)。葛云峰等[18]借助PIV技術(shù)獲取高速遠(yuǎn)程滑坡中任意時(shí)刻碎屑顆粒流的運(yùn)動(dòng)信息，揭示了滑體中后部速度波動(dòng)及碰撞次數(shù)能量均低于前部顆粒運(yùn)動(dòng)演化規(guī)律。目前的研究較少涉及到降雨入滲過程和內(nèi)排土場(chǎng)淺層滑坡的全過程，以及從滑體的運(yùn)動(dòng)特征角度來探究淺層滑坡的破壞特征。

因此，筆者以內(nèi)蒙古元寶山露天煤礦內(nèi)排土場(chǎng)為原型，依據(jù)相似理論進(jìn)行物理模擬試驗(yàn)，獲得內(nèi)排土場(chǎng)降雨入滲與淺層滑坡破壞規(guī)律，在此基礎(chǔ)之上，推導(dǎo)出內(nèi)排土場(chǎng)淺層滑坡穩(wěn)定性分析方程及淺層滑坡的啟動(dòng)判據(jù)，并且借助離散元方法(MatDEM)的快速GPU矩陣計(jì)算模擬出降雨入滲過程和淺層滑坡全過程，掌握滑坡過程中的位移，速度和能量轉(zhuǎn)換規(guī)律，最終獲得降雨條件下內(nèi)排土場(chǎng)的淺層破壞與運(yùn)動(dòng)特征，成果可為內(nèi)排土場(chǎng)淺層滑坡機(jī)理研究提供參考。

1 內(nèi)排土場(chǎng)降雨入滲特征

1.1 研究區(qū)概況

元寶山露天煤礦位于內(nèi)蒙古赤峰市，季風(fēng)性氣候，夏季雨量集中，一般為230～270 mm，占全年降水量總量的68%，并常以大雨或暴雨的形式出現(xiàn)，容易引發(fā)內(nèi)排土場(chǎng)滑坡和泥石流災(zāi)害。

圖1為研究區(qū)照片，礦區(qū)地層以泥質(zhì)粉砂巖、礫巖、砂礫巖和泥巖為主，地層傾角3°～14°，內(nèi)排土場(chǎng)土料主要為煤礦剝離的黃褐色細(xì)粒土質(zhì)砂，主要由亞黏土、洪積圓礫以及泥礫組成，作為廢棄松散土料堆積在現(xiàn)有地質(zhì)體上，內(nèi)排土場(chǎng)土料性質(zhì)相近。

圖1 研究區(qū)實(shí)物(鏡頭西南朝向)Fig.1 Physical map of the study area(The lens faces southwest)

根據(jù)相似理論，以研究區(qū)南幫內(nèi)排土場(chǎng)的典型臺(tái)階(紅圈部分)為研究對(duì)象，坡角34°，表土覆蓋厚12 m，按1∶60幾何比例進(jìn)行縮尺，建立降雨型內(nèi)排土場(chǎng)滑坡物理模擬試驗(yàn)?zāi)Ｐ?，如圖2,3所示。試驗(yàn)裝置由透明鋼化有機(jī)玻璃拼接而成，其尺寸為100 cm×24 cm×50 cm，主要包括：模型槽，降雨系統(tǒng)，測(cè)量系統(tǒng)和數(shù)據(jù)采集系統(tǒng)等。

圖2 實(shí)驗(yàn)裝置Fig.2 Diagram of experimental setup

圖3 模型試驗(yàn)Fig.3 Model test

本次試驗(yàn)是結(jié)合研究區(qū)域?qū)嶋H降雨特征，取0.04 mm/min的降雨強(qiáng)度，持續(xù)降雨30 min，觀察降雨過程中雨水的入滲和內(nèi)排土場(chǎng)失穩(wěn)破壞特征。

推導(dǎo)出模型內(nèi)土料主要參數(shù)的相似系數(shù)，如式(1)～(3)所示。試驗(yàn)土料取自內(nèi)排土場(chǎng)現(xiàn)場(chǎng)，依據(jù)土體顆粒級(jí)配相似原則[19]，通過篩分法獲得初始模型土料，進(jìn)行物理力學(xué)參數(shù)測(cè)試，根據(jù)測(cè)試結(jié)果和相似系數(shù)進(jìn)行反復(fù)調(diào)整土料級(jí)配，最終獲得土料的物理力學(xué)參數(shù)對(duì)比結(jié)果見表1。

表1 土料的物理力學(xué)參數(shù)Table 1 Physical and mechanical parameters of waste dump

Cρ=Cφ=Cγ=Cθ=Cg=1

(1)

Ch=1

(2)

(3)

式中，Cρ,Cφ,Cγ,Cθ,Cg分別為土的密度、內(nèi)摩擦角、容重、含水率和重力加速度的相似比；Ch為幾何尺寸相似比；Ck為滲透系數(shù)相似比；Cq為降雨強(qiáng)度相似比。

1.2 試驗(yàn)結(jié)果

從連續(xù)記錄的試驗(yàn)過程中可以看出，降雨過程中內(nèi)排土場(chǎng)整體穩(wěn)定性較好。降雨初期，在雨滴沖擊作用下，坡面形成大小不一的濺蝕凹槽，加上雨水逐漸濕潤(rùn)坡面土并填補(bǔ)坡表裂隙，使得降在斜坡表面上的水短時(shí)間不能入滲到土體內(nèi)部，而是在坡頂積水成灘，在斜坡面匯集于濺蝕凹槽內(nèi)，如圖4(a)所示。當(dāng)濺蝕凹槽內(nèi)的水達(dá)到一定量時(shí)，溢出凹槽并連通幾個(gè)凹槽開始形成小徑流，或沒有形成有效的徑流而是呈現(xiàn)斜坡面流。隨著降雨的持續(xù)，浸潤(rùn)峰不斷向下運(yùn)移，導(dǎo)致坡面無黏性土層飽和之后多表現(xiàn)為侵蝕破壞， 坡體則已經(jīng)開始發(fā)生淺層破壞，且破壞規(guī)模逐漸增大，如圖4(b)所示。從圖4(c)淺層滑坡結(jié)束后的坡面侵蝕情況來看坡體已發(fā)生整體失穩(wěn)并破壞。

圖4 浸潤(rùn)峰運(yùn)移時(shí)的坡表面侵蝕情況Fig.4 Slope erosion during the migration of wetting peak

此外，王樂、蘇永華等[20-21]分別研究了不同降雨雨型和間歇性降雨對(duì)邊坡穩(wěn)定性的影響，發(fā)現(xiàn)也存在類似的降雨入滲特征，即呈平行于坡表面的淺層破壞特征。持續(xù)性降雨將會(huì)導(dǎo)致淺層滑坡現(xiàn)象的間隙性重復(fù)出現(xiàn)，坡表面破壞程度也逐漸加大。面流形成的縱向與橫向破壞同時(shí)進(jìn)行，說明隨著降雨的持續(xù)，表土層呈達(dá)到飽和—淺層滑坡—新表土層達(dá)到飽和的模式循環(huán)出現(xiàn)。

為了進(jìn)一步探究?jī)?nèi)排土場(chǎng)淺層滑坡機(jī)理，通過降雨入滲改進(jìn)模型進(jìn)行淺層滑坡的啟動(dòng)判據(jù)分析，并借助數(shù)值模擬來搭建上述模型試驗(yàn)降雨入滲至淺層滑坡啟動(dòng)時(shí)刻的狀態(tài)，以此來重點(diǎn)研究淺層滑坡的運(yùn)動(dòng)特征。

1.3 改進(jìn)的降雨入滲模型

試驗(yàn)測(cè)得降雨條件下內(nèi)排土場(chǎng)坡面處體積含水率曲線，如圖5所示。結(jié)合有機(jī)玻璃觀察雨水的實(shí)際入滲情況，可以得出降雨初期坡表面土壤逐漸從非飽和向飽和過渡，且從坡面土較長(zhǎng)時(shí)間處于濕潤(rùn)狀態(tài)可知，坡面土在降雨初期的雨水附著率要小于雨水入滲率。當(dāng)坡表面出現(xiàn)飽和區(qū)時(shí)，飽和區(qū)與干土區(qū)之間存在一定厚度的過渡區(qū)。

圖5 體積含水率實(shí)測(cè)曲線Fig.5 Measured curve of volumetric water content

由于雨水入滲過程中飽和區(qū)和干土區(qū)之間過渡區(qū)的存在，在充分考慮了過渡區(qū)厚度和累積入滲量的基礎(chǔ)上，對(duì)未考慮過渡區(qū)的經(jīng)典Green-Ampt入滲模型進(jìn)行改進(jìn)，獲得了適用于內(nèi)排土場(chǎng)降雨入滲實(shí)際特征的計(jì)算模型，計(jì)算簡(jiǎn)圖如圖6所示。

注：H為積水深度，cm；L為雨水入滲深度，cm；b為過渡區(qū)厚度，cm；Z為縱向深度，cm；θ為體積含水率；θi為初始體積含水率；θs為飽和體積含水率；t0為傳感器數(shù)據(jù)發(fā)生變化時(shí)刻；t1為傳感器數(shù)據(jù)趨于穩(wěn)定時(shí)刻。圖6 改進(jìn)前后的Green-Ampt模型Fig.6 Green-Ampt model before and after improvement

如圖6(b)所示，雨水累積入滲總量與入滲深度關(guān)系是一個(gè)分段函數(shù)：當(dāng)入滲深度L小于過渡區(qū)厚度b時(shí)，土體未出現(xiàn)飽和區(qū)，根據(jù)實(shí)測(cè)數(shù)據(jù)將累積入滲量I與入滲深度L進(jìn)行擬合可以得到I1段函數(shù)；當(dāng)入滲深度大于過渡區(qū)厚度b時(shí)，土體中出現(xiàn)飽和區(qū)，此時(shí)累積入滲量與入滲深度之間的關(guān)系為I2段函數(shù)。具體計(jì)算公式見式(4)和式(5):

(4)

(5)

式中，θ′為過渡區(qū)內(nèi)累積入滲量；dt為I2段的某微小時(shí)間段；dI為dt時(shí)間段內(nèi)的總累積入滲量；i為單位時(shí)間內(nèi)的累積入滲量；L1為飽和區(qū)厚度，cm；K為過渡區(qū)內(nèi)平均滲透系數(shù)，cm/min；Sf為過渡區(qū)內(nèi)平均水土吸力水頭，cm。

由于由實(shí)測(cè)數(shù)據(jù)可以得到累積入滲量、浸潤(rùn)峰深度和時(shí)間3者之間的關(guān)系(I1段函數(shù))，因此，只對(duì)I2段函數(shù)計(jì)算式進(jìn)行推導(dǎo)，聯(lián)立式(4)中I2段函數(shù)和式(5)可得

(6)

式中，dL為dt時(shí)間段末首時(shí)刻的入滲深度差值。

因此，I2段函數(shù)入滲深度與時(shí)間之間的關(guān)系為

(7)

考慮到地表積水深度H很小時(shí)、或者入滲時(shí)間t較長(zhǎng)而導(dǎo)致L較大時(shí)H可以忽略，根據(jù)入滲深度L與過渡區(qū)厚度b之間的關(guān)系，得到入滲時(shí)間和入滲深度的分段函數(shù)：

(8)

通過試驗(yàn)獲取入滲模型所需參數(shù)見表2。

表2 入滲模型參數(shù)Table 2 Parameters of infiltration model

將參數(shù)代入式(8)，整理得出雨水入滲深度和時(shí)間之間的關(guān)系

(9)

2 考慮入滲特征的淺層穩(wěn)定性分析

如前文所述，淺層滑坡多呈近直線型。因此，在考慮雨水入滲特征的基礎(chǔ)上，推導(dǎo)內(nèi)排土場(chǎng)穩(wěn)定性計(jì)算方程，圖7為計(jì)算簡(jiǎn)圖。

注：w為浸潤(rùn)峰以上某一條形土的質(zhì)量；ua為孔隙氣壓力，暫不考慮氣體影響，即為大氣壓力ua=0；uw為孔隙水壓力；ua-uw=-uw為土體浸潤(rùn)峰處的基質(zhì)吸力Ψm；τm為破壞面的下滑力；σn為破壞面的正應(yīng)力；β為土體坡角圖7 穩(wěn)定性計(jì)算簡(jiǎn)圖Fig.7 Simplified diagram of stability calculation

邊坡穩(wěn)定系數(shù)由浸潤(rùn)峰以上整個(gè)土體的抗滑力和滑動(dòng)力之比所得，浸潤(rùn)峰處的抗滑力采用非飽和強(qiáng)度公式進(jìn)行求解，其上部分土體的滑動(dòng)力是通過計(jì)算上部土體總重沿平行于浸潤(rùn)峰方向的分力求解。根據(jù)非飽和土摩爾-庫侖準(zhǔn)則和極限平衡法得到邊坡穩(wěn)定系數(shù)計(jì)算公式為

(10)

式中，τf為非飽和土抗剪強(qiáng)度；c′和φ′分別為土的有效黏聚力和內(nèi)摩擦角；φb為抗剪強(qiáng)度隨基質(zhì)吸力變化的吸力摩擦角。

式(10)充分考慮了非飽和土的強(qiáng)度特性，計(jì)算的時(shí)候考慮了浸潤(rùn)層厚度內(nèi)不同深度因含水率不同而不同的土的容重，通過對(duì)降雨入滲過程中不同時(shí)刻實(shí)測(cè)的體積含水率變化數(shù)據(jù)進(jìn)行多項(xiàng)式擬合推算出不同深度的土體容重。

γ(z)=[θ(z)+ρd]g?z∈(L-b,L)

(11)

式中，z為浸潤(rùn)層厚度內(nèi)某一深度；γ(z)為z深度的土體容重；θ(z)為z深度的土體體積含水率；g為重力加速度。

浸潤(rùn)峰處的正應(yīng)力以及浸潤(rùn)峰以上土體的滑動(dòng)力計(jì)算時(shí)分2種情況考慮：第1種是當(dāng)土體在降雨初期時(shí)，土體表層未出現(xiàn)飽和層；另一種是土體表層出現(xiàn)飽和層，計(jì)算公式為

第1階段(L≤b)：

(12)

第2階段(L≥b)：

(13)

式中，γt為土的飽和容重。

整理得邊坡在浸潤(rùn)峰處的穩(wěn)定系數(shù)，第1階段(L≤b)：

(14)

第2階段(L≥b)：

(15)

根據(jù)體積含水率的實(shí)測(cè)數(shù)據(jù)對(duì)過渡區(qū)土體重度與入滲深度之間的關(guān)系進(jìn)行擬合，擬合的關(guān)系式為

γ(L)=16.593-0.320L+0.164L2-0.011L3

(16)

擬合系數(shù)R2=0.977 1，說明擬合效果好，測(cè)定土的非飽和抗剪強(qiáng)度參數(shù)見表3。

表3 土的非飽和抗剪強(qiáng)度參數(shù)Table 3 Unsaturated shear strength parameters of soil

將參數(shù)代入式(14)和(15)，計(jì)算得出邊坡穩(wěn)定系數(shù)隨雨水入滲深度之間的關(guān)系，如圖8所示。

圖8 邊坡穩(wěn)定系數(shù)與雨水入滲深度的關(guān)系Fig.8 Relationship between the slope stability factor and the depth of rainwater infiltration

根據(jù)圖8的計(jì)算結(jié)果可知，排土臺(tái)階在雨水入滲過程中，隨著入滲深度的增加，斜坡表面的穩(wěn)定性逐漸降低。當(dāng)浸潤(rùn)峰為2 cm時(shí)，土體穩(wěn)定系數(shù)為1.37，邊坡表面只出現(xiàn)濺蝕凹槽的擴(kuò)展破壞而未發(fā)生淺層滑坡的現(xiàn)象，如圖4(a)所示，說明失穩(wěn)破壞土體并不全為浸潤(rùn)峰以上土體；當(dāng)浸潤(rùn)峰為10 cm時(shí)，邊坡穩(wěn)定系數(shù)為0.8，坡體淺層滑動(dòng)已發(fā)生，其中坡表飽和層的無黏性土多表現(xiàn)為侵蝕破壞，如圖4(b)所示。當(dāng)浸潤(rùn)峰在2 cm附近時(shí)，浸潤(rùn)峰以上邊坡的穩(wěn)定系數(shù)下降速率較快，說明此范圍雨水的入滲深度對(duì)邊坡穩(wěn)定性影響較大；當(dāng)浸潤(rùn)峰在6 cm附近時(shí)，穩(wěn)定系數(shù)下降速率較慢，說明此范圍雨水的入滲深度對(duì)邊坡穩(wěn)定性影響較小。為了進(jìn)一步探究邊坡穩(wěn)定系數(shù)對(duì)雨水入滲深度的敏感性程度，以圖8中曲線斜率為縱坐標(biāo)，雨水入滲深度為橫坐標(biāo)繪制出圖9，結(jié)果顯示邊坡穩(wěn)定系數(shù)對(duì)雨水入滲深度的變化可劃分為高敏感、敏感削弱和低敏感3個(gè)階段。

圖9 邊坡穩(wěn)定系數(shù)變化速率隨雨水入滲深度之間的關(guān)系Fig.9 Relationship between the slope stability factor and the time of rainwater infiltration

由式(14)可知，當(dāng)浸潤(rùn)峰為4.8 cm時(shí)，邊坡穩(wěn)定系數(shù)為1，邊坡即將發(fā)生淺層滑動(dòng)。對(duì)于黏性土邊坡，邊坡滑坡將呈現(xiàn)整體滑動(dòng)，既滑動(dòng)的各層黏性土體趨于同步；對(duì)于內(nèi)排土場(chǎng)無黏性土邊坡，浸潤(rùn)峰以上土體也同樣屬于整體失穩(wěn)。只是無黏性土趨于飽和時(shí)，土粒間的束縛會(huì)進(jìn)一步減弱，則越趨于飽和的土體越易呈顆粒流滑動(dòng)，從而表現(xiàn)為侵蝕現(xiàn)象[22]。此外，越近飽和區(qū)土體越先啟動(dòng)且越易滑動(dòng)，隨著趨于飽和土體的先行滑動(dòng)，暴露出的下方土體也逐步趨于新的飽和，在物理模擬試驗(yàn)中觀測(cè)出的無黏性土邊坡循環(huán)破壞模式也證實(shí)了這一點(diǎn)。因此，本節(jié)提出的淺層穩(wěn)定性分析方法具有適用性。

為了進(jìn)一步研究?jī)?nèi)排土場(chǎng)淺層滑坡的運(yùn)動(dòng)特征，即滑坡內(nèi)部的速度場(chǎng)、能量轉(zhuǎn)移等特征，利用MatDEM軟件對(duì)降雨入滲過程和淺層滑坡運(yùn)動(dòng)過程進(jìn)行模擬。

3 降雨入滲與淺層滑坡數(shù)值模擬

MatDEM[23]軟件是用于模擬顆粒構(gòu)件的離散元方法軟件。該軟件通過一定粒徑范圍的顆粒緊密隨機(jī)堆載且壓實(shí)來為之后的模型切割做準(zhǔn)備，再對(duì)指定范圍的顆粒群體附加各項(xiàng)參數(shù)，使其在微觀上以法向彈簧和切向彈簧的線彈性彈簧形式達(dá)到模擬特定巖土體的目的[24]。此外，還可以結(jié)合理論或經(jīng)驗(yàn)公式來模擬出物理實(shí)驗(yàn)中雨水入滲并軟化土壤的過程。然后，施加重力進(jìn)行力平衡迭代，將原先設(shè)置的任意不平衡初始條件迭代到力平衡的穩(wěn)定狀態(tài)。迭代平衡過程是空氣阻尼與摩擦碰撞耗散動(dòng)能的過程。數(shù)值模擬的詳細(xì)方法和技術(shù)方面信息可參考其他論文[25-26]。

3.1 模型中宏、微觀參數(shù)

土料的實(shí)際宏觀參數(shù)包括彈性模量(E)、泊松比(ν)、天然密度(ρ)、抗壓強(qiáng)度(Cu)，拉伸強(qiáng)度(Tu)和內(nèi)摩擦因數(shù)(μi)，其中密度和內(nèi)摩擦因數(shù)可直接獲得，彈性模量和泊松比可根據(jù)經(jīng)驗(yàn)獲取[27]，因排土場(chǎng)土料為散體材料，拉伸強(qiáng)度近似為0，抗壓強(qiáng)度(Cu)可通過式(7)[28]獲得

(17)

最終得到實(shí)際宏觀參數(shù)見表4。

表4 土料的實(shí)際宏觀參數(shù)Table 4 Actual macroscopic parameters afterconversion of soil material

在使用MatDEM軟件時(shí)，模型中的土料需要對(duì)應(yīng)的宏、微觀參數(shù)，其中微觀參數(shù)，包括法向剛度(Kn)、剪切剛度(Kτ)、抗剪強(qiáng)度(Fso)、斷裂位移(Xb)和摩擦因數(shù)(μp)。模型中土料的宏、微觀參數(shù)可先通過相關(guān)公式來快速確定，然后經(jīng)訓(xùn)練優(yōu)化后獲得[23,29]，見表5，6。

表5 模型中需對(duì)應(yīng)的宏參數(shù)Table 5 Corresponding macro parameters in the model

表6 模型中需對(duì)應(yīng)的微參數(shù)Table 6 Corresponding micro-parameters in the model

3.2 幾何模型搭建

MatDEM軟件不僅可以模擬出降雨入滲至淺層滑坡的啟動(dòng)時(shí)刻，還可以模擬出淺層滑坡土體的具體破壞過程。根據(jù)內(nèi)排土場(chǎng)剖面模型尺寸搭建約252 534個(gè)球單元的二維初始模型，其中單元平均直徑約為0.001 1 m，平均質(zhì)量約為2.267 mg。根據(jù)模型試驗(yàn)結(jié)果，得出最終包括約80 687個(gè)球單元的二維幾何模型，如圖10所示。建立內(nèi)排土場(chǎng)外形并將訓(xùn)練后的微參數(shù)賦值給干土層，即圖10藍(lán)色部分。當(dāng)雨水入滲至淺層滑坡啟動(dòng)時(shí)刻，坡表土體的內(nèi)摩擦因數(shù)降低至0.425。

圖10 二維土坡模型Fig.10 Two-dimensional soil slope model diagram

3.3 降雨入滲模型

參考前人對(duì)多場(chǎng)耦合DEM模型的研究[30-32]，考慮了水分的不均勻分布和傳遞，提出了降雨入滲并軟化土壤的離散元模型。離散土粒之間發(fā)生機(jī)械和水力相互作用，使土壤軟化，土粒的內(nèi)摩擦因數(shù)隨其含水量而變化。在宏觀層面的模擬過程中，以坡表濕潤(rùn)向形成飽和區(qū)與過渡區(qū)的趨勢(shì)進(jìn)行。

在本模型中，每個(gè)離散土粒都是大量土壤顆粒、孔隙和孔隙水的組合，而不是單個(gè)黏土礦物或微觀孔隙。2個(gè)土粒之間的聯(lián)系是土壤的一部分與另一部分之間的聯(lián)系，以及允許水流動(dòng)的路徑，如圖11(a)所示。除了半徑、質(zhì)量和速度等基本屬性外，還為土粒分配了一個(gè)參數(shù)，該參數(shù)表示一定體積土壤中的平均含水量，可由式(18)獲得：

Dω=Mω

(18)

式中，Dω為數(shù)值模型中單個(gè)土粒的平均含水量；M為單個(gè)土粒的平均質(zhì)量；ω為土壤含水率。

暴露在空氣中的土粒，定義為坡面土粒。降雨時(shí)，坡面土粒的水含量增加并變得大于相鄰?fù)亮５乃俊Ｈ缓笏μ荻闰?qū)動(dòng)水從富水土粒流向缺水土粒。水流量qw與含水量梯度正相關(guān)，計(jì)算公式為

qw=KAdω/r12

(19)

式中，K為滲透系數(shù),2.08×10-7m/s；A為滲透路徑等效面積，取2個(gè)土粒大圓的平均面積，兩粒大圓半徑分別為R1和R2，如圖11(b)所示；dω和r12為2個(gè)土粒之間的含水量差及其球心之間的距離。

圖11 降雨入滲和水分運(yùn)移Fig.11 Rainfall infiltration and moisture transport

雨水附著上坡面土粒的過程是雨水入滲擠壓和替代坡面土中空氣的水-氣二相流過程。坡面土粒雨水附著流量iω可借助Green-Ampt提出干土在表面積水條件下入滲的式(20)[33]來計(jì)算替代：

(20)

式中，KS為飽和滲透系數(shù)；Iq為降雨強(qiáng)度。

當(dāng)土粒的含水率增加時(shí)，它的內(nèi)摩擦因數(shù)會(huì)減低。一般來說，內(nèi)摩擦角響應(yīng)于含水量的變化呈非線性衰減趨勢(shì)，符合一階指數(shù)衰減。一些研究人員提出了一個(gè)實(shí)用模型[34]，具體計(jì)算見式(21)：

φ=26exp(-ω/18.4)+21.9

(21)

坡面土趨于飽和及降雨入滲并軟化土壤的過程中單次迭代的主要程序如下：

(1)搜索坡面土粒并升高其含水量。坡面土粒暴露在空氣中，每次迭代時(shí)它們會(huì)按式(20)進(jìn)行增加含水量，這模擬了坡面土逐漸趨于飽和。

(2)計(jì)算相鄰元素之間的水分轉(zhuǎn)移。首先，計(jì)算模型中的水力梯度。然后可以從式(18)和(19)中獲得水分轉(zhuǎn)移量，從而模擬實(shí)現(xiàn)雨水入滲過程。

(3)分析雨水軟化。當(dāng)含水量增多時(shí)，土壤軟化，可以從式(21)計(jì)算出土粒因含水量上升而降低的內(nèi)摩擦因數(shù)。在微觀尺度上，單元參數(shù)需要重新計(jì)算。

(4)平衡模型。隨著單元參數(shù)的變化，模型可能變得不平衡。因此，需要更新元素位置和連接狀態(tài)。如果2個(gè)單元之間的法向力或切向力超過相應(yīng)的限制，則粒子間的連接會(huì)斷裂，將不會(huì)有水分轉(zhuǎn)移。

重復(fù)步驟(1)～(4)，隨著模擬的進(jìn)行，坡面土的水分逐漸被附著，雨水逐漸入滲，坡面土壤趨于飽和，直到含水量上升到預(yù)設(shè)的最終值，相應(yīng)的土壤內(nèi)摩擦因數(shù)減低到最低值。

4 模擬結(jié)果與分析

4.1 降雨入滲模型

圖12模擬了降雨使坡面土體的含水量增加，雨水進(jìn)而逐漸入滲的過程。坡面土濕潤(rùn)后，較長(zhǎng)時(shí)間內(nèi)其含水量沒明顯上升，但雨水入滲深度逐漸加深(圖12(a)～(d))。當(dāng)雨水入滲到一定深度后，坡面土的含水量開始有明顯上升(圖12(e)～(f))。這與圖5所展示的坡面土體積含水率實(shí)測(cè)曲線基本相符，兩者均體現(xiàn)出前期坡面土的雨水附著率小于雨水入滲率。從式(19)，(20)可知，在沒有入滲的情況下，坡面土含水量將進(jìn)行線性增加，直至達(dá)到飽和含水量。當(dāng)坡面土濕潤(rùn)到一定程度時(shí)，坡面處顆粒間開始存有含水量差并加大，水流量也相應(yīng)增大，這會(huì)導(dǎo)致坡面處的雨水附著率小于水流量率(圖12(a)～(d))。隨著雨水入滲深度的逐漸加大，導(dǎo)致水分運(yùn)移路徑加大，從而確保了坡面處顆粒水分的積累以及維持了顆粒間較小的含水量差。此時(shí)坡面土含水量才開始有較為明顯的上升(圖12(e)～(f))。

圖12 降雨入滲至淺層滑坡啟動(dòng)過程Fig.12 Diagram of the starting process of rainfall infiltration to shallow landslide

坡體含水率的漸變分布與改進(jìn)后Green-Ampt模型里飽和區(qū)、過渡區(qū)和干土區(qū)的分布基本相符。隨著雨水的不斷入滲，坡體出現(xiàn)明顯的近飽和區(qū)，土體強(qiáng)度逐漸降低，此時(shí)邊坡呈現(xiàn)整體失穩(wěn)的特征，圖12(f)中的坡體有較小的滑動(dòng)，說明此時(shí)的雨水入滲軟化狀態(tài)已經(jīng)啟動(dòng)了淺層滑坡。

4.2 淺層滑坡的運(yùn)動(dòng)過程

數(shù)值模型經(jīng)過力平衡迭代到穩(wěn)定狀態(tài)。圖13顯示了滑坡在不同時(shí)刻的變形破壞狀態(tài)。在初始階段，淺層滑坡土層的斜坡表面先產(chǎn)生變形(圖13(a))，淺層滑坡土層后開始整體滑動(dòng)變形(圖13(b))，分別對(duì)應(yīng)室內(nèi)降雨入滲試驗(yàn)的徑流侵蝕剛連通和淺土層開始平行滑坡時(shí)間。隨后，位于底部的淺層滑坡土層抵達(dá)坡腳，滑坡舌也逐漸成形(圖13(c)～(e))。在淺層滑坡土層的滑動(dòng)帶動(dòng)下，干土層的接觸表面也開始部分滑動(dòng)，直至整個(gè)滑坡體經(jīng)過坡腳減速并靜態(tài)沉積在后部平臺(tái)上(圖13(f)～(h))。淺層滑坡模擬結(jié)果與物理模擬試驗(yàn)中淺層滑坡結(jié)束后的坡面侵蝕情況基本吻合(圖13(h)，圖4(c))。此外，持續(xù)的降雨會(huì)在新斜坡處產(chǎn)生新的浸潤(rùn)土層，加之原先坡頂未滑動(dòng)土體受雨水的進(jìn)一步入滲，故淺層滑坡很容易被降雨重新激活。

圖13 滑坡運(yùn)動(dòng)過程Fig.13 Landslide movement process

4.3 速度

圖14顯示了整個(gè)滑坡體在不同時(shí)刻的速度分布。當(dāng)產(chǎn)生明顯位移，滑坡體向下移動(dòng)，邊坡整體失穩(wěn)。在早期階段，淺層滑坡土層的表面首先能量積聚(圖14(a))，隨著能量積聚不斷增加，淺層滑坡土層的坡面較快運(yùn)動(dòng)(圖14(b))，即雨水持續(xù)入滲過程中徑流侵蝕已連通。然后，隨著滑坡體的運(yùn)動(dòng)，滑坡體的重力勢(shì)能部分轉(zhuǎn)化為動(dòng)能，滑坡體的速度繼續(xù)增加(圖14(c)～(e))?；伦畲笏俣?約1.58 m/s)在經(jīng)過0.633 5 s后出現(xiàn)(圖14(e))。此后，由于滑坡體內(nèi)部顆粒的碰撞和摩擦，速度逐漸降低(圖14(f)～(g))。1.077 s時(shí)，滑坡體相對(duì)靜止，沉積在坡腳處及后方平臺(tái)上(圖14(h))。此外由圖14可知，從坡面到滑動(dòng)面的滑動(dòng)速度逐漸減小，換言之，滑動(dòng)面附近的土體比位于滑體上部的土體移動(dòng)得慢[18]。這與內(nèi)排土場(chǎng)無黏性土邊坡越趨于飽和的坡表土越先啟動(dòng)且越易滑動(dòng)相符合。

圖14 不同時(shí)刻滑坡的速度場(chǎng)Fig.14 Velocity field of the landslide at different times

圖15給出了滑體的運(yùn)動(dòng)速度變化曲線?！扒?.25%”曲線代表速度在前6.25%的土體平均速度，依此類推。由于滑體速度向內(nèi)部遞減，“前6.25%”曲線反映滑坡表層土體的平均滑動(dòng)速度。最大速度反映表面跳躍土體的速度。因此，圖15在一定程度上反映了淺層滑坡體由底部向表層的速度變化規(guī)律。從圖15中可以看到，滑體速度在前0.38 s快速增加至表層約1.5 m/s，平均約1.1 m/s；滑體的平均速度在0.38～0.45 s時(shí)達(dá)到最大，此時(shí)段表層的平均速度約為整個(gè)滑體平均速度的3.5倍，揭示了滑坡速度由表層向下遞減；而在0.93 s時(shí)，表層速度約為整個(gè)滑體平均速度的5倍，單元速度向下快速減小，這與模型試驗(yàn)中的觀察到的現(xiàn)象一致；0.56 s后，最大速度出現(xiàn)強(qiáng)烈跳動(dòng)，反映少量土體的跳躍運(yùn)動(dòng)。總體上看，滑坡初始時(shí)運(yùn)動(dòng)速度迅速增加，并保持短時(shí)的快速運(yùn)動(dòng)，然后逐漸減小，最后逐漸趨于穩(wěn)定。

圖15 滑坡體單元平均速度變化曲線Fig.15 Variation curves of average velocity of landslide body unit

4.4 位移

滑坡在不同時(shí)刻的位移變化如圖16所示，模擬滑坡的最大位移約為0.7 m。在早期階段，滑坡體開始表層滑動(dòng)且初始移動(dòng)緩慢，然后，隨著滑坡體滑動(dòng)范圍加大及重力勢(shì)能的轉(zhuǎn)換，位移逐漸增加?；瑒?dòng)過程中，位移從滑動(dòng)面向坡面逐漸增大，滑坡中部和上部的土體比下部土體滑動(dòng)得更遠(yuǎn)，呈現(xiàn)類似“紊流”特征。隨后，由于空氣阻尼和顆粒之間發(fā)生的摩擦碰撞，導(dǎo)致滑坡體減速。最后，位移逐漸減小，直至靜止。在整個(gè)滑動(dòng)過程中，滑體前緣的位移始終比后緣大。因此，可以添加類似聚丙烯酰胺等化工材料到坡表土來增加表層土黏性，進(jìn)而束縛住表層土，這會(huì)在一定程度上降低淺層滑坡中的坡面侵蝕規(guī)模[35-36]。

圖16 滑坡在不同時(shí)刻的位移變化Fig.16 Displacement changes of landslides at different moments

4.5 能量轉(zhuǎn)換

圖17給出了各種形式的能量演化過程。在整個(gè)滑動(dòng)過程中，重力勢(shì)能主要轉(zhuǎn)化為熱能和動(dòng)能，無彈性勢(shì)能的轉(zhuǎn)化。0.271 5 s后，整個(gè)滑體逐漸脫離母體(圖13(d))。隨著滑體向下運(yùn)動(dòng)，重力勢(shì)能不斷減小，動(dòng)能逐漸增大(圖17(a))。

由于空氣阻尼的存在，產(chǎn)生少許阻尼熱，滑動(dòng)帶各顆粒之間的摩擦產(chǎn)生大量的摩擦熱(圖17(b))。因此，滑體中的熱量逐漸增加，機(jī)械能逐漸減少?？梢杂^察到，動(dòng)能最初保持穩(wěn)定，然后逐漸增加，最后逐漸降低。1.077 0 s后，當(dāng)滑坡運(yùn)動(dòng)停止時(shí)，動(dòng)能接近0。同時(shí)，熱量顯著增加，表明動(dòng)能已經(jīng)轉(zhuǎn)化為熱能，說明摩擦是整個(gè)滑坡運(yùn)動(dòng)過程中能量耗散的主要來源。因此，在條件允許下，應(yīng)以較緩坡度對(duì)內(nèi)排土場(chǎng)邊坡進(jìn)行堆積，這會(huì)有效地在淺層滑坡過程中降低滑坡體的動(dòng)能積攢且增加滑坡體間的摩擦損耗[37-38]。

圖17 能量、熱量時(shí)間歷史曲線Fig.17 Energy,thermal time history curves

5 結(jié) 論

(1)降雨入滲條件下，內(nèi)排土場(chǎng)坡面侵蝕由濺蝕凹槽開始，整體呈現(xiàn)為溯源侵蝕和淺層破壞特征，隨著降雨的持續(xù)，呈達(dá)到飽和—淺層滑坡—新表土層達(dá)到飽和的模式循環(huán)出現(xiàn)。

(2)在考慮實(shí)際非飽和浸潤(rùn)層的基礎(chǔ)上，獲得了改進(jìn)的Green-Ampt入滲模型，并依此推導(dǎo)出內(nèi)排土場(chǎng)穩(wěn)定性分析方程。

(3)根據(jù)雨水入滲深度，內(nèi)排土場(chǎng)穩(wěn)定系數(shù)的變化可劃分為高敏感、敏感削弱和低敏感3個(gè)階段。

(4)滑坡體呈浸潤(rùn)峰之上土體的淺層破壞特征，內(nèi)排土場(chǎng)滑坡體各層均呈現(xiàn)滑坡初始時(shí)運(yùn)動(dòng)速度迅速增加，并保持短時(shí)的快速運(yùn)動(dòng)，然后逐漸減小，最后逐漸趨于穩(wěn)定的規(guī)律，但滑坡表層、前緣土體的運(yùn)動(dòng)速度比底層、后緣高。

(5)滑動(dòng)期間由勢(shì)能轉(zhuǎn)換來的大部分動(dòng)能被滑帶附近的摩擦所消散，持續(xù)的降雨會(huì)產(chǎn)生新的飽和層和浸潤(rùn)層進(jìn)而使淺層滑坡重復(fù)發(fā)生。