基于廣義Hoek-Brown準則的多參數強度折減法在邊坡穩定性分析中的應用

劉昱君， 劉明揚， 杜文杰

(1.中建國際投資(湖北)有限公司， 武漢 430071； 2.中國建筑國際集團有限公司， 香港 999077； 3.中國科學院 武漢巖土力學研究所,巖土力學與工程國家重點實驗室， 武漢 430071； 4.中國科學院大學， 北京 100049)

一直以來，道路交通、水利水電等工程在國家發展戰略中占有重要地位。中國中西部多山，地質條件復雜多變，工程施工過程中需要對山體進行開挖，這就對原山體造成了一定擾動，巖體穩定性因此降低。很容易在外力作用下(二次開挖、極端惡劣天氣等)發生山體滑坡、泥石流等自然災害。保障這些工程的安全重點在于保障道路、鐵路沿線開挖山體以及水利水電工程附近山體的穩定，這些都可以歸結為邊坡穩定性問題[1-3]。

近年來，隨著巖土體的本構模型和有限元理論的發展，有限元強度折減法已被廣泛應用于邊坡穩定性分析和其他相關工程問題上，成為邊坡穩定性分析應用最廣的方法之一。強度折減法的發展得益于子技術的發展：失效準則已不再局限于單一的計算不收斂準則，目前已經發展到塑性區貫通、塑性應變能突變、最大水平和豎直位移突變等多種失穩判據可以選用。同時，三維有限元模型強度折減法得到了很大的發展，這對于解決巖土工程中很多復雜邊坡問題、提高計算精度是一個全新的思路。不過，有限元強度折減法在非連續介質，比如支擋結構與巖土介質共同作用方面，計算精度的問題依然有待提高。

1999年，Griffith首次采用有限元強度折減法進行邊坡穩定性分析，求解得到的安全系數與傳統極限平衡分析法得到的安全系數在數值上十分接近，表明用有限元強度折減法進行邊坡穩定性分析是可以有效指導實踐的[4-5]。當然，隨著有限元的發展和廣泛的分析應用，現有強度折減法的問題也暴露了出來，僅僅對強度參數黏聚力c、內摩擦角φ進行折減的方案忽略了邊坡強度參數和變形參數之間存在的對應關系；另外對抗剪強度參數黏聚力c、內摩擦角φ進行折減模擬得到的結果往往與實際邊坡破壞相差較大，如果以塑性區貫通為失效準則，模型底部出現一部分不合理的塑性區，甚至有時會先于破壞面貫通，另外實際邊坡頂部的拉應力區也無法在模型中體現出來，這與實際邊坡破壞是有差別的。于是，在抗剪強度折減法的基礎上有學者相繼提出了對其他參數，主要是變形參數進行折減的方案，這些學者的研究成果進一步推動了強度折減法的完善和發展，更好地與工程實際相結合并正確指導工程實際[6-10]。

鄭宏等[11]指出如果只考慮對強度參數的折減，塑性區將集中在邊坡深部，模型底部在約束的影響下發生低應力破壞而出現塑性區，這與實際工程情況不符；袁維等[12]證明了彈性模量和泊松比對安全系數和塑性區分布有一定影響，對模型中下部的剪切塑性區影響較大，抗拉強度對邊坡的安全系數和塑性區分布的影響較大，因此有必要對變形參數進行折減；施建勇等[13]指出，只對抗剪強度參數進行折減并沒有從土體的本構關系上全面反映強度折減法的內涵，對此提出了在抗剪強度參數折減的基礎上，同時調整初始切線模量以滿足土體的變形性質。

現針對有限元強度折減法研究成果及其存在的不足，基于Hoek-Brown準則，考慮巖體彈性模量的沿深度遞增關系，提出一種針對彈性模量進行折減的方案；并結合邊坡算例進行穩定性計算驗證本文方案的合理性。

1 基于Hoek-Brown準則的折減方法探討

1.1 廣義Hoek-Brown準則

20世紀80年代，Hoek和Brown基于巖體性狀理論，結合巖體工程分類理論，統計并分析了大量現場試驗資料和三軸試驗資料，總結出一個巖體破壞的半經驗公式，Hoek-Brown準則[14-15]其極限主應力表達式為

(1)

(2)

(3)

(4)

式中：σ1為巖體破壞時的最大主應力；σ3為作用的最小主應力；σc為完整巖石單軸抗壓強度；m為反映巖石的軟硬程度的系數；s為反映巖體的破碎程度的系數；m與s均為由巖體工程分類和三軸試驗數據總結得到的與巖體破損程度有關的經驗參數；mb表示與巖體特征相關的經驗參數；mt表示巖石量綱為1的經驗參數；a表示與巖體特征相關的經驗參數；D是擾動系數；GSI為地質強度指標。

1.2 多參數協調折減方案的理論推導

假定巖土材料充滿自重應力場作用下的半無限空間，各方向主應力可以由下式得出

(5)

(6)

σc=ε0E

(7)

式中：ν為泊松比；h為深度；γ為巖體重度。文獻[16]提出了完整巖石抗壓強度與彈性模量的關系，由式(7)表達，其中ε0為應變參數，將式(5)～式(7)代入式(1)后得

(8)

文獻[17]中利用國內外關于地應力的實測資料，結合實驗室三軸試驗數據，通過數值方法得到彈性模量沿深度的變化規律：10 km埋深以上的土體彈性模量是沿深度線性增加的，因此假定埋深與巖土體彈性模量滿足

E=E0+ah

(9)

將式(9)變形為h=(E-E0)/a代入式(8)，得

(10)

化簡式(10)得

(11)

求解二元二次方程并化簡方程根，舍去負根，得

(12)

化簡得到Ei的表達式為

(13)

式(13)中除了νi、Ei外均是系數(參數取值參考文獻[16-18])，因此Ei可以認為是νi的函數。基于以上推導，建立了彈性模量與泊松比之間的協調折減關系。

1.3 算法

1.2節中所提出的彈性模量和泊松比關系式計算安全系數的建議算法如下：

1)確定巖體的參數c、φ、m、s、E0、a。

2)根據已有參數求得β=sinφ/(1-2ν)。

3)取折減系數Ki，黏聚力、內摩擦角按ci=c/Ki、φi=arctan(tanφ/Ki)進行折減。

4)對泊松比按νi=(1-sinφ/β)/2進行折減。

5)對彈性模量按式(13)進行折減。

6)得到的ci、φi、νi、Ei代入模型計算，若邊坡達到失穩破壞，輸出Ki為安全系數；否則調整Ki重復第3)～6)步直至邊坡發生失穩破壞。

2 非線性強度折減法算例分析

2.1 邊坡算例介紹

邊坡算例如圖1所示，模擬過程中使用彈性和塑性兩種巖土材料，采用非關聯流動法則和Mohr-Coulomb屈服準則，邊坡尺寸如圖1所示。對該邊坡算例用本文所提方法進行穩定性分析，為避免深部塑性區率先貫通使計算終止，模擬時以塑性區貫通和計算不收斂作為復合失穩判據。該邊坡材料參數見表1，式(13)中各參數取值參考文獻[18]。

圖1 邊坡模型

表1 邊坡模型材料參數

按照彈性力學規定，對于邊坡這樣截面尺寸相比長度很小的實體，可以簡化為平面應變問題。假定邊坡所受外力不隨縱深方向變化，位移和應變都發生在截面內。

對邊坡模型沿深度分層，以深度800 m處作為分界線，模型上部巖體為理想彈塑性材料，下部為彈性材料(圖1)，彈性模量隨深度逐層增大，給出如下理由：

1)模型分層考慮實際邊坡深部土體的壓實，設置合理的彈性模量變化來模擬彈塑性變形，其塑性區的發展、應力的分布符合實際情況。

2)在式(13)的推導過程中涉及彈性模量隨深度變化的關系式，因此對邊坡模型沿深度分層以保證推導過程中各彈性模量保持一致。

3)考慮模型分層，塑性區下部的彈性材料單元可以產生一定的垂直變形和水平變形，減小了由于邊界效應在邊坡下部出現的塑性區。

2.2 折減方案對比設計

為了證明多參數協調折減的合理性，根據兩種折減辦法設計實驗方案進行比對：①抗剪強度折減法，僅考慮折減黏聚力c和內摩擦角φ；②基于本文提出的泊松比和彈性模量折減關系，同時考慮黏聚力c、內摩擦角φ、泊松比ν和彈性模量E的多參數協調折減方法。針對以上折減方法設計如下方案：

方案1：根據抗剪強度折減法，首先取Ki=1.2，對黏聚力c、內摩擦角φ兩個參數按傳統抗剪強度折減法進行折減。分別取度折減系數分別為Ki=1.4、Ki=1.6、Ki=1.8、Ki=2.0、Ki=2.2、Ki=2.4、Ki=2.6、Ki=2.8、Ki=3.0的邊坡材料本構模型。折減系數與對應調整的輸入參數見表2。

表2 傳統抗剪強度折減法折減后的參數

方案2：根據本文提出的多參數折減方案，對照流程圖，首先取Ki=1.2，對黏聚力c、內摩擦角φ、泊松比ν和彈性模量E折減。同上，定義強度折減系數分別為Ki=1.4、Ki=1.6、Ki=1.8、Ki=2.0、Ki=2.2、Ki=2.4、Ki=2.6、Ki=2.8、Ki=3.0的邊坡材料本構模型，調整后的參數c、φ、ν和E見表3。

表3 本文折減方案折減后的參數

2.3 計算結果與分析

2.3.1 基于折減方案1的計算結果與分析

輸入表2中各折減系數對應的參數進行求解，讀入不同折減參數下的結果文件，在后處理中顯示邊坡塑性應變云圖，查看不同折減系數對應的模型塑性區分布，如圖2所示。

圖2 傳統抗剪強度折減方案折減系數對應的累積塑性應變分布

折減系數Ki=3.0時，迭代計算不收斂，模擬終止。為找到準確的折減系數，對折減系數在2.8～3.0區間內細分，細化折減參數折減后的黏聚力和摩擦角，見表4。

表4 抗剪強度折減法細化折減后的參數

各折減系數對應邊坡模型的塑性區分布如圖3所示，可以看到，當折減系數取Ki=2.92時，折減的參數使模型塑性區完全貫通，此時邊坡將發生失穩破壞，塑性區上部的土體將沿貫通的塑性區滑動，因此方案1安全系數取為2.92。

由圖2和圖3可以看出，除了滑動面的塑性區貫通，模型深部塑性應變明顯，甚至先于滑動面貫通，這明顯不符合實際工程情況。鄭宏等在文獻[11]中指出，坡底的塑性區是由于低應力水平作用下，由于邊界約束效應引起的破壞，與彈塑性材料的選取有關。

2.3.2 基于折減方案2的計算結果與分析

按照方案2的折減方法，對模型輸入表3參數進行穩定性分析。在迭代計算收斂的情況下，在后處理中查看不同折減系數下模型塑性區分布，如圖4所示。

圖4 本文方案折減系數對應的累積塑性應變分布

折減系數Ki=3.0時，模型計算不收斂，為找到使邊坡發生失穩時的折減系數，在2.8～3.0區間內對折減系數細分，對應細化折減后的參數見表5。

表5 本文方案細化折減后的參數

細分后各折減系數對應模型的塑性區分布如圖5所示，可以看到，當折減系數Ki=2.86時，輸入的計算參數可以使邊坡模型塑性區完全貫通，沿貫通的塑性區邊坡發生失穩破壞，因此方案2所得折減系數取為2.86。

圖5 本文方案細化折減系數對應的累積塑性應變分布

2.3.3 小結

對比兩種方案的計算結果，方案2求解得到的安全系數偏小；方案1模型破壞過程中底部塑性區明顯，方案2有效消除了底部不合理塑性區；對于實際邊坡破壞時坡頂出現的拉應力區，兩種方案都沒有體現出來，對比見表6。

表6 兩種方案的對比

3 結論

基于當前有限元強度折減法研究成果及其存在的不足，基于Hoek-Brown準則，考慮巖體彈性模量的沿深度遞增關系，提出了一種針對彈性模量進行折減的強度折減法；結合邊坡算例進行穩定性計算，通過ANSYS建立有限元模型，以塑性區貫通為邊坡失穩判據，模擬巖石邊坡的漸進劣化失穩破壞，最后將結果與抗剪強度折減法所得到的安全系數進行比較，得出以下結論：

1)在強度折減計算過程中，黏聚力c、內摩擦角φ、泊松比ν和彈性模量E的折減不是任意的，而是應該按照一定的協調關系進行同步折減。

2)本文提出的多參數協調折減方案得到的安全系數略小于傳統抗剪強度折減法計算的安全系數，穩定性評價較為嚴格，消除了傳統方法折減后模型底部出現的不合理塑性區，表明本文提出的多參數協調折減方案相比傳統抗剪強度折減法貼近工程實際，可以有效指導邊坡防治。