張 晗高 楊李 濱李 軍吳偉樂

1.長安大學地質工程與測繪學院,陜西 西安 710054;

2.中國地質科學院地質力學研究所,北京 100081;

3.自然資源部活動構造與地質安全重點實驗室,北京 100081

0 引言

全球氣候變化造成了極端天氣事件的增多,進而導致特大型地質災害頻發(fā)。強降雨條件已成為誘發(fā)滑坡及復合型地質災害的最主要原因之一(孫廣忠,1988;殷躍平,1998;李濱等,2016;高楊等,2017;閆金凱等,2020)。復合型地質災害是指包含兩種或兩種以上基本地質災害類型,并可引發(fā)堰塞壩和涌浪等次生效應的鏈式地質災害(Varnes,1978;Yin et al.,2009;殷躍平等,2010,2021;Hungr et al.,2014)。在中國西南山區(qū),高位崩滑體失穩(wěn)啟動后,經過勢動能轉化不斷碰撞解體形成碎屑流,往往在水動力作用下極易向泥石流或水石流轉化,形成復合型地質災害,對人類的生命和財產安全造成巨大的威脅。針對水動力作用造成滑坡遠程運動理論,國內外學者開展了大量研究?；滤樾剂髟谒畡恿l件下向泥石流的轉換現象已經被大量地質工作者通過野外調查和試驗證明(Plafker and Ericksen,1978;Voight et al.,1983;崔鵬,1991;崔鵬和關君蔚,1993;Iverson,1997;Takarada et al.,1999;余國安,2022)。固體和流體不同的相互作用不僅可以劃分滑坡或泥石流的類型,還使這類災害具有獨特的破壞性。碎屑流-泥石流復合型滑坡的流變試驗和數值研究表明不同液體體積分數、液體性質、固液速度差和固體顆粒尺寸是影響固液相互作用,從而影響滑坡運動距離的關鍵因素(Gao et al.,2022a)。流體存在會增加孔隙水壓力,有效降低顆粒間摩擦并減少能量損耗,提高了滑坡的運動特性。對比流體拖曳力,具有黏性和不可壓縮性的液體存在使得運動滑體的邊界層保有較高的孔隙水壓力,在高速剪切、振動作用下形成超孔隙水壓力,導致滑體摩擦阻力下降,增加滑坡運動距離(Hutchinson and Bhandari,1971;Evans et al.,2001;Sassa et al.,2004;殷躍平等,2010;Iverson and George,2016)。Strom and Abdrakhmatov(2018)在進行大量滑坡調查后發(fā)現,滑坡高速運動造成的飽水下墊層液化是其側向擴離運動及覆蓋范圍增大的重要原因。碎屑流向泥石流轉化形成的復合型滑坡在世界范圍內多次發(fā)生,給人類生命財產安全帶來巨大威脅,例如:1970年 秘 魯Yugay滑 坡(Plafker and Ericksen,1978),2002年 俄 羅 斯Kolka滑 坡(Evans et al.,2009),2006年菲律賓Leyte滑坡(Evans et al.,2007),2008年四川文家溝滑坡(余斌等,2010),2014年重慶咸池水庫滑坡(衛(wèi)童瑤,2021;Gao et al.,2022a;Li et al.,2022),2019年貴州雞場滑坡(Gao et al.,2020;李壯等,2020;高浩源,2021)。

為開展高位滑坡后破壞成災范圍和成災風險的定量評估,進行科學避險,數值模擬技術成為反演分析滑坡運動過程和動力學機制的最為高效的手段之一。目前研究兩相復合型滑坡的數值模擬方法包括歐拉-歐拉模型(Reynolds,1976;Iverson and Denlinger,2001;Pitman and Le,2005;Pailha and Pouliquen,2009;George and Iverson,2011;Bouchut et al.,2015,2016)、拉格朗日-拉格朗日模型(DEM-SPH等;Tan and Chen,2017;徐文 杰,2020)、歐 拉-拉 格 朗 日 模 型(SPHDCDEM、DEM-CFD等;Shan and Zhao,2014;Jing et al.,2019)。歐拉-拉格朗日模型較其他模型能夠更好地還原固液兩相流的運動過程,但由于現有數值模型沒有充分考慮拖曳效應及技術限制,還未實現液體和固體顆粒的共同運動。

2014年9月2日晚,重慶市奉節(jié)縣竹園鎮(zhèn)無山村所在山體發(fā)生高位滑坡(以下簡稱無山坪滑坡),滑坡體內房屋全部倒塌,共1504間房屋被毀壞。由于及時的降雨警報,無山村所有居民在山體滑坡發(fā)生前20分鐘被疏散?；卵貛r層層面產生整體位移,前緣滑體大面積滑塌,滑塌的巖土體混和暴雨形成泥石流,沿溝道沖入并堵塞岔河,屬典型的滑坡轉化為泥石流的復合型滑坡。文中采用基于兩相耦合模型和SPH算法技術結合的滑坡后破壞數值模擬平臺(LPF3D),針對無山坪滑坡,采用現場調查、遙感影像對比和數值模擬等方法對該滑坡的水力學特征進行了初步探討,分析了滑坡基本情況和成災機理,模擬并重現了滑坡運動堆積過程,為類似固液耦合復合型滑坡的預測和反演提供幫助。

1 研究區(qū)地質環(huán)境條件

研究區(qū)位于重慶市東北部奉節(jié)縣竹園鎮(zhèn),距離奉節(jié)主縣城約80 km,中心地理坐標為109°01′17″E,30°29′19″N。研究區(qū)滑坡前后緣公路交通條件較好,與鄉(xiāng)鎮(zhèn)、奉節(jié)縣城相連?；掳l(fā)生后公路損壞嚴重。地形地貌屬構造剝蝕斜坡、中低山地貌,地形總體呈南高北低?；潞缶壸罡唿c海拔為1330 m,前緣最低點位于岔河河谷,海拔約800 m,相對高差550 m;地形坡度一般為15°～35°,多陡坎。山體表面發(fā)育多條沖溝,滑體失穩(wěn)下滑后,沿運動路徑中的沖溝滑動,并在前緣沿東西向岔河內堆積。

研究區(qū)位于新華夏系第三沉降帶四川盆地東端,地處川東坳褶帶、大巴山南緣弧形褶皺帶及鄂湘黔隆褶帶接合部。地質構造以褶皺為主,斷裂欠發(fā)育?；绿幱诩t巖向斜軸部近西端(圖1)。由于早期的構造作用和風化作用,部分巖體破碎,有利于雨水的滲透和地下水的儲存。

圖1 研究區(qū)地質構造圖Fig.1 Geological structural map of the study area

滑坡區(qū)域覆蓋第四系松散堆積物,主要由黏性土、砂泥巖碎塊石和少量生活及建筑垃圾組成。基巖為侏羅系中統(tǒng)新田溝組(J2x)的砂巖和泥巖,上部以砂巖為主,中部砂、泥巖互層或混雜沉積成巖,下部砂巖較泥巖發(fā)育。從無山坪滑坡的后緣到前緣,地層的傾角逐漸減小到零;在滑坡的前端甚至出現倒置的地層。該區(qū)砂巖強度較高,容易沿著巖層的傾角發(fā)生緩慢變形,從而在斜坡表面逐漸形成許多卸荷裂縫,這些裂縫加速了地表水和雨水的滲入(殷躍平和胡瑞林,2004;劉新榮等,2008;王志儉等,2008;劉小紅等,2015)。砂巖和泥巖之間的薄弱界面和廣泛拉張裂縫的發(fā)展有利于雨水的滲透,促進了滑坡形成(Li et al.,2022)。

研究區(qū)多年年平均降雨量由川鄂交界處的2000～2100 mm及北部邊緣的1500～1600 mm向長江河谷降低為1100 mm左右,其中3—8月降雨量占全年的68%以上。滑坡發(fā)生一周內,重慶市奉節(jié)縣竹園鎮(zhèn)遭到百年一遇最大集中降雨,奉節(jié)縣氣象局數據顯示奉節(jié)縣竹園鎮(zhèn)累計降雨量可達266 mm。短時強降雨是包括無山坪滑坡在內的許多地質災害的主要誘發(fā)因素。

2 無山坪滑坡運動堆積特征

2014年8月31至9月1日,重慶市奉節(jié)縣竹園鎮(zhèn)累計降雨量達到266 mm,數小時內誘使無山坪滑坡發(fā)生整體位移,并出現大面積滑塌,滑塌的巖土體同地表徑流混和形成泥石流,沿溝槽沖入岔河,堵塞岔河形成兩處堰塞湖。無山坪滑坡目前仍殘余大量潛在不穩(wěn)定體,平臺前部與右側共發(fā)育有五條溝道,坡體內發(fā)育有大量拉張裂隙,若遇強降雨,滑坡極有可能發(fā)生再次下滑,威脅滑坡下方岔河兩岸居民及基礎設施安全。

通過現場調查、衛(wèi)星和無人機影像分析(圖2),滑坡失穩(wěn)后發(fā)生解體,沿N20°W和N8°E兩個方向高位剪出,約4×106m3滑體殘留在滑源區(qū),2.83×106m3滑體向下游流態(tài)化運動堆積。強降雨在研究區(qū)內造成了豐富的地表徑流,滑坡失穩(wěn)啟動后,滑體同豐富的地表徑流混合形成泥石流復合災害,泥石流沿著無山坪平臺前部的五條溝道向下游流動,匯入坡腳處東西向的主河道。根據無山坪滑坡的堆積特征,從總體上將該滑坡-泥石流分為滑源區(qū)、流通區(qū)和堆積區(qū)(圖3)。

圖2 無山坪滑坡滑前滑后遙感對比圖Fig.2 Image comparison before and after the Wushanping landslide

(1)滑源區(qū)

滑源區(qū)的平面形態(tài)呈“舌型”,順坡向滑坡體總體長約500 m,前源寬約820 m,后緣寬約360 m,平均厚度約19 m。平面面積約39.2×104m2,體積約6.83×106m3?；逻\動過程中,滑體沿與巖層之間的滑面旋轉和滑動,導致后緣垂直錯位超過30 m,滑后可見砂、泥巖互層的基巖出露,與水平夾角約30°(圖4a)?；磪^(qū)中部地形較高,在滑動過程中形成了許多長10～20 m的裂縫,伴隨著不均勻的下沉。隨處可見無山村被毀道路與房屋(圖4b)。

(2)流通區(qū)

流通區(qū)主要位于高程1000 m以下的溝道處,斜坡基巖主要以第四系堆積物和砂泥巖為主。滑源區(qū)巖體失穩(wěn)后,在巨大的勢能和雨水拖曳的作用下,快速順溝道向下運動(溝道位置見圖3a)。流通區(qū)順坡向長約400 m,寬約700 m,面積約28×104m2。滑坡剪出口位于一個陡峭斜坡的自由面,坡度為40°～50°,距離坡腳的高度為150～200 m(圖4c、4d)。

圖3 無山坪滑坡平剖面圖Fig.3 Profile and plan of the Wushanping landslide

(3)堆積區(qū)

滑前在堆積區(qū)有一條梅溪河的支流岔河流過山谷,方向為160°。由于暴雨的緣故,滑體在陡坡上運動解體轉化為泥石流,大量流固混合物沿溝道沖向河流。最終形成順滑向670 m、順河長250 m、平均厚度14 m、體積約為2.3×106m3的堰塞體,將岔河堵塞。此外,大量滑坡碎屑繼續(xù)向下流動約1000 m。泥石流和堰塞湖淹沒了原河谷兩邊的一些房屋(圖4e、4f)。

圖4 滑源區(qū)、流通區(qū)及堆積區(qū)現場調查照片Fig.4 Site photos of the slide source area, propagation area and accumulation area

3 固液耦合計算方法

針對此類兩相流復合型滑坡,文中基于SPH方法和兩相耦合動力學理論研發(fā)的LPF3D(Landslides post failure 3D)數值計算方法,能夠處理具有自由表面、變形邊界、運動交界以及極大變形的問題,具有計算速度快和計算精度高的適用性。流體采用基于連續(xù)介質的SPH方法計算,SPH粒子承載著質量、密度、速度、加速度等物理量,粒子間通過核函數相互作用,適用于大變形情況。固體顆粒體基于顆粒動力學模型將傳統(tǒng)SPH方法改造成適用于離散顆粒的SPH方法,采用理想彈塑性本構進行求解,使得SPH粒子不僅承載質量、密度、速度、加速度等各種物理量,同時承載著顆粒的粒徑、體積分數等顆粒屬性(Cui et al.,2021;Chen and Yan,2021)。

3.1 控制方程

不同狀態(tài)下滑體材料的控制方程由質量和動量守恒原理描述,見公式(1)和(2)。

連續(xù)方程:

運動方程:

公式中:α和β分別是x、y兩個方向上的分量;ρ為材料密度;vα為速度在x方向上的分量;xα和xβ分別為x在x和y兩個方向上分量;t為時間;f代表其他外力(如重力、相間作用力);d/dt為方程全導數?？倧埩喀姚力峦ǔ７譃楦飨蛲詨毫和剪應力τ兩部分。表示為:

其中δαβ是克羅內克函數,當α=β時,δαβ=1;當α≠β時,δαβ=0。

3.2 計算本構

3.2.1 流體本構

該本構模型以計算純水流體和濃縮流體為主,流體密度ρ與壓強P關系為:

其中P0=ρ0c20/γ,為初始壓強;ρ0為流體初始密度;γ與流體的可壓縮性相關,γ=7;c0為初始聲速,為保證流體可壓縮性,一般取c0=(10～40)vmax,vmax為流體最大速度。

對于牛頓流體,剪應力ταβ(公式5)與剪應變率εαβ(公式6)成比例,比例系數為黏性系數μ;對于非牛頓流體,黏性系數μ為剪應變率函數。

其中·v為速度的散度,其他參數同公式(1)與公式(3)。

3.2.2 固體本構

靜壓力Ps是使用平均應力的標準定義從顆粒構成方程中直接計算的,表示為:

其中σxx、σyy和σzz是應力張量σγγ在x、y和z方向上的分量。當顆粒處于準靜態(tài)狀態(tài)時,顆粒整體幾乎不產生大變形,主要以彈性變形為主,在此階段,根據線彈性模型(胡克定律)計算顆粒的應力-應變關系。

其中σ·αβ為增量 形 式 的 應 力 分 量,G為剪 切 模量,K為彈性體積模量,E為彈性模量(楊氏模量),υ為泊松比,e·αβ為偏剪切應變率張量,ε·αβ為應變率張量。

ταβ為偏差應力分量;為第二不變量;為應變率張量,其定義與摩擦系數μ(I)相對應,這也意味著體積分數對慣性常數的單調性相關。

此處, 實驗和數值模擬表明在函數μ(I)中,在非常低的慣性常數I(準靜態(tài))下,最小值μp逐漸增加到I增加時的有限值μ2。I0是一個常數,公式(4)中的參數取決于材料特性。例如,典型值時

μp=tan(21°),μ2=tan(33°),I0=0.28(MiDi, 2004)。

慣性常數表示慣性時間尺度(d2ρp/Ps)0.5與宏觀變形時間尺度之間的比率;d表示顆粒粒徑;ρp表示顆粒密度。

3.3 相間作用

大量研究表明,兩相流運動中顆粒體和流體的相間作用對滑體運動范圍至關重要(Davies,1990;Iverson et al.,2010;Tayyebi et al.,2022;Gao et al.,2022b)。計算固液相間拖曳力時,采用Gidaspow Gidaspow,提出的公式(1994),即對于密相的計算采用Ergun方程以及對于稀相的計算采用Wen-Yu方程(Ergun,1952;Wen and Yu,1966):

其中β為流體與固體間的動量傳遞系數;CD為曳力系數;φP為固體內摩擦角;μf為流體黏度;dP為固體顆粒粒徑;ρf為流體密度;vf、vP分別代表流體和固體的速度;αf為流體體積分數。

曳力系數CD為:

相對雷諾數ReP定義為:

為消除兩個方程間的不連續(xù)性,引入松弛因子φfP對過渡區(qū)域中的動量交換系數進行光滑,

因此,動量交換系數β可以表示為:

可得作用于單位質量顆粒上的曳力R′fP為:

其中參數含義同公式(16)。

3.4 邊界條件

3.4.1 法向邊界力

SPH作為無網格方法,不能像網格方法那樣直接將界面力施加在邊界點上。邊界力求解采用強洪夫等(2011)提出的罰函數方法求解顆粒與邊界之間的法向邊界力。

3.4.2 切向邊界力

滑體與基底的切向力根據滑體材料性質分別選擇,流體和顆粒體分別采用庫倫摩擦模型和層流黏滯模型(Hungr, 1995)。

固體-庫侖摩擦模型:

公式中:T為基底剪切阻力;γ為重度;Hi為流體厚度;θ為運動路徑坡角;ac=v2/R為離心加速度,大小取決于運動路徑的曲率半徑R,v為運動速度;g是重力加速度,取9.8 m·s-2;ru為孔隙水壓力系數;φ為摩擦角。

流體-層流模型:

公式中:T為基底剪切阻力,A為面積,v為流體速度,μ為動態(tài)黏滯系數,H為流體深度。

4 數值模擬及結果分析

4.1 模型建立及參數選取

為定量描述研究區(qū)地形地貌特征,采用1∶2000的三維等高線數據建立研究區(qū)地質模型,在滑源區(qū)建立滑體與流體模型。此次模擬不考慮前期滑體失穩(wěn)過程,假設滑坡由降雨持續(xù)影響并突然爆發(fā),目的是為了模擬流體與巖土體顆粒的運動堆積過程。

滑體材料物理力學參數采用實際參數。滑坡巖性主要為侏羅系中風化砂巖,巖石的物理力學參數為ρ=2240 kg/m3,粒徑為0.1 m。流體采用泥漿參數,ρ=1200 kg/m3,黏滯系數為0.2 Pa·s。為探討水動力對固液兩相流滑坡后破壞過程的影響,此次模擬共設計四種不同工況(表1)。①工況Ⅰ:純固體顆粒下滑,不考慮孔隙水壓力和相間作用力。②工況Ⅱ:考慮孔隙水壓力對下滑運動過程中固體顆粒的影響。③工況Ⅲ:考慮運動過程中流體對固體顆粒持續(xù)拖曳的影響。④工況Ⅳ:同時考慮流體孔隙水壓力與拖曳效應的影響。

表1 無山坪滑坡LPF模擬參數Table 1 LPF simulation parameters of the wushanping landslide

4.2 結果分析

分別取顆粒運動停止時刻為計算終止時間。不同工況下無山坪滑坡模擬結果如圖5所示。工況Ⅰ:純固體顆粒下滑,大部分運動顆粒停留在溝道內或溝道出口處,均未抵達下方河道,最大運動速度為23 m/s,最遠運動距離為900 m。工況Ⅱ:滑體底部孔隙水壓力存在降低了摩擦,使固體顆粒運動距離較工況Ⅰ更遠,大部分顆粒停留在溝槽出口位置,最大運動速度為26 m/s,最遠運動距離為1150 m。工況Ⅲ:僅考慮流體拖曳作用時,固體顆粒與流體耦合向下流動,在各溝道出口處均形成了不同程度的沖積扇,具有明顯的流態(tài)化特征,最大運動速度為31.5 m/s,最遠運動距離為1250 m。工況Ⅳ:孔隙水壓力減阻和流體拖曳力增程共同作用時,滑坡的運動距離更遠,致災范圍更廣,與無山坪滑坡災后實際堆積距離、堆積范圍最為接近,最大運動速度為34 m/s,最遠運動距離為1300 m。

圖5 四種工況下的堆積結果圖Fig.5 Accumulations of the Wupingshan landslide under four working conditions

通過對固液耦合作用的對比,結合堆積模擬結果,可以發(fā)現工況Ⅳ得到的計算結果與實際最為接近。因此選用工況Ⅳ代表無山坪滑坡泥石流后破壞過程,并對其運動堆積特性進行論述。

文中選取0～40 s和180 s六個關鍵節(jié)點來分析無山坪滑坡破壞后的運動過程。工況Ⅳ下無山坪滑坡泥石流隨時間變化的運動過程如圖6所示,棕色顆粒表示固體顆粒,藍色顆粒表示流體顆粒。整個運動過程中固液整體的運動速度如圖7所示,在180 s時,顆粒的平均速度已無法影響滑坡的堆積狀態(tài),故180 s時運動已停止。從圖中可以看出,在0～10 s之間滑坡高位啟動,固體顆粒與流體混合順多處溝道快速向下流動。10～20 s,滑坡前緣抵達下方處,大約在16 s左右,速度達到最高值34 m/s。20～40 s,由于泥漿的拖曳作用,后部固體顆粒沿著溝道繼續(xù)運動不斷匯入堆積區(qū),堆積厚度不斷增大。在40 s左右,大部分顆粒運動到河道位置,并在此處淤積,平均速度明顯下降。在180 s左右,整體顆粒接近于靜止運動,計算停止。

圖6 工況Ⅳ下運動過程圖Fig.6 Diagrams showing the fluid-solid coupled movement of the Wushanping landslide under working condition Ⅳ

圖7 工況Ⅳ下運動速度圖Fig.7 Velocity diagrams of the Wupingshan landslide under working condition Ⅳ

滑坡運動過程中固體堆積厚度變化等值線圖如圖8所示,圖中紅色線代表真實滑坡,紅黃藍三種色系體現固體顆粒不同堆積厚度變化情況?；逻\動全過程,巖土體材料混合流體順多條溝道向下流動,逐漸在各溝道處堆積。顆粒停止運動時,滑坡前緣在河道處呈流態(tài)化堆積,最大運動距離達1300 m,最大堆積厚度約21.5 m,堆積形態(tài)與現場調查結果較為一致。

圖8 工況Ⅳ下滑坡堆積厚度圖Fig.8 Diagrams showing the deposition thickness with time of the Wupingshan landslide under working condition Ⅳ

流體與固體顆粒間的速度差是決定拖曳力大小的關鍵因素,為進一步了解無山坪滑坡運動過程中拖曳力及孔隙水對固體顆粒運動的影響,選?、谔枩系捞巸蓚€速度對比曲線圖作為對比。四種工況下滑坡前緣固體顆粒的速度對比曲線如圖9a所示,水和相間作用力的影響下,固體顆粒的運動速度有明顯的增大,且孔隙水作用大于相間作用。同時有流體參與運動時,滑坡運動時間更長,運動距離更遠。工況Ⅲ、工況Ⅳ下滑坡前緣固體顆粒及其周圍流體的速度對比曲線如圖9b所示,滑動和加速階段,工況Ⅲ和工況Ⅳ的流體速度相同,但固體顆粒速度有明顯不同,故滑坡啟動和加速主要受孔隙水壓力的影響?；聹p速階段,兩種工況下固體顆粒速度隨流體速度減小,且變化趨勢基本一致,故相間作用力在滑坡減速階段起主要作用。

圖9 無山坪滑坡速度曲線圖Fig.9 Velocity change of the front-edge granules under different working conditions

5 結論

文中用自主研發(fā)的LPF3D軟件,以重慶奉節(jié)無山坪滑坡為例,模擬再現了無山坪滑坡運動堆積過程,探討分析了遠程滑坡后破壞運動過程中的水動力作用,認為流體會在孔隙水壓力減阻和拖曳力增程兩方面影響高位滑坡的遠程運動,取得了以下認識。

(1)通過四種工況對無山坪滑坡泥石流的模擬結果分析,孔隙水壓力和拖曳力的共同作用會使滑坡的運動速度、運動距離增大。表明水動力在滑坡運動過程中的動力學效應主要分為液化效應和拖曳效應兩種,兩種效應的共同作用會使碎屑流轉化為泥石流,導致遠程成災。

(2)基于SPH的多相耦合計算模型,將流體狀態(tài)方程、彈塑性本構方程和相間作用力進行搭建結合,模擬結果較好地再現了無山坪滑坡的運動堆積。表明基于多相耦合理論的LPF3D計算方法能夠更真實地還原滑坡后破壞過程,對復合型滑坡的模擬具有較好適用性。該方法對于計算自由表面流、大變形和損傷破壞問題具有較大優(yōu)勢,避免了網格法在計算時遇到的界面追蹤困難、網格扭曲和纏繞、計算量大等問題。

(3)模擬結果顯示:無山坪滑坡全程最大運動速度為34 m/s,最大堆積厚度為21.5 m,最大堆積面積為0.12 km2,最遠運動距離為1300 m。在前期強降雨影響下,無山坪滑坡滑體材料充分飽水,強度降低。同時降雨入滲補給使滑面處孔隙水無法即時排出,導致孔隙水壓力增加,摩擦系數降低,滑坡高速啟動,并在運動過程中與流體混合形成泥石流。流固兩相間拖曳力的存在進一步促進滑坡顆粒流化運動,提高了其運動速度,使得滑坡遠程成災。