基于SPH深度積分模型的滑坡數(shù)值分析

陶言祺，李家鋼

(1.大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室, 遼寧 大連 116024;2.中海油研究總院, 北京 100081)

基于SPH深度積分模型的滑坡數(shù)值分析

陶言祺1，李家鋼2

(1.大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室, 遼寧 大連 116024;2.中海油研究總院, 北京 100081)

滑坡會(huì)對(duì)工程結(jié)構(gòu)物造成嚴(yán)重破壞,滑移速度和距離是量化和分析滑坡的兩個(gè)重要參數(shù)。針對(duì)滑移速度和距離預(yù)測(cè)一直比較復(fù)雜的現(xiàn)狀,建立考慮土體固結(jié)和侵蝕效應(yīng)的控制方程,選用摩擦流變模型,采用SPH深度積分算法,對(duì)滑坡進(jìn)行了模擬研究。對(duì)比不同坡度、接觸摩擦系數(shù)和侵蝕率條件下的滑移體速度、長度的時(shí)程曲線,分析不同條件下滑移體的最終沉積形態(tài),整理了最大滑移距離和速度,討論其變化規(guī)律。研究成果可為滑坡災(zāi)害預(yù)警提供技術(shù)參考。

滑移速度；SPH深度積分法；摩擦流變模型；侵蝕效應(yīng)

滑坡是一種常見的地質(zhì)災(zāi)害,廣泛分布在山區(qū)、河谷地帶,危害很大。國內(nèi)外對(duì)于滑坡的研究主要集中在誘發(fā)機(jī)制、失穩(wěn)分析、災(zāi)害評(píng)估等方面。由于地震、暴雨、人類活動(dòng)影響等造成的滑坡每年都造成巨大的經(jīng)濟(jì)損失和人員傷亡。因此,研究滑坡對(duì)于工程設(shè)計(jì)和防災(zāi)減災(zāi)具有重要意義。

針對(duì)滑坡[1]問題,國內(nèi)外學(xué)者多采用SPH法[2]、離散單元法[3]、無網(wǎng)格伽遼金法[4]、細(xì)胞自動(dòng)機(jī)法[5]、非連續(xù)變形分析方法[6]等方法進(jìn)行研究,并取得了一定研究成果,但是針對(duì)滑移速度和距離的計(jì)算一直是比較復(fù)雜的問題[7]。Pastor等[8]采用SPH深度積分法,考慮土骨架與孔隙流體的動(dòng)力耦合和侵蝕效應(yīng),建立了耦合孔隙水壓力的深度積分模型。依據(jù)Biot-Zienkiewicz方程的u-p格式進(jìn)行積分,將流動(dòng)-固結(jié)方程簡(jiǎn)化成耦合孔壓的深度積分方程。二維深度積分模型能夠簡(jiǎn)便、準(zhǔn)確地計(jì)算任意時(shí)刻滑移體的滑移速度、高度和滑移距離等。Pastor等用SPH深度積分法對(duì)山體滑坡、泥石流等問題進(jìn)行了一系列研究和驗(yàn)證,對(duì)香港青山滑坡和加拿大Frank滑坡進(jìn)行了模擬[9-10],結(jié)果顯示SPH深度積分法能夠較為精確地模擬陸地滑坡問題,同時(shí)計(jì)算效率高。

綜上所述,本文基于二維SPH深度積分法,建立了滑坡模型,并考慮土骨架和孔隙流體的動(dòng)力耦合,主要分析了坡度、滑移體與坡床間的接觸摩擦系數(shù)以及侵蝕率對(duì)滑移距離和速度的影響。

1 基本理論

滑坡過程中需要考慮滑移體的自重作用、摩擦作用,且被侵蝕、挾帶的土體的慣性阻力會(huì)增加滑移體底部的剪切阻力[11]。

1.1 SPH深度積分模型

SPH深度積分法通過以下公式得到離散化的深度積分方程。函數(shù)的近似逼近：

(1)

式中：φ(xj)是空間變量xj的函數(shù)；插值核函數(shù)Wij=W(xj-xi，hsml)，常采用三次樣條曲線插值格式[12]，hsml為光滑長度；mj為質(zhì)量；ρj為密度。

深度積分公式：

(2)

式中：φi(x3)是沿深度方向變量x3的函數(shù)，[z,z+h]為x3的積分區(qū)間。

在控制方程的簡(jiǎn)化過程中,假定流體各項(xiàng)同性、孔隙流體不可壓縮,在土骨架中的流動(dòng)滿足達(dá)西定律,忽略孔隙水相對(duì)于土骨架的加速度;滑移體為飽和土體,其自由面為排水面;滑移速度可以分解成滑動(dòng)速度和滲流速度,土體固結(jié)由滲流速度決定。

(1) 準(zhǔn)拉格朗日格式控制方程

質(zhì)量守恒方程:

(3)

動(dòng)量守恒方程：

(4)

豎直向固結(jié)方程：

(5)

式中：hi為各滑移條塊的高度；mj=Ωjhj，為虛擬質(zhì)量，Ωj各滑移條塊底面積；vij=vi-vj，vj為各滑移條塊的深度平均速度；er體積增長率；gradWij為插值核函數(shù)Wij的梯度值；Pi為孔隙水壓力；σi為應(yīng)力張量，以拉應(yīng)力為正；b為體加速度；Nb為坡面空間梯度值；tb為作用于滑移體底面的應(yīng)力；P1i為超孔隙水壓力；cv為土體的滲透固結(jié)系數(shù)。

(2) 摩擦流變模型

基于SPH深度積分法,Pastor引入一系列流變模型,并進(jìn)行了修正[9，13]。成果表明,土體沿坡面滑動(dòng)過程中,需要考慮摩擦作用,并對(duì)摩擦進(jìn)行線性動(dòng)量修正。根據(jù)Coulomb摩擦理論,滑移體與坡面接觸位置，摩擦引起的包含超孔隙水壓力項(xiàng)的剪應(yīng)力為：

(6)

(3) 侵蝕準(zhǔn)則

Pastor等引入了Hungr侵蝕定律[11]模擬侵蝕：

(7)

式中：Er為侵蝕率；V0和Vf分別為滑移體的初始和最終體積;l為流徑的長度。

1.2 SPH深度積分法的求解方法

SPH深度積分法采用顯示四階Runge-Kutta法求解。滑移條塊受自重、摩擦作用、被侵蝕土體的慣性阻力,定義總作用為：N=Nin+Nex。

滑移條塊間的相互作用Nin：

(8)

自重、摩擦作用、侵蝕作用之和Nex為：

(9)

為保證解的收斂性,時(shí)間步長[12]設(shè)定為：

(10)

式中：Δl為質(zhì)點(diǎn)間的間距；g為重力加速度；h為滑移體最大高度。

滑移速度和距離的計(jì)算公式為：

(11)

(12)

2 SPH深度積分法模擬滑坡

滑移體沿坡度為θ的斜坡滑動(dòng),滑移體密度ρ為1 800 kg/m3,接觸摩擦系數(shù)tanφb為0.1～0.29[10],甚至更小[14-15];侵蝕率Er為0～1×10-4/m[16]；滲透固結(jié)系數(shù)cv設(shè)為2.71×10-4m2/s[16]；重力加速度取為9.8 m/s2。

滑坡示意圖如圖1所示,水平方向?yàn)閤方向,z坐標(biāo)豎直向上。坡面采用整體坐標(biāo)系定義,斜坡段水平長度為12 000 m,水平緩沖段長度為5 000 m。土條高度h采用相對(duì)坐標(biāo),即為相對(duì)坡面的高度?；企w的初始形狀通過多組(x，h)坐標(biāo)定義為對(duì)稱型拋物線,沿坡面長度為280 m,最大高度為11 m。對(duì)滑移體進(jìn)行剖分,其中滑移體沿x方向等分150份。定義滑移體的最大高度為hmax、前后端的水平距離為滑移體長度L,滑移體前端的滑移距離為D、前端速度為vf。

圖1 滑坡示意圖

3 數(shù)值計(jì)算與對(duì)比分析

依據(jù)上述理論,對(duì)不同工況下的陸地滑坡進(jìn)行數(shù)值模擬,討論相關(guān)參數(shù)的變化對(duì)滑坡的影響。

3.1 坡度對(duì)滑坡的影響

在接觸摩擦系數(shù)tanφb=0.1、侵蝕率Er=0時(shí),分別針對(duì)坡度θ=1°、5°、10°、13°、14°、15°六種工況進(jìn)行模擬,結(jié)果見圖2、圖3。

圖2 坡度對(duì)滑坡的影響

圖3 滑移體最終沉積形態(tài)

結(jié)果表明,坡度為1°、5°、10°三種工況下,僅滑移體的前端發(fā)生小距離的滑動(dòng),后端沒有滑動(dòng),滑坡很快停止,沉積形態(tài)與初始形態(tài)基本一致;坡度為13°、14°、15°三種工況下，滑移體發(fā)生了完全滑動(dòng),同時(shí)隨著坡度的增大,最大滑移距離和速度均增大,但在緩沖段沉積高度和長度非常接近。坡度為13°時(shí),沉積物頂部波動(dòng)較明顯。

3.2 接觸摩擦系數(shù)對(duì)滑坡的影響

在坡度θ=15°、侵蝕率Er=0時(shí),分別設(shè)定初始接觸摩擦系數(shù)tanφb=0.1、0.105、0.11、0.115、0.12,對(duì)五種工況進(jìn)行計(jì)算,結(jié)果如圖4、圖5所示。由圖4(a)可見,滑移體在斜坡段上加速滑動(dòng)時(shí),最大速度和加速度隨接觸摩擦系數(shù)的增大而減小。一旦進(jìn)入緩沖段減速時(shí),加速度基本相同。由圖4(b)和圖5可知,前四種工況,在進(jìn)入緩沖段之前,滑移體長度隨時(shí)間的變化規(guī)律基本一致,最終沉積形態(tài)也基本一致。當(dāng)tanφb=0.12時(shí),滑移體滑動(dòng)約670 m后在斜坡段上停止,未進(jìn)入緩沖段,滑移體后端堆積非常明顯。綜合分析表明,接觸摩擦系數(shù)較為敏感,其參數(shù)變化對(duì)滑坡的最大滑移速度和距離影響很大。

圖4 摩擦系數(shù)對(duì)滑坡的影響

圖5 滑移體最終沉積形態(tài)

3.3 侵蝕率對(duì)滑坡的影響

滑移體侵蝕坡面,并挾帶土體共同滑移。接觸摩擦系數(shù)tanφb=0.1、坡度θ=15°時(shí)，分別針對(duì)侵蝕率Er=0、10-7/m、10-6/m、10-5/m、10-4/m五種工況進(jìn)行計(jì)算,圖6給出了vf、L隨時(shí)間的變化曲線。當(dāng)侵蝕率為0、10-7/m、10-6/m時(shí),vf-t、L-t曲線基本重合,前端滑移距離也非常接近,說明侵蝕率小于10-6/m時(shí),侵蝕效果不明顯。由圖6、圖7可見，當(dāng)侵蝕率增大至10-5/m時(shí),滑移體最大沉積高度和長度均比前三種工況大。計(jì)算結(jié)果顯示侵蝕率為10-4/m時(shí),滑移體最終沉積高度為10.12 m,接近初始高度;沉積長度為1 019 m,約為初始長度的4倍,滑移體體量明顯增大,說明滑移體對(duì)坡床發(fā)生了顯著侵蝕。

圖6 侵蝕率對(duì)滑坡的影響

圖7 滑移體最終沉積形態(tài)

由于被侵蝕土體的初始動(dòng)量為零,根據(jù)動(dòng)量守恒可知,滑移體的加速度隨著體量的增大而減小。因此,侵蝕率越大,滑移體的最大滑移距離和速度越小。

4 結(jié) 論

通過SPH深度積分法對(duì)不同工況的滑坡進(jìn)行了數(shù)值模擬,并討論了最大滑移距離和速度及最終沉積形態(tài)的影響因素。分析表明：

(1) 當(dāng)坡度低于10°時(shí),未發(fā)生整體滑動(dòng)。隨著坡度的增長,滑移體發(fā)生整體滑動(dòng),同時(shí)隨著坡度的增大,最大滑移距離和速度均增大,但坡度變化對(duì)最終沉積形態(tài)影響較?。?/p>

(2) 接觸摩擦系數(shù)比較敏感,其參數(shù)變化對(duì)滑坡的最大滑移距離和速度影響很大；

(3) 當(dāng)侵蝕率達(dá)到10-4/m以上時(shí),侵蝕效應(yīng)明顯,滑移體的最大滑移距離和速度減小,但滑移體體量明顯增大；

(4) 建立的SPH深度積分模型能夠考慮滑移體滑移過程中土體固結(jié)和侵蝕效應(yīng),計(jì)算滑移速度和距離簡(jiǎn)便,適于模擬滑坡。

[1] 郝明輝，許 強(qiáng)，楊 磊，等.滑坡-碎屑流物理模型試驗(yàn)及運(yùn)動(dòng)機(jī)制探討[J].巖土力學(xué)，2014，35(S1):127-132.

[2] 黃 雨，郝 亮，謝 攀，等.土體流動(dòng)大變形的SPH數(shù)值模擬[J].巖土工程學(xué)報(bào)，2009，31(10):1520-1524.

[3] 申 通，王運(yùn)生，吳龍科.重慶小南?；滦纬蓹C(jī)制離散元模擬分析[J].巖土力學(xué)，2014，35(S2):667-675.

[4] Arimoto S, Murakami A. Characteristics of localized behavior of saturated soil by EFG[C]//Proceedings of 1st International Workshop on New Frontiers in Computational Geomechanics, Banff, 2002：169-172.

[5] Baxter G W, Behringer R P. Cellular automata models of granular flow[J]. Physical Review A, 1990,42(2):1017-1020.

[6] Sasaki T, Ohnishi Y, Yohinaka R. Discontinues deformation analysis and its application to rock mechanics problems[J]. Journal of Geotechnical Engineering, JSCE, 1994，Ⅲ-27:11-20.

[7] 秦 云，姜清輝，郭慧黎.滑坡速度預(yù)測(cè)的計(jì)算方法探討[J].巖土力學(xué)，2008,29(S1):373-378.

[8] Pastor M, Haddad B, Sorbino G, et al. A depth-integrated coupled SPH model for flow-like landslides and related phenomena[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2009,33:143-172.

[9] Pastor M, Blanc T, Pastor M J. A depth-integrated visco-plastic model for dilatant saturated cohesive-frictional fluidized mixtures: Application to fast catastrophic landslides[J]. Journal of Non-Newtonian Fluid Mechanics, 2009,158:142-153.

[10] Thomas Blanc. Numerical simulation of debris flows with the 2D-SPH depth integrated model[D]. Vienna, University of Natural Resources and Applied Life Sciences, 2008.

[11] Pirulli M, Pastor M. Numerical study on the entrainment of bed material into rapid landslides[J]. Geotechnique, 2012,62(11):959-972.

[12] Monaghan J J, Gingold R A. Shock simulation by the particle method SPH[J]. Journal of Computational Physics, 1983,52:374-389.

[13] Pastor M, Quecedo M, Gonzalez E, et al. A simple approximation to bottom friction for Bingham fluid depth integrated models[J]. Journal of Hydraulic Engineering ASCE, 2004,130(2):149-155.

[14] McDougall S, Hungr O. A model for the analysis of rapid landslide motion across three-dimensional terrain[J]. Canadian Geotechnical Journal, 2004,41:1084-1097.

[15] McDougall S, Hungr O. Dynamic modeling of entrainment in rapid landslide[J]. Canadian Geotechnical Journal, 2005,42:1437-1448.

[16] Pastor M. The “SafeLand” Project: Living with landslide risk in Europe, Landslide runout: Guidelines for using a simple propagation code[R]. Caminos, Fundación Agustín de Betancourt ETS de Ingenieros de Caminos, 2012.

Numerical Analysis on Landslides Based on A SPH Depth Integral Model

TAO Yanqi1, LI Jiagang2

(1.StateKeyLaboratoryofCoastalandOffshoreEngineering,DalianUniversityofTechnology,Dalian,Liaoning116024,China;2.GeneralResearchInstituteofChinaNationalOffshoreOilCorporation,Beijing100081,China)

Landslides can cause serious damages to the structure of the engineering, sliding velocity and distance are two important parameters in the quantification and analysis of landslides. However, the existing prediction method of sliding velocity and distance is complicated. In view of this situation, a governing equation was established considering soil consolidation and erosion effect. The SPH depth integral method was used to simulate landslides by the selected frictional rheological model. Comparing the time history curve of velocity, length of sliding body under the different conditions of slope angle, contact frictional coefficient and erosion rate, the final deposit shapes of the sliding body were analyzed, the maximum sliding distance and velocity were summarized and the underlying influence of related physical parameters were discussed. These research results can provide some technical reference for the disaster warning of landslides.

sliding velocity; SPH depth integral model; frictional rheological model; erosion effect

10.3969/j.issn.1672-1144.2015.05.022

2015-05-03

2015-06-01

國家自然科學(xué)基金(50909014)；國家重大科技專項(xiàng)“南海深水油氣開發(fā)示范工程”子課題“南海北部陸坡區(qū)地質(zhì)災(zāi)害風(fēng)險(xiǎn)評(píng)價(jià)預(yù)測(cè)研究”(2011ZX05056-001-02)；中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)資金資助

陶言祺(1989—)，男，安徽蚌埠人，碩士，主要從事水利規(guī)劃和水工設(shè)計(jì)工作。E-mail:1056047556@qq.com

李家鋼(1964—)，男，遼寧大連人，工程師，碩士，主要從事海洋地質(zhì)災(zāi)害防治研究。E-mail:lijg2@cnooc.com.cn

P642

A

1672—1144(2015)05—0112—05