滑坡堆積體反粒序現象的離散元數值分析*

楊小龍 王 剛

(①山地城鎮建設與新技術教育部重點實驗室,重慶 400045,中國)(②重慶大學,土木工程學院,重慶 400045,中國)

0 引 言

關于大型滑坡的反粒序堆積現象,最早是由Heim(1932)在對Goldau滑坡進行調查研究時提出的。所謂反粒序堆積現象,它是指滑坡體在移動過程中逐漸形成大粒徑顆粒在上、小粒徑顆粒在下的反粒序堆積結構的現象。在這種堆積結構中,位于滑坡堆積體上部的硬殼層主要以體積較大的塊石、碎石等粗顆粒為主,細顆粒含量較少; 位于中部的主體層,其粒徑分布范圍極廣,以厘米級以下的細顆粒為主,同時夾帶少量的粗顆粒; 位于底部的基底層平均粒徑最小。在Heim(1932)之后,學者們對國內外眾多的大型滑坡進行了現場調查研究,發現了大量存在反粒序堆積現象的地質證據。1986年,Cruden et al.(1986)在對1903年發生于加拿大亞伯達省的Frank滑坡進行詳細的現場調查研究后發現,Frank滑坡堆積體中普遍存在反粒序堆積現象,Cruden et al.(1986) 還通過滑坡堆積體的典型剖面對顆粒粒徑進行了較為詳細的統計分析,對反粒序堆積結構進行了定量分析。Boultbee et al.(2006)借助WipFrag圖片分析方法,對2002年發生于加拿大不列顛哥倫比亞省的Zymoetz River滑坡堆積體中的反粒序堆積現象進行了定量描述。王玉峰等(2012)通過對文家溝滑坡、牛圈溝滑坡以及謝家店子滑坡這3個大型滑坡進行現場勘察,對這3個大型滑坡堆積體中的反粒序堆積結構進行了量化,并探討了反粒序堆積結構的形成機理。此外,學者們還在大量的大型滑坡堆積中發現了反粒序堆積結構(Hewitt,1988; Luzio et al.,2004;Strom,2004; Evans et al.,2006; 鄭光等,2019)。

反粒序是大型滑坡的一種典型堆積現象,但是到目前為止,對于大型滑坡堆積體中反粒序堆積現象的形成機理仍頗有爭議?？偨Y已有研究成果,反粒序堆積現象的形成原因主要有以下兩種觀點。

一是顆粒遷移。Bagnold(1954)將快速移動過程中呈離散流狀態的滑坡體視為牛頓流體,針對滑坡堆積體的反粒序堆積現象提出了顆粒間離散應力的概念,認為顆粒在移動過程中的碰撞會產生離散應力,該離散應力的大小與顆粒粒徑以及剪切應變速率有關,剪切應變速率越大或顆粒粒徑越大,顆粒越容易向滑坡體上部移動,從而在滑坡堆積體中形成大粒徑顆粒在上、小粒徑顆粒在下的反粒序堆積結構。Savage et al.(1988)結合物理模型試驗從定量化的角度比較合理的解釋了在移動過程中呈密集流狀態的滑坡體的顆粒分選過程,提出了“隨機振動篩分”機制,該機制是一種重力作用下的、與顆粒粒徑大小相關的孔隙填充機制?；麦w內的小顆粒逐漸向下方隨機產生的孔隙內轉移,從而形成反粒序堆積結構。Zhang et al.(2011)通過對東河口大型滑坡進行現場調查,研究了滑坡體在滑動過程中顆粒流動的特征以及顆粒尺寸、密度等沿運動路徑的變化,認為堆積體中形成反粒序堆積結構的原因是由于力學效應(擠壓驅逐機制)和幾何效應(隨機振動篩分機制)的共同作用。

二是動力破碎。Strom(1994)通過對發生在帕米爾高原以及天山山脈的多處大型滑坡事件進行研究,對各滑坡堆積體中呈現的成層性現象進行了分析,認為滑坡堆積體的反粒序堆積現象是滑坡體在運動過程中的動力破碎造成的。此外,Strom還通過對高加索山脈Ardon河岸的3處高速滑坡堆積體中反粒序堆積結構進行研究,進一步證實了動力破碎對滑坡體在滑動過程中發生反粒序堆積現象的影響。Imre et al.(2010)探究了滑坡體在滑動過程中的動力破碎現象,認為滑坡體在移動過程中基本保持了初始滑坡體的層序關系,但由于滑坡體內不同深度處巖體受到的作用力不同,使得不同深度處巖體的動力破碎程度不同,從而使滑坡堆積體形成反粒序結構。申智好等(2021)基于顆粒離散元數值分析方法,通過調控滑體組成顆粒的級配分布實現對巖塊碎屑化效應的表征,為辨識巖塊破碎效應對碎屑流滑坡的運移規律及相關災害防治提供了參考。

在對大型滑坡進行研究的早期,反粒序堆積被認為是貫穿整個滑坡堆積體厚度的現象,但近年來,一些學者在對滑坡堆積體進行更為系統的現場調查和粒徑組成分析后,認為反粒序堆積現象僅發生在滑坡堆積體的上表面附近,而位于滑坡體下方的不同粒徑的顆粒仍主要呈混合狀態(Dunning,2004; Crosta et al.,2007; Dufresne et al.,2016)。本文基于顆粒流程序,探究滑坡體顆粒分形維數、滑坡體積以及加速區長度對滑坡堆積體中各粒徑顆粒層序排列的影響規律,以此來探討大型滑坡堆積體出現反粒序堆積現象的原因。

1 數值模型的建立和驗證

為了探究滑坡體在滑動過程中的顆粒分選過程,林小龍(2019)通過一系列物理模型試驗研究了分形維數對反粒序堆積現象的影響,本文根據其試驗裝置建立的二維數值模型如圖1所示。由圖1可知,整個數值模型由傾斜滑動面、水平滑動面、初始滑坡體以及擋板構成,可將模型劃分為物源區、加速運動區和減速堆積區3部分,坐標原點O設在傾斜滑動面與水平滑動面的交點處。物源區初始滑坡體的寬度為W0,高度為H0,滑動面傾角為θ,加速區長度為L0,滑坡堆積體前緣到坐標原點O的距離為L。

顆粒流程序由顆粒(Ball)與墻體(Wall)這兩種最基本的實體構建模型,交替使用運動定律和力-位移定理來分別不斷更新顆粒的位置和接觸力來模擬散粒體的運動,非常適合模擬滑坡的運動過程(Lin et al.,2015; Feng et al.,2017; Borykov et al.,2019; 陳骎等,2020; 劉倫杰等,2020; 潘清等,2020; 李坤等,2021)。本文采用顆粒流計算程序PFC2D通過模擬圖1所示的滑坡模型來探究滑坡體顆粒分形維數(Sammis et al.,1987; Tyler et al.,1992; Crosta et al.,2007)、滑坡體積以及加速區長度對滑坡堆積體中各粒徑顆粒層序排列的影響規律。

由于是二維模型,構成滑坡堆積體的各顆粒實際上是具有單位厚度的半徑為r的圓盤,因此這里所講的堆積體的體積實際上是堆積體剖面的面積。數值模擬方案如表1所列,包含了實驗尺度的數值模型方案(1～6mm)和現場尺度的數值模擬方案(2～12m),所有數值模擬方案的滑動面傾角均為40°。

表1 數值模擬計算方案

顆粒質量-粒徑分布的分形模型由Tyler et al.(1992)提出,具體計算方法為:

(1)

式中:M(r

(2)

由上式可知,在雙對數坐標下,如果ln(M(r

圖2 不同分形維數模擬試樣數據

由于物理模型試驗(林小龍,2019)中采用的是無黏性顆粒材料,因此,在顆粒流程序中選用線性接觸模型來模擬各單元之間的接觸行為,并通過摩擦和阻尼來耗散數值模型中的能量,線性接觸模型示意圖如圖1所示。由圖可知,線性接觸模型由法向和切向的線性彈簧(kn和ks)以及法向和切向的阻尼器(βn和βs)組成。切向接觸力與法向接觸力間的大小關系遵循庫侖摩擦定律,顆粒間的摩擦系數為μc,顆粒與基底的摩擦系數為μw。各參數的大小通過參數標定得到,具體的標定過程為:通過不斷調整顆粒摩擦系數、基底摩擦系數以及阻尼等參數,使得數值模型中顆粒堆積體的自然休止角與實際相等,如圖3所示。最終得到材料的力學參數如表2所列。

表2 數值模型參數

圖3 顆粒流程序參數標定

為了確保顆粒流數值模型的可靠性,首先對林小龍(2019)的物理模型試驗進行模擬。林小龍(2019)以高速遠程滑坡的分形維數分布特征為依據,通過自行設計的模型試驗裝置,研究了不同分形維數顆粒流堆積體的反粒序特征,結果表明:隨著顆粒流分形維數的增大,細顆粒含量的增加,顆粒流運動過程中不同粒徑顆粒的分選作用減弱,導致堆積體中反粒序的發育程度降低。圖4為本文借助顆粒流程序得到的數值模擬結果與林小龍(2019)做的物理模型試驗結果的對比圖。由圖可知,對于數值模擬與物理模型試驗得到的堆積體,兩者在前緣移動距離、堆積體長度以及堆積體最大高度等形體參數上都很接近。因此,可以證明本文的數值模型的計算結果是可靠的。

圖4 數值模擬結果與物理模型試驗結果對比

2 結果分析

2.1 分形維數的影響

圖5給出了3個大顆粒在滑坡體運動過程中移動示意圖。由圖可知,這3個大顆粒一開始位于初始滑坡體的底部,滑坡體開始滑動后,在t=0.4s時這3個大顆粒仍位于滑坡體的下方,但當t=0.8s時,這3個大顆粒已經移動到了滑坡體的表面,揭示了大顆粒在滑坡體滑動過程中向滑坡體上方移動的過程。

圖5 顆粒移動示意圖

圖6給出了顆粒粒徑分形維數對滑坡體滑動過程中各粒級顆粒平均高度的影響。各粒級顆粒平均高度通過下式進行計算:

(3)

圖6 分形維數D對各粒級顆粒平均高度的影響

式中:hd為某一粒級顆粒距滑動面的平均高度,在本文的數值模型中顆粒粒級從小到大分為d1到d6共6個粒級;hi為某一粒級第i個顆粒距滑動面的高度;n為該粒級顆粒的顆?？倐€數。

由圖6可得到以下結果。(1)D=1.5。當在t=0.2s時,d6粒級的顆粒就已經率先開始分離。當t=0.4s時,d1和d6粒級的顆粒已經分別從滑坡體內部移動到了滑坡體的最下方和最上方。當t=0.6時,d2粒級的顆粒已經完成分離; 同時,d3、d4和d5粒級的顆粒也開始分離。當t=1.0s時,d3、d4和d5粒級的顆粒完成分離,此時,整個滑坡體已經形成反粒序堆積結構。(2)D=3.5。當t<0.8s時,d6粒級的顆粒一直在滑坡體的最上方,而d1～d5粒級的顆粒在0.6s前仍處于混合狀態。當t=0.8s時,d1～d3粒級的顆粒完成分離,移動到了滑坡體的下部區域。當t>0.8s時,d5粒級顆粒移動到了滑坡體的最上方,而d6粒級顆粒卻開始從滑坡體最上方向滑坡體內部移動,最終與d3粒級顆?；旌显谝黄?滑坡體沒能形成反粒序堆積結構。

2.2 滑坡體積的影響

圖7給出了滑坡體積對滑坡體滑動過程中各粒級顆粒平均高度的影響。由圖可知,當滑坡體的顆粒粒徑分形維數D=1.5,加速區長度L0=1000m時,隨著滑坡體積的增加,滑坡堆積體的反粒序堆積現象越來越不明顯。(1)滑坡體積V=44388m3。當t=8.0s時,d1粒級的顆粒就已經率先開始分離,移動到滑坡體的最下方,并且d5和d6粒級顆粒一起移動到了滑坡體的最上方,但d5和d6粒級顆粒之間還未開始分離,d2、d3和d4粒級的顆粒之間也未分離。當t=12s時,d5和d6粒級顆粒之間開始分離,d2、d3和d4粒級的顆粒之間也開始分離。當t=20s時,所有粒級的顆粒完成分離,并且各粒級顆粒按照小粒徑顆粒在下大粒徑顆粒在上的順序排列,已經形成反粒序堆積結構。(2)滑坡體積V=177553m3。當t=12s時,d1和d6粒級的顆粒就已經率先完成分離,分別移動到滑坡體的最下方和最上方。當t=32s時,d2粒級的顆粒完成分離,其平均高度只大于d1粒級顆粒。當t=36s時,d3粒級的顆粒完成分離,其平均高度只大于d1和d2粒級顆粒。但是,直到滑坡體完成整個堆積過程,d4和d5粒級的顆粒仍未分離,滑坡體只完成了部分粒級顆粒的分選過程,沒能形成整個滑坡體的反粒序堆積結構。(3)滑坡體積V=537024m3?；麦w中各粒級顆粒的分離程度顯著降低,直到t=32s時,d1粒級顆粒才完成分離,移動到了整個滑坡體的最下方。但直到滑坡體完成整個堆積過程,滑坡體中d5和d6粒級顆粒之間以及d2、d3和d4粒級顆粒之間仍未分離,因此在堆積體中未能形成明顯的反粒序堆積結構。計算結果表明滑坡體積越大,反粒序現象越不明顯,這與Dufresne et al.(2016)通過現場調查得出的結論一致,即反粒序堆積現象僅發生在滑坡堆積體的上表面附近,而位于滑坡體下方的不同粒徑的顆粒仍主要呈混合狀態。

圖7 滑坡體積V對各粒級顆粒平均高度的影響

2.3 加速區長度的影響

圖8給出了加速區長度L0對滑坡體滑動過程中各粒級顆粒平均高度的影響。由圖可知,當滑坡體的顆粒粒徑分形維數D=1.5、滑動面傾角θ=40°、滑坡體積V=177553m3時,隨著加速區長度L0的增加,滑坡堆積體的反粒序堆積現象越來越明顯。(1)加速區長度L0=1000m。當t=12s時,d1和d6粒級的顆粒就已經率先完成分離,分別移動到滑坡體的最下方和最上方。當t=32s時,d2粒級的顆粒完成分離,其平均高度只大于d1粒級顆粒。當t=36s時,d3粒級的顆粒完成分離,其平均高度只大于d1和d2粒級顆粒。但是,直到滑坡體完成整個堆積過程,d4和d5粒級的顆粒仍未分離。因此,滑坡體積V=177553m3的滑坡體只完成了部分粒級顆粒的顆粒分選,仍有d4和d5粒級的顆粒仍未分離,沒能形成整個滑坡體的反粒序堆積結構。(2)加速區長度L0=2000m。當在t=8.0s時,d6粒級的顆粒就已經率先完成分離,移動到滑坡體的最上方。當t=56s時,所有粒級的顆粒均完成分離,并且各粒級顆粒按照小粒徑顆粒在下大粒徑顆粒在上的順序排列,已經形成典型的反粒序堆積結構。但值得注意的是,d4和d5粒級顆粒的平均高度很接近,表明分離程度并不高。(3)加速區長度L0=4000m。當t=24s時,d1和d6粒級的顆粒就已經率先完成分離,分別移動到滑坡體的最下方和最上方。當t=48s時,所有粒級顆粒完成分離,形成小粒徑顆粒在下、大粒徑顆粒在上的典型反粒序堆積結構。

圖8 加速區長度L對各粒級顆粒平均高度的影響

3 反粒序形成機制分析

上述數值模擬結果表明,顆粒粒徑分形維數、滑坡體積以及加速區長度對滑坡體在滑動過程中的顆粒分選程度均有較大影響。反粒序的形成主要通過篩分作用和剪切作用來完成,篩分作用使小顆粒向滑坡體下方移動,剪切作用使大顆粒向滑坡體上方移動。

由不同粒徑大小顆粒組成的滑坡體在滑動過程中會發生剪切,使滑坡體發生膨脹,顆粒間將不斷形成大大小小的孔隙,位于孔隙上方的小顆粒在篩分作用下不斷地通過孔隙而向下移動,如圖9所示。當滑坡體大顆粒粒徑較大且含量較高(即分形維數較小)時,篩分作用更顯著,反粒序現象也就越明顯。篩分作用還與顆粒受到的上覆壓力大小有關,當顆粒受到的上覆壓力較大時,顆粒間的孔隙形成過程將會受到一定程度的抑制,從而使得反粒序程度減小。此外,這種篩分作用還與加速區長度有關,加速區長度越大,滑坡體就能得到越充分的“伸展”,從而使滑坡體底部顆粒受到的上覆壓力減小,增強顆粒分選程度; 同時加速區長度越大,滑坡體也就有更多的時間來完成反粒序的分選過程。

圖9 顆粒分選過程示意圖

滑坡體在滑動過程中會發生剪切變形,滑坡體內的各顆粒會因為受到剪切作用而發生滾動,當大顆粒受到的剪力達到某一臨界值后,大顆?？梢暂p易地通過滾動越過位于其前方的小顆粒,向滑坡體上方移動,最終在滑坡體內形成大顆粒在上、小顆粒在下地反粒序堆積結構(Dasgupta et al.,2011),其示意圖如圖9所示。這種在剪切作用下的滾動過程也與顆粒受到的上覆壓力大小有關,當顆粒受到的上覆壓力較大時,該過程將會受到一定程度的抑制,從而使得反粒序程度減小。

4 雞尾山滑坡模擬

2009年6月5日14:51,重慶市武隆縣鐵礦鄉雞尾山山體,體積約5×106m3的完整性很高的厚層二疊系茅口組灰巖山體經漫長的蠕滑后,由于山體前部“關鍵塊體”被剪斷而突然發生大規模崩滑破壞,滑坡體在躍下超過50m高的陡坎后,獲得巨大的動能,與基巖發生碰撞并迅速解體,發生高速遠程滑坡,滑坡體沿溝谷向下游運動,在溝道里形成平均厚約30m,縱向長度約2200m的堆積區(許強等,2009,2016),最終造成74人死亡。下面結合該滑坡實例進一步對滑坡堆積體的反粒序堆積現象進行分析。

圖10給出了雞尾山滑坡沿其滑動方向的剖面圖。由圖可知,(1)滑源區滑坡體的海拔最高點位于滑坡體的最左側,海拔約為1500m,最低點海拔約為1257m; (2)位于滑源區右下方有一小部分松散土體,該部分土體在滑坡發生后將被滑源區滑坡體鏟刮而沿滑動面向下移動; (3)滑坡體沿滑動面的水平移動距離約為1716m。

圖10 雞尾山滑坡剖面圖(改自Yin et al.(2011))

由于滑源區滑坡體的主要組成部分為棲霞組灰巖,其物理力學參數如表3所列(Zhang et al.,2019),為了便于計算,數值模型中滑源區滑坡體均采用灰巖。采用單軸壓縮試驗來對數值模型中滑坡體材料的細觀參數進行標定,顆粒間的接觸模型采用平行黏接模型。常用參數標定方法為試錯法,即多次調整平行黏接模型的各參數大小,直到通過單軸壓縮試驗得到的試樣的應力-應變曲線與實際一致。最終得到的細觀參數如表4所列?；磪^滑坡體右下方的松散土體,由于不知道其物理力學參數,所以本文借鑒Zhang et al.(2019)計算時確定的細觀參數。對松散土體采用接觸黏結模型,顆粒的密度ρ=1500kg·m-3,法向剛度kn=1.0×107N·m-1和切向剛度ks=1.0×107N·m-1,法向黏接強度TF=0.4×106MPa,切向黏接強度SF=0.4×106MPa,摩擦系數μc=0.4。

表3 雞尾山滑坡巖體力學參數

表4 顆粒流程序細觀參數

圖11給出了雞尾山發生滑坡前后的對比圖。由圖可知,整個滑坡過程可以分為失穩、鏟刮、堆積3個階段,持續時間為150s。在滑坡體滑動過程中,滑源區巖體下方的松散土體被滑源區巖體鏟刮并最終堆積在滑坡堆積體的最前方。模擬計算結果與實際結果一致。

圖11 雞尾山滑坡數值模型圖

圖12給出了各粒徑小組的顆粒的平均高度在滑坡體滑動過程中的變化曲線,由圖可知得到以下幾點認識。(1)在滑坡體初始模型中,由于顆粒為隨機生成,所以各粒徑小組的顆粒處于混亂狀態,此時各粒徑小組顆粒的相對簡易滑坡面的平均高度hd為:hd1=57.0m,hd2=57.5m,hd3=57.0m,hd4=58.2m,hd5=56.6m,hd6=61.2m?？梢钥吹?d5小組的顆粒位于滑坡最下方; 其次是d1和d3小組,這兩個小組的顆粒的平均高度恰好相等; 位于d1和d3小組上方的小組分別為d2、d4和d6小組,這3個小組按從下往上的順序排列。在初始模型中,各小組顆粒除了d1和d3小組,其余小組顆粒間的平均高度差在0.5m以上,彼此間的分離程序較大。(2)在t=0～18s這段時間,d1小組顆粒的平均高度逐漸減小。當t=18s時,d1小組顆粒的平均高度比其余小組的都小,表明d1小組的顆粒已經在顆粒的分選過程中移動到了滑坡體的最下方,并且在此后的滑動過程中,d1小組的顆粒一直位于滑坡體的最下方。(3)在t=0～50s這段時間,d2～d5 這4個小組顆粒的平均高度差值逐漸減小。這是因為滑坡體在越過陡坎時會發生破碎,使各粒徑小組顆粒間的相對位置發生變化,使其處于更加混亂的狀態,導致各小組顆粒間的平均高度差值減小。(4)在t=50～58s這段時間,d2和d3兩個小組的顆粒在顆粒分選過程中向滑坡體下方滑動,與d1小組顆粒按照小粒徑顆粒在下、大粒徑顆粒在上的反粒序堆積結構,且在此后的滑動過程中一直維持著反粒序堆積結構。(5)在t=58～88s這段時間,d4、d5和d6這3個小組顆粒繼續在動力篩分作用和剪切作用的共同工作用下進行著顆粒分選過程,當t=88s時,d5小組的顆粒完成顆粒分選,與d4和d6小組的顆粒完成分離。但d4和d6小組的顆粒仍沒有完成分離,以混合的狀態位于滑坡體的最上方。(6)在t=88～150s這段時間,滑坡體逐漸堆積,在t=150s時完成堆積過程,d1、d2、d3和d5這4個小組的層序關系沒有發生改變,一直維持著反粒序堆積結構。但是,直到滑坡體完成堆積,d4和d6小組的顆粒仍未完成分離。整個滑坡堆積體沒能形成嚴格意義上的反粒序堆積結構,只形成了部分反粒序堆積結構。

圖12 各粒徑小組顆粒平均高度變化曲線

5 結 論

本文根據已有的物理模型試驗,在顆粒流程序中建立數值模型,研究了滑坡體中顆粒粒徑的分形維數、滑坡體積以及加速區長度對滑坡體在滑動過程中的顆粒分選程度的影響,以探究滑坡堆積體出現反粒序堆積現象的原因?？梢缘贸鲆韵鲁醪浇Y論:

(1)滑坡體堆積體的反粒序堆積現象是動力破碎以及顆粒分選的結果,顆粒分選是在動力篩分作用和剪切作用的共同作用下完成的?；麦w在滑動過程中因發生顆粒破碎而形成不同粒徑大小的顆粒,動力篩分作用使小粒徑顆粒通過大粒徑顆粒間的孔隙向滑坡體下方移動,剪切作用使大顆粒越過小顆粒向上移動,最終使不同粒徑的顆粒完成顆粒分選過程,形成反粒序堆積結構。

(2)顆粒粒徑分形維數、滑坡體積以及加速區長度均對滑坡體內的顆粒分選程度有較大影響,顆粒分選程度隨分形維數的減小而增強,隨滑坡體積的增大而減弱,隨加速區長度的增加而增強。

(3)數值模擬再現了雞尾山滑坡堆積體中存在的反粒序堆積現象?；磪^巖體由于碰撞以及受到的拉力超過其抗拉強度等原因而發生破碎,形成粒徑大小不一的顆粒,并在滑動過程中完成顆粒分選,形成大粒徑顆粒在上、小粒徑顆粒在下的反粒序堆積結構。