基于有限元的邊坡穩(wěn)定可靠度改進(jìn)一次二階矩法研究

汪安南，林潮寧，李同春，2，齊慧君，邵天成

(1.河海大學(xué)水利水電學(xué)院，江蘇 南京 210024；2.水安全與水科學(xué)協(xié)同創(chuàng)新中心，江蘇 南京 210024)

1 研究背景

邊坡失穩(wěn)引起的滑坡災(zāi)害已成為全球范圍內(nèi)主要地質(zhì)災(zāi)害之一[1]。由大型邊坡失穩(wěn)和滑坡引起的工程事故歷史上屢見不鮮[2]，且其造成的后果也是不可估計的，因此邊坡的抗滑穩(wěn)定安全度研究成為國內(nèi)外工程界十分關(guān)心的研究課題[3]。目前，在邊坡穩(wěn)定性分析方面的研究大致可以分為2類：一類是由極限平衡法和數(shù)值計算方法所組成的確定性分析方法[4]；另一類是考慮參數(shù)隨機(jī)性、變異性的可靠性分析方法[5]。相比于確定性分析，可靠性分析方法考慮了參數(shù)不確定性對邊坡穩(wěn)定性的影響，更接近實際情況。

常用的邊坡可靠性分析方法主要有蒙特卡洛法、響應(yīng)面法及一次二階矩法等。李育超等[6]將蒙特卡洛法與有限元法相結(jié)合，利用隨機(jī)走步法不斷更新試算滑動面，確定了邊坡臨界滑動面及其對應(yīng)的最小安全系數(shù)；何淑軍等[7]采用蒙特卡洛法對夏呀河四級滑坡進(jìn)行可靠性分析，定量地表達(dá)了該滑坡的安全程度；周罕等[8]采用響應(yīng)面法和強(qiáng)度折減法相結(jié)合的分析方法，把安全系數(shù)和破壞概率作為綜合評價邊坡穩(wěn)定性的指標(biāo)，從不同角度對邊坡可靠性進(jìn)行分析；陳昌富等[9]基于響應(yīng)面法建立了一種高效的邊坡可靠指標(biāo)和失效概率近似計算方法，大大降低了構(gòu)造響應(yīng)面函數(shù)的計算工作量；朱彥鵬等[10]基于改進(jìn)一次二階矩法建立了錨定板擋土墻的系統(tǒng)可靠度分析模型，分析了參數(shù)在不同變異系數(shù)下對系統(tǒng)可靠度指標(biāo)的影響和土性參數(shù)敏感性。通常，當(dāng)功能函數(shù)為隱式時，僅蒙特卡洛法和響應(yīng)面法可用，但蒙特卡洛法的模擬次數(shù)多，計算效率低；響應(yīng)面法通過不同的響應(yīng)面函數(shù)重構(gòu)隱式功能函數(shù)，但采用傳統(tǒng)響應(yīng)面形式處理復(fù)雜問題計算精度有限，采用機(jī)器學(xué)習(xí)算法構(gòu)建響應(yīng)面函數(shù)涉及超參數(shù)選取等問題，計算復(fù)雜度高[11]。

考慮上述方法在處理邊坡隱式極限狀態(tài)函數(shù)時的問題，本文提出一種基于有限元的邊坡穩(wěn)定可靠度改進(jìn)一次二階矩法，定義極限狀態(tài)函數(shù)為滑體極限狀態(tài)下抗滑力與下滑力的差值，考慮黏聚力及內(nèi)摩擦角作為隨機(jī)變量并且服從獨(dú)立正態(tài)分布，求出極限狀態(tài)方程中隨機(jī)變量的偏導(dǎo)數(shù)，基于改進(jìn)一次二階矩法編寫了求解邊坡失穩(wěn)概率和可靠指標(biāo)的計算程序，通過算例驗證了所提方法的適用性，并對某庫岸邊坡進(jìn)行了可靠性分析。

2 顯式極限狀態(tài)函數(shù)的建立及改進(jìn)一次二階矩法原理

2.1 邊坡顯式極限狀態(tài)函數(shù)

有限元迭代解法[3]是一種以剛體極限平衡基本思想為基礎(chǔ)的直接迭代解法，是一種在滑動面已知情況下邊坡穩(wěn)定分析的有效方法，不僅可計算出邊坡沿已知滑動面的安全系數(shù)，且能同時給出邊坡臨界失穩(wěn)時的位移場與應(yīng)力場。

用有限元迭代解法分析時，邊坡穩(wěn)定的非線性特性主要反映在滑動面上，滑動面上的強(qiáng)度準(zhǔn)則一般采用摩爾-庫倫準(zhǔn)則?；瑒用娴膯卧问街饕捎梦墨I(xiàn)[12]中提出的常規(guī)的矩形(平面)或立方體(空間)的統(tǒng)一模式。當(dāng)邊坡處于臨界穩(wěn)定狀態(tài)時，滑動面上的滑動力和阻滑力滿足

(1)

式中，m為危巖體內(nèi)滑動面單元個數(shù)；Wi為第i個滑動面單元上的法向力；Pi為第i個滑動面單元上的切向力；Ai為第i個滑動面單元上扣除受拉區(qū)后的面積；ci為第i個滑動面單元上的黏聚力；φi為第i個滑動面單元上的內(nèi)摩擦角；K為假定安全系數(shù)；K*為真實的抗滑穩(wěn)定安全系數(shù)。

由于力的矢量特征，安全系數(shù)為在邊坡潛在滑動面上總抗滑力矢在整體下滑趨勢方向上投影的代數(shù)和與總下滑力矢在此方向上投影代數(shù)和的比值。當(dāng)滑體處于極限平衡狀態(tài)時，安全系數(shù)即為極限抗滑力與下滑力的代數(shù)和比值。將邊坡滑動面上的滑動力達(dá)到抗滑力作為邊坡極限狀態(tài)的標(biāo)志，故將抗力函數(shù)R和作用函數(shù)S相等的情況作為極限狀態(tài)，即

(2)

式中，S(g)為作用函數(shù)；R(g)為抗力函數(shù)；g為函數(shù)變量；Z0為功能函數(shù)。對隨機(jī)變量求偏導(dǎo)數(shù)得

(3)

式中，Z為功能函數(shù)；fi為第i個滑動面單元上的內(nèi)摩擦系數(shù)；X為隨機(jī)變量的向量表示，X1i、X2i分別為X的第1個、第2個隨機(jī)變量。

2.2 改進(jìn)一次二階矩法的基本原理及計算步驟

(4)

(5)

根據(jù)隨機(jī)變量線性組合的性質(zhì)及可靠指標(biāo)的定義，可得

(6)

對式(3)、(4)進(jìn)行變換，得設(shè)計驗算點(diǎn)在原始X空間中的坐標(biāo)為

xi=μXi+βσXicosθXi

(7)

3 算例驗證

3.1 計算模型

含有1個確定性滑面的邊坡算例是邊坡穩(wěn)定分析中非常經(jīng)典的例子[13]，其二維模型結(jié)構(gòu)清晰，受力狀態(tài)簡單，可以得到其力學(xué)上的解析解。單一接觸面巖質(zhì)滑坡體幾何尺寸見圖1。巖土體抗剪強(qiáng)度參數(shù)視為隨機(jī)變量，服從正態(tài)分布，具體數(shù)值見表1。

圖1 單一接觸面巖質(zhì)滑坡體截面(單位：m)

表1 巖土體物理力學(xué)參數(shù)

相關(guān)研究顯示，變量之間的相關(guān)性對可靠指標(biāo)的影響較小，可以認(rèn)為是相互獨(dú)立的[14-15]，暫不考慮黏聚力與內(nèi)摩擦角的相關(guān)性。目前研究中[16-17]，對于巖土體黏聚力及內(nèi)摩擦角變異系數(shù)的取值均在0.1～0.3范圍內(nèi)，本算例巖質(zhì)材料的黏聚力及內(nèi)摩擦角的變異系數(shù)取值為0.15。邊坡顯式極限狀態(tài)函數(shù)中滑動面上各單元的切向力及法向力采用非線性有限元方法求取，綜合考慮時間和計算精度，建立有限元計算模型，見圖2。共560個單元，609個節(jié)點(diǎn)。底部基巖的兩側(cè)和底部邊界均為固定約束，上部滑動塊體為自由邊界。

圖2 有限元計算模型

3.2 計算過程與分析

荷載考慮邊坡自身重力，采用非線性有限元計算滑動面上單元切向力、法向力，建立邊坡滑動面的顯式極限狀態(tài)函數(shù)；依據(jù)驗算點(diǎn)法的邊坡可靠性分析流程，進(jìn)行可靠指標(biāo)的迭代計算。按照本文的邊坡可靠指標(biāo)計算流程，共進(jìn)行2輪可靠指標(biāo)的迭代計算即滿足精度要求，所以該邊坡的可靠指標(biāo)為1.493 9，其破壞概率為6.76%。

通常，蒙特卡洛法計算結(jié)果可作為基準(zhǔn)解[18]。為驗證本文所提方法計算結(jié)果的準(zhǔn)確性，結(jié)合蒙特卡洛法與有限元強(qiáng)度折減法，經(jīng)過100 000次的蒙特卡洛模擬，得到該邊坡的失效概率為6.18%，計算耗時約為6.0×105s(166.67 h)，總抗滑力投影與下滑力投影差值分布見圖3。從圖3可知，本文方法計算結(jié)果與蒙特卡洛法計算結(jié)果相比，誤差在3%以內(nèi)，說明本文方法計算出來的失效概率與蒙特卡洛法的計算結(jié)果具有可比性，計算結(jié)果合理可行。

圖3 邊坡總抗滑力投影與下滑力投影差值分布

同時，采用其他目前常用的可靠性分析方法對該邊坡算例進(jìn)行計算，結(jié)果見表2?；诒疚慕⒌娘@式極限狀態(tài)函數(shù)，采用均值一次二階矩法對該邊坡算例進(jìn)行可靠性分析，失效概率為12.04%，此結(jié)果與蒙特卡洛法計算結(jié)果相比，誤差為94.82%。采用基于二次多項式響應(yīng)面函數(shù)的改進(jìn)一次二階矩法對該邊坡算例進(jìn)行可靠性分析，失效概率為5.21%，此結(jié)果與蒙特卡洛法計算結(jié)果相比，誤差為15.70%。說明本文所提方法相比于均值一次二階矩法、響應(yīng)面法具有更高的計算精度。此外，相比于響應(yīng)面法的較高計算復(fù)雜度，本文方法計算過程簡單易于操作，能考慮滑坡體、巖基力學(xué)特性的不均勻性及滑動面上的應(yīng)力分布影響。

表2 可靠性計算結(jié)果

4 工程應(yīng)用

4.1 工程概況

某庫岸邊坡位于右岸壩前，距大壩900～1 700 m，坡高650～700 m，為花崗巖巖質(zhì)邊坡，陡傾岸內(nèi)及岸外構(gòu)造發(fā)育。根據(jù)地質(zhì)勘探結(jié)果，該庫岸巖體分為散體結(jié)構(gòu)、碎裂結(jié)構(gòu)、塊裂結(jié)構(gòu)和鑲嵌一次塊狀結(jié)構(gòu)4種。其中，散體結(jié)構(gòu)水平厚度一般小于30 m，碎裂結(jié)構(gòu)一般在40～60 m，塊裂結(jié)構(gòu)一般在50～80 m。此外，岸坡平臺還分布幾條拉裂陷落帶，其中規(guī)模較大的有LF54和LF55等。3號梁庫岸巖體結(jié)構(gòu)見圖4。

圖4 庫岸巖體結(jié)構(gòu)

4.2 計算模型

本文選取3號梁地質(zhì)剖面，底部全長1.582 km，高差0.788 km。在此基礎(chǔ)上，建立了1組整體穩(wěn)定分析模型和3組局部穩(wěn)定性分析模型，見圖5、6。為考慮該變形體是否存在整體失穩(wěn)的可能，對庫岸沿變形和非變形界限的抗滑穩(wěn)定性進(jìn)行分析，選取了3號梁地質(zhì)剖面，以地質(zhì)勘查推測的變形范圍界限和碎裂巖體下限為假定滑動面。對于局部穩(wěn)定分析，指定滑動面以頂部平臺出露的LF54為后緣拉開面，以巖體結(jié)構(gòu)面為假定滑面，以不同埋深深度分別進(jìn)行分析，局部計算模型見圖6。

圖5 整體穩(wěn)定計算模型

圖6 局部穩(wěn)定計算模型

由于本文給出的極限狀態(tài)函數(shù)含有邊坡滑動面各單元的切向力及法向力，采用非線性有限元法計算，根據(jù)工程地質(zhì)剖面圖，建立相應(yīng)的數(shù)值模型，見圖7。綜合考慮計算時間和精度，有限元計算模型共1 616個單元，1 736個節(jié)點(diǎn)；邊界條件為下部固定，左右兩側(cè)水平約束，上部為自由邊界。

圖7 穩(wěn)定性分析模型

覆蓋層除外，計算模型巖體參數(shù)取值見表3?？紤]滑動面抗剪強(qiáng)度參數(shù)指標(biāo)c、f為隨機(jī)變量進(jìn)行可靠度分析，結(jié)合同類巖土工程經(jīng)驗，可以認(rèn)為各滑動面強(qiáng)度參數(shù)服從正態(tài)分布，統(tǒng)計特性指標(biāo)見表4。由于邊坡上建筑物荷載條件未明確，故未考慮坡體加載的影響。

表3 巖體力學(xué)參數(shù)取值

表4 結(jié)構(gòu)面參數(shù)統(tǒng)計特征值

4.3 計算過程與分析

采用非線性有限元法計算滑動面單元的切向力和法向力，建立滑動面的顯式極限狀態(tài)函數(shù)；依據(jù)驗算點(diǎn)法的邊坡可靠度分析流程，對整體及局部共4種假定滑動面進(jìn)行可靠指標(biāo)的迭代計算，直到滿足精度為止；為對邊坡的安全狀態(tài)有更全面的認(rèn)識，采用巖體及各滑動面強(qiáng)度參數(shù)的平均值作為標(biāo)準(zhǔn)值，采用有限元強(qiáng)度折減法進(jìn)行安全系數(shù)計算，計算結(jié)果見表5。

從表5可知，整體及局部的安全系數(shù)均高于現(xiàn)行邊坡規(guī)范抗滑穩(wěn)定安全系數(shù)控制下限值1.15，然而局部的失穩(wěn)概率較高，需適當(dāng)?shù)牟扇」こ檀胧┮越档惋L(fēng)險。該計算結(jié)果與實際監(jiān)測中邊坡局部出現(xiàn)拉裂、位錯及坡面坍塌等變形破壞但整體保持穩(wěn)定的現(xiàn)象互相印證。安全系數(shù)與可靠性分析對邊坡局部狀態(tài)判斷出現(xiàn)矛盾的原因是：①表征邊坡穩(wěn)定性的安全系數(shù)沒有考慮材料參數(shù)的變異性，低估了邊坡局部的失穩(wěn)風(fēng)險；②散體結(jié)構(gòu)帶和碎裂結(jié)構(gòu)帶淺表部材料參數(shù)雖然均值較高，但是變異性大(滑動面巖土體的變異系數(shù)都達(dá)到0.3)，導(dǎo)致邊坡局部失穩(wěn)概率較大。

表5 各計算模型可靠指標(biāo)計算結(jié)果

因此，應(yīng)采取分級放坡[8]，加強(qiáng)對邊坡的支護(hù)處理。為保持邊坡的長久穩(wěn)定，應(yīng)植草對坡面進(jìn)行防沖刷處理，對外露的坡面采用格構(gòu)式等形式進(jìn)行護(hù)坡。在做好坡上截水、坡體排水、坡面防護(hù)等綜合治理的同時，應(yīng)加強(qiáng)對邊坡的監(jiān)測及數(shù)據(jù)分析，做好災(zāi)害的主動防范，提前預(yù)報。

5 結(jié) 語

針對蒙特卡洛法與響應(yīng)面法在處理邊坡隱式極限狀態(tài)函數(shù)時的問題，本文提出了一種基于有限元的邊坡穩(wěn)定可靠度改進(jìn)一次二階矩法，建立了邊坡顯式極限狀態(tài)函數(shù)，并應(yīng)用于某庫岸邊坡可靠性分析中，得出以下結(jié)論：

(1)與已有的可靠性分析不同，本文定義極限狀態(tài)函數(shù)為滑體極限狀態(tài)下抗滑力與下滑力的差值，考慮黏聚力及內(nèi)摩擦角作為隨機(jī)變量并且服從獨(dú)立正態(tài)分布，求出極限狀態(tài)方程中隨機(jī)變量的偏導(dǎo)數(shù)，在降低計算復(fù)雜度，提升計算效率的同時保證計算精度，計算精度上優(yōu)于響應(yīng)面法，計算效率上優(yōu)于蒙特卡洛法。

(2)采用本文提出的邊坡可靠度分析方法，多角度綜合分析了某庫岸邊坡的穩(wěn)定性。由于滑坡穩(wěn)定性分析不可能考慮全部的影響因素，因此除了巖土體的抗剪強(qiáng)度參數(shù)黏聚力、內(nèi)摩擦角外，考慮其他參數(shù)影響的邊坡可靠性分析需要進(jìn)一步研究。