曾 逸

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

計算分兩步進(jìn)行,在第一步中,采用土體的真實(shí)強(qiáng)度參數(shù)模擬整個加載過程;在第二步中,控制溫度場隨時步的增加而變化,從而帶動強(qiáng)度參數(shù)不斷折減。在每一時步,考察該時刻對應(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)力在單元中重分配。若計算收斂,則時步自動增加一個 △t,強(qiáng)度參數(shù)折減到t+△t時刻對應(yīng)的值,繼續(xù)計算該時刻的應(yīng)力和變形;若計算不收斂,則減少增量步再次迭代。但這個求解過程,都是由 ABAQUS自動完成的,無需重新編制程序或人為干預(yù)計算過程。

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

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

在邊坡穩(wěn)定分析中,采用不同的屈服準(zhǔn)則,往往也會使計算結(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)流動法則進(jìn)行有限元計算。

另一方面,不同研究者對極限狀態(tài)的確定是不相同的。有的認(rèn)為計算不收斂便是極限狀態(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)致不同的計算結(jié)果,需要加以考察。

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

不收斂時刻 t并不能直接代入式(4)來確定安全系數(shù)。根據(jù)分析問題和采用控制條件的不同,不收斂時 δmax/H值的范圍是不同的。δmax/H可能大于 0.1,可以判斷為邊坡已經(jīng)破壞,只是由于采用單元類型的數(shù)值穩(wěn)定性和程序極強(qiáng)的非線性計算功能而繼續(xù)計算。也可能 δmax/H小于0.05,說明邊坡變形還沒有突變,可能是模型或控制條件不當(dāng)導(dǎo)致程序提前不收斂。只有當(dāng)不收斂時對應(yīng)的 δmax/H在 0.05～0.1之間時,才可把不收斂時刻 t代入式(4)來確定安全系數(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給出無量綱位移 δmax/H與時步 t的關(guān)系。取 δmax/H急劇增大,且其值位于 0.05～0.1之間時,對應(yīng)的強(qiáng)度折減系數(shù)作為邊坡整體的穩(wěn)定系數(shù),可以取 t=0.3115,對應(yīng) δmax/H=0.05,安全系數(shù) Fs=1.390。

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

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

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

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

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

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

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

圖 5給出無量綱位移 δmax/H與時步 t的關(guān)系,可以得出當(dāng) t=0.304時達(dá)到極限狀態(tài),對應(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之間 ,可以直接采用計算不收斂作為判斷標(biāo)準(zhǔn)。

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

5 結(jié)語

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

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

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

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

〔3〕鄭穎人,趙尚毅.巖土工程極限分析有限元法及其應(yīng)用[J].土木工程學(xué)報,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é)報,2003,23(3):1～8.

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

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