彭曉鋼 李 嘉 李有志 江 建

深圳市天?。瘓F(tuán)）股份有限公司 廣東 深圳 518034

準(zhǔn)確判斷邊坡的穩(wěn)定性對路基、邊坡、基坑等工程非常重要，但是常規(guī)有限元分析方法難以直接評價邊坡穩(wěn)定性，Zienkiewicz等[1]提出了強(qiáng)度折減法，該方法定義土體所能提供的最大抗剪強(qiáng)度和土體內(nèi)實際剪應(yīng)力之比為抗剪強(qiáng)度折減系數(shù)，計算時先取較小的折減系數(shù)Ftrial，保證開始時邊坡處于穩(wěn)定狀態(tài)，然后逐漸增大折減系數(shù)Ftrial，直至Ftrial增大到某一值時邊坡處于極限狀態(tài)，即荷載引起的土體中實際剪應(yīng)力恰好等于按照實際強(qiáng)度指標(biāo)折減后的抗剪強(qiáng)度指標(biāo)。Duncan[2]將邊坡安全系數(shù)定義為使邊坡剛好達(dá)到臨界破壞狀態(tài)時對土的剪切強(qiáng)度進(jìn)行折減的程度，這一概念也就是傳統(tǒng)意義上的邊坡穩(wěn)定安全系數(shù)。

強(qiáng)度折減法提出伊始受到計算條件的限制普及難度較大，Griffiths等[3]利用Fortran語言編寫有限元程序首次實現(xiàn)了邊坡強(qiáng)度折減法分析。隨著計算機(jī)硬件和數(shù)值軟件的發(fā)展，強(qiáng)度折減法愈加受到重視，很多巖土工程商業(yè)軟件中已經(jīng)可以較為方便地應(yīng)用強(qiáng)度折減法，如美國Itasca公司開發(fā)的FLAC、荷蘭Technical University of Delft開發(fā)的Plaxis、瑞士Zace Services公司開發(fā)的Z-SOIL.PC、韓國POSCO公司開發(fā)的Midas/GTS模塊等。

1 強(qiáng)度折減法在Abaqus軟件中的實現(xiàn)

Abaqus是當(dāng)前應(yīng)用最為廣泛的非線性有限元軟件之一，該軟件沒有提供強(qiáng)度折減法，但可以通過場變量的設(shè)置簡便地實現(xiàn)[4-6]。

強(qiáng)度折減法的基本原理實質(zhì)就是強(qiáng)度參數(shù)c（黏聚力）、φ（內(nèi)摩擦角）的逐漸降低導(dǎo)致土體強(qiáng)度不足以承擔(dān)當(dāng)前應(yīng)力，在有限元中即反映為強(qiáng)度參數(shù)c、φ的逐漸降低導(dǎo)致某單元應(yīng)力超過了屈服面產(chǎn)生塑性變形，超出屈服面的應(yīng)力將通過塑性變形傳遞到周圍土體單元，當(dāng)應(yīng)力超過屈服面的單元形成貫通的曲面時土體將失去穩(wěn)定。

在Abaqus軟件中，可以通過設(shè)置強(qiáng)度參數(shù)隨場變量變化來實現(xiàn)強(qiáng)度折減，場變量沒有具體的物理意義，使計算分析更少受到人為假設(shè)的影響。在Abaqus軟件中，場變量Ffield隨時間增量步線性變化，材料參數(shù)和場變量之間有如下關(guān)系：

其中，c、φ是土體所能提供的抗剪強(qiáng)度指標(biāo)，c'、φ'是某時間增量步下強(qiáng)度折減后土體的抗剪強(qiáng)度指標(biāo)。

對比式（1）和式（2），土體到達(dá)極限狀態(tài)時的抗剪強(qiáng)度折減系數(shù)Ftrial和場變量Ffield有如下關(guān)系：

在Abaqus軟件計算過程中，場變量隨著增量步時間t而變化，通過設(shè)置材料參數(shù)隨場變量變化，實現(xiàn)材料強(qiáng)度折減過程。在Abaqus軟件中可實現(xiàn)增量步時間t的自動變化，進(jìn)而實現(xiàn)強(qiáng)度的自動折減，并根據(jù)設(shè)定的破壞判據(jù)找到所對應(yīng)的增量步時間，然后可得安全系數(shù)，并可得到相應(yīng)的滑動面。

2 Abaqus軟件中強(qiáng)度折減法的失穩(wěn)判定

正確判斷失穩(wěn)狀態(tài)影響到強(qiáng)度折減有限元法得到的折減系數(shù)，進(jìn)而影響安全系數(shù)的計算結(jié)果，強(qiáng)度折減有限元法的失穩(wěn)判據(jù)主要有3類[7-8]，即：數(shù)值迭代不收斂判據(jù)、特征部位位移突變判據(jù)、廣義塑性應(yīng)變或等效塑性應(yīng)變貫通判據(jù)。

2.1 分析模型

為便于討論比較，以文獻(xiàn)[3]中的案例為對象，用Abaqus軟件進(jìn)行強(qiáng)度折減有限元分析，對上述失穩(wěn)判據(jù)進(jìn)行對比。

Example1為均質(zhì)邊坡，坡高H＝10 m，坡角β＝30°，土體重度γ＝20 kN/m3，黏聚力c＝5 kPa，內(nèi)摩擦角φ＝30°。

建立邊坡有限元模型，幾何尺寸和文獻(xiàn)[3]中一致，坡頂和坡腳到模型邊界均取為12 m，坡腳以下土體厚度取10 m，假定彈性模量E＝100 MPa，為消除因網(wǎng)格劃分造成的計算結(jié)果誤差，采用和文獻(xiàn)[3]一致的網(wǎng)格劃分，采用縮減積分平面應(yīng)變單元CPE4R。約束模型左右方向邊界的水平方向位移，底端邊界的水平方向和豎直方向位移。有限元模型如圖1所示。

圖1 土坡強(qiáng)度折減有限元模型

對于實際計算邊坡安全系數(shù)的工程問題，強(qiáng)度折減系數(shù)Ftrial≥1，即根據(jù)土體實測強(qiáng)度指標(biāo)進(jìn)行強(qiáng)度折減，計算得到實際邊坡的安全系數(shù)。但對于有限元計算，不考慮實際意義，強(qiáng)度折減系數(shù)可以小于1，即折減系數(shù)起到擴(kuò)大強(qiáng)度指標(biāo)的作用，折減后強(qiáng)度指標(biāo)反而更趨于安全。

計算分2步進(jìn)行：第1步對邊坡施加豎直向下的重力，為使重力加載步順利完成，重力加載后多數(shù)區(qū)域處于彈性狀態(tài)，并使強(qiáng)度折減計算平穩(wěn)過渡，采用Ftrial＝0.5時的參數(shù)，即黏聚力c＝10.0 kPa，內(nèi)摩擦角φ＝49.0°；第2步對強(qiáng)度參數(shù)進(jìn)行折減，設(shè)置最終狀態(tài)為Ftrial＝2.0，此時對應(yīng)的土體強(qiáng)度參數(shù)為黏聚力c＝2.5 kPa，內(nèi)摩擦角φ＝16.7°。

經(jīng)過計算，可得不同泊松比v時的塑性區(qū)分布，如圖2所示。

圖2 不同泊松比v時的塑性區(qū)分布

由圖2可知，泊松比取值對塑性區(qū)的分布有較大影響，泊松比取較小值時比取較大值時的塑性區(qū)范圍大，泊松比取較大值，直至計算因不收斂而終止時的塑性區(qū)范圍依然很小，但當(dāng)泊松比取更小值時甚至?xí)a(chǎn)生塑性區(qū)貫通的現(xiàn)象，因此采用廣義塑性應(yīng)變或等效塑性應(yīng)變貫通判據(jù)不能準(zhǔn)確判斷土坡失穩(wěn)。

有限元最后的總位移等值線包含了第1步中施加重力荷載引起的位移，無法判斷滑動面的位置，通過增量步之間的位移可反映邊坡失穩(wěn)趨勢，計算采用最后2個增量步之間的位移進(jìn)行判斷，計算得到最后一步的位移增量等值線云圖如圖3所示，根據(jù)圖3可判斷滑動面的位置，和圖4所示的文獻(xiàn)[3]中滑動面形狀對比也極為接近。

2.2 特征部位位移突變對邊坡失穩(wěn)判定的影響

分別選取坡頂點(diǎn)和坡腳點(diǎn)作為特征點(diǎn)，將計算得到的強(qiáng)度折減系數(shù)Ftrial同特征點(diǎn)豎向位移U之間的關(guān)系繪制成曲線，如圖5所示。

圖3 最后一步的位移增量等值線云圖

圖4 文獻(xiàn)[3]中的滑動面形狀

圖5 坡頂和坡腳特征點(diǎn)Ftrial-U關(guān)系曲線

由圖5可知，隨著強(qiáng)度折減系數(shù)增大，土體強(qiáng)度指標(biāo)降低，直至計算因不收斂而終止，坡腳和坡頂特征點(diǎn)豎向位移均呈現(xiàn)較小泊松比時位移較大的規(guī)律，符合彈性理論。其中坡腳特征點(diǎn)豎向位移沒有明顯突變，較大泊松比時豎向位移還有略微減小的趨勢；坡頂特征點(diǎn)位移突變較明顯，然而位移突變點(diǎn)的強(qiáng)度折減系數(shù)差異較大，選取曲線水平段結(jié)束點(diǎn)作為突變點(diǎn)，強(qiáng)度折減系數(shù)Ftrial如表1所示。

表1 坡頂特征點(diǎn)位移突變點(diǎn)強(qiáng)度折減系數(shù)Ftrial

由表1可知，塑性區(qū)的分布會影響到坡頂特征點(diǎn)位移突變點(diǎn)的強(qiáng)度折減系數(shù)Ftrial，隨泊松比增大，塑性區(qū)范圍減小，由此判斷的強(qiáng)度折減系數(shù)Ftrial呈增大趨勢，且該方法需根據(jù)經(jīng)驗判定位移突變點(diǎn)，沒有明確的指標(biāo)作為判別依據(jù)，有一定的經(jīng)驗誤差。

2.3 數(shù)值迭代不收斂對邊坡失穩(wěn)判定的影響

將數(shù)值計算不收斂導(dǎo)致計算終止時的強(qiáng)度折減系數(shù)Ftrial匯總?cè)绫?所示。

表2 不收斂終止時的強(qiáng)度折減系數(shù)Ftrial

由表2可知，隨泊松比的增大，雖然塑性區(qū)的分布發(fā)生了相應(yīng)的變化，但計算因不收斂而終止時的強(qiáng)度折減系數(shù)Ftrial的差異很小，而且該方法無需根據(jù)經(jīng)驗判定位移的突變點(diǎn)。

有限元中引起計算不收斂的因素有很多，有些學(xué)者質(zhì)疑以數(shù)值計算不收斂作為邊坡失穩(wěn)破壞依據(jù)是否存在其他人為導(dǎo)致的結(jié)果偏差，可通過對相同模型進(jìn)行不同的設(shè)置并進(jìn)行對比來說明這一問題。

該問題討論前首先要進(jìn)行假設(shè)，即計算采用的有限元模型是正確的，如果模型建得不對或采用的有限元程序本身有缺陷，由此而引起的有限元數(shù)值計算不收斂本身就是毫無意義的。

上述模型計算采用的計算步長為1，初始增量步為0.1，最小增量步為1×10－5，其中最小增量步是判斷收斂與否的關(guān)鍵，將最小增量步進(jìn)一步減小，分別設(shè)置為1×10－6和1×10－7，將計算結(jié)果進(jìn)行對比如表3所示。

表3 收斂標(biāo)準(zhǔn)對計算結(jié)果的影響

Abaqus軟件中默認(rèn)的最小增量步為1×10－5，在將最小增量步進(jìn)一步調(diào)小至1×10－6和1×10－7后，計算因不收斂而終止時的強(qiáng)度折減系數(shù)變化不大，說明默認(rèn)的增量步已可以保證足夠的計算精度，計算結(jié)果也表明有限元收斂條件的設(shè)置對結(jié)果影響不大，因此可以把有限元計算是否收斂作為邊坡失穩(wěn)的判據(jù)。

3 結(jié)語

利用Abaqus軟件進(jìn)行有限元強(qiáng)度折減計算的程序可靠，經(jīng)過對數(shù)值迭代不收斂判據(jù)、特征部位位移突變判據(jù)、廣義塑性應(yīng)變或等效塑性應(yīng)變貫通3類判據(jù)的計算分析，得到如下結(jié)論：

1）泊松比取值對塑性區(qū)的分布有較大影響，泊松比取較小值時比取較大值時的塑性區(qū)范圍大，泊松比取較大值，直至計算因不收斂而終止時的塑性區(qū)范圍依然很小，但當(dāng)泊松比取更小值時甚至?xí)a(chǎn)生塑性區(qū)貫通的現(xiàn)象，因此采用廣義塑性應(yīng)變或等效塑性應(yīng)變貫通判據(jù)不能準(zhǔn)確判斷土坡失穩(wěn)。

2）塑性區(qū)的分布會影響到坡頂特征點(diǎn)位移突變點(diǎn)的強(qiáng)度折減系數(shù)Ftrial，隨泊松比增大，塑性區(qū)范圍減小，由此判斷的強(qiáng)度折減系數(shù)Ftrial呈增大趨勢，且該方法需根據(jù)經(jīng)驗判定位移突變點(diǎn)，沒有明確的指標(biāo)作為判別依據(jù)，有一定的經(jīng)驗誤差。

3）數(shù)值計算是否收斂與邊坡是否失穩(wěn)存在對應(yīng)關(guān)系，計算因不收斂而終止時的強(qiáng)度折減系數(shù)變化不大，說明默認(rèn)的增量步已可以保證足夠的計算精度，計算結(jié)果也表明有限元收斂條件的設(shè)置對結(jié)果影響不大，建議利用有限元計算收斂作為失穩(wěn)判定的依據(jù)。