基于DEM的碎屑流運動特性數值模擬

趙 川,付成華

(西華大學能源與動力工程學院,四川 成都 610039)

基于DEM的碎屑流運動特性數值模擬

趙 川,付成華

(西華大學能源與動力工程學院,四川 成都 610039)

采用顆粒離散元法建立了碎屑流三維數值模型,分別從整體和局部對碎屑流運動特性進行深入分析。結果表明:碎屑流的運動過程可分為啟動加速、高速滑動和減速堆積3個階段;邊坡底面平整時,堆積體中大顆粒堆積在表面,而細小顆粒位于中下部位;碎屑流運動速度和能量損失受顆粒間的接觸碰撞影響大;碎屑流不同部位運動速度最大值相同,達到速度峰值的時間與顆粒的初始水平位置有關,與深度無關;堆積體分區明顯,堆積體中的顆粒按初始水平位置依次排布;坡腳位置堆積深度最大,往兩側逐漸減小;碎屑流在運動初期和末期橫向運動明顯,擴大了災害影響范圍。

碎屑流;三維顆粒離散元法;碰撞接觸;堆積體;數值試驗

碎屑流是一種滑坡巖體在初始狀態或者高速運動過程中轉化為碎屑物質的運動形式,其某些運動特性與流體相似。碎屑流是一種常見的自然災害,有運動速度快、攜帶巖體量大和破壞力強等特點,對人民生命財產安全構成巨大威脅。我國地形條件復雜,高山丘陵廣泛分布,發生地質災害的頻率很高,尤其是我國西南地區在遭受了汶川地震后,邊坡表層巖體破碎嚴重,更是碎屑流災害的高發地帶[1]。黃潤秋[2]分析了20世紀以來中國的大型滑坡,將其分為巖質滑坡、土層滑坡和松散堆積層滑坡3類。本文研究的碎屑流不考慮顆粒單元間的黏結作用,屬于松散堆積層滑坡。

20世紀以來,國內外學者對顆粒碎屑物質進行了大量研究，《Science》于2005年把顆粒物質列為100個科學難題之一[3-4];Iverson等[5]基于復雜的三維地形,對沙粒崩塌(granular avalanche)進行了室內物理實驗,并采用連續介質模型進行模擬,模擬結果與實驗結果基本吻合;Hussin等[6]以MassMov2D為工具,通過反分析確定了某滑坡帶的計算參數,對潛在的滑坡進行了預測;王玉峰等[7]對汶川地震觸發的高速遠程碎屑流堆積反粒序特征進行了分析,揭示了碎屑流堆積體內部所特有的堆積規律;張明等[8]對目前碎屑流的研究成果進行了系統的歸納和總結。總的看來,已有的碎屑流研究對象主要針對室內小尺寸模型或已經發生的具體滑坡。具體的滑坡案例受地形條件影響大,室內模型難以反映真實高速碎屑流的某些運動特性,其研究結論難以反映一般規律[9-12]。

為了解決上述碎屑流研究存在的問題,本文采用三維顆粒離散元法(discrete element method，DEM)建立了標準的碎屑流數值模型,選用含有棱角的顆粒單元,并通過漏斗模型試驗標定計算參數,對碎屑流運動特性進行了深入分析探討。

1 碎屑流數值模型建立及參數選取

1.1 離散元模型構建

碎屑流數值模型主要由碎屑堆積槽、滑動面和堆積面組成。滑動面和堆積面構成一個邊坡系統，坡高為10 m,坡比為1∶2,堆積面長為10 m;碎屑堆積槽用于放置顆粒物質,槽外形為三棱柱體,尺寸為1.5 m×3 m×5 m,如圖1所示。

圖1 碎屑流離散元模型尺寸(單位:m)

由于真實巖體碎屑外形并非標準球體,而是含有一定棱角的球體,因此模型試驗中構建了含有棱角的顆粒單元,如圖2所示,其中圖2(a)為真實碎屑流塊體,圖2(b)為構建的離散元塊體。顆粒大小按高斯分布隨機生成,直至充滿堆積槽,顆粒總數為3 707,總填充體積22.5 m3,最終顆粒模型試樣如圖3所示。

圖2 數值計算顆粒單元外形

圖3 碎屑流離散元模型

滑坡碎屑流運動過程中,巖塊之間存在大量的碰撞和接觸,為準確描述DEM中顆粒單元運動過程的力學行為,計算采用Hertz-Mindlin無滑動接觸模型,主要包括法向力Fn、法向阻尼力Fnd、切向力Ft和切向阻尼力Ftd4個參數,計算公式[13-14]如下:

(1)

(2)

Ft=-Stδt

(3)

(4)

式中:R*為等效顆粒半徑,m;Y*為等效楊氏模量,Pa;δn、δt分別為兩個顆粒接觸面法向與切向重疊量,m;m*為等效質量,kg;vn、vt分別為兩個顆粒接觸面法向和切向相對速度,m/s;β為阻尼力系數;Sn、St分別為兩個顆粒接觸面法向與切向剛度。

1.2 數值計算參數標定

碎屑流的運動受摩擦系數的影響很大,為校核相關摩擦系數,對試驗選用的顆粒采用漏斗模型試驗進行標定,如圖4所示。由于本試驗中碎屑顆粒間不存在黏結力,其天然休止角即為松散堆積體的內摩擦角φ。數值試驗中滑動面為平面,為模擬邊坡表面覆蓋的松散堆積層[15],基底摩擦系數取0.42,顆粒間摩擦系數取0.4。具體方法是首先設置好基底摩擦系數和顆粒間摩擦系數,通過不斷調整滾動摩擦系數重復計算,將堆積體的天然休止角為22°(arctan 0.4=arctanH/R=21.8°,其中H為堆積體中心高度,R為堆積體底半徑)時所采用的參數作為碎屑流模型的參數,最終得到的計算參數見表1。

圖4 離散元顆粒漏斗數值模型試驗

表1 DEM模型相關參數取值

圖5 不同時刻碎屑流運動速度分布(單位：m/s)

2 數值試驗結果及分析

2.1 整體運動特性

為分析巖體碎屑流不同時刻的運動狀態,圖5給出了碎屑流的運動速度分布。在時間t=0 s時,撤掉試樣前端擋板,在重力作用下,顆粒開始沿著坡面向下運動。由于顆粒并不是標準的球形顆粒,而是含有棱角的顆粒單元,所以整個運動過程中主要以平面滑動為主,僅在顆粒流表面存在滾動現象。運動過程中,前端顆粒由于所受約束小,運動速度較大,最大顆粒速度超過7 m/s。t=4 s時,前端顆粒運動到堆積面,并與底面發生碰撞,動能損失較大,速度逐漸減小,并堆積在坡腳位置。最終形成的堆積角在10°左右,僅為初始坡角26.6°的38%。

圖6為整個運動過程顆粒物質平均速度隨時間的變化,10 s時碎屑的平均速度僅為0.002 m/s,運動基本停止。根據曲線變化趨勢,可將碎屑流的運動過程劃分為啟動加速、高速滑動和減速堆積3個運動階段,分別對應的時間段為0～4 s、4～6 s和6～10 s,加速和減速階段所用的時間大體相當,均占總耗時的40%左右,而高速滑動僅占總耗時的20%。

圖6 碎屑流平均速度的變化過程

另外,從最后的堆積體顆粒分布可以發現,大顆粒位于堆積體表面,細小顆粒位于堆積體中下部位。究其原因,是由于在運動過程中,顆粒間相互碰撞摩擦,顆粒間存在一定大小的空隙,細小顆粒沿著空隙下落到堆積體底部,從而形成大顆粒物質在表面、細小顆粒在中下部的現象,這與實際碎屑流形成的堆積序列相符。

模擬時不考慮巖體間的黏結力,僅考慮相互碰撞作用。圖7為碎屑流運動過程中顆粒單元相互接觸數量隨時間的變化,初始時刻堆積較為緊密,總的接觸數量大約有11 500;隨著運動速度的增大接觸數量逐漸減小,在t=5.4 時接觸數量減小到最小值1 841,說明此時顆粒最為稀疏,顆粒間相互作用最少;此后隨著顆粒流前緣堆積體逐漸形成,接觸數量又逐漸增加,最后可達13 500左右,略大于初始時刻的數量。對比圖6和圖7可見,碎屑流的平均速度隨時間的變化曲線與單元接觸數量隨時間的變化曲線的變化趨勢正好相反,且速度的最大值正好對應接觸數的最小值,均發生在t=5.4 s時,說明顆粒越稀疏,運動速度越大,顆粒越密集意味著更多的碰撞,而碰撞將導致碎屑流內部能量大量耗散,故運動速度也隨之減小。

圖7 顆粒單元接觸數量的變化過程

圖8為單元的運動跡線和速度,可以發現大部分單元的速度分布基本是連續變化的,表現出一定的流體性質,但少數表面不連續的紅色線條顯示出存在顆粒單元在某個時段內脫離碎屑流運動,而單獨以較大的速度在表面發生跳動,這是由于表面部分高速運動的顆粒之間發生偏心碰撞和回彈造成的。

圖8 碎屑流DEM模型單元運動跡線及速度(單位:m/s)

2.2 局部運動特性分析

為分析碎屑流不同部位在運動過程中的差異,將模型試樣分成5個區域,每部分顆粒總數占總顆粒數的20%左右,并分別編號A、B、C、D、E,如圖9所示,其中A區代表后部,B區代表中上部,C區代表中下部,D區代表前端上部,E區代表前端下部。圖10為5個區域內顆粒的運動速度隨時間的變化過程,可以發現,位于前端的C和E區的顆粒最先啟動,在前0.5 s內,曲線斜率很大,說明速度增長很

圖9 顆粒分區

圖10 不同分區顆粒的運動速度

快,這是由于前端顆粒在初始時刻存在掉落現象,隨后曲線斜率保持穩定,在t=4.5 s時速度達到最大值;而位于中部的B和D區顆粒啟動稍晚,在加速運動過程中,速度增長保持穩定,在t=6.5 s左右速度達到峰值;位于最后端的A區則最后啟動,且初始時刻速度增加幅度不大,這是因為中部的B區和D區顆粒阻擋和延遲了A區顆粒的運動,隨后速度逐漸增加,并在t=8 s時速度達到最大值。綜合分析可見,5個區域的顆粒速度在不同時刻達到最大值,且最大值基本相同,約為4.5 m/s;B區顆粒速度稍小,為4.4 m/s。

圖11為5個區域內顆粒在水平方向x、豎直方向y、橫向z的運動速度分量,可以發現x方向的速度最大,最大峰值速度約4 m/s;y方向的速度次之,最大速度約2 m/s;z方向的速度最小,最大也僅有0.04 m/s。各部分顆粒的x和y方向速度幾乎同時達到最大值,且曲線變化趨勢也基本相同,其中C和E區顆粒在前0.5 s內y方向速度快速增大,這是顆粒突然釋放且前緣部分受到阻力很小造成的(圖11(b))。z方向顆粒運動速度一直較小,具體分析曲線變化趨勢可以發現:在碎屑流運動初期,位于前緣的C和E區顆粒橫向運動較為明顯;運動中期各部分顆粒z方向運動速度均很小,此階段內顆粒堆積密度最小,相互之間的碰撞擠壓也最少;在最后的堆積期間,A、C和D區的顆粒橫向運動速度相對大一些。總體來說,在運動初期和最后堆積期,顆粒沿運動方向速度較小時,z方向的速度相對較大;而當顆粒沿運動方向高速運動時,z方向速度相對較小(圖11(c))。

圖12為碎屑流最終的堆積狀態,A區顆粒完全堆積在坡面上,B和D區大部分堆積在坡面上,極少部分堆積在坡腳以下; C和E區的顆粒完全堆積在底面上,堆積長度較大,越往前堆積厚度越小;坡腳位置的堆積厚度最大,最大堆積厚度可達1 m。

圖12 碎屑流最終堆積狀態

3 結 論

a. 根據碎屑流速度時間曲線,可將其運動過程劃分為3個運動階段:啟動加速段、高速滑動段和減速堆積段。

b. 碎屑流在平整坡面上運動時,存在大量接觸和碰撞,小顆粒不斷通過空隙向下運動,最終造成了明顯的堆積分布序列:大顆粒位于堆積體表面,小顆粒位于堆積體底部,在坡腳位置堆積深度最大;碎屑流的速度隨時間的變化曲線與單元接觸數量隨時間變化曲線的變化趨勢恰好相反,且速度的最大值正好對應接觸數量的最小值。

c. 不同區域內顆粒速度最大值基本相同,但達到最大值的時刻有所不同,這與顆粒的初始水平位置相關,而與其深度無關;各區域內顆粒水平和豎直方向運動速度較大,且變化趨勢基本同步,而橫向速度相對較小;堆積體分區明顯,堆積體中的顆粒按初始水平位置依次排布。

d. 碎屑流沿運動方向高速運動時,橫向速度很小;在碎屑流運動初期和堆積末期,水平和豎直方向運動速度較小,橫向速度較大,碎屑流橫向運動更為明顯,擴大了災害橫向影響范圍。

[1] 鄭穎人,陳祖煜,王恭先,等.邊坡與滑坡工程治理[M].北京:人民交通出版社,2007.

[2] 黃潤秋.20世紀以來中國的大型滑坡及其發生機制[J].巖石力學與工程學報,2007,26(3):433-454.(HUANG Runqiu.Large-scale landslides and their sliding mechanisms in China since the 20th century[J].Chinese Journal of Rock Mechanics and Engineering,2007,26(3): 433-433.(in Chinese))

[3] HEWITT K.Catastrophic landslide deposits in the Karakoram Himalaya[J].Science,1988,242:64-67.

[4] 孫其誠,王光謙.顆粒物質力學導論[M].北京:科學出版社,2009.

[5] IVERSON R M,LOGAN M,DENLINGER R P.Granular avalanches across irregular three-dimensional terrain: 2.experimental tests[J].Journal of Geophysical Research: Earth Surface,2004,109(1):337-357.

[6] HUSSIN H Y,CIUREAN R,FRIGERIO S,et al.Assessing the effect of mitigation measures on landslide hazard using 2D numerical runout modeling[J].Landslide Science for a Safer Geoenvironment,2014,2(6):679-684.

[7] 王玉峰,程謙恭,朱圻.汶川地震觸發高速遠程滑坡-碎屑流堆積反粒序特征及機制分析[J].巖石力學與工程學報,2012,31(6): 1089-1106.(WANG Yufeng,CHENG Qiangong,ZHU Qi.Inverse grading analysis of deposit from rock avalanches triggered by Wenchuan Earthquake[J].Chinese Journal of Rock Mechanics and Engineering,2012,31(6): 1089-1106.(in Chinese))

[8] 張明,殷躍平,吳樹仁,等.高速遠程碎屑流運動機理研究發展現狀與展望[J].工程地質學報,2010,18(6):805-817.(ZHANG Ming,YIN Yueping,WU Shuren,et al.Development status and prospects of studies on kinematics of long runout rock avalanches[J].Journal of Engineering Geology,2010,18(6): 805-817.(in Chinese))

[9] BANTON J,VILLARD P,JONGMANS D,et al.Two-dimensional discrete element models of debris avalanches: parameterization and the reproducibility of experimental results[J].Journal of Geophysical Research: Earth Surface,2009,114(4):60-71.

[10] LO Chaiming,LIN Minglang,TANG Chaolung,et al.A kinematic model of the Hsiaolin landslide calibrated to the morphology of the landslide deposit[J].Engineering Geology,2011,123(1):22-39.

[11] HARALD T,WANG Y,CHIOU M C,et al.Flow-obstacle interaction in rapid granular avalanches: DEM simulation and comparison with experiment[J].Granular Matter,2009,11(4):209-220.

[12] 畢鈺璋,付躍升,何思明,等.牛眠溝地震滑坡碎屑化全過程離散元模擬[J].中國地質災害與防治學報,2015,26(3):17-25.(BI Yuzhang,FU Yuesheng,HE Siming,et al.Simulation of the whole process of Niumiangou creek rock avalanche triggered by the earthquake using a distinct element method[J].The Chinese Journal of Geological Hazard and Control,2015,26(3):17-25.(in Chinese))

[13] 徐泳,孫其誠,張凌,等.顆粒離散法研究進展[J].力學進展,2003,33(2):251-260.(XU Yong,SUN Qicheng,ZHANG Ling,et al.Advances in discrete element methods for particulate materials[J].Advances in Mechanics,2003,33 (2): 251-260.(in Chinese))

[14] 王國強,郝萬軍,王繼新.離散單元法及其在EDEM上的實踐[M].西安: 西北工業大學出版社,2010.

[15] 錢家歡,殷宗澤.土工原理與計算[M].2版.北京:中國水利水電出版社,1996.

Numerical simulation of movement characteristics of debris flow based on DEM

ZHAO Chuan, FU Chenghua

(SchoolofEnergyandPowerEngineering,XihuaUniversity,Chengdu610039,China)

A three-dimensional numerical model for debris flow was established with the particle discrete element method (DEM), and the movement characteristics of debris flow were analyzed from the full and local aspects. The results show that the movement process of debris flow can be divided into three stages, including starting acceleration, high-speed sliding, and reduced-speed deposit. When the bottom of the slope is flat, large particles are found on the surface of the deposition body, but small particles are in the middle and low areas. The velocity and energy loss of debris flow is significantly influenced by the contacts and collisions between particles. The maximum speeds of different parts of the debris flow are almost the same, and the time required to reach the speed peak is related to the initially horizontal position of particles, but not related to the deposit depth. Area division is distinct in the deposition body, and the particles of the deposition body are arranged according to their initially horizontal positions. The maximum deposit height is at the slope toe, and the deposit height gradually decreases towards both sides. Significant lateral movement occurs in the debris flow in the early and late stages, which can expand the scope of a disaster.

debris flow； 3-D particle discrete element； numerical experiment； contact collision； deposition body

西華大學流體及動力機械教育部重點實驗室(重點研究基地)開放課題(szjj2015-039)

趙川(1989—),男,碩士研究生,主要從事巖土工程穩定性研究。E-mail:zhaochuanvip@163.com

付成華(1978—),女,副教授,博士,主要從事水利水電工程研究。E-mail:fuchh_xhu@163.com

10.3880/j.issn.1006-7647.2017.02.008

P642.21

：A

：1006-7647(2017)02-0043-05

2016-02-26 編輯：鄭孝宇)