彭曉鋼 李 嘉 李有志 江 建
深圳市天?。瘓F(tuán))股份有限公司 廣東 深圳 518034
準(zhǔn)確判斷邊坡的穩(wěn)定性對路基、邊坡、基坑等工程非常重要,但是常規(guī)有限元分析方法難以直接評價(jià)邊坡穩(wěn)定性,Zienkiewicz等[1]提出了強(qiáng)度折減法,該方法定義土體所能提供的最大抗剪強(qiáng)度和土體內(nèi)實(shí)際剪應(yīng)力之比為抗剪強(qiáng)度折減系數(shù),計(jì)算時(shí)先取較小的折減系數(shù)Ftrial,保證開始時(shí)邊坡處于穩(wěn)定狀態(tài),然后逐漸增大折減系數(shù)Ftrial,直至Ftrial增大到某一值時(shí)邊坡處于極限狀態(tài),即荷載引起的土體中實(shí)際剪應(yīng)力恰好等于按照實(shí)際強(qiáng)度指標(biāo)折減后的抗剪強(qiáng)度指標(biāo)。Duncan[2]將邊坡安全系數(shù)定義為使邊坡剛好達(dá)到臨界破壞狀態(tài)時(shí)對土的剪切強(qiáng)度進(jìn)行折減的程度,這一概念也就是傳統(tǒng)意義上的邊坡穩(wěn)定安全系數(shù)。
強(qiáng)度折減法提出伊始受到計(jì)算條件的限制普及難度較大,Griffiths等[3]利用Fortran語言編寫有限元程序首次實(shí)現(xiàn)了邊坡強(qiáng)度折減法分析。隨著計(jì)算機(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模塊等。
Abaqus是當(dāng)前應(yīng)用最為廣泛的非線性有限元軟件之一,該軟件沒有提供強(qiáng)度折減法,但可以通過場變量的設(shè)置簡便地實(shí)現(xiàn)[4-6]。
強(qiáng)度折減法的基本原理實(shí)質(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)力超過屈服面的單元形成貫通的曲面時(shí)土體將失去穩(wěn)定。
在Abaqus軟件中,可以通過設(shè)置強(qiáng)度參數(shù)隨場變量變化來實(shí)現(xiàn)強(qiáng)度折減,場變量沒有具體的物理意義,使計(jì)算分析更少受到人為假設(shè)的影響。在Abaqus軟件中,場變量Ffield隨時(shí)間增量步線性變化,材料參數(shù)和場變量之間有如下關(guān)系:
其中,c、φ是土體所能提供的抗剪強(qiáng)度指標(biāo),c'、φ'是某時(shí)間增量步下強(qiáng)度折減后土體的抗剪強(qiáng)度指標(biāo)。
對比式(1)和式(2),土體到達(dá)極限狀態(tài)時(shí)的抗剪強(qiáng)度折減系數(shù)Ftrial和場變量Ffield有如下關(guān)系:
在Abaqus軟件計(jì)算過程中,場變量隨著增量步時(shí)間t而變化,通過設(shè)置材料參數(shù)隨場變量變化,實(shí)現(xiàn)材料強(qiáng)度折減過程。在Abaqus軟件中可實(shí)現(xiàn)增量步時(shí)間t的自動變化,進(jìn)而實(shí)現(xiàn)強(qiáng)度的自動折減,并根據(jù)設(shè)定的破壞判據(jù)找到所對應(yīng)的增量步時(shí)間,然后可得安全系數(shù),并可得到相應(yīng)的滑動面。
正確判斷失穩(wěn)狀態(tài)影響到強(qiáng)度折減有限元法得到的折減系數(shù),進(jìn)而影響安全系數(shù)的計(jì)算結(jié)果,強(qiáng)度折減有限元法的失穩(wěn)判據(jù)主要有3類[7-8],即:數(shù)值迭代不收斂判據(jù)、特征部位位移突變判據(jù)、廣義塑性應(yīng)變或等效塑性應(yīng)變貫通判據(jù)。
為便于討論比較,以文獻(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)格劃分造成的計(jì)算結(jié)果誤差,采用和文獻(xiàn)[3]一致的網(wǎng)格劃分,采用縮減積分平面應(yīng)變單元CPE4R。約束模型左右方向邊界的水平方向位移,底端邊界的水平方向和豎直方向位移。有限元模型如圖1所示。
圖1 土坡強(qiáng)度折減有限元模型
對于實(shí)際計(jì)算邊坡安全系數(shù)的工程問題,強(qiáng)度折減系數(shù)Ftrial≥1,即根據(jù)土體實(shí)測強(qiáng)度指標(biāo)進(jìn)行強(qiáng)度折減,計(jì)算得到實(shí)際邊坡的安全系數(shù)。但對于有限元計(jì)算,不考慮實(shí)際意義,強(qiáng)度折減系數(shù)可以小于1,即折減系數(shù)起到擴(kuò)大強(qiáng)度指標(biāo)的作用,折減后強(qiáng)度指標(biāo)反而更趨于安全。
計(jì)算分2步進(jìn)行:第1步對邊坡施加豎直向下的重力,為使重力加載步順利完成,重力加載后多數(shù)區(qū)域處于彈性狀態(tài),并使強(qiáng)度折減計(jì)算平穩(wěn)過渡,采用Ftrial=0.5時(shí)的參數(shù),即黏聚力c=10.0 kPa,內(nèi)摩擦角φ=49.0°;第2步對強(qiáng)度參數(shù)進(jìn)行折減,設(shè)置最終狀態(tài)為Ftrial=2.0,此時(shí)對應(yīng)的土體強(qiáng)度參數(shù)為黏聚力c=2.5 kPa,內(nèi)摩擦角φ=16.7°。
經(jīng)過計(jì)算,可得不同泊松比v時(shí)的塑性區(qū)分布,如圖2所示。
圖2 不同泊松比v時(shí)的塑性區(qū)分布
由圖2可知,泊松比取值對塑性區(qū)的分布有較大影響,泊松比取較小值時(shí)比取較大值時(shí)的塑性區(qū)范圍大,泊松比取較大值,直至計(jì)算因不收斂而終止時(shí)的塑性區(qū)范圍依然很小,但當(dāng)泊松比取更小值時(shí)甚至?xí)a(chǎn)生塑性區(qū)貫通的現(xiàn)象,因此采用廣義塑性應(yīng)變或等效塑性應(yīng)變貫通判據(jù)不能準(zhǔn)確判斷土坡失穩(wěn)。
有限元最后的總位移等值線包含了第1步中施加重力荷載引起的位移,無法判斷滑動面的位置,通過增量步之間的位移可反映邊坡失穩(wěn)趨勢,計(jì)算采用最后2個(gè)增量步之間的位移進(jìn)行判斷,計(jì)算得到最后一步的位移增量等值線云圖如圖3所示,根據(jù)圖3可判斷滑動面的位置,和圖4所示的文獻(xiàn)[3]中滑動面形狀對比也極為接近。
分別選取坡頂點(diǎn)和坡腳點(diǎn)作為特征點(diǎn),將計(jì)算得到的強(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)降低,直至計(jì)算因不收斂而終止,坡腳和坡頂特征點(diǎn)豎向位移均呈現(xiàn)較小泊松比時(shí)位移較大的規(guī)律,符合彈性理論。其中坡腳特征點(diǎn)豎向位移沒有明顯突變,較大泊松比時(shí)豎向位移還有略微減小的趨勢;坡頂特征點(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)驗(yàn)判定位移突變點(diǎn),沒有明確的指標(biāo)作為判別依據(jù),有一定的經(jīng)驗(yàn)誤差。
將數(shù)值計(jì)算不收斂導(dǎo)致計(jì)算終止時(shí)的強(qiáng)度折減系數(shù)Ftrial匯總?cè)绫?所示。
表2 不收斂終止時(shí)的強(qiáng)度折減系數(shù)Ftrial
由表2可知,隨泊松比的增大,雖然塑性區(qū)的分布發(fā)生了相應(yīng)的變化,但計(jì)算因不收斂而終止時(shí)的強(qiáng)度折減系數(shù)Ftrial的差異很小,而且該方法無需根據(jù)經(jīng)驗(yàn)判定位移的突變點(diǎn)。
有限元中引起計(jì)算不收斂的因素有很多,有些學(xué)者質(zhì)疑以數(shù)值計(jì)算不收斂作為邊坡失穩(wěn)破壞依據(jù)是否存在其他人為導(dǎo)致的結(jié)果偏差,可通過對相同模型進(jìn)行不同的設(shè)置并進(jìn)行對比來說明這一問題。
該問題討論前首先要進(jìn)行假設(shè),即計(jì)算采用的有限元模型是正確的,如果模型建得不對或采用的有限元程序本身有缺陷,由此而引起的有限元數(shù)值計(jì)算不收斂本身就是毫無意義的。
上述模型計(jì)算采用的計(jì)算步長為1,初始增量步為0.1,最小增量步為1×10-5,其中最小增量步是判斷收斂與否的關(guān)鍵,將最小增量步進(jìn)一步減小,分別設(shè)置為1×10-6和1×10-7,將計(jì)算結(jié)果進(jìn)行對比如表3所示。
表3 收斂標(biāo)準(zhǔn)對計(jì)算結(jié)果的影響
Abaqus軟件中默認(rèn)的最小增量步為1×10-5,在將最小增量步進(jìn)一步調(diào)小至1×10-6和1×10-7后,計(jì)算因不收斂而終止時(shí)的強(qiáng)度折減系數(shù)變化不大,說明默認(rèn)的增量步已可以保證足夠的計(jì)算精度,計(jì)算結(jié)果也表明有限元收斂條件的設(shè)置對結(jié)果影響不大,因此可以把有限元計(jì)算是否收斂作為邊坡失穩(wěn)的判據(jù)。
利用Abaqus軟件進(jìn)行有限元強(qiáng)度折減計(jì)算的程序可靠,經(jīng)過對數(shù)值迭代不收斂判據(jù)、特征部位位移突變判據(jù)、廣義塑性應(yīng)變或等效塑性應(yīng)變貫通3類判據(jù)的計(jì)算分析,得到如下結(jié)論:
1)泊松比取值對塑性區(qū)的分布有較大影響,泊松比取較小值時(shí)比取較大值時(shí)的塑性區(qū)范圍大,泊松比取較大值,直至計(jì)算因不收斂而終止時(shí)的塑性區(qū)范圍依然很小,但當(dāng)泊松比取更小值時(shí)甚至?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)驗(yàn)判定位移突變點(diǎn),沒有明確的指標(biāo)作為判別依據(jù),有一定的經(jīng)驗(yàn)誤差。
3)數(shù)值計(jì)算是否收斂與邊坡是否失穩(wěn)存在對應(yīng)關(guān)系,計(jì)算因不收斂而終止時(shí)的強(qiáng)度折減系數(shù)變化不大,說明默認(rèn)的增量步已可以保證足夠的計(jì)算精度,計(jì)算結(jié)果也表明有限元收斂條件的設(shè)置對結(jié)果影響不大,建議利用有限元計(jì)算收斂作為失穩(wěn)判定的依據(jù)。