土質(zhì)滑坡失穩(wěn)、運動及沖擊壓力物質(zhì)點法模擬研究*

王 升 曾 鵬 李天斌 孫小平

(地質(zhì)災害防治與地質(zhì)環(huán)境保護國家重點實驗室(成都理工大學)，成都 610059，中國)

0 引 言

滑坡災害是全世界發(fā)生最為頻繁的災害之一，備受國際關(guān)注。其中流動型滑坡往往具有速度快、距離長的特點，對生命財產(chǎn)構(gòu)成極大威脅。然而目前對于此類滑坡的預測尚面臨諸多挑戰(zhàn)。已有諸多學者(Wartman et al.，2005；Lin et al.，2006；Tohari et al.，2007；杜婷婷等，2018；李楠等，2019；孫萍等，2019；趙曉彥等，2019)采用物理模型試驗來研究滑坡的失穩(wěn)、運動和沖擊過程。但室內(nèi)試驗往往受到尺度效應的影響難以復現(xiàn)滑坡的原位狀態(tài)，而現(xiàn)場試驗研究又成本巨大，因此通過數(shù)值分析方法深入研究滑坡失穩(wěn)、運動與沖擊過程對滑坡風險評價和管理具有重要意義。

目前國內(nèi)外用于模擬滑坡問題的數(shù)值方法主要分為連續(xù)介質(zhì)力學方法和非連續(xù)介質(zhì)力學方法兩大類。傳統(tǒng)的連續(xù)介質(zhì)力學方法(比如有限元方法，F(xiàn)EM)在模擬大變形問題時，由于網(wǎng)格畸形，使得計算意外結(jié)束。非連續(xù)介質(zhì)力學方法(比如離散元方法，DEM)雖然對網(wǎng)格的依賴性較低，但是計算效率較低。因此，很多無網(wǎng)格方法相繼出現(xiàn)，例如光滑粒子流體動力學方法(smoothed particle hydrodynamics，SPH)、粒子有限元方法(particle finite element method，PFEM)、物質(zhì)點法(material point method，MPM)等。其中：物質(zhì)點法能準確、高效地捕捉材料的大變形運動過程。相比SPH，物質(zhì)點法對模型邊界條件的施加更加簡單。相比PFEM，物質(zhì)點法的計算代價更低(Soga et al.，2016；王斌等，2017)。

物質(zhì)點法是Sulsky et al.(1994)于1994年提出來的。20多年來，物質(zhì)點法在模擬滑坡問題上得到了較廣泛的應用。王斌等(2017)將物質(zhì)點法應用于邊坡的穩(wěn)定性評價中，研究表明：物質(zhì)點法能穩(wěn)定、高效地評價邊坡的初始破壞以及誘發(fā)的眾多次生災害。史卜濤等(2016)將物質(zhì)點強度折減法運用到邊坡的穩(wěn)定性分析中，為邊坡穩(wěn)定性分析提供了嶄新的思路。謝艷芳等(2018)基于物質(zhì)點法對新磨村滑坡進行了二維的動力演化分析。孫玉進等(2018)、張巍等(2017)、Shi et al.(2019)采用物質(zhì)點法構(gòu)建了深圳滑坡的二維模型，從不同角度對該滑坡的運動過程和堆積特征等進行了模擬研究。Yerro et al.(2016)采用物質(zhì)點法分析了由孔隙水壓力增加引起的邊坡穩(wěn)定問題。Andersen et al.(2010)利用物質(zhì)點法對邊坡失穩(wěn)過程及其對建筑穩(wěn)定性的影響進行了分析。Xu et al.(2018)采用三維物質(zhì)點法對中國紅石巖滑坡的破壞和運動進行了模擬分析。在地表運動物質(zhì)沖擊壓力計算方面，Armanini(1997)提出泥石流沖擊壓力由動壓力和平均靜止壓力兩部分組成，獲得了廣泛應用。Ceccato et al.(2017)在物質(zhì)點法的框架下，基于線性脈沖原理(Jóhannesson et al.，2009)，提出了土質(zhì)滑坡沖擊壓力的計算公式。Li et al.(2018)通過建立物質(zhì)點法剛?cè)峤佑|模型，分析了泥石流沿山坡流下沖擊擋土墻的理想情況。Feng et al.(2019)采用耦合SPH和FEM模擬了深圳滑坡的三維運動過程和對建筑物的沖擊作用。

滑坡災害的風險評價需要了解滑坡從穩(wěn)定、失穩(wěn)到運動、沖擊和停積的全過程狀態(tài)和特征?，F(xiàn)有滑坡數(shù)值模擬研究大多關(guān)注滑坡的失穩(wěn)和運動過程，或假定失穩(wěn)條件下進行滑坡運動和沖擊特征數(shù)值模擬研究。同時，這些研究多以特殊的滑坡幾何模型為基礎(chǔ)，其工程適用性尚不明確。本文擬在物質(zhì)點法(MPM)框架下對深圳滑坡失穩(wěn)、運動和沖擊特征開展數(shù)值模擬，并與實際觀察到的滑坡滑面、運動及堆積特征進行逐一對比，提供一套再現(xiàn)滑坡穩(wěn)定—失穩(wěn)—運動—沖擊—停積全生命周期的數(shù)值模擬方法，從而為單體滑坡風險定量評估提供物理力學模型基礎(chǔ)。

1 物質(zhì)點法基本理論

1.1 物質(zhì)點法基本原理

物質(zhì)點法結(jié)合拉格朗日法和歐拉法的雙重優(yōu)勢，將連續(xù)體材料離散為一系列帶有質(zhì)量、應變等所有信息的物質(zhì)點。背景網(wǎng)格不攜帶材料信息，僅用于求解控制方程。如圖1(Soga et al.，2016)所示，在每一個時間步中，物質(zhì)點的所有信息通過形函數(shù)映射給背景網(wǎng)格結(jié)點，求解后，網(wǎng)格結(jié)點又將所有計算后的信息映射回各物質(zhì)點，從而得到這些物質(zhì)點在下一個時刻的初始信息。在每個時間步結(jié)束后，物質(zhì)點法會丟棄掉畸形網(wǎng)格，因而在模擬大變形上具有顯著優(yōu)勢。

1.2 物質(zhì)點法控制方程及弱形式

1.2.1 更新拉格朗日型控制方程

(1)

邊界條件：

(2)

初始條件：

ui(X，0)=ui0(X)

(3)

ui，t(X，0)=ui，t0(X)

(4)

1.2.2 更新拉格朗日式的弱形式

(5)

1.3 物質(zhì)點離散

物質(zhì)點法將材料進行離散，用物質(zhì)點來代替相應的離散體，且采用背景網(wǎng)格進行計算。而背景網(wǎng)格形函數(shù)的導數(shù)在網(wǎng)格之間不連續(xù)，當物質(zhì)點跨越背景網(wǎng)格時會產(chǎn)生數(shù)值噪聲。Bardenhagen et al.(2004)提出了基于特征函數(shù)χp(x)的廣義插值物質(zhì)點法(generalized interpolation material point method，GIMP)，從而有效降低了因物質(zhì)點跨越背景網(wǎng)格而產(chǎn)生的數(shù)值噪聲，進一步提高了計算精度。

χp(x)用于確定質(zhì)點所占據(jù)的空間區(qū)域，χp相當于質(zhì)點在空間所有質(zhì)點中所占據(jù)的體積份數(shù)，故滿足：

(6)

物質(zhì)點的體積可以表示為：

(7)

式中：Ωp為物質(zhì)點p在現(xiàn)時構(gòu)形中的所占據(jù)的區(qū)域。

(8)

(9)

(10)

從而，可將弱形式方程(5)表示成如下形式：

(11)

利用背景網(wǎng)格將虛位移近似為：

(12)

將式(12)代入式(11)，得到

(13)

式中：

(14)

(15)

(16)

(17)

(18)

式中：NI(x)為背景網(wǎng)格結(jié)點I的形函數(shù)。

1.4 本構(gòu)模型

基巖采用彈性模型，土體選用Drucker-Prager屈服準則下的理想彈塑性模型。在彈性階段，應力-應變關(guān)系為線性，即：

σij=λεkkδij+2μεij

(19)

式中：εkk為體積應變；δij為柯式δ(Kronecher-δ)；εij為應變張量；λ、μ為拉梅常數(shù)。彈性模量E、泊松比v滿足：

(20)

(21)

在剪切方向上，采用Drucker-Prager屈服面，即：

(22)

式中：J2為偏應力張量的第二不變量；qφ為摩擦系數(shù)；I1為應力張量的第一不變量；kc為純剪切狀態(tài)時的屈服應力。材料常數(shù)qφ和kc與內(nèi)摩擦角φ和黏聚力c的關(guān)系為：

(23)

(24)

2 深圳滑坡物質(zhì)點法模型構(gòu)建

2.1 滑坡基本概況

2015年12月20日深圳光明新區(qū)發(fā)生特大滑坡，摧毀33棟建筑物，共造成77人死亡和失蹤。滑坡主滑方向N20°W，運動距離為1100im。滑坡平面圖如圖2所示。深圳滑坡主要由底部花崗巖基巖和上部人工堆填體組成。其中堆填體可細分為第一層建筑渣土和第二層采石場棄土。圖3為深圳滑坡A-A′剖面圖(Yin et al.，2016)，滑坡整體堆積較緩，堆積坡角為6.1°?；露逊e區(qū)最大堆積厚度為21.3im，平均堆積厚度為12.9im。從圖3中可以看出在滑坡堆積區(qū)范圍內(nèi)存在建筑物1#～4#，本文將以這4排建筑物所處的位置，進行滑坡運動物質(zhì)沖擊壓力和沖擊荷載計算。

2.2 物質(zhì)點法模型構(gòu)建

2.2.1 物質(zhì)點模型

本文基于MATLAB開發(fā)了物質(zhì)點法的前處理程序。背景計算網(wǎng)格采用八節(jié)點六面體單元，每個網(wǎng)格含有8個物質(zhì)點，網(wǎng)格尺寸為2im，物質(zhì)點間距為1im?？傆嬌晌镔|(zhì)點127i118個，其中72i996個物質(zhì)點代表基巖(圖4中藍色點)，54i122個物質(zhì)點代表土體(圖4中紅色點)。

2.2.2 計算參數(shù)

模型計算參數(shù)的取值采用Yin et al.(2016)對深圳滑坡開展的環(huán)剪試驗結(jié)果，如表1(Yin et al.，2016)所示。根據(jù)有效應力原理：

表1 深圳滑坡物質(zhì)點法輸入?yún)?shù)(Yin et al.，2016)Table 1 Parameters used in MPM(Yin et al.，2016)

τ=c′+(σ-u)tanφ′

(25)

式中：σ為總應力；u為孔隙水壓力；c′為有效黏聚力；φ′為有效內(nèi)摩擦角。土體采用理想彈塑性模型，故峰值抗剪強度與殘余抗剪強度取值相同。因此，c′和φ′取Yin et al.(2016)環(huán)剪試驗結(jié)果0ikPa和24°。

考慮地下水的作用，將孔隙水壓力系數(shù)表示為：

ru=u/σ

(26)

那么式(25)可改寫為：

τ=0+σ(1-ru)tanφ′

(27)

采用高楊(2018)給出的深圳滑坡孔隙水壓力系數(shù)ru=0.75。按照式(28)，對φ′進行等效替換(Hungr et al.，2009)，得到等效內(nèi)摩擦角φe。

τ=σ(1-ru)tanφ′=σtanφe

(28)

3 深圳滑坡模擬結(jié)果分析

核心計算程序采用開源代碼MPM3D-F90(張雄等，2013)。首先采用動態(tài)松弛技術(shù)(Oakley et al.，1995)生成邊坡初始應力場，時間為10is。在物質(zhì)點法中，重力在初始時刻加載，會產(chǎn)生動力響應，故而需要在邊坡中生成初始應力場，以滿足靜力平衡。在初始應力場的基礎(chǔ)上，修改土體黏聚力c為0(孫玉進等，2015)，進行深圳滑坡動力學模擬，動力計算總時間為70is，時間步長為0.849×10-4s。

3.1 滑面識別

采用經(jīng)典的K均值聚類算法(Huang et al.，2013)進行滑面識別。對于一組觀察值(d1，d2，…，dN)，且每個觀察值都是一個m維向量，K均值聚類算法的目的就是通過不斷地使方差函數(shù)(J)最小化從而將N個觀察值聚成K類。

(29)

式中：ck(k=1，…K)是聚類中心。K均值聚類算法在每個迭代計算周期內(nèi)分為以下兩步：

將物質(zhì)點的相對位移作為觀察值(d1，d2，…，dN)，即采用一維聚類方法，N為物質(zhì)點總數(shù)。滑面提取結(jié)果如圖5所示，與實際滑面對比信息見表2。

表2 模擬滑面與實際滑面對比信息一覽表Table 2 Information of simulated slip surface and observed slip surface for comparison

計算結(jié)果表明，基于物質(zhì)點法模擬提取的滑面與實際滑面吻合較好。模擬滑面的最大深度與實際滑面的最大深度基本一致。模擬滑面的后緣—剪出口水平距離以及高差均與實際滑面相差較小。實際滑面在中部及剪出口處均為近水平狀，而模擬滑面稍有起伏。實際滑裂面后緣較陡，而模擬的滑裂面后緣較緩。值得注意的是，實際滑面是Yin et al.(2016)通過6個鉆孔數(shù)據(jù)推測而來的，且相鄰鉆孔的最大距離達到了150im，因而實際滑面并不能完全定量反映深圳滑坡的真實情況，僅能作為參考。

3.2 滑坡運動特征分析

3.2.1 運動過程

深圳滑坡的動態(tài)數(shù)值模擬過程如圖6所示。坡趾最開始產(chǎn)生最大剪切應變帶，并不斷延伸至貫通整個滑裂面(第13is左右)。其后剪出口物質(zhì)點開始加速運動，進而促進滑坡整體加速運動(13～20is)。后緣物質(zhì)成階梯狀不斷往下潰散，經(jīng)高速運動后逐漸停積并形成堆積體(20～60is)。模擬結(jié)果表明，后緣物質(zhì)點運動距離較短，約為300im，中部物質(zhì)點運動距離約為350im，剪出口處物質(zhì)點運動距離最長，達500im。

3.2.2 運動速度

圖4中特征物質(zhì)點A、B、C的運動速度模擬結(jié)果如圖7所示。特征物質(zhì)點運動速度整體上呈現(xiàn)先增大后減小的趨勢，不同位置的物質(zhì)點運動速度差別較大。滑坡運動物質(zhì)的峰值速度從前緣(物質(zhì)點C)到中部(物質(zhì)點B)再到后緣(物質(zhì)點A)呈現(xiàn)遞減的規(guī)律，這表明在滑坡的運動過程中，前緣的運動物質(zhì)速度最快，中部的運動物質(zhì)速度較快，而后緣的運動物質(zhì)速度相對較慢?；虑熬壙拷舫隹谔幍奈镔|(zhì)點C的速度增加最快，第20is達到最大速度16.7im·s-1。而中部物質(zhì)點B出現(xiàn)峰值速度12.3im·s-1的時間晚于物質(zhì)點C，反映了中部物質(zhì)點的速度遞增較慢。后緣靠近滑裂面處的物質(zhì)點A的運動速度在第10is達到峰值速度7.9im·s-1之后呈現(xiàn)震蕩變化，反映了深圳滑坡后緣物質(zhì)的復雜運動特性。A點峰值速度出現(xiàn)的時間早于C點與B點，這表明滑坡后緣物質(zhì)階梯狀潰散的快速性，同時體現(xiàn)了滑坡運動物質(zhì)能量傳遞的過程。

3.3 滑坡堆積特征分析

深圳滑坡堆積特征的模擬結(jié)果如圖8和表3所示?？傮w而言，物質(zhì)點法模擬的深圳滑坡堆積形態(tài)接近實際的滑坡堆積形態(tài)。模擬的堆積坡角5.2°與實際堆積坡角6.1°基本一致。模擬的平均堆積厚度為10.3im，實際的平均堆積厚度為12.9im，兩者基本吻合。模擬的最大堆積厚度(25.2im)與實際最大堆積厚度(21.3im)也基本一致。但模擬結(jié)果的運動距離(1260im)比實際運動距離(1100im)略大，主要由于本文將三維滑坡運動過程簡化為二維模型，未考慮滑坡運動物質(zhì)的橫向擴展特征。通過與實際堆積對比，初步驗證了物質(zhì)點法在模擬滑坡大變形運動過程及其停積特征上的適用性。

表3 深圳滑坡模擬堆積與實際堆積對比信息一覽表Table 3 Information of simulated accumulation and actual accumulation for comparison

3.4 運動物質(zhì)沖擊壓力與沖擊荷載計算

3.4.1 計算方法

Armanini(1997)認為泥石流對建(構(gòu))筑結(jié)構(gòu)造成的總沖擊壓力理論上可表達為動壓力與平均靜止壓力之和，其本質(zhì)是將泥石流考慮為不可壓縮流體，采用動量守恒定理推導而來的。本文在考慮土質(zhì)滑坡失穩(wěn)、運動的基礎(chǔ)上，為建立滑坡運動物質(zhì)對其運移路徑上任意點的沖擊壓力計算理論，將滑坡運動物質(zhì)考慮為不可壓縮流體，采用Armanini(1997)的沖擊壓力(Pa)計算公式：

Pa=Pd+Ps

(30)

Pd=ρv2

(31)

Ps=0.5ρgh

(32)

式中：Pd為動壓力；Ps為平均靜止壓力；ρ為滑坡運動物質(zhì)的密度；v為沖擊力作用面上的所有物質(zhì)點的平均速度；h為物質(zhì)點到達某點時的堆積厚度；g為重力加速度。

取單位寬度，相對應地，滑坡運動物質(zhì)的總沖擊荷載為沖擊力與靜荷載之和：

Fa=Pdh+Psh=ρhv2+0.5ρgh2

(33)

3.4.2 計算結(jié)果

針對圖3所示的建筑物所處的位置1#～4#進行沖擊壓力的計算，結(jié)果如圖9和表4所示。結(jié)果表明滑坡運動物質(zhì)在點1#～4#所產(chǎn)生的沖擊壓力都表現(xiàn)為先快速增加至峰值再逐漸降低至穩(wěn)定水平，反映了到達點1#～4#的物質(zhì)點平均運動速度為先增加再遞減。其中在點1#處產(chǎn)生的峰值沖擊壓力明顯高于其他位置處的峰值沖擊壓力，這表明物質(zhì)點沖擊1#處的時間段正好是整個滑坡運動速度達到最大的時間段。沖擊壓力曲線也反映了滑坡運動物質(zhì)對運移路徑上某確定位置的沖擊壓力分為動壓力和靜壓力，動壓力占主導地位，當滑坡運動停止時，沖擊壓力退化為靜壓力。

表4 在某點處產(chǎn)生的沖擊壓力計算結(jié)果Table 4 Computed impact pressure at certain locations

峰值沖擊壓力隨計算點與滑坡體之間距離的增加而減小。對距離滑坡后緣滑裂面最近的點1#(660im)造成的峰值沖擊壓力達1667.9ikPa，對距離滑坡后緣滑裂面最遠的點4#(983im)造成的峰值沖擊壓力為490.8ikPa。根據(jù)Zanchetta et al.(2004)的研究結(jié)果，當框架結(jié)構(gòu)的建筑物所受沖擊壓力超過100ikPa時，會遭受嚴重損毀。本文的計算結(jié)果在一定程度上解釋了深圳滑坡導致建筑物大面積損毀和淤埋的原因。

沖擊荷載的計算結(jié)果如圖10和表5所示，滑坡運動物質(zhì)在點1#～4#產(chǎn)生的沖擊荷載均表現(xiàn)為先增加后降低至穩(wěn)定水平的特點，這是由沖擊壓力曲線的特點決定的。在點1#處產(chǎn)生的峰值沖擊荷載的大小是其穩(wěn)定沖擊荷載大小的5倍，這再次表明，沖擊荷載由沖擊力和靜荷載兩部分組成，在滑坡運動過程中，沖擊荷載受運動速度的影響很大，當滑坡運動停止時，沖擊荷載退化為靜荷載。

表5 在某點處產(chǎn)生的沖擊荷載計算結(jié)果Table 5 Computed impact load at certain locations

4 結(jié) 論

本文將二維物質(zhì)點法與Drucker-Prager屈服準則結(jié)合，建立了土質(zhì)滑坡失穩(wěn)、運動、堆積及沖擊壓力的數(shù)值模擬方法，并成功運用到深圳滑坡中。模擬得到了深圳滑坡的滑面、運動距離、運動速度、堆積厚度和沖擊壓力等關(guān)鍵信息。綜合本文的研究，主要得出以下幾點結(jié)論：

(1)物質(zhì)點法不僅可以模擬滑坡運動過程與停積特征，同樣可以模擬滑坡失穩(wěn)過程并提取滑裂面，是再現(xiàn)滑坡穩(wěn)定—失穩(wěn)—運動—沖擊—停積全生命周期的可靠數(shù)值模擬方法。

(2)采用聚類算法提取的深圳滑坡滑面與實際滑面吻合較好，驗證了此類算法在基于物質(zhì)點法的滑坡失穩(wěn)過程數(shù)值模擬中的適用性。

(3)物質(zhì)點法可再現(xiàn)深圳滑坡的漸進式潰散失穩(wěn)特征。模擬結(jié)果顯示滑坡前緣物質(zhì)點運動距離最大，其運動速度明顯高于中部和后緣的物質(zhì)點。

(4)物質(zhì)點法模擬的深圳滑坡運動距離(1260im)比實際運動距離(1100im)略大，整體堆積特征與現(xiàn)場實際觀察結(jié)果吻合較好。

(5)計算得到的滑坡運動物質(zhì)對其運移路徑上4個點的沖擊壓力與沖擊荷載整體上呈現(xiàn)出先快速增大至峰值，然后逐漸衰減為靜壓力或靜荷載的特征，且峰值沖擊壓力與峰值沖擊荷載隨計算點的位置與滑坡體距離的增長而減小。該計算結(jié)果可為滑坡風險定量評估提供重要數(shù)據(jù)支撐。