堆載誘發(fā)型土質(zhì)滑坡運(yùn)動(dòng)特征物質(zhì)點(diǎn)法模擬

吳方東，張 巍，史卜濤，施 斌，張 云，鄭 帥

(南京大學(xué)地球科學(xué)與工程學(xué)院，江蘇 南京 210046)

堆載誘發(fā)型土質(zhì)滑坡運(yùn)動(dòng)特征物質(zhì)點(diǎn)法模擬

吳方東，張 巍，史卜濤，施 斌，張 云，鄭 帥

(南京大學(xué)地球科學(xué)與工程學(xué)院，江蘇 南京 210046)

坡頂堆載是人類工程活動(dòng)誘發(fā)滑坡的主因之一。物質(zhì)點(diǎn)法(MPM)屬于一種無(wú)網(wǎng)格數(shù)值計(jì)算方法，它能夠有效模擬滑坡大變形全過(guò)程物質(zhì)行為與運(yùn)動(dòng)特征。文章基于線性形函數(shù)離散方法、MUSL求解格式及Drucker-Prager屈服準(zhǔn)則，建立了可用于滑坡全過(guò)程模擬的單套單相物質(zhì)點(diǎn)模型；通過(guò)對(duì)比干燥鋁棒堆積物模擬砂堆失穩(wěn)過(guò)程的基準(zhǔn)試驗(yàn)結(jié)果，對(duì)模型有效性進(jìn)行了驗(yàn)證。對(duì)堆載誘發(fā)型土質(zhì)滑坡典型工況進(jìn)行了物質(zhì)點(diǎn)法全過(guò)程模擬，獲得了滑坡全過(guò)程中典型時(shí)刻坡體形態(tài)、塑性應(yīng)變分布以及控制點(diǎn)滑速演化趨勢(shì)。結(jié)果表明：算例堆載誘發(fā)型土質(zhì)滑坡屬推移型滑坡，具有漸進(jìn)性破壞特征，可分為坡頂壓縮、局部蠕滑、加速滑動(dòng)與減速滑動(dòng)等四個(gè)階段。參數(shù)分析結(jié)果亦表明，堆載誘發(fā)型土質(zhì)滑坡前緣物質(zhì)運(yùn)動(dòng)特征量均與堆載量間存在強(qiáng)正相關(guān)性、而與土體黏聚力及內(nèi)摩擦角存在強(qiáng)負(fù)相關(guān)性。統(tǒng)計(jì)29種典型工況，分別建立了峰值滑動(dòng)加速度、最大滑速、最大滑距及坡體最大動(dòng)能等運(yùn)動(dòng)特征量與堆載量、土體黏聚力及內(nèi)摩擦角之間的線性回歸方程，可用于堆載誘發(fā)型土質(zhì)滑坡致災(zāi)行為預(yù)測(cè)。

土質(zhì)滑坡；堆載；物質(zhì)點(diǎn)法；運(yùn)動(dòng)參數(shù)

斜坡棄土堆載、坡頂興建房屋或構(gòu)筑物等人類工程活動(dòng)是誘發(fā)滑坡的主因之一[1～3]。國(guó)內(nèi)外學(xué)者已對(duì)堆載誘發(fā)型滑坡開展了一系列研究。試驗(yàn)方面，胡田飛等[4]基于室內(nèi)土質(zhì)邊坡模型，對(duì)分級(jí)加載作用下滑坡演化過(guò)程進(jìn)行了觀察，發(fā)現(xiàn)堆載誘發(fā)型滑坡加速滑動(dòng)階段與劇滑啟動(dòng)的過(guò)渡時(shí)間極短，滑動(dòng)面貫通后前緣坡體會(huì)迅速解體剪出。Li等[5]在室內(nèi)對(duì)一個(gè)砂土邊坡模型進(jìn)行了加載破壞試驗(yàn)，并嘗試使用離散元方法對(duì)斜坡破壞行為進(jìn)行模擬。此外，Zhu等[6]還采用分布式光纖傳感技術(shù)，表征土質(zhì)邊坡模型加載破壞試驗(yàn)過(guò)程中的應(yīng)變場(chǎng)信息。理論方面，陳春利等[7]基于有限元模擬，分析了人工堆載黃土滑坡變形機(jī)制與破壞機(jī)理，結(jié)果表明，邊坡堆載改變了坡體內(nèi)的應(yīng)力分布，致使坡體中、下部土體剪應(yīng)力逐漸接近其抗剪強(qiáng)度，邊坡穩(wěn)定性隨之降低直至失穩(wěn)破壞。王洪兵等[8]采用有限差分法分析云南寶騰高速公路路塹邊坡，發(fā)現(xiàn)坡體中部及坡腳堆載對(duì)邊坡穩(wěn)定有利，坡頂堆載對(duì)坡體穩(wěn)定性最為不利。此外，董夫錢等[9]以浙江上三公路路塹1號(hào)滑坡為例，采用有限元方法分析發(fā)現(xiàn)堆載除直接改變邊坡穩(wěn)定性之外，還將間接影響到坡體內(nèi)部地下水的滲流條件，并進(jìn)一步降低坡體穩(wěn)定性。

值得指出，既有研究多集中于堆載誘發(fā)型滑坡失穩(wěn)機(jī)制的探討，而研究斜坡從變形演化發(fā)展到失穩(wěn)運(yùn)動(dòng)的全過(guò)程，能更全面地表征滑坡的地質(zhì)力學(xué)行為及其動(dòng)力演化機(jī)制[10～12]。數(shù)值模擬是研究滑坡的有力工具，但傳統(tǒng)的極限平衡計(jì)算方法只能對(duì)斜坡體的穩(wěn)定性進(jìn)行分析[12]；有限元方法在處理大變形計(jì)算時(shí)則因網(wǎng)格畸變而出現(xiàn)數(shù)值困難[13～14]，這些方法都難以表征滑坡全過(guò)程中的物質(zhì)運(yùn)動(dòng)特征。

物質(zhì)點(diǎn)法屬于一種無(wú)網(wǎng)格法，其在數(shù)值計(jì)算中不需要生成網(wǎng)格，而按任意分布的坐標(biāo)點(diǎn)來(lái)構(gòu)造插值函數(shù)實(shí)現(xiàn)離散的控制方程，可模擬各種復(fù)雜形態(tài)流場(chǎng)[14]，因而適用于滑坡大變形及超大變形模擬。該方法還可導(dǎo)入各種土體本構(gòu)模型，考慮流固耦合[15～16]，而與SPH等其他無(wú)網(wǎng)格法相比較，物質(zhì)點(diǎn)法具有更高的計(jì)算效率與數(shù)值穩(wěn)定性[15]。目前國(guó)內(nèi)外使用物質(zhì)點(diǎn)法模擬土質(zhì)滑坡運(yùn)動(dòng)全過(guò)程已取得一系列進(jìn)展[17～21]，已被用于模擬降雨入滲誘發(fā)非飽和土滑坡等復(fù)雜物質(zhì)運(yùn)動(dòng)問(wèn)題[22]。本文采用物質(zhì)點(diǎn)法模型模擬堆載誘發(fā)型土質(zhì)滑坡的典型工況，重現(xiàn)了滑坡運(yùn)動(dòng)全過(guò)程，獲取了典型時(shí)刻坡體的滑速、滑距、塑性應(yīng)變等物質(zhì)運(yùn)動(dòng)及變形特征。

1 物質(zhì)點(diǎn)模型

1.1基本描述

物質(zhì)點(diǎn)法采用拉格朗日質(zhì)點(diǎn)系和歐拉網(wǎng)格雙重描述的方法。如圖1，離散的質(zhì)點(diǎn)系攜帶介質(zhì)的所有物質(zhì)信息，計(jì)算網(wǎng)格僅用于動(dòng)量方程的求解與空間導(dǎo)數(shù)求解。

圖1 物質(zhì)點(diǎn)法示意圖Fig.1 Schematic diagram of the material point method

不考慮熱交換，物質(zhì)點(diǎn)法控制方程包括虛功方程、動(dòng)量方程、幾何方程、本構(gòu)方程以及邊界條件方程等，具體參見文[19]、[21]原理部分。

1.2物質(zhì)點(diǎn)離散

采用具有線性形函數(shù)的背景網(wǎng)格進(jìn)行離散，連續(xù)體的密度可近似表示為：

式中：np——質(zhì)點(diǎn)總數(shù)；

mp——質(zhì)點(diǎn)p的質(zhì)量；

δ——Dirac Delta函數(shù)；

xi——空間坐標(biāo)；

xip——質(zhì)點(diǎn)p的坐標(biāo)。

據(jù)文[13]可知，動(dòng)量方程和給定面力邊界條件的等效積分弱形式(即動(dòng)量方程)為式(2)：

引入式(1)可將虛功方程式(2)簡(jiǎn)化為求和的形式：

式中，用帶有下標(biāo)p的量來(lái)表示物質(zhì)點(diǎn)攜帶的相應(yīng)物理量(下同)，同時(shí)引入假想邊界層厚度h將虛功方程式左端最后一項(xiàng)轉(zhuǎn)化為體積分。

用物質(zhì)點(diǎn)法求解時(shí)，質(zhì)點(diǎn)與背景網(wǎng)格不發(fā)生相對(duì)移動(dòng)，故可以通過(guò)建立在背景網(wǎng)格結(jié)點(diǎn)上的有限元形函數(shù)NI(xi)實(shí)現(xiàn)質(zhì)點(diǎn)和背景網(wǎng)格結(jié)點(diǎn)之間的信息映射，即：

其中，xip、uip、uip，j、δuip分別為質(zhì)點(diǎn)p的坐標(biāo)、位移、位移的導(dǎo)數(shù)以及虛位移，δuiI為結(jié)點(diǎn)的虛位移，用帶有下標(biāo)I的量來(lái)表示網(wǎng)格結(jié)點(diǎn)的變量。

式(5)～(7)結(jié)合式(3)，同時(shí)考慮結(jié)點(diǎn)虛位移δuiI在本質(zhì)邊界Γu上為零，且在其它結(jié)點(diǎn)處具有任意性，得到背景網(wǎng)格結(jié)點(diǎn)的運(yùn)動(dòng)方程：

?

其中

為第I個(gè)結(jié)點(diǎn)在i方向上的動(dòng)量；

為背景網(wǎng)格的質(zhì)量矩陣；

為結(jié)點(diǎn)內(nèi)力；

為結(jié)點(diǎn)外力。

當(dāng)采用集中質(zhì)量矩陣時(shí)，

式(9)簡(jiǎn)化為：

則背景網(wǎng)格結(jié)點(diǎn)的運(yùn)動(dòng)方程為：

?

1.3求解格式

求解動(dòng)量方程時(shí)，需要將某時(shí)刻質(zhì)點(diǎn)的質(zhì)量、動(dòng)量等信息映射到背景網(wǎng)格，以計(jì)算背景網(wǎng)格結(jié)點(diǎn)的質(zhì)量和動(dòng)量等。動(dòng)量方程采用顯式時(shí)間積分求解。網(wǎng)格結(jié)點(diǎn)的速度變化量可以用于計(jì)算質(zhì)點(diǎn)的應(yīng)變?cè)隽?，從而?duì)質(zhì)點(diǎn)進(jìn)行應(yīng)力更新。

根據(jù)應(yīng)力更新時(shí)所采用的結(jié)點(diǎn)速度的不同，物質(zhì)點(diǎn)法可分為USF、USL和MUSL等不同的求解格式。在時(shí)間步開始時(shí)進(jìn)行應(yīng)力更新的格式稱為USF格式；在時(shí)間步結(jié)束時(shí)進(jìn)行應(yīng)力更新的格式稱為USL格式；在時(shí)間步結(jié)束時(shí)的質(zhì)點(diǎn)動(dòng)量映射到背景網(wǎng)格后計(jì)算結(jié)點(diǎn)速度，并用此時(shí)的結(jié)點(diǎn)速度更新應(yīng)力，稱之為MUSL格式。本文模型采用MUSL格式，它具有較好的能量守恒特性[13]。

1.4本構(gòu)模型

物質(zhì)點(diǎn)法可分別采用單套、兩套與三套物質(zhì)點(diǎn)，并分別基于單相、兩相與三相理論表征土體力學(xué)行為[22]，文[21]據(jù)此整理了單套單相、單套兩相、單套三相、兩套兩相、兩套三相與三套三相等6種土體物質(zhì)點(diǎn)表征方法。為簡(jiǎn)化算法，提高計(jì)算效率，本文采用單套單相物質(zhì)點(diǎn)表征模型，該模型適用于飽和土體，包括水飽和土與氣飽和土(干燥土體)。

考慮到邊坡初始應(yīng)力場(chǎng)對(duì)坡體穩(wěn)定性的影響，通常采用線彈性模型確定施加重力荷載后邊坡的初始應(yīng)力場(chǎng)。對(duì)于各向同性線彈性本構(gòu)模型，焦曼應(yīng)力率和變形率之間的關(guān)系為

邊坡開始產(chǎn)生塑性變形直至失穩(wěn)出現(xiàn)大變形階段，則應(yīng)采用彈塑性本構(gòu)模型。Mohr-Coulomb強(qiáng)度準(zhǔn)則雖廣泛應(yīng)用于巖土領(lǐng)域，但考慮其屈服面是由六個(gè)平面圍成的錐形表面，容易導(dǎo)致數(shù)值求解困難。本文采用的Drucker-Prager屈服準(zhǔn)則，其屈服面為一光滑圓錐面。屈服函數(shù)fs為：

式中：τ——等效剪應(yīng)力；

σm——球應(yīng)力；

qφ——摩擦系數(shù)；

kφ——純剪切狀態(tài)時(shí)的屈服應(yīng)力。

qφ，kφ與材料的黏聚力c和摩擦角φ的關(guān)系由下式確定：

加號(hào)表示DP屈服面在π平面上內(nèi)接MC屈服面，減號(hào)則表示外接。

2 模型驗(yàn)證

Bui等采用干燥鋁棒堆積物進(jìn)行試驗(yàn)，模擬了干燥砂堆的失穩(wěn)過(guò)程[23]。鋁棒被排列在一個(gè)長(zhǎng)方體區(qū)域內(nèi)，其長(zhǎng)、高、寬分別為0.20 m、0.10 m、0.05 m。材料摩擦角約19.8°，泊松比0.3，平均體積模量0.7 MPa。

建立以上鋁棒堆積物的單套單相物質(zhì)點(diǎn)模型，示意如圖2。

圖2 物質(zhì)點(diǎn)模型(80 000個(gè)物質(zhì)點(diǎn))Fig.2 MPM model (80 000 particles)

X=0 m,X=0.2 m取對(duì)稱邊界條件。Y=0 m,Y=0.05 m取對(duì)稱邊界條件。Z=0 m取固定邊界條件；Z=0.1 m取自由邊界條件。計(jì)算網(wǎng)格間距為1 mm，物質(zhì)點(diǎn)間距為0.5 mm，物質(zhì)點(diǎn)總數(shù)為80 000個(gè)。典型計(jì)算時(shí)刻鋁棒堆積物失穩(wěn)形態(tài)與塑性應(yīng)變?cè)茍D見圖3。

圖3 失穩(wěn)形態(tài)與塑性應(yīng)變?cè)茍DFig.3 Instability morphology and plastic strain contour map

圖4為試驗(yàn)過(guò)程中鋁棒堆積物的最終失穩(wěn)形態(tài)。圖5為分別采用SPH模型[23]和本文MPM模型所得最終失穩(wěn)形態(tài)界面與試驗(yàn)失穩(wěn)形態(tài)界面的比較結(jié)果。對(duì)比三者結(jié)果發(fā)現(xiàn)，堆積物滑動(dòng)跡線與穩(wěn)定后的坡面跡線吻合程度較高，驗(yàn)證了前文所建立的單套單相物質(zhì)點(diǎn)模型的合理性與有效性。

圖4 鋁棒堆積物失穩(wěn)最終形態(tài)(Bui，2008)Fig.4 Final form of the aluminum rods debris

圖5 模擬與實(shí)驗(yàn)結(jié)果比較Fig.5 Comparison of the results of simulation and experiment

3 土質(zhì)滑坡物質(zhì)點(diǎn)法模擬

建立一坡高25 m、坡角45°、地基厚10 m的理想均勻土質(zhì)邊坡模型。采用規(guī)則的八結(jié)點(diǎn)六面體單元，網(wǎng)格間距1 m，質(zhì)點(diǎn)間距0.5 m，每個(gè)計(jì)算網(wǎng)格內(nèi)在水平和豎直方向各有2個(gè)物質(zhì)點(diǎn)，垂直紙面方向有1個(gè)物質(zhì)點(diǎn)。

根據(jù)背景網(wǎng)格三維坐標(biāo)值確定每個(gè)物質(zhì)點(diǎn)的幾何位置。物質(zhì)點(diǎn)的質(zhì)量為其所在網(wǎng)格區(qū)域所代表土體單元的質(zhì)量。如圖6所示，邊坡離散體設(shè)置時(shí)，由于其坡面線切割六面體背景網(wǎng)格，故而在坡面附近的物質(zhì)點(diǎn)將偏離背景網(wǎng)格的體心，需按照三角形區(qū)域的形心布置，坡面物質(zhì)點(diǎn)的質(zhì)量也是三角形區(qū)域代表的土體質(zhì)量。

根據(jù)以上方法所建立的物質(zhì)點(diǎn)離散模型如圖7所示。模型底部采用固定邊界，兩側(cè)施加對(duì)稱邊界條件。此外，通過(guò)約束三維模型的y方向位移來(lái)模擬平面應(yīng)變狀態(tài)，即模型的前后面亦施加對(duì)稱邊界條件。在坡肩、坡中與坡趾等位置坡面上分別設(shè)置A、B、C三個(gè)計(jì)算監(jiān)測(cè)點(diǎn)。

土體的主要物理力學(xué)參數(shù)取值：密度ρ=2 000 kg/m3，彈性模量E=70 MPa，泊松比u=0.3，黏聚力c=20 kPa，內(nèi)摩擦角φ=25°。

圖6 物質(zhì)點(diǎn)及網(wǎng)格布置示意圖Fig.6 Schematic diagram of the material point and grids layout

圖7 邊坡物質(zhì)點(diǎn)模型(10 975個(gè)物質(zhì)點(diǎn))Fig.7 Slope MPM model (10 975 particles)

4 滑坡物質(zhì)運(yùn)動(dòng)特征

4.1荷載設(shè)置

在坡頂設(shè)置高度為5 m的堆載體,分布范圍自坡肩A點(diǎn)向內(nèi)20 m。為了消除參數(shù)分析過(guò)程中由于堆載體體積變化而對(duì)坡體運(yùn)動(dòng)可能產(chǎn)生的影響，文章通過(guò)改變堆載體的重力加速度的方式來(lái)模擬不同的坡頂荷載。另外，為了減小計(jì)算過(guò)程中的數(shù)值震蕩，施加重力荷載時(shí)均采用線性加載的方式。此外，考慮到邊坡初始應(yīng)力場(chǎng)對(duì)坡體穩(wěn)定性的影響，0～3 s時(shí)刻間采用線彈性本構(gòu)模型，3 s后采用Drucker-Prager彈塑性本構(gòu)[19]，如圖8所示。圖9為重力施加完成后邊坡內(nèi)部初始應(yīng)力場(chǎng)分布。

圖8 加載模式Fig.8 Loading mode

圖9 初始應(yīng)力場(chǎng)Fig.9 Initial stress field

4.2滑坡全過(guò)程分析

基于圖9的初始應(yīng)力場(chǎng)進(jìn)行加載，并對(duì)堆載體賦3倍的重力加速度值,即相當(dāng)于294 kPa的坡頂荷載，計(jì)算時(shí)間為15 s。圖10繪出了滑坡各典型時(shí)刻坡體位移形態(tài)與塑性應(yīng)變區(qū)演化過(guò)程。

圖11繪出了坡體失穩(wěn)滑動(dòng)前后A、B、C三個(gè)測(cè)點(diǎn)的水平滑速曲線。由圖可見，從第3秒坡體失穩(wěn)后3個(gè)監(jiān)測(cè)點(diǎn)均出現(xiàn)加速段，隨后在第6～7秒內(nèi)依次轉(zhuǎn)變?yōu)闇p速段，計(jì)算時(shí)刻10 s左右滑動(dòng)停止，其中，坡腳C點(diǎn)在起滑后迅速達(dá)到0.45 m/s滑速，隨后基本以勻速滑動(dòng)至4 s時(shí)刻，才繼續(xù)加速滑動(dòng)。3個(gè)測(cè)點(diǎn)中，坡中B點(diǎn)最大滑速最大，達(dá)到5.8 m/s；坡肩A點(diǎn)最大滑速達(dá)到5 m/s；坡角C點(diǎn)最大滑速最小，僅1.5 m/s。分別采用式(20)、(21)指數(shù)衰減函數(shù)對(duì)A點(diǎn)滑速曲線的上升段和下降段進(jìn)行分段擬合；而B、C兩點(diǎn)滑速曲線則可由式(22)、(23)表示的Gauss函數(shù)進(jìn)行全段擬合，擬合函數(shù)的表達(dá)式具體如下：

圖10 滑坡位移形態(tài)與塑性應(yīng)變區(qū)演化Fig.10 Landslide displacement morphology and plastic strain zone evolution

圖11 滑速演化趨勢(shì)Fig.11 Slide speed evolution trends

結(jié)合圖10、11分析可知，本算例堆載誘發(fā)型土質(zhì)滑坡全過(guò)程可劃分為以下4個(gè)階段:

(1)Ⅰ坡頂壓縮

如圖10(b)所示，坡頂土體在堆載的作用下發(fā)生壓縮變形，堆載體兩側(cè)邊緣以及坡腳產(chǎn)生應(yīng)力集中，開始出現(xiàn)塑性應(yīng)變。

(2)Ⅱ局部蠕滑

如圖10(c)、(d)所示，荷載的持續(xù)作用使得弧形塑性區(qū)(潛在滑帶)首先從坡腳開始形成并逐漸向坡體內(nèi)部發(fā)育?；w前緣開始鼓脹隆起，坡腳處土體的剪應(yīng)力超過(guò)了抗剪強(qiáng)度并開始向前蠕滑。

(3)Ⅲ加速滑動(dòng)

如圖10(e)所示，持續(xù)的蠕滑使得弧形塑性區(qū)完全貫通，滑體前緣進(jìn)一步鼓脹形成滑舌，最終沿滑帶從坡腳處剪出產(chǎn)生整體滑動(dòng)，滑速迅速增大進(jìn)入巨滑階段，如圖11所示。

(4)Ⅳ減速滑動(dòng)

如圖11所示，從6～7 s階段，坡體先后由加速滑動(dòng)階段轉(zhuǎn)入減速滑動(dòng)階段，滑坡體繼續(xù)前移，滑舌下方出現(xiàn)高塑性應(yīng)變區(qū)，滑坡前緣最大滑距達(dá)到17 m，如圖10(f)所示。

可見，本算例斜坡在荷載作用初期以坡頂壓縮變形為主，隨坡體內(nèi)塑性應(yīng)變的累積，逐漸向臨空面發(fā)展，沿弧形滑帶發(fā)生整體滑動(dòng)，表現(xiàn)出推移式滑坡的變形破壞特點(diǎn)，即以壓剪破壞為主[24]，這與文[4]所觀察到的試驗(yàn)現(xiàn)象相吻合。邊坡破壞一般表現(xiàn)為滑帶形成、發(fā)育直至整體貫通后滑動(dòng)的漸進(jìn)破壞過(guò)程[24-25]，算例堆載誘發(fā)型滑坡除具有這種階段性演化特征之外，還表現(xiàn)出劇滑前的蠕動(dòng)變形時(shí)間短，滑帶形成后迅速解體破壞的特點(diǎn)[4]。

5 參數(shù)分析

高速運(yùn)動(dòng)的滑體攜帶巨大的動(dòng)能，具有極高的沖擊力。因此，運(yùn)動(dòng)過(guò)程中動(dòng)能大小可以用來(lái)評(píng)估滑坡的運(yùn)動(dòng)災(zāi)害性[26]。通過(guò)改變不同坡頂荷載和土體物理力學(xué)參數(shù)，分析了滑坡動(dòng)能的變化情況，結(jié)果如圖12所示。圖12(a)為考慮不同坡頂荷載影響，且保持黏聚力為20 kPa，內(nèi)摩擦角為25°的各工況；圖12(b)為考慮不同土體內(nèi)摩擦角影響，且保持坡頂堆載為294 kPa，土體黏聚力20 kPa的各工況；圖12(c)為考慮不同土體黏聚力影響，且保持坡頂堆載為294 kPa，土體內(nèi)摩擦角為25°的各工況。

從圖12可見，滑坡動(dòng)能峰值隨坡頂堆載增大而接近線性增加，隨內(nèi)摩擦角以及黏聚力的增大而接近線性的減小，且動(dòng)能變化曲線基本都呈現(xiàn)出上升段與下降段基本對(duì)稱的特征。圖12(a)表明，隨坡頂荷載的增加，系統(tǒng)動(dòng)能峰值隨之增加，而達(dá)到峰值所需要的時(shí)間略有減小，坡體達(dá)到變形穩(wěn)定的時(shí)間隨之縮短。圖12(b)表明，內(nèi)摩擦角的增大使得系統(tǒng)動(dòng)能峰值明顯減小。在高摩擦角情況下，坡體滑動(dòng)達(dá)到動(dòng)能峰值的時(shí)間也有較明顯的推遲，滑坡體需要更長(zhǎng)的時(shí)間達(dá)到變形穩(wěn)定階段。圖12(c)表明，黏聚力的變化只對(duì)系統(tǒng)動(dòng)能峰值有影響，即使在高黏聚力狀態(tài)下，滑坡達(dá)到變形穩(wěn)定所需要的時(shí)間也基本不變。需指出， Drucker-Prager屈服準(zhǔn)則未考慮材料的軟化特性，因此在時(shí)間尺度上，計(jì)算時(shí)間與實(shí)際滑坡時(shí)間尺度的對(duì)應(yīng)關(guān)系需另文研究。

總體而言，坡頂荷載的增加與土體強(qiáng)度的降低，將使得滑坡動(dòng)能峰值顯著增加。同時(shí)，高荷載、低摩擦角還將縮短坡體動(dòng)能達(dá)到峰值的時(shí)間，使得滑坡的突發(fā)性變得更加顯著。

滑速、滑距以及動(dòng)能等滑坡運(yùn)動(dòng)特征參數(shù)對(duì)于評(píng)估滑坡致災(zāi)范圍和致災(zāi)程度具有重要意義。為了進(jìn)一步研究上覆荷載、黏聚力和內(nèi)摩擦角對(duì)滑速、滑距以及坡體動(dòng)能等參數(shù)的影響程度，建立了29種不同堆載量、土體黏聚力與內(nèi)摩擦角組合工況下邊坡的物質(zhì)點(diǎn)模型，對(duì)滑坡前緣C點(diǎn)的峰值加速度、滑速、最大滑距以及坡體最大動(dòng)能進(jìn)行了統(tǒng)計(jì)(表1)。

圖12 不同條件下系統(tǒng)動(dòng)能變化曲線Fig.12 System kinetic energy variation curves under different conditions

黏聚力/kPa內(nèi)摩擦角/(°)總荷載/kPa滑速峰值/(m·s-1)加速度峰值/(m·s-2)坡體滑距/m動(dòng)能峰值/kJ202519613505212153820252451570681525592025294186082173596202534321110721463415252942551102349291625294236110224654172529421510421436718252942011002041041925294192094193841212529418208717335322252941760831631122325294166073152878242529415406414265425252941520581324202625294149058122180272529414405211196928252941310451017622925294120038915633025294099032813712020294430167286814202129439612926602620222943031322353052023294261109214677202429422110419411220262941600721631052027294150062142607202829412005012215520292940980381117902030294072020101382

對(duì)以上計(jì)算結(jié)果進(jìn)行回歸統(tǒng)計(jì)分析，并依次用函數(shù)式(24)～(27)建立C點(diǎn)速度峰值v、加速度峰值a、滑距s以及最大動(dòng)能k與黏聚力c、內(nèi)摩擦角φ和坡頂堆載p三者的線性表達(dá)式，四式的相關(guān)性系數(shù)分別達(dá)到0.92、0.98、0.99和0.99，表明三個(gè)影響因素與滑速、加速度和滑距之間存在強(qiáng)線性相關(guān)性。

觀察上式可以發(fā)現(xiàn)，黏聚力和內(nèi)摩擦角與以上4個(gè)滑坡運(yùn)動(dòng)特征量呈負(fù)線性相關(guān)性，即二者減小都會(huì)使得峰值加速度、最大滑速、滑距及動(dòng)能顯著增加，且內(nèi)摩擦角的影響程度遠(yuǎn)甚于黏聚力。而堆載量則與以上4個(gè)滑坡運(yùn)動(dòng)特征量呈正相關(guān)性，其增大將導(dǎo)致峰值加速度、最大滑速、最大滑距及最大動(dòng)能增大。

6 結(jié)論

基于物質(zhì)點(diǎn)法基本理論，采用數(shù)值算例，模擬了堆載誘發(fā)型土質(zhì)滑坡運(yùn)動(dòng)全過(guò)程，分析了物質(zhì)運(yùn)動(dòng)特征，得出以下結(jié)論：

(1)物質(zhì)點(diǎn)法是一種能夠有效模擬滑坡大變形的無(wú)網(wǎng)格數(shù)值方法。該方法可對(duì)滑坡運(yùn)動(dòng)全過(guò)程進(jìn)行重現(xiàn)，可提取出任意計(jì)算時(shí)刻坡體的滑速、滑距、塑性應(yīng)變等物質(zhì)運(yùn)動(dòng)及變形特征。

(2)算例堆載誘發(fā)型土質(zhì)滑坡屬推移式滑坡，其力學(xué)機(jī)制為壓剪破壞，以沿弧形滑帶的整體滑動(dòng)為主，并具有漸進(jìn)式破壞特征，可劃分為坡頂壓縮、局部蠕滑、加速滑動(dòng)與減速滑動(dòng)等4個(gè)階段。

(3)堆載誘發(fā)型土質(zhì)滑坡運(yùn)動(dòng)特征與上覆堆載量和土體物理力學(xué)性質(zhì)密切相關(guān)。上覆堆載量越大，土體內(nèi)摩擦角越小，滑坡越趨向于在短時(shí)間內(nèi)達(dá)到較高滑速，突發(fā)性更為明顯，體現(xiàn)出“突然滑動(dòng)，迅速穩(wěn)定”的運(yùn)動(dòng)特征，而土體黏聚力對(duì)此影響則相對(duì)較小。

(4)統(tǒng)計(jì)了29種典型工況下不同堆載量與土體物理力學(xué)指標(biāo)數(shù)據(jù)，建立了滑坡峰值加速度、最大滑速、最大滑距與最大動(dòng)能的線性回歸方程，可用于此類滑坡致災(zāi)行為預(yù)測(cè)。

[1] 張茂省, 李同錄. 黃土滑坡誘發(fā)因素及其形成機(jī)理研究[J]. 工程地質(zhì)學(xué)報(bào), 2011, 19(4): 530-540. [ZHANG M S, LI T L. Triggering factors and forming mechanism of loess landslides[J]. Journal of Engineering Geology, 2011, 19(4): 530-540. (in Chinese)]

[3] 劉傳正. 南昆鐵路八渡滑坡成因機(jī)理新認(rèn)識(shí)[J]. 水文地質(zhì)工程地質(zhì), 2007, 34(5): 1-5. [LIU C Z. A new discussion about genesis and failure mechanism of Badu landslides in Nanning-Kunming railway[J]. Hydrogeology amp; Engineering Geology, 2007,34(5): 1-5. (in Chinese)]

[4] 胡田飛. 堆載誘發(fā)型滑坡變形演化機(jī)理的模型試驗(yàn)研究[J]. 石家莊鐵道大學(xué)學(xué)報(bào)(自然科學(xué)版), 2016,29(1):38-42. [HU T F. Model test study of deformation mechanism of landslide induced by loading[J].Journal of Shijiazhuang Railway Institute,2016, 29(1):38-42.(in Chinese)]

[5] LI N, CHENG Y M. Laboratory and 3-D distinct element analysis of the failure mechanism of a slope under external surcharge[J]. Natural Hazards and Earth System Sciences, 2015, 15:35-43.

[6] ZHU H H, SHI B, ZHANG J,etal. Distributed fiber optic monitoring and stability[J]. Journal of Mountain Science, 2014, 11(4): 979-989.

[7] 陳春利, 李同錄, 賀凱, 等. 人工堆載誘發(fā)黃土滑坡失穩(wěn)機(jī)制分析[J]. 中國(guó)地質(zhì)災(zāi)害與防治學(xué)報(bào), 2014, 25(1): 1-5. [CHEN C L, LI T L, HE K,etal. Mechanism of loess landslide induced by loading[J]. The Chinese Journal of Geological Hazard and Control, 2014,25(1): 1-5. (in Chinese)]

[8] 王洪兵, 姚磊華, 文海. 基于FLAC3D進(jìn)行堆載作用下邊坡穩(wěn)定分析[J]. 勘察科學(xué)技術(shù), 2012 (4): 5-10. [WANG H B, YAO L H, WEN H. Stability analysis of slope under preloading based on FLAC3D[J].Site Investigation Science and Technology, 2012 (4): 5-10. (in Chinese)]

[9] 董夫錢, 繆志順, 呂慶, 等. 公路堆載誘發(fā)型滑坡穩(wěn)定性分析[J]. 巖石力學(xué)與工程學(xué)報(bào), 2004 (增刊): 4517-4520. [DONG F Q, MIAO Z S, LV Q,etal. Mechanism analysis of landslide induced by accumulation loading of highway[J]. Chinese Journal of Rock Mechanics and Engineering, 2004 (Sup1): 4517-4520. (in Chinese)]

[10] 張偉鋒, 黃潤(rùn)秋, 裴向軍. 大光包滑坡運(yùn)動(dòng)特征及其過(guò)程分析[J]. 工程地質(zhì)學(xué)報(bào), 2015, 23(5): 866-885. [ZHANG W F, HUANG R Q, PEI X J. Analysis on kinematics characteristics and movement process of Daguangbao landslide[J].Journal of Engineering Geology, 2015, 23(5): 866-885. (in Chinese)]

[11] 羅鋒, 柴波, 方恒, 等. 順層巖質(zhì)滑坡運(yùn)動(dòng)過(guò)程數(shù)值模擬研究[J]. 水文地質(zhì)工程地質(zhì), 2016, 43(3): 124-130. [LUO F,CHAI B,FANG H,etal. Numerical simulation of the bedding rock landslide motion process[J]. Hydrogeology amp; Engineering Geology, 2016, 43(3): 124-130. (in Chinese)]

[12] 殷坤龍, 姜清輝, 汪洋. 滑坡運(yùn)動(dòng)過(guò)程仿真分析[J]. 地球科學(xué)(中國(guó)地質(zhì)大學(xué)學(xué)報(bào)), 2002,27(5): 632-636. [YIN K L, JIANG Q H, WANG Y. Simulation of landslide movement process by discontinuous deformation analysis[J]. Earth Science(Journal of China University of Geosciences), 2002, 27(5):632-636. (in Chinese)]

[13] 張雄, 廉艷平, 劉巖,等. 物質(zhì)點(diǎn)法[M]. 北京：清華大學(xué)出版社, 2013. [ZHANG X, LIAN Y P, LIU Y,etal. Material point method[M]. Beijing: Tsinghua University Press, 2013. (in Chinese)]

[14] 張雄, 劉巖, 馬上. 無(wú)網(wǎng)格法的理論及應(yīng)用[J]. 力學(xué)進(jìn)展, 2009,39(1): 1-36. [ZHANG X, LIU Y, MA S. Meshfree methods and their applications[J]. Advances in Mechanics, 2009,39(1) : 1-36. (in Chinese)]

[15] York II A R, Sulsky D, Schreyer H L. The material point method for simulation of thin membranes[J]. International Journal for Numerical Methods in Engineering,1999,44: 1429-1456.

[16] York II A R, Sulsky D, Schreyer H L. Fluid membrane interaction based on the material point method[J]. International Journal for Numerical Methods in Engineering, 2000,48: 901-924.

[17] Beuth L, Wi?cknoski Z, Vermeer P. Solution of quasi-static large-strain problem by the material point method[J]. International Journal for Numerical and Analytical method in Geomechanics, 2011,35(13): 1451-1465.

[18] Andersen S, Andersen L. Modelling of landslides with the material-point method[J]. Computational Geosciences, 2010, 14(1): 137-147.

[19] 史卜濤, 張?jiān)? 張巍. 邊坡穩(wěn)定性分析的物質(zhì)點(diǎn)強(qiáng)度折減法[J]. 巖土工程學(xué)報(bào), 2016, 38(9): 1678-1684. [SHI B T, ZHANG Y, ZHANG W. Strength reduction material point method for slope stability[J]. Chinese Journal of Geotechnical Engineering, 2016,38(9): 1678-1684. (in Chinese)]

[20] 孫玉進(jìn), 宋二祥. 大位移滑坡形態(tài)的物質(zhì)點(diǎn)法模擬[J]. 巖土工程學(xué)報(bào), 2015, 37(7): 1218-1225. [SUN Y J, SONG E X. Simulation of large-displacement landslide by material point method[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(7): 1218-1325. (in Chinese)]

[21] 張巍,史卜濤,施斌,等. 土質(zhì)滑坡運(yùn)動(dòng)全過(guò)程的物質(zhì)點(diǎn)法模擬及其應(yīng)用[J]. 工程地質(zhì)學(xué)報(bào), 2017, 25(3)：815-823. [ZHANG W, SHI B T, SHI B,etal. Run-out process simulation of soil landslide using material point mehtod and its application[J]. Journal of Engineering Geology, 2017，25(3)：815-823. (in Chinese)]

[22] Yerro A, Alonso E, Pinyol N. The material point method for unsaturated soils[J]. Géotechnique,2015, 61(3): 201-217.

[23] BUI H H, FUKAGAWA R, SAKO K,etal. Lagrangian meshfree particles method (SPH) for large deformation and failure flows of geomaterial using elastic-plastic soil constitutive model[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2008, 32(12):1537-1570.

[24] 王寶亮, 彭盛恩, 陳洪凱. 推移式滑坡形成機(jī)制的力學(xué)演繹[J]. 地質(zhì)災(zāi)害與環(huán)境保護(hù), 2010 (2): 74-77. [WANG B L,PENG S E,CHEN H K. Mechanical deduction of formation mechanism for the thrust load caused landslide[J].Journal of Geological Hazards and Environment Preservation, 2010 (2): 74-77.(in Chinese)]

[25] 盧應(yīng)發(fā), 黃學(xué)斌, 劉德富. 推移式滑坡漸進(jìn)破壞機(jī)制及穩(wěn)定性分析[J]. 巖石力學(xué)與工程學(xué)報(bào), 2016, 35(2): 333-345. [LU Y F,HUANG X B,LIU D F. Mechanism and stability analyses of progressive failure of thrust-type landslides[J].Chinese Journal of Rock Mechanics and Engineering, 2016, 35(2): 333-345. (in Chinese)]

[26] 王效寧, 杜永廉, 巫錫勇. 巖質(zhì)順層滑坡的運(yùn)動(dòng)災(zāi)害性模擬試驗(yàn)研究[J]. 水文地質(zhì)工程地質(zhì), 1988, 15(6): 1-5. [WANG X N, DU Y L, WU X Y. Simulation and experimental study on movement severe of bedding rock landslide[J].Hydrogeology amp; Engineering Geology, 1988, 15(6):1-5. (in Chinese)]

責(zé)任編輯

：汪美華

Run-outcharacteristicsimulationofasurcharge-inducedsoillandslideusingthematerialpointmethod

WU Fangdong, ZHANG Wei, SHI Butao, SHI Bin, ZHANG Yun, ZHENG Shuai

(SchoolofEarthSciencesandEngineering,NanjingUniversity,Nanjing,Jiangsu210046,China)

Spoil ground surcharge is one of the main factors inducing landslides in hilly grounds triggered by human engineering activities. The material point method (MPM) belongs to one of the meshless numerical analysis method, which can effectively simulate the material behavior and the run-out characteristics of landslides. Based on the discrete method using the linear shape function, the MUSL solving format and the Drucker-Prager yield criterion, we develop a single-layer and single-phase MPM model to simulate the run-out process of landslides. In comparison with the benchmark experiment results of the process of a sand pile losing its stability and simulated with accumulated dry aluminum bars, the proposed MPM model is verified. The run-out process of the typical scenario of the surcharged-induced soil landslide is simulated using MPM. The related results during the representative moments of the run-out process, including the morphology of slope mass, the plastic strain distribution and sliding velocity evolution trends of reference points, are obtained. The results of the numerical example show that the surcharge-induced soil landslide belongs to the thrust-type landslide, and is of the progressive failure. The whole run-out process can be divided into 4 stages, namely, the slope crest compression, local creep sliding, accelerated sliding and decelerated sliding. Concerning the sliding front, the results of parametric analysis also show that a strong positive correlation property exists between both of the kinetic representative parameters and the surcharge amount, and a strong negative correlation property exits between the kinetic representative parameters and the cohesion or internal friction angle of the soil. It is noted that the kinetic representative parameters include the maximal sliding acceleration, velocity, distance and kinetic energy of the slope mass. Based on 29 typical scenarios, the linear regression equations of all of the above 4 kinetic representative parameters, denoted by the surcharge amount, cohesion and internal friction angel of the soil, are established to predict the disastrous behavior of the surcharge-induced soil landslides.

soil landslide; surcharge; material point method; motion parameters

10.16030/j.cnki.issn.1000-3665.2017.06.19

P642.22

A

1000-3665(2017)06-0126-09

2016-12-27;

2017-04-05

國(guó)家自然科學(xué)基金重點(diǎn)項(xiàng)目(41230636)；江蘇省自然科學(xué)基金項(xiàng)目(BK20160366)；蘇州市科技計(jì)劃項(xiàng)目(SYG201613)

吳方東(1993-)，男，碩士研究生，主要從事地質(zhì)工程大變形數(shù)值模擬研究。E-mail: njuwfd@163.com

張巍(1974-)，男，副教授，從事工程地質(zhì)數(shù)值方法研究。E-mail: wzhang@nju.edu.cn