劉 鑫 王 宇 李典慶

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

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

0 引 言

巖土工程中不確定性不可避免，其中一個重要的不確定性是土體參數(shù)內(nèi)在的空間變異性( Phoon et al.，1999; Wang et al.，2016b; 趙晶等，2018) 。由于邊坡總是沿著土體中的最軟弱區(qū)域滑動，土體參數(shù)的空間變異性與邊坡失穩(wěn)密切關(guān)聯(lián)。大量研究已表明土體參數(shù)的空間變異性對邊坡失效概率有較大的影響，這些研究大多采用隨機(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) 。然而，由于極限平衡方法或有限元方法的限制，它們難以模擬土體的大變形破壞過程。實(shí)際上，邊坡失穩(wěn)是一個涉及啟動、滑動和堆積3 個階段的動態(tài)過程。但極限平衡方法僅滿足靜態(tài)平衡方程( 陳祖煜，2003) ，有限元方法在模擬大變形時可能遇到網(wǎng)格破壞問題導(dǎo)致計(jì)算不收斂。邊坡失穩(wěn)的全過程將決定邊坡的失事后果，因而有必要研究考慮土體參數(shù)的空間變異性情形下邊坡的大變形破壞過程及其特征。

近年來，物質(zhì)點(diǎn)法( Material Point Method，MPM) ( Sulsky et al.，1994; 張雄等，2013) 被廣泛應(yīng)用于巖土領(lǐng)域的大變形問題中，如邊坡漸進(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)法與有限元法類似，是一種基于連續(xù)介質(zhì)力學(xué)的數(shù)值模擬工具。物質(zhì)點(diǎn)法將問題域同時劃分為歐拉背景網(wǎng)格和拉格朗日物質(zhì)點(diǎn)，背景網(wǎng)格用于求解動量方程，拉格朗日物質(zhì)點(diǎn)攜帶材料和場變量( 如速度、加速度) 信息，使物質(zhì)點(diǎn)可以自由變形，避免了有限元里重新劃分網(wǎng)格帶來的問題( Yerro et al.，2016) 。這使得物質(zhì)點(diǎn)法適于分析邊坡穩(wěn)定性分析和大變形破壞全過程( 史卜濤等，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ù)的空間變異性情形下的邊坡大變形破壞模式。該方法同時利用極限平衡方法簡單、高效的優(yōu)勢和物質(zhì)點(diǎn)方法模擬土體大變形破壞的能力，能夠?qū)崿F(xiàn)高效的邊坡可靠度分析和風(fēng)險評估。

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

隨機(jī)物質(zhì)點(diǎn)法( Wang et al.，2016a) 一般采用直接蒙特卡洛模擬產(chǎn)生N 個隨機(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ì)算耗時較長，直接蒙特卡洛模擬的計(jì)算總耗時可能過高。例如，對于邊坡失效概率Pf=0.001，至少需要產(chǎn)生10 000 個隨機(jī)樣本。若采用隨機(jī)物質(zhì)點(diǎn)法完整分析所有樣本，考慮每次分析耗時15 min，不考慮并行情況下計(jì)算總耗時將高達(dá)15萬分鐘，即2500 h。這對于工程實(shí)踐是難以接受的。

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

本文提出的隨機(jī)極限平衡－物質(zhì)點(diǎn)法將結(jié)合隨機(jī)極限平衡方法的計(jì)算效率優(yōu)勢和隨機(jī)物質(zhì)點(diǎn)方法模擬大變形破壞的能力，用于高效的邊坡可靠度分析與風(fēng)險評估( Liu et al.，2019) 。圖1 繪制了本文方法的框架，包含4 個部分。第1 部分是采用傳統(tǒng)極限平衡方法和物質(zhì)點(diǎn)方法分別建立確定性分析的邊坡模型，并保證兩個模型的一致性，即擁有相同的幾何形狀、土體參數(shù)等。第2 部分采用隨機(jī)場理論( 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 個隨機(jī)場樣本，采用隨機(jī)極限平衡方法進(jìn)行不確定性分析，計(jì)算所有隨機(jī)場樣本的安全系數(shù)值。由于隨機(jī)極限平衡方法計(jì)算簡便、耗時短，可以輕易地完成N 個隨機(jī)場樣本的不確定性分析。接下來將所有隨機(jī)場樣本根據(jù)安全系數(shù)值從小到大排序，后續(xù)計(jì)算過程被劃分為k 個迭代步。對于第k=1 次迭代，從排序后的樣本序列中取第1 個隨機(jī)場樣本采用隨機(jī)物質(zhì)點(diǎn)分析，判斷是否發(fā)生大變形破壞，若發(fā)生破壞則令失效樣本總數(shù)目加1，采用式( 1) 計(jì)算失效概率。若不滿足終止條件，則迭代步k=k+1。對于第k 次迭代( k＞1) ，從排序后的樣本序列中取第k 個隨機(jī)場樣本采用隨機(jī)物質(zhì)點(diǎn)方法分析，判斷是否發(fā)生大變形破壞，若發(fā)生破壞則更新失效樣本數(shù)目，并采用式(1) 計(jì)算失效概率。迭代計(jì)算的終止條件定義為當(dāng)失效概率收斂到一個確定的值，即迭代過程中計(jì)算的失效概率經(jīng)過kt個迭代步保持幾乎不變。根據(jù)經(jīng)驗(yàn)，kt可取為1%～10%的樣本總數(shù)目N。

2 邊坡大變形破壞模式

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

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

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

首先采用直接蒙特卡洛模擬產(chǎn)生40 000 個隨機(jī)場樣本，采用文中的隨機(jī)極限平衡－物質(zhì)點(diǎn)方法分析每個隨機(jī)場樣本是否發(fā)生大變形破壞( 閾值為最大相對位移超過1 m，見Liu et al. ( 2019) ) 。最終得到一共有520 個樣本發(fā)生破壞，對應(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較為吻合。使用本文方法，無需采用隨機(jī)物質(zhì)點(diǎn)方法分析全部40 000 個隨機(jī)場樣本，實(shí)際僅耗費(fèi)2512 個隨機(jī)場樣本即達(dá)到收斂，計(jì)算消耗僅約為原直接蒙特卡洛模擬的6%，大大提高了計(jì)算效率。

對于本算例，本文方法通過隨機(jī)物質(zhì)點(diǎn)方法分析可得到不同的破壞模式，根據(jù)滑動深度不同，可劃分為淺層、深層、中層破壞模式，另外還有漸進(jìn)破壞模式涉及到多重破壞模式同時或漸進(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)度分布圖。對于本算例，淺層破壞模式被定義為僅發(fā)生在上層土體。該破壞模式土體沿著上下層土的交界面滑動，順著坡面下滑并堆積在坡腳處。圖3 顯示該破壞模式發(fā)生位移約為30 m。極限平衡方法得到的臨界滑動面( 即對應(yīng)安全系數(shù)最小的滑動面) 與物質(zhì)點(diǎn)法中實(shí)際發(fā)生的滑動體形態(tài)較為一致。如圖3b 所示，上層土靠近坡面部分的不排水抗剪強(qiáng)度較低，低于左側(cè)土體，顯著低于下層土。這種情況容易產(chǎn)生淺層破壞。

2.2 深層破壞模式

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

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

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

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

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

圖8 繪制了漸進(jìn)破壞模式對應(yīng)的不排水抗剪強(qiáng)度分布。它的土體不排水抗剪強(qiáng)度分布與中層破壞模式的類似，但上層土的強(qiáng)度較低，上層土在第1個中層破壞模式發(fā)生后失去支撐，導(dǎo)致了第2 個淺層破壞模式的發(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 種典型破壞模式，并分析了對應(yīng)不排水抗剪強(qiáng)度與邊坡破壞模式之間的關(guān)系。得到如下結(jié)論:

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

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

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

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

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