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

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

堆載誘發(fā)型土質滑坡運動特征物質點法模擬

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

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

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

土質滑坡；堆載；物質點法；運動參數(shù)

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

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

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

1 物質點模型

1.1基本描述

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

圖1 物質點法示意圖Fig.1 Schematic diagram of the material point method

不考慮熱交換，物質點法控制方程包括虛功方程、動量方程、幾何方程、本構方程以及邊界條件方程等，具體參見文[19]、[21]原理部分。

1.2物質點離散

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

式中：np——質點總數(shù)；

mp——質點p的質量；

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

xi——空間坐標；

xip——質點p的坐標。

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

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

式中，用帶有下標p的量來表示物質點攜帶的相應物理量(下同)，同時引入假想邊界層厚度h將虛功方程式左端最后一項轉化為體積分。

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

其中，xip、uip、uip，j、δuip分別為質點p的坐標、位移、位移的導數(shù)以及虛位移，δuiI為結點的虛位移，用帶有下標I的量來表示網(wǎng)格結點的變量。

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

?

其中

為第I個結點在i方向上的動量；

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

為結點內力；

為結點外力。

當采用集中質量矩陣時，

式(9)簡化為：

則背景網(wǎng)格結點的運動方程為：

?

1.3求解格式

求解動量方程時，需要將某時刻質點的質量、動量等信息映射到背景網(wǎng)格，以計算背景網(wǎng)格結點的質量和動量等。動量方程采用顯式時間積分求解。網(wǎng)格結點的速度變化量可以用于計算質點的應變增量，從而對質點進行應力更新。

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

1.4本構模型

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

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

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

式中：τ——等效剪應力；

σm——球應力；

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

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

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

加號表示DP屈服面在π平面上內接MC屈服面，減號則表示外接。

2 模型驗證

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

建立以上鋁棒堆積物的單套單相物質點模型，示意如圖2。

圖2 物質點模型(80 000個物質點)Fig.2 MPM model (80 000 particles)

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

圖3 失穩(wěn)形態(tài)與塑性應變云圖Fig.3 Instability morphology and plastic strain contour map

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

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

圖5 模擬與實驗結果比較Fig.5 Comparison of the results of simulation and experiment

3 土質滑坡物質點法模擬

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

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

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

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

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

圖7 邊坡物質點模型(10 975個物質點)Fig.7 Slope MPM model (10 975 particles)

4 滑坡物質運動特征

4.1荷載設置

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

圖8 加載模式Fig.8 Loading mode

圖9 初始應力場Fig.9 Initial stress field

4.2滑坡全過程分析

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

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

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

圖11 滑速演化趨勢Fig.11 Slide speed evolution trends

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

(1)Ⅰ坡頂壓縮

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

(2)Ⅱ局部蠕滑

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

(3)Ⅲ加速滑動

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

(4)Ⅳ減速滑動

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

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

5 參數(shù)分析

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

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

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

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

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

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

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

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

6 結論

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

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

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

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

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

[1] 張茂省, 李同錄. 黃土滑坡誘發(fā)因素及其形成機理研究[J]. 工程地質學報, 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]. 水文地質工程地質, 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]. 石家莊鐵道大學學報(自然科學版), 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]. 中國地質災害與防治學報, 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進行堆載作用下邊坡穩(wěn)定分析[J]. 勘察科學技術, 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]. 巖石力學與工程學報, 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] 張偉鋒, 黃潤秋, 裴向軍. 大光包滑坡運動特征及其過程分析[J]. 工程地質學報, 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] 羅鋒, 柴波, 方恒, 等. 順層巖質滑坡運動過程數(shù)值模擬研究[J]. 水文地質工程地質, 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] 殷坤龍, 姜清輝, 汪洋. 滑坡運動過程仿真分析[J]. 地球科學(中國地質大學學報), 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] 張雄, 廉艷平, 劉巖,等. 物質點法[M]. 北京：清華大學出版社, 2013. [ZHANG X, LIAN Y P, LIU Y,etal. Material point method[M]. Beijing: Tsinghua University Press, 2013. (in Chinese)]

[14] 張雄, 劉巖, 馬上. 無網(wǎng)格法的理論及應用[J]. 力學進展, 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] 史卜濤, 張云, 張巍. 邊坡穩(wěn)定性分析的物質點強度折減法[J]. 巖土工程學報, 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] 孫玉進, 宋二祥. 大位移滑坡形態(tài)的物質點法模擬[J]. 巖土工程學報, 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] 張巍,史卜濤,施斌,等. 土質滑坡運動全過程的物質點法模擬及其應用[J]. 工程地質學報, 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]. 地質災害與環(huán)境保護, 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] 盧應發(fā), 黃學斌, 劉德富. 推移式滑坡漸進破壞機制及穩(wěn)定性分析[J]. 巖石力學與工程學報, 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] 王效寧, 杜永廉, 巫錫勇. 巖質順層滑坡的運動災害性模擬試驗研究[J]. 水文地質工程地質, 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)]

責任編輯

：汪美華

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

國家自然科學基金重點項目(41230636)；江蘇省自然科學基金項目(BK20160366)；蘇州市科技計劃項目(SYG201613)

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

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