劉 鑫 王 宇 李典慶

( ①武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室 武漢 430072)

( ②香港城市大學(xué)建筑學(xué)及土木工程學(xué)系 香港 999077)

0 引 言

巖土工程中不確定性不可避免，其中一個(gè)重要的不確定性是土體參數(shù)內(nèi)在的空間變異性( Phoon et al.，1999; Wang et al.，2016b; 趙晶等，2018) 。由于邊坡總是沿著土體中的最軟弱區(qū)域滑動(dòng)，土體參數(shù)的空間變異性與邊坡失穩(wěn)密切關(guān)聯(lián)。大量研究已表明土體參數(shù)的空間變異性對(duì)邊坡失效概率有較大的影響，這些研究大多采用隨機(jī)極限平衡方法( Random Limit Equilibrium，RLEM) ( EI Ramly et al.，2002; Cho，2007; 楊 繼 紅 等，2007; Ji et al.，2012; 蔣水華等，2014; 楊智勇等，2019) 或隨機(jī)有限元方法( Random Finite Element Method，RFEM)( Griffiths et al.，2004，2009; Hicks et al.，2010; 肖特等，2014; Li et al.，2016a) 。然而，由于極限平衡方法或有限元方法的限制，它們難以模擬土體的大變形破壞過(guò)程。實(shí)際上，邊坡失穩(wěn)是一個(gè)涉及啟動(dòng)、滑動(dòng)和堆積3 個(gè)階段的動(dòng)態(tài)過(guò)程。但極限平衡方法僅滿(mǎn)足靜態(tài)平衡方程( 陳祖煜，2003) ，有限元方法在模擬大變形時(shí)可能遇到網(wǎng)格破壞問(wèn)題導(dǎo)致計(jì)算不收斂。邊坡失穩(wěn)的全過(guò)程將決定邊坡的失事后果，因而有必要研究考慮土體參數(shù)的空間變異性情形下邊坡的大變形破壞過(guò)程及其特征。

近年來(lái)，物質(zhì)點(diǎn)法( Material Point Method，MPM) ( Sulsky et al.，1994; 張雄等，2013) 被廣泛應(yīng)用于巖土領(lǐng)域的大變形問(wèn)題中，如邊坡漸進(jìn)破壞( Zabala et al.，2011; Yerro et al.，2016) 、多相介質(zhì)失穩(wěn)( Bandara et al.，2015; Yerro et al.，2015; Wang et al.，2018) 、管涌( Yerro et al.，2017) 、靜力觸探試驗(yàn)( Ceccato et al.，2016) 等。物質(zhì)點(diǎn)法與有限元法類(lèi)似，是一種基于連續(xù)介質(zhì)力學(xué)的數(shù)值模擬工具。物質(zhì)點(diǎn)法將問(wèn)題域同時(shí)劃分為歐拉背景網(wǎng)格和拉格朗日物質(zhì)點(diǎn)，背景網(wǎng)格用于求解動(dòng)量方程，拉格朗日物質(zhì)點(diǎn)攜帶材料和場(chǎng)變量( 如速度、加速度) 信息，使物質(zhì)點(diǎn)可以自由變形，避免了有限元里重新劃分網(wǎng)格帶來(lái)的問(wèn)題( Yerro et al.，2016) 。這使得物質(zhì)點(diǎn)法適于分析邊坡穩(wěn)定性分析和大變形破壞全過(guò)程( 史卜濤等，2016; 王斌等，2017; 張巍等，2017) 。為了考慮土體參數(shù)的空間變異性，已有學(xué)者提出了結(jié)合直接蒙特卡洛模擬的隨機(jī)物質(zhì)點(diǎn)法( Random Material Point Method，RMPM) ( Wang et al.，2016a) ，但計(jì)算消耗較大，阻礙了它在實(shí)際工程中的應(yīng)用。

為此，本文采用一種隨機(jī)極限平衡－物質(zhì)點(diǎn)法分析考慮土體參數(shù)的空間變異性情形下的邊坡大變形破壞模式。該方法同時(shí)利用極限平衡方法簡(jiǎn)單、高效的優(yōu)勢(shì)和物質(zhì)點(diǎn)方法模擬土體大變形破壞的能力，能夠?qū)崿F(xiàn)高效的邊坡可靠度分析和風(fēng)險(xiǎn)評(píng)估。

1 隨機(jī)極限平衡－物質(zhì)點(diǎn)法

隨機(jī)物質(zhì)點(diǎn)法( Wang et al.，2016a) 一般采用直接蒙特卡洛模擬產(chǎn)生N 個(gè)隨機(jī)樣本進(jìn)行不確定分析，邊坡失效概率Pf的表達(dá)式為:

式中，Nf，RMPM為失效樣本總數(shù)，即發(fā)生土體大變形破壞的隨機(jī)樣本數(shù)目; N 為隨機(jī)樣本總數(shù)，建議至少取N=10/Pf保證失效概率估計(jì)的變異系數(shù)小于0.3( Ang et al.，2007) 。

然而，邊坡失穩(wěn)通常是小概率事件，直接蒙特卡洛模擬需要產(chǎn)生大量的樣本，考慮到每次物質(zhì)點(diǎn)法分析本身計(jì)算耗時(shí)較長(zhǎng)，直接蒙特卡洛模擬的計(jì)算總耗時(shí)可能過(guò)高。例如，對(duì)于邊坡失效概率Pf=0.001，至少需要產(chǎn)生10 000 個(gè)隨機(jī)樣本。若采用隨機(jī)物質(zhì)點(diǎn)法完整分析所有樣本，考慮每次分析耗時(shí)15 min，不考慮并行情況下計(jì)算總耗時(shí)將高達(dá)15萬(wàn)分鐘，即2500 h。這對(duì)于工程實(shí)踐是難以接受的。

事實(shí)上，蒙特卡洛模擬的N 個(gè)樣本中僅有Nf，RMPM個(gè)失效樣本，其他的未失效樣本占大部分，對(duì)應(yīng)在物質(zhì)點(diǎn)法分析中沒(méi)有發(fā)生大變形破壞。所以沒(méi)有必要采用物質(zhì)點(diǎn)法分析這些穩(wěn)定的樣本，這樣可以大大提高計(jì)算效率。

本文提出的隨機(jī)極限平衡－物質(zhì)點(diǎn)法將結(jié)合隨機(jī)極限平衡方法的計(jì)算效率優(yōu)勢(shì)和隨機(jī)物質(zhì)點(diǎn)方法模擬大變形破壞的能力，用于高效的邊坡可靠度分析與風(fēng)險(xiǎn)評(píng)估( Liu et al.，2019) 。圖1 繪制了本文方法的框架，包含4 個(gè)部分。第1 部分是采用傳統(tǒng)極限平衡方法和物質(zhì)點(diǎn)方法分別建立確定性分析的邊坡模型，并保證兩個(gè)模型的一致性，即擁有相同的幾何形狀、土體參數(shù)等。第2 部分采用隨機(jī)場(chǎng)理論( Vanmarcke，2010) 作為不確定性模型模擬空間變異的土體，許多文獻(xiàn)( 蔣水華等，2014; 肖特等，2014) 已有介紹，這里不再贅述。第3 部分是采用蒙特卡洛模擬進(jìn)行不確定性分析。第4 部分為結(jié)果后處理，處理隨機(jī)極限平衡方法和隨機(jī)物質(zhì)點(diǎn)方法得到的結(jié)果，如計(jì)算邊坡失效概率，研究邊坡大變形破壞模式等。

圖1 隨機(jī)極限平衡－物質(zhì)點(diǎn)法框架Fig. 1 Illustration of the proposed framework

首先利用蒙特卡洛模擬產(chǎn)生N 個(gè)隨機(jī)場(chǎng)樣本，采用隨機(jī)極限平衡方法進(jìn)行不確定性分析，計(jì)算所有隨機(jī)場(chǎng)樣本的安全系數(shù)值。由于隨機(jī)極限平衡方法計(jì)算簡(jiǎn)便、耗時(shí)短，可以輕易地完成N 個(gè)隨機(jī)場(chǎng)樣本的不確定性分析。接下來(lái)將所有隨機(jī)場(chǎng)樣本根據(jù)安全系數(shù)值從小到大排序，后續(xù)計(jì)算過(guò)程被劃分為k 個(gè)迭代步。對(duì)于第k=1 次迭代，從排序后的樣本序列中取第1 個(gè)隨機(jī)場(chǎng)樣本采用隨機(jī)物質(zhì)點(diǎn)分析，判斷是否發(fā)生大變形破壞，若發(fā)生破壞則令失效樣本總數(shù)目加1，采用式( 1) 計(jì)算失效概率。若不滿(mǎn)足終止條件，則迭代步k=k+1。對(duì)于第k 次迭代( k＞1) ，從排序后的樣本序列中取第k 個(gè)隨機(jī)場(chǎng)樣本采用隨機(jī)物質(zhì)點(diǎn)方法分析，判斷是否發(fā)生大變形破壞，若發(fā)生破壞則更新失效樣本數(shù)目，并采用式(1) 計(jì)算失效概率。迭代計(jì)算的終止條件定義為當(dāng)失效概率收斂到一個(gè)確定的值，即迭代過(guò)程中計(jì)算的失效概率經(jīng)過(guò)kt個(gè)迭代步保持幾乎不變。根據(jù)經(jīng)驗(yàn)，kt可取為1%～10%的樣本總數(shù)目N。

2 邊坡大變形破壞模式

為了驗(yàn)證本文提出的隨機(jī)極限平衡－物質(zhì)點(diǎn)法，采用一個(gè)雙層土坡算例分析。圖2 繪制了該土坡的幾何形狀。上層土高18 m，下層土高10 m，坡比為3︰4，水平距離被延長(zhǎng)到112 m 用于堆積滑動(dòng)體。已有研究采用隨機(jī)有限元方法對(duì)該土坡進(jìn)行可靠度分析( Li et al.，2016b) ，為了保持一致，土體參數(shù)取自該文獻(xiàn)。兩層土的重度均等于19 kN·m－3。兩層土的不排水強(qiáng)度均服從對(duì)數(shù)正態(tài)分布。上層土不排水強(qiáng)度的均值為80 kPa，變異系數(shù)為0.30; 下層土不排水強(qiáng)度的均值為120 kPa，變異系數(shù)為0.30。采用指數(shù)型隨機(jī)場(chǎng)模擬土體不排水強(qiáng)度的空間變異性，有關(guān)細(xì)節(jié)可參考文獻(xiàn)( Li et al.，2015) 。水平相關(guān)長(zhǎng)度為24 m，豎直相關(guān)長(zhǎng)度為2.4 m。圖2 繪制了隨機(jī)場(chǎng)網(wǎng)格，網(wǎng)格大小為1 m×1 m，采用喬列斯基分解離散隨機(jī)場(chǎng)。

圖2 二層土坡幾何形狀Fig. 2 Geometry of a two-layer soil slope

對(duì)于邊坡穩(wěn)定分析，本文與前人一樣假設(shè)滑動(dòng)面為圓弧形并采用簡(jiǎn)化畢肖普法計(jì)算邊坡安全系數(shù)。對(duì)于物質(zhì)點(diǎn)方法分析，整個(gè)區(qū)域被離散為1 m×1 m 的網(wǎng)格，每個(gè)網(wǎng)格均勻分布4 個(gè)物質(zhì)點(diǎn)。采用考慮線(xiàn)性軟化的Drucker-Prager 模型作為本構(gòu)模擬土體 的 彈 塑 性( Bandara et al.，2015; Wang et al.，2018) 。土體的殘余強(qiáng)度取原始強(qiáng)度的50%，硬化模量為－50 kPa( 即軟化) 。物質(zhì)點(diǎn)法總模擬時(shí)長(zhǎng)為20 s，每個(gè)時(shí)間步長(zhǎng)約為0.000 75 s，重力加速度等于9.81 m·s－2。

首先采用直接蒙特卡洛模擬產(chǎn)生40 000 個(gè)隨機(jī)場(chǎng)樣本，采用文中的隨機(jī)極限平衡－物質(zhì)點(diǎn)方法分析每個(gè)隨機(jī)場(chǎng)樣本是否發(fā)生大變形破壞( 閾值為最大相對(duì)位移超過(guò)1 m，見(jiàn)Liu et al. ( 2019) ) 。最終得到一共有520 個(gè)樣本發(fā)生破壞，對(duì)應(yīng)隨機(jī)物質(zhì)點(diǎn)方法的失效概率等于Pf，RLEM=520/40000=1.30×10－2。該失效概率與隨機(jī)極限平衡方法得到的1.12×10－2和隨機(jī)有限元方法( Li et al.，2016b) 得到的1.06×10－2較為吻合。使用本文方法，無(wú)需采用隨機(jī)物質(zhì)點(diǎn)方法分析全部40 000 個(gè)隨機(jī)場(chǎng)樣本，實(shí)際僅耗費(fèi)2512 個(gè)隨機(jī)場(chǎng)樣本即達(dá)到收斂，計(jì)算消耗僅約為原直接蒙特卡洛模擬的6%，大大提高了計(jì)算效率。

對(duì)于本算例，本文方法通過(guò)隨機(jī)物質(zhì)點(diǎn)方法分析可得到不同的破壞模式，根據(jù)滑動(dòng)深度不同，可劃分為淺層、深層、中層破壞模式，另外還有漸進(jìn)破壞模式涉及到多重破壞模式同時(shí)或漸進(jìn)發(fā)生。

2.1 淺層破壞模式

圖3 邊坡淺層破壞模式位移與不排水抗剪強(qiáng)度分布Fig. 3 Displacement and undrained shear strength distributions of a shallow slope failure mode

圖3a、圖3b 分別繪制了一種典型邊坡淺層破壞模式的位移與不排水抗剪強(qiáng)度分布圖。對(duì)于本算例，淺層破壞模式被定義為僅發(fā)生在上層土體。該破壞模式土體沿著上下層土的交界面滑動(dòng)，順著坡面下滑并堆積在坡腳處。圖3 顯示該破壞模式發(fā)生位移約為30 m。極限平衡方法得到的臨界滑動(dòng)面( 即對(duì)應(yīng)安全系數(shù)最小的滑動(dòng)面) 與物質(zhì)點(diǎn)法中實(shí)際發(fā)生的滑動(dòng)體形態(tài)較為一致。如圖3b 所示，上層土靠近坡面部分的不排水抗剪強(qiáng)度較低，低于左側(cè)土體，顯著低于下層土。這種情況容易產(chǎn)生淺層破壞。

2.2 深層破壞模式

圖4a 繪制了一種典型的深層破壞模式。深層破壞模式被定義為滑動(dòng)帶底部低于坡腳高度的破壞模式。該破壞模式滑動(dòng)體積一般大于淺層破壞模式。該滑動(dòng)體整體沿著某一個(gè)圓心轉(zhuǎn)動(dòng)發(fā)生失穩(wěn)，滑動(dòng)體形態(tài)保持較為完整，故被細(xì)分為徑向深層破壞模式。圖4b 繪制了對(duì)應(yīng)的不排水強(qiáng)度分布。下層土底部及坡腳附近的不排水抗剪強(qiáng)度均較低，且下層土的軟弱區(qū)與上層土的軟弱區(qū)相連形成了弧形軟弱帶，可能促使徑向深層滑動(dòng)模式的發(fā)生。

圖5a 繪制了另一種典型的深層破壞模式。與圖4a 繪制的徑向深層破壞模式不同，圖5a 中滑動(dòng)體主要沿著邊坡底部水平切向滑動(dòng)，轉(zhuǎn)動(dòng)較小。且滑動(dòng)體被分為兩個(gè)三角形的滑動(dòng)體，兩個(gè)滑動(dòng)體之間存在明顯的錯(cuò)動(dòng)帶。這可能是由于滑動(dòng)體前后的滑動(dòng)速度不一致造成滑動(dòng)體沿著某個(gè)傾斜的軟弱帶發(fā)生錯(cuò)動(dòng)。這種破壞模式被分類(lèi)為切向深層破壞模式。該破壞模式的滑動(dòng)距離約為24 m，大于圖4a中的徑向深層滑動(dòng)模式。從圖5b 中可以看出下層土的底部存在連續(xù)的水平軟弱帶，使得滑動(dòng)體更易于沿著底部水平運(yùn)動(dòng)。值得注意的是，該破壞模式的形態(tài)不符合圓弧滑動(dòng)面假設(shè)，說(shuō)明圓弧滑動(dòng)面的假設(shè)在考慮土體參數(shù)的空間變異性情形時(shí)可能不成立。

圖4 邊坡徑向深層破壞模式位移與不排水抗剪強(qiáng)度分布Fig. 4 Displacement and undrained shear strength distributions of a radial deep slope failure mode

圖5 邊坡切向深層破壞模式位移與不排水抗剪強(qiáng)度分布Fig. 5 Displacement and undrained shear strength distributions of a tangential deep slope failure mode

2.3 中層破壞模式

圖6 邊坡中層破壞模式位移與不排水抗剪強(qiáng)度分布Fig. 6 Displacement and undrained shear strength distributions of an intermediate slope failure mode

中層破壞模式的滑動(dòng)深度介于淺層和深層破壞模式之間，通常沿著坡腳切向滑動(dòng)。圖6a 繪制了一種典型的中層破壞模式的位移分布。與淺層破壞模式相比，中層破壞模式的滑動(dòng)體積更大; 與深層破壞模式相比，滑動(dòng)體的重力勢(shì)能轉(zhuǎn)換為動(dòng)能后無(wú)需抵抗重力做功，故中層破壞模式的滑動(dòng)距離通常大于淺層和深層破壞模式。圖6b 繪制了該破壞模式的不排水抗剪強(qiáng)度分布。可以看出，坡腳高度之上的軟弱區(qū)與上層土的軟弱區(qū)相連，坡腳以下土體強(qiáng)度相對(duì)較高，從而導(dǎo)致了這種破壞模式。

2.4 漸進(jìn)破壞模式

除此之外，本研究還發(fā)現(xiàn)了漸進(jìn)破壞模式，即多種破壞模式同時(shí)或漸進(jìn)發(fā)生。圖7 繪制了一種典型的漸進(jìn)破壞模式。圖7a 和圖7b 顯示了時(shí)間為第4 s 和第8 s 時(shí)一個(gè)中層破壞模式在滑動(dòng)，滑動(dòng)形態(tài)與極限平衡方法得到的臨界滑動(dòng)面一致。圖7c表明在第10 s 時(shí)一個(gè)新的淺層破壞模式逐漸形成，最終破壞形態(tài)見(jiàn)圖7d，最終滑動(dòng)距離達(dá)到40 m。由于漸進(jìn)破壞模式涉及到多種破壞模式的發(fā)生，多個(gè)模式的滑動(dòng)體積更大，滑動(dòng)距離更遠(yuǎn)，造成的后果通常更為嚴(yán)重。需要指出的是，由于極限平衡方法或有限元方法難以模擬土體的大變形破壞全過(guò)程，它們難以得到后續(xù)的破壞模式，也就無(wú)法用于分析漸進(jìn)破壞模式。說(shuō)明了物質(zhì)點(diǎn)法在模擬邊坡大變形破壞全過(guò)程的優(yōu)勢(shì)。

圖8 繪制了漸進(jìn)破壞模式對(duì)應(yīng)的不排水抗剪強(qiáng)度分布。它的土體不排水抗剪強(qiáng)度分布與中層破壞模式的類(lèi)似，但上層土的強(qiáng)度較低，上層土在第1個(gè)中層破壞模式發(fā)生后失去支撐，導(dǎo)致了第2 個(gè)淺層破壞模式的發(fā)生。

圖7 邊坡漸進(jìn)破壞模式位移演化Fig. 7 The displacement evolution of a progressive slope failure mode

3 結(jié) 論

本文采用隨機(jī)極限平衡－物質(zhì)點(diǎn)法分析了考慮土體參數(shù)空間變異性的邊坡大變形破壞模式，得到了包括淺層、中層、深層和漸進(jìn)一共4 種典型破壞模式，并分析了對(duì)應(yīng)不排水抗剪強(qiáng)度與邊坡破壞模式之間的關(guān)系。得到如下結(jié)論:

(1) 邊坡的破壞模式與土體強(qiáng)度參數(shù)的空間分布密切相關(guān)。不同的土體強(qiáng)度參數(shù)空間分布可能引起不同的破壞模式，導(dǎo)致不同的失事后果。因此，在進(jìn)行邊坡風(fēng)險(xiǎn)評(píng)估與加固時(shí)，應(yīng)重視勘察工作與勘察數(shù)據(jù)的利用，盡可能地量化土體參數(shù)的空間變異性。從而針對(duì)不同的可能破壞模式，做出更有效的加固措施。

( 2) 與極限平衡方法和有限元方法相比，物質(zhì)點(diǎn)法具備土體大變形分析的能力，既可用于分析邊坡是否失穩(wěn)，又具備模擬邊坡土體大變形破壞的全過(guò)程的能力。本研究提出的隨機(jī)極限平衡－物質(zhì)點(diǎn)法可以同時(shí)利用極限平衡方法的計(jì)算效率優(yōu)勢(shì)和物質(zhì)點(diǎn)方法的土體大變形分析能力，可以實(shí)現(xiàn)高效地邊坡風(fēng)險(xiǎn)評(píng)估。

( 3) 在考慮土體參數(shù)的空間變異性時(shí)，多數(shù)情況下極限平衡方法得到的圓弧滑動(dòng)面與物質(zhì)點(diǎn)法得到的滑動(dòng)體形態(tài)較為一致，但在某些情況下圓弧滑動(dòng)面假設(shè)可能失效( 如存在較大連通區(qū)域的軟弱帶時(shí)) 。此外，極限平衡方法難以用于考慮漸進(jìn)破壞模式，即多個(gè)破壞模式同時(shí)或漸進(jìn)發(fā)生情形。

( 4) 本文所采用的物質(zhì)點(diǎn)方法僅分析單相土體，未來(lái)仍需拓展到多相用于考慮土－水的耦合作用，考慮孔隙水壓力等更復(fù)雜工況。本文提出的隨機(jī)極限平衡－物質(zhì)點(diǎn)法仍然可用于邊坡可靠度分析和風(fēng)險(xiǎn)評(píng)估。

( 5) 可以看出與傳統(tǒng)有限元方法相比，物質(zhì)點(diǎn)法適于邊坡和其他巖土問(wèn)題的大變形模擬和破壞分析，可以利用現(xiàn)有的土體本構(gòu)模型，但物質(zhì)點(diǎn)法在巖土領(lǐng)域的研究尚處于初期階段，存在一些問(wèn)題有待完善，如考慮土－水的多相耦合作用、隱式物質(zhì)點(diǎn)法、應(yīng)力波動(dòng)和與有限元相比計(jì)算精度不高等問(wèn)題。