溫控參數(shù)有限元法在邊坡穩(wěn)定分析中的應(yīng)用

曾 逸

(雙流縣水務(wù)局,四川雙流,610200)

1 引言

目前邊坡工程中常用的穩(wěn)定分析方法是基于極限平衡原理的傳統(tǒng)方法。傳統(tǒng)穩(wěn)定分析方法積累了豐富的經(jīng)驗(yàn),成果可靠。但是由于傳統(tǒng)方法做了很多假設(shè)和簡(jiǎn)化,例如需要假設(shè)滑裂面的位置、不考慮條分土體內(nèi)部的應(yīng)力應(yīng)變關(guān)系等,不能得出滑坡體真實(shí)的應(yīng)力和變形。

近年來(lái),不少研究者嘗試采用有限元方法來(lái)進(jìn)行邊坡的穩(wěn)定性分析。有限元法能全面滿足靜力平衡條件、應(yīng)變相容條件和非線性本構(gòu)關(guān)系,可以不受邊坡形狀不規(guī)則、材料不均勻的限制,不需預(yù)先假定滑裂面,能模擬滑坡的可能運(yùn)動(dòng)方向,確定破壞區(qū)的位置和破壞發(fā)展情況。

由于許多邊坡通常并不受到表面分布載荷的作用,目前大多數(shù)學(xué)者是通過(guò)土體抗剪強(qiáng)度參數(shù)的折減而超載使邊坡達(dá)到極限狀態(tài)的。即假定實(shí)際的土體強(qiáng)度參數(shù)(內(nèi)聚力及內(nèi)摩擦角)折減某一倍數(shù) Fs后,作為一組新的材料參數(shù)代入進(jìn)行有限元計(jì)算,并判斷邊坡在完成全部加載過(guò)程之前是否達(dá)到極限狀態(tài)。如果經(jīng)過(guò)某個(gè) Fs值的折減能夠使邊坡達(dá)到極限狀態(tài),而經(jīng)過(guò)比該值稍小的Fs值的折減不能使邊坡達(dá)到極限狀態(tài),此時(shí)的 Fs即為安全系數(shù)。通過(guò)這種方法還可以同時(shí)得到臨界滑裂面的位置。

在通常的有限元分析程序中,材料參數(shù)是需要在輸入文件中給定的。為了按照上述方法確定折減系數(shù),需要修改輸入文件中的 c和f值進(jìn)行反復(fù)的試算。顯然,為了求得安全系數(shù),計(jì)算量大,即使采用二分法,一般也需要計(jì)算 10次左右[1]。

另一種作法是,首先按照實(shí)際的強(qiáng)度參數(shù)指標(biāo)進(jìn)行模擬計(jì)算,完成全部的加載過(guò)程后,再分步降低土體的強(qiáng)度參數(shù)值進(jìn)行迭代計(jì)算,直至邊坡進(jìn)入極限狀態(tài),由此得出相應(yīng)的安全系數(shù)[2]。這種作法應(yīng)該更加符合實(shí)際情況,而且整個(gè)分析過(guò)程是一次完成,不用修改輸入文件反復(fù)試算。但是,文獻(xiàn)[2]所采用的方法需要調(diào)整屈服面并修正應(yīng)力,從而計(jì)算體系的不平衡力。這些算法都需要重新編制有限元程序才能實(shí)現(xiàn)。本文利用大多數(shù)通用有限元軟件(如 ABAQUS)中都具有的可使材料參數(shù)隨溫度變化的功能,用溫度場(chǎng)控制強(qiáng)度參數(shù)變化來(lái)實(shí)現(xiàn)參數(shù)的連續(xù)折減。這樣作不僅能最大限度地利用 ABAQUS本身所具有的計(jì)算精度高、非線性求解能力強(qiáng)、可靠性強(qiáng)、運(yùn)算速度快、使用方便等優(yōu)點(diǎn),而且整個(gè)分析過(guò)程也是一次完成,大大簡(jiǎn)化了計(jì)算工作量。

2 溫控參數(shù)有限元法的基本原理

2.1 抗剪強(qiáng)度的折減

為了使邊坡達(dá)到極限狀態(tài),需要對(duì)內(nèi)聚力、摩擦角和膨脹角分別按式(1)進(jìn)行折減。采用有效應(yīng)力指標(biāo),

其中,c′f、φ′f和 ψ′f分別為折減 Fs倍后的有效內(nèi)聚力、摩擦角和膨脹角。按照前述方法確定的臨界折減系數(shù) Fs,即為邊坡穩(wěn)定的安全系數(shù)。

對(duì)土坡強(qiáng)度折減的程度,即土的實(shí)際強(qiáng)度與極限狀態(tài)時(shí)所采用強(qiáng)度的比值,具有強(qiáng)度儲(chǔ)備的物理意義,和傳統(tǒng)極限平衡方法是一致的。

2.2 溫控參數(shù)有限元法的基本原理

傳統(tǒng)的離散試驗(yàn)算法是強(qiáng)度參數(shù)不隨時(shí)步變化,每一組試驗(yàn)參數(shù)都要進(jìn)行一個(gè)完整的時(shí)步計(jì)算。顯然,要找出一個(gè)使邊坡剛好達(dá)到極限狀態(tài)的 Fs是比較繁瑣的。而在 ABAQUS程序中,可以利用其現(xiàn)成的材料參數(shù)隨溫度場(chǎng)變量的變化而變化的功能,定義材料強(qiáng)度指標(biāo)隨著溫度場(chǎng)的變化而變化。此溫度場(chǎng)只是一個(gè)變量場(chǎng),不代表真實(shí)溫度,只起到帶動(dòng)材料參數(shù)變化的作用。如果給定其熱膨脹系數(shù)為 0,那么溫度變化不會(huì)給結(jié)構(gòu)帶來(lái)應(yīng)力和變形上的變化。

對(duì)于安全系數(shù)在 1.0～10.0之間的一般情況,當(dāng)溫度場(chǎng)變量 θ由 0增加到 100時(shí),定義 c′f和 tanφ′f隨著 θ線性折減到 0.1c′f和 0.1tanφ′f,即 1/Fs由 1.0折減到 0.1。其數(shù)學(xué)關(guān)系式為:

靜力分析中的時(shí)步不代表真實(shí)時(shí)間,只代表“載荷”的變化過(guò)程。當(dāng) t由 0增加到 1.0時(shí),定義溫度場(chǎng) θ隨時(shí)步 t由 0線性增加到 100。即

把式(3)代入式(2),從而得出安全系數(shù) Fs和時(shí)步 t之間關(guān)系:

這樣,通過(guò)遞推,實(shí)現(xiàn)了材料強(qiáng)度參數(shù)與時(shí)步t的一一對(duì)應(yīng),并隨著 t的增加而線性折減。對(duì)于安全系數(shù)小于 1.0的情況也可以采用同樣的方法,只要增大 1/Fs的值即可。

計(jì)算分兩步進(jìn)行,在第一步中,采用土體的真實(shí)強(qiáng)度參數(shù)模擬整個(gè)加載過(guò)程;在第二步中,控制溫度場(chǎng)隨時(shí)步的增加而變化,從而帶動(dòng)強(qiáng)度參數(shù)不斷折減。在每一時(shí)步,考察該時(shí)刻對(duì)應(yīng)的強(qiáng)度參數(shù)下每一點(diǎn)的應(yīng)力狀態(tài),與破壞準(zhǔn)則相比較。如果某一點(diǎn)上的應(yīng)力位于屈服面內(nèi),則認(rèn)為該點(diǎn)處于彈性響應(yīng),如果應(yīng)力位于屈服面上,則按照彈塑性理論,將屈服應(yīng)力在單元中重分配。若計(jì)算收斂,則時(shí)步自動(dòng)增加一個(gè) △t,強(qiáng)度參數(shù)折減到t+△t時(shí)刻對(duì)應(yīng)的值,繼續(xù)計(jì)算該時(shí)刻的應(yīng)力和變形;若計(jì)算不收斂,則減少增量步再次迭代。但這個(gè)求解過(guò)程,都是由 ABAQUS自動(dòng)完成的,無(wú)需重新編制程序或人為干預(yù)計(jì)算過(guò)程。

3 屈服準(zhǔn)則和極限狀態(tài)判據(jù)

這表示邊坡變形突變,達(dá)到極限狀態(tài)。如果δmax/H大于 0.1,則做出 δmax/H～ t關(guān)系曲線,找出 δmax/H位于 0.05～0.1之間且其值突變時(shí)刻對(duì)應(yīng)的 Fs作為安全系數(shù)。此時(shí),邊坡變形急劇增加,Fs差別甚小。這一土坡破壞標(biāo)準(zhǔn)便于數(shù)值計(jì)算,避免采用不收斂作為破壞標(biāo)準(zhǔn)可能導(dǎo)致結(jié)果偏大或者偏小的缺點(diǎn)。

在邊坡穩(wěn)定分析中,采用不同的屈服準(zhǔn)則,往往也會(huì)使計(jì)算結(jié)果有較大的差異。采用外接圓的Duncker-Prager模型,得到的安全系數(shù)比傳統(tǒng)結(jié)算結(jié)果大 25%左右[3]。本文采用巖土工程中廣泛采用的 Mohr-Coulomb準(zhǔn)則,大量的工程實(shí)踐證明了該準(zhǔn)則的精度和可靠性。

為了與采用圓弧滑裂面進(jìn)行極限平衡分析的經(jīng)典算例進(jìn)行比較,本文采用膨脹角為 0的非關(guān)聯(lián)流動(dòng)法則進(jìn)行有限元計(jì)算。

另一方面,不同研究者對(duì)極限狀態(tài)的確定是不相同的。有的認(rèn)為計(jì)算不收斂便是極限狀態(tài)[3～5],或者以一定迭代次數(shù)作為極限狀態(tài)的判定準(zhǔn)則[6]。此外還有人把塑性區(qū)的連通看成是達(dá)到了極限狀態(tài)[7],或者以邊坡體內(nèi)某一幅值的廣義剪應(yīng)變從邊坡底部下方向坡頂上方貫通作為極限狀態(tài)的標(biāo)志[8]。但是,不同的收斂準(zhǔn)則或不同的迭代方法可能導(dǎo)致不同的計(jì)算結(jié)果,需要加以考察。

采用有限元程序進(jìn)行計(jì)算,有時(shí)邊坡變形甚至大于坡高,計(jì)算仍收斂;有時(shí)邊坡幾乎沒(méi)有變形,計(jì)算已經(jīng)不收斂。如果按照通常的計(jì)算不收斂作為判斷標(biāo)準(zhǔn),肯定會(huì)有很大的誤差。本文針對(duì)通用有限元軟件 ABAQUS,采用四節(jié)點(diǎn)的四邊形單元進(jìn)行計(jì)算。考察無(wú)量綱位移 δmax/H與時(shí)步 t的關(guān)系,其中 δmax為邊坡體的最大節(jié)點(diǎn)位移,H為邊坡高度。通過(guò)以下兩個(gè)算例發(fā)現(xiàn),當(dāng) t大于某一數(shù)值后,δmax/H會(huì)快速增長(zhǎng),并導(dǎo)致程序不收斂而中止計(jì)算。

不收斂時(shí)刻 t并不能直接代入式(4)來(lái)確定安全系數(shù)。根據(jù)分析問(wèn)題和采用控制條件的不同,不收斂時(shí) δmax/H值的范圍是不同的。δmax/H可能大于 0.1,可以判斷為邊坡已經(jīng)破壞,只是由于采用單元類型的數(shù)值穩(wěn)定性和程序極強(qiáng)的非線性計(jì)算功能而繼續(xù)計(jì)算。也可能 δmax/H小于0.05,說(shuō)明邊坡變形還沒(méi)有突變,可能是模型或控制條件不當(dāng)導(dǎo)致程序提前不收斂。只有當(dāng)不收斂時(shí)對(duì)應(yīng)的 δmax/H在 0.05～0.1之間時(shí),才可把不收斂時(shí)刻 t代入式(4)來(lái)確定安全系數(shù) Fs,因?yàn)?/p>

4 算例分析

4.1 均質(zhì)土坡的安全系數(shù)

土坡幾何模型及四邊形單元網(wǎng)格劃分如圖 1所示,坡高 20m,坡度 1∶2,右邊界距坡趾和左邊界距坡頂均為 40m,底邊界距坡底 10m。材料參數(shù):E=100MPa,ν=0.3,γ=20kN/m3,c=20kPa,φ=20°。基底采用剛性邊界,兩側(cè)為水平約束,上部邊界為自由邊界。

圖 1 均質(zhì)土坡的有限元網(wǎng)格

圖 2給出無(wú)量綱位移 δmax/H與時(shí)步 t的關(guān)系。取 δmax/H急劇增大,且其值位于 0.05～0.1之間時(shí),對(duì)應(yīng)的強(qiáng)度折減系數(shù)作為邊坡整體的穩(wěn)定系數(shù),可以取 t=0.3115,對(duì)應(yīng) δmax/H=0.05,安全系數(shù) Fs=1.390。

圖 2 無(wú)量綱位移與時(shí)步t的關(guān)系

圖 3(a)和圖 3(b)對(duì)應(yīng)于極限狀態(tài) Fs=1.390時(shí)節(jié)點(diǎn)位移矢量和變形網(wǎng)格。由位移場(chǎng)和變形場(chǎng)的分布情況判斷,這種判定準(zhǔn)則是合理的。對(duì)于給定的問(wèn)題,所得到的有限元數(shù)值解與 Bishop和 Morgenstern的極限平衡解 Fs=1.380、Grif-fiths和 Lame的有限元數(shù)值解 Fs=1.4都很接近[5]。這種確定安全系數(shù)的方法避免了目前采用數(shù)值計(jì)算不收斂作為判斷標(biāo)準(zhǔn)的人為不確定性。事實(shí)上,由于 ABAQUS的非線性計(jì)算能力很強(qiáng),本算例在 t=0.335時(shí)候才開始不收斂。如果按照計(jì)算不收斂作為判斷標(biāo)準(zhǔn),安全系數(shù) Fs=1.432,對(duì)應(yīng)的 δmax/H=2.0,邊坡已經(jīng)破壞。

圖 3 對(duì)應(yīng)于Fs=1.390的網(wǎng)格變形圖(a)節(jié)點(diǎn)位移矢量;(b)變形網(wǎng)格

4.2 非均質(zhì)土坡的安全系數(shù)

本例來(lái)源于澳大利亞計(jì)算機(jī)應(yīng)用協(xié)會(huì)(ACADS)的一個(gè)考題[9],幾何模型和材料分區(qū)如圖 4所示,材料特性如表 1所示。

圖 4 非均質(zhì)邊坡的幾何模型

表 1 各土層的強(qiáng)度參數(shù)指標(biāo)

圖 5給出無(wú)量綱位移 δmax/H與時(shí)步 t的關(guān)系,可以得出當(dāng) t=0.304時(shí)達(dá)到極限狀態(tài),對(duì)應(yīng)的 δmax/H=0.08,安全系數(shù) Fs=1.377。這與推薦的裁判答案 Fs=1.39誤差 1%。此例收斂于 t=0.304,δmax/H=0.08,位于 0.05 ～ 0.1之間 ,可以直接采用計(jì)算不收斂作為判斷標(biāo)準(zhǔn)。

圖 5 無(wú)量綱位移與時(shí)步t的關(guān)系

5 結(jié)語(yǔ)

在 ABAQUS程序中,采用溫控參數(shù)法來(lái)計(jì)算邊坡的穩(wěn)定安全系數(shù),不僅能最大限度地利用ABAQUS本身所具有的計(jì)算精度高、非線性求解能力強(qiáng)、可靠性強(qiáng)、運(yùn)算速度快、使用方便等優(yōu)點(diǎn),而且不用另外編制程序,整個(gè)分析過(guò)程也是一次完成的,大大簡(jiǎn)化了計(jì)算工作量。

本文通過(guò)兩個(gè)算例,得出取 δmax/H在 0.05～0.1之間時(shí),對(duì)應(yīng)的 Fs作為整體安全系數(shù)這一土坡破壞標(biāo)準(zhǔn),便于數(shù)值計(jì)算,避免采用不收斂作為破壞標(biāo)準(zhǔn)的缺點(diǎn)。

〔1〕程 曄,趙明華,曹文貴.基樁下溶洞頂板穩(wěn)定性評(píng)價(jià)的強(qiáng)度折減有限元法[J].巖土工程學(xué)報(bào),2005,27(1):38～41.

〔2〕宋二祥,高 翔,邱 玥.基坑土釘支護(hù)安全系數(shù)的強(qiáng)度參數(shù)折減有限元方法[J].巖土工程學(xué)報(bào),2005,27(3):258～263.

〔3〕鄭穎人,趙尚毅.巖土工程極限分析有限元法及其應(yīng)用[J].土木工程學(xué)報(bào),2005,38(1):91～98.

〔4〕Dawson E M,Roth W H,Drescher A.Slope stability analysis by strength reduction[J].Geotechnique,1999,49(6):835～840.

〔5〕Griffiths D V,Lame P A.Slope stability analysis by finite elements[J].Geotechnique,1999,49(3):387～403.

〔6〕Ugai K.A method of calculation of total factor of safety of slopes by elasto-plastic FEM[J].Soils and Foundations,1989.29(2):190～195.

〔7〕欒茂田,武亞軍,年廷凱.強(qiáng)度折減有限元法中邊坡失穩(wěn)的塑性區(qū)判據(jù)及其應(yīng)用[J].防災(zāi)減災(zāi)工程學(xué)報(bào),2003,23(3):1～8.

〔8〕連鎮(zhèn)營(yíng),韓國(guó)城,孔憲京.強(qiáng)度折減法研究開挖邊坡的穩(wěn)定性[J].巖土工程學(xué)報(bào),2001,23(4):406～411.

〔9〕陳祖煜.土質(zhì)邊坡穩(wěn)定分析:原理·方法·程序[M].北京:中國(guó)水利水電出版社,2003.