基于Linear Cohesion模型的含水堆積體滑坡過程數(shù)值模擬

高道路， 王林峰

(1.重慶交通大學(xué) 土木工程學(xué)院， 重慶 400074； 2.重慶交通大學(xué) 河海學(xué)院， 重慶 400074)

1 研究背景

近年來，世界各地發(fā)生滑坡等地質(zhì)災(zāi)害的頻率越來越高，造成了大量人員傷亡和財產(chǎn)損失[1-3]?；掳磶r體性質(zhì)可劃分為巖質(zhì)滑坡、土層滑坡和松散堆積層滑坡3類，其中，松散堆積層滑坡由于巖體結(jié)構(gòu)組成極為復(fù)雜，如何準確地計算其運動特征是當前的研究熱點和難點[4]。

目前，滑坡災(zāi)害的研究主要包括室內(nèi)試驗和數(shù)值分析，小規(guī)?；驴梢酝ㄟ^室內(nèi)試驗進行模擬并得到較為可靠的結(jié)果，但數(shù)值分析法可以不受尺度的限制，能夠模擬更加真實尺寸的滑坡運動過程[5-6]。Zhou等[7]采用離散元法中的干顆粒流模型研究了固體顆粒沿斜坡通道的接觸行為，確定了顆粒流滑動速度分布和相應(yīng)的剪切速率。劉思杰等[8]采用三維顆粒離散元對不同粒徑的碎屑流運動進行分析，發(fā)現(xiàn)粒徑組成越不均勻，滑坡運動速度和破壞力越大，滑動距離也越遠。Lo等[9]利用離散元程序PFC對某滑坡的運動和土體破壞的運動過程進行了系統(tǒng)分析，發(fā)現(xiàn)不同摩擦系數(shù)和粘結(jié)強度對滑坡堆積范圍的影響。Mead等[10]采用數(shù)值分析法模擬了在不規(guī)則地形上干燥細沙堆積體的運動過程，數(shù)值試驗結(jié)果與室內(nèi)物理試驗結(jié)果基本一致，證明了采用數(shù)值方法模擬堆積體滑坡是合理可行的。

大量相關(guān)文獻研究成果表明，采用數(shù)值試驗分析計算滑坡的運動過程可以得到較為準確的結(jié)果，特別是基于非連續(xù)介質(zhì)理論的離散元法應(yīng)用尤為廣泛?？傮w來說，已有大量學(xué)者采用數(shù)值分析法模擬各類滑坡體的運動，并取得了豐碩成果，但在滑坡運動模擬過程中直接考慮水的黏附作用的相關(guān)研究還很少。

因此，本文在傳統(tǒng)離散單元法軟球模型的基礎(chǔ)上，直接引入顆粒單元間的液橋力，提出采用Linear Cohesion接觸模型模擬含水濕潤狀態(tài)下顆粒堆積體的滑動過程，并進一步研究不同液橋力大小對堆積體滑坡各項運動特征的影響。

2 接觸模型簡介

離散單元法(Discrete element method，以下簡稱DEM)是Cundall等[11]在20世紀70年代首次提出的描述節(jié)理巖體和非連續(xù)介質(zhì)力學(xué)行為的有效數(shù)值方法。DEM認為節(jié)理巖體由離散的塊體和塊體之間的節(jié)理面組成，不僅可以考慮巖體內(nèi)部結(jié)構(gòu)的大位移、旋轉(zhuǎn)和滑動，甚至還可以計算塊體的分離過程，能客觀真實地反映裂隙巖體非線性大變形的特點，特別適用于裂隙巖體的應(yīng)力變形分析[12]。與傳統(tǒng)的有限元方法相比，DEM方法提供了豐富的顆粒單元接觸模型，可以很好地解決裂隙巖體的大變形問題。

本文主要計算松散堆積體在干燥和含水濕潤兩種狀態(tài)下的滑動過程。其中，針對松散堆積體在干燥狀態(tài)下的失穩(wěn)、滑動和堆積整個過程，接觸模型采用了常用的Hertz-Mindlin滑動模型，如圖1所示，為兩顆粒單元i和j的接觸示意，圖1中所示的kn、kt、kr分別為法向剛度、切向剛度和滾動剛度，dn、dt、dr分別為法向阻尼、切向阻尼和滾動阻尼。

圖1 Hertz-Mindlin顆粒接觸模型

對于模擬含水濕潤狀態(tài)下松散堆積體滑動過程，由于兩相互接觸的顆粒間形成了液橋，因此，需引入液橋力來模擬。本文采用Linear Cohesion接觸模型模擬含水濕潤狀態(tài)下松散堆積體滑動過程，該模型是在傳統(tǒng)Hertz-Mindlin軟球無滑動接觸模型的基礎(chǔ)上增加一個法向黏聚力。此時，堆積體顆粒在滑動過程中受3種力的作用，包括重力、顆粒間接觸力和液橋力。除此之外，互相接觸的顆粒單元還受到切向力矩和滾動摩擦力矩的作用，規(guī)定第i個顆粒的運動方程[13-14]如下:

(1)

(2)

式中:m為顆粒的質(zhì)量，kg；I為顆粒的轉(zhuǎn)動慣量,kg/m2；g為重力加速度,m/s2；ni為與小球i接觸的顆?？倲?shù)量；V為運動速度,m/s；ω為角速度,rad/s；t為時間,s；Fn,ij和Ft,ij分別為法向和切向作用力,N；Fcoh,ij為兩顆粒間的液橋力,N；Tt,ij和Tr,ij分別為切向力矩和滾動摩擦力矩,N·m。Linear Cohesion接觸模型模中的液橋力Fcoh定義如下：

Fcoh=kA

(3)

式中:A為兩顆粒的接觸面積m2;k為顆粒間的粘結(jié)能量密度,J/m3。

3 模型建立與計算條件

利用Auto-CAD建立了三維松散堆積滑坡數(shù)值試驗?zāi)Ｐ停⑵鋵?dǎo)入離散元程序EDEM2017中進行計算，建立的滑坡數(shù)值模型如圖2所示，包括底板、側(cè)板和堆積體物源。底板水平長度為20 m，坡長16.77 m，水槽寬度為1.5 m，滑坡模型的坡角設(shè)置為27°(1∶2)，邊坡垂直高度為7.5 m，斜槽水平投影長度為15 m，松散堆積物由1 470個顆粒組成，顆粒形狀參照實際巖屑進行建模，如圖3所示。堆積體物源總質(zhì)量為5 687 kg，堆積物總體積約為3.2 m3。

圖2 堆積體滑坡數(shù)值模型

圖3 松散堆積體的顆粒單元

此外，采用Linear cohesion模型模擬含水松散堆積體滑坡運動時，需要輸入相應(yīng)的能量密度k，由于缺乏相關(guān)的資料作為參考，本文通過反復(fù)迭代試算，并討論不同能量密度取值對含水松散堆積體滑坡運動特性的影響，最終采用的能量密度從1 000 J/m3開始遞增至6 000 J/m3。其余計算參數(shù)還包含顆粒密度、泊松比和摩擦系數(shù)等，參考了類似的顆粒離散元數(shù)值試驗并進行優(yōu)化[15]，詳細參數(shù)取值見表1。

4 結(jié)果與分析

在本次松散堆積體滑坡數(shù)值試驗中，計算研究了干燥顆粒和含水濕潤顆粒的運動過程和堆積特點，包含能量密度k取0(干燥狀態(tài))、1 000、2 000、3 000、4 000、5 000和6 000 J/m3等7種計算條件，運動過程中不同時刻的滑動速度如圖4～8所示，由于圖幅較多，僅給出部分結(jié)果。

表1 離散元模型相關(guān)參數(shù)

分析可知，當能量密度k取0時，即原堆積體是干燥的，此時松散顆粒堆積體滑動時具有良好的流動性，滑槽上顆粒的速度在滑動過程中不斷變化。在t=3.1 s時，前顆粒到達水平面，最大速度超過7 m/s，最終當t=8 s時運動基本停止，顆粒流堆積在斜槽和平面上，運動的距離較長(圖4)。此外，為了分析能量密度k對松散堆積體滑坡運動的影響，計算了松散濕顆粒在考慮能量密度的情況下的滑動過程(圖5～8)。分析可知，隨著能量密度k的逐漸增大，顆粒堆積體滑動時的流動性越來越小，如圖6(k=2 000 J/m3)所示，堆積體在斜坡上運動時前段解體為小顆粒群運動，而后段滑動時基本保持結(jié)構(gòu)完整，最終當t=8 s時運動基本停止，但此時堆積體主要堆積在平面上，極少數(shù)顆粒堆積在斜坡上。能量密度越大，堆積體滑動到坡腳位置所消耗的時間也越長，如k取0，為3.1 s，而取1 000和2 000 J/m3時延遲到3.4和3.7 s。

當能量密度k超過4 000 J/m3時，松散濕顆粒堆積體在斜槽上滑動時(t=0～4 s)均能保持初始狀態(tài)的結(jié)構(gòu)完整，這是由于足夠的能量密度k賦予了顆粒間較大的液橋力，可以抵消外界一定強度的作用力，當所受力不超過液橋力時，顆粒粘結(jié)結(jié)構(gòu)依舊能保持完整。如圖7中當k=4 000 J/m3時，僅當松散堆積體沖撞到底部平面時，由于慣性作用，前端顆粒速度遠大于后端顆粒，堆積體前端開始出現(xiàn)拉裂隙，最前端的顆粒逐漸團脫離母體運動。而圖8中當k=6 000 J/m3時，整個堆積體滑動過程中結(jié)構(gòu)一直保持初始的狀態(tài)，即便在水平面上運動時顆粒粘結(jié)依舊未被破壞。

以上分析表明，較小的能量密度可以使松散濕顆粒堆積物在滑動過程中發(fā)生崩解，形成數(shù)個破碎的巖屑顆粒群單獨運動，并形成較遠的滑動距離。

圖9給出了不同能量密度參數(shù)時松散堆積體整體平均運動速度和動能隨時間的變化情況。

圖9表明，隨著堆積體顆粒間的能量密度逐漸增大，滑動達到最大平均速度和動能峰值的時刻越來越早，數(shù)值越來越大，但是當k≥5 000 J/m3時，堆積體滑動的各項運動特征差別非常小。

圖4 不同時刻干燥堆積體滑動速度分布(k=0)

圖5 不同時刻含水堆積體滑動速度分布(k=1 000 J/m3)

圖6 不同時刻含水堆積體滑動速度分布(k=2 000 J/m3)

圖7 不同時刻含水堆積體滑動速度分布(k=4 000 J/m3)

圖8 不同時刻含水堆積體滑動速度分布(k=6 000 J/m3)

例如，當k取5 000和6 000 J/m3時，堆積體運動的最大速度分別為6.69和6.71 m/s，最大動能分別為127和128 kJ，差別很小，對應(yīng)的時刻均在t=3.8 s。而當k取2 000 J/m3時，堆積體運動的最大速度和最大動能分別為4.65 m/s和63.5 kJ，對應(yīng)的時刻為t=4.4 s。

進一步分析發(fā)現(xiàn)，當k≥4 000 J/m3時，堆積體滑動速度曲線峰值前段幾乎呈線性增加，這是由于此時的堆積體結(jié)構(gòu)在坡面保持整體運動。動能曲線圖9(b)與9(a)中的速度曲線趨勢基本一致，因為速度與動能之間的關(guān)系為Ek=1/2(mv2)。此外，當k不超過3 000 J/m3時，松散堆積體取不同能量密度時滑動最大速度和最大動能差異不大。

圖10為不同能量密度松散堆積體x、y、z方向的滑動速度。在x方向(圖10(a))，沿滑動運動方向，變化趨勢與圖9(a)中體積速度的變化趨勢非常相似，同樣，能量密度越大，堆積體達到x方向的最大平均速度越早，最大平均速度越大。在y方向(圖10(b))，速度負值表示垂直向下的滑動運動，當能量密度k<3 000 J/m3時，堆積體y方向的最大平均速度差別不大，約在1.8 m/s左右，而當k≥5 000 J/m3時，y方向的最大平均速度幾乎可以達到3 m/s。此外，在z方向，該速度代表濕顆粒沉積物的橫向運動，所有情況下z方向的最大平均速度約為2～3 cm/s，集中在4～6 s，由于側(cè)板的阻擋作用，其大小比x和y方向小兩個數(shù)量級，如圖10(c)所示。

圖9 不同k值松散堆積體滑動平均速度和動能變化

圖10 不同k值松散堆積體滑坡x、y、z方向平均速度變化

表2給出了松散堆積體滑動后的堆積長度統(tǒng)計結(jié)果，結(jié)合圖4～8可知，隨著堆積體能量密度的逐漸增加，斜坡上和平面上的堆積長度均有逐漸減小的趨勢，如k=1 000 J/m3時，斜坡上的堆積長度和平面上的堆積長度分別為4.72和7.59 m，當能量密度增加到3 000 J/m3，堆積長度分別減小到0.96和4.73 m，當能量密度增加到5 000 J/m3，堆積長度進一步減小到0和5.30 m。以上分析說明，能量密度較小的情況下，濕顆粒堆積體流動性更強，在實際滑坡運動過程中，滑坡災(zāi)害的影響范圍更廣。

表2 含水堆積體滑坡滑動方向堆積長度

圖11為不同能量密度堆積體滑動過程中顆粒接觸數(shù)量的變化情況。結(jié)果表明，當能量密度不超過3 000 J/m3時，顆粒的接觸數(shù)量在第1階段(t=0～4 s)逐漸減少，然后逐漸增加，最后趨于穩(wěn)定并保持不變。而當能量密度k≥4 000 J/m3時，顆粒的接觸數(shù)量在第1階段(t=0～4 s)保持不變，在t=4 s時出現(xiàn)陡然增加，這是由于堆積體剛好撞擊到底部平板上，然后接觸數(shù)量逐漸減小，最后趨于穩(wěn)定并保持不變。

圖11 松散堆積體滑動過程顆粒接觸數(shù)量變化

5 結(jié) 論

本文在傳統(tǒng)離散單元法軟球模型的基礎(chǔ)上，提出了一種采用Linear Cohesion接觸模型模擬含水濕潤狀態(tài)下的松散堆積體滑坡運動過程的分析方法，該模擬可以有效地考慮相互接觸顆粒單元之間的液橋力，并引入了能量密度參數(shù)k進行表征。通過建立含水堆積體滑坡離散元數(shù)值試驗?zāi)Ｐ?，研究了不同能量密度值對堆積體滑坡各項運動特征的影響，得到以下結(jié)論：

(1)隨著堆積體顆粒間能量密度逐漸增大，堆積體滑動時的流動性越來越小，堆積體滑動到坡腳位置所消耗的時間也越長，達到最大平均速度和動能峰值的時刻越來越早，數(shù)值越來越大。當k≥5 000 J/m3時，堆積體滑動的各項運動特征差別非常小。

(2)較小能量密度的松散濕顆粒堆積物在滑動過程中極易發(fā)生崩解，形成的巖屑顆粒群單獨運動，流動性更強，在斜坡和平面上的堆積長度更長，此類滑坡災(zāi)害的影響范圍更廣。

(3)能量密度對含水堆積體滑坡過程中顆粒的碰撞接觸影響較大，當能量密度k≤3 000 J/m3時，顆粒的接觸數(shù)量先減少然后增加，最后趨于穩(wěn)定；當能量密度k≥4 000 J/m3時，接觸數(shù)量先保持不變，之后陡然增加，最終趨于穩(wěn)定。

以上結(jié)論表明Linear Cohesion接觸模型可以較好地反映含水濕顆粒間的相互作用，可用于模擬含水松散堆積體滑坡的整個運動過程，為今后此類滑坡災(zāi)害的數(shù)值分析預(yù)測和災(zāi)害評估提供一種新的研究途徑。能量密度k是Linear Cohesion接觸模型模擬含水堆積體顆粒運動的關(guān)鍵參數(shù)，在今后的研究中有必要對堆積體的含水量與膠結(jié)能量密度的關(guān)系進行更深入的探討，進而可以根據(jù)不同含水條件選取相應(yīng)的能量密度值進行含水松散堆積體滑坡的數(shù)值模擬計算。