基于DEM的碎屑流運動特性數(shù)值模擬

趙 川,付成華

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

基于DEM的碎屑流運動特性數(shù)值模擬

趙 川,付成華

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

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

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

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

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

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

1 碎屑流數(shù)值模型建立及參數(shù)選取

1.1 離散元模型構(gòu)建

碎屑流數(shù)值模型主要由碎屑堆積槽、滑動面和堆積面組成?；瑒用婧投逊e面構(gòu)成一個邊坡系統(tǒng)，坡高為10 m,坡比為1∶2,堆積面長為10 m;碎屑堆積槽用于放置顆粒物質(zhì),槽外形為三棱柱體,尺寸為1.5 m×3 m×5 m,如圖1所示。

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

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

圖2 數(shù)值計算顆粒單元外形

圖3 碎屑流離散元模型

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

(1)

(2)

Ft=-Stδt

(3)

(4)

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

1.2 數(shù)值計算參數(shù)標(biāo)定

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

圖4 離散元顆粒漏斗數(shù)值模型試驗

表1 DEM模型相關(guān)參數(shù)取值

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

2 數(shù)值試驗結(jié)果及分析

2.1 整體運動特性

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

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

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

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

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

圖7 顆粒單元接觸數(shù)量的變化過程

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

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

2.2 局部運動特性分析

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

圖9 顆粒分區(qū)

圖10 不同分區(qū)顆粒的運動速度

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

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

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

圖12 碎屑流最終堆積狀態(tài)

3 結(jié) 論

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

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

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

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

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

[2] 黃潤秋.20世紀(jì)以來中國的大型滑坡及其發(fā)生機(jī)制[J].巖石力學(xué)與工程學(xué)報,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] 孫其誠,王光謙.顆粒物質(zhì)力學(xué)導(dǎo)論[M].北京:科學(xué)出版社,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] 王玉峰,程謙恭,朱圻.汶川地震觸發(fā)高速遠(yuǎn)程滑坡-碎屑流堆積反粒序特征及機(jī)制分析[J].巖石力學(xué)與工程學(xué)報,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] 張明,殷躍平,吳樹仁,等.高速遠(yuǎn)程碎屑流運動機(jī)理研究發(fā)展現(xiàn)狀與展望[J].工程地質(zhì)學(xué)報,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].中國地質(zhì)災(zāi)害與防治學(xué)報,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ìn)展[J].力學(xué)進(jìn)展,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].西安: 西北工業(yè)大學(xué)出版社,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

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

趙川(1989—),男,碩士研究生,主要從事巖土工程穩(wěn)定性研究。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 編輯：鄭孝宇)