溫控參數有限元法在邊坡穩定分析中的應用

曾 逸

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

1 引言

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

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

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

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

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

2 溫控參數有限元法的基本原理

2.1 抗剪強度的折減

為了使邊坡達到極限狀態,需要對內聚力、摩擦角和膨脹角分別按式(1)進行折減。采用有效應力指標,

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

對土坡強度折減的程度,即土的實際強度與極限狀態時所采用強度的比值,具有強度儲備的物理意義,和傳統極限平衡方法是一致的。

2.2 溫控參數有限元法的基本原理

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

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

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

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

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

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

3 屈服準則和極限狀態判據

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

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

為了與采用圓弧滑裂面進行極限平衡分析的經典算例進行比較,本文采用膨脹角為 0的非關聯流動法則進行有限元計算。

另一方面,不同研究者對極限狀態的確定是不相同的。有的認為計算不收斂便是極限狀態[3～5],或者以一定迭代次數作為極限狀態的判定準則[6]。此外還有人把塑性區的連通看成是達到了極限狀態[7],或者以邊坡體內某一幅值的廣義剪應變從邊坡底部下方向坡頂上方貫通作為極限狀態的標志[8]。但是,不同的收斂準則或不同的迭代方法可能導致不同的計算結果,需要加以考察。

采用有限元程序進行計算,有時邊坡變形甚至大于坡高,計算仍收斂;有時邊坡幾乎沒有變形,計算已經不收斂。如果按照通常的計算不收斂作為判斷標準,肯定會有很大的誤差。本文針對通用有限元軟件 ABAQUS,采用四節點的四邊形單元進行計算。考察無量綱位移 δmax/H與時步 t的關系,其中 δmax為邊坡體的最大節點位移,H為邊坡高度。通過以下兩個算例發現,當 t大于某一數值后,δmax/H會快速增長,并導致程序不收斂而中止計算。

不收斂時刻 t并不能直接代入式(4)來確定安全系數。根據分析問題和采用控制條件的不同,不收斂時 δmax/H值的范圍是不同的。δmax/H可能大于 0.1,可以判斷為邊坡已經破壞,只是由于采用單元類型的數值穩定性和程序極強的非線性計算功能而繼續計算。也可能 δmax/H小于0.05,說明邊坡變形還沒有突變,可能是模型或控制條件不當導致程序提前不收斂。只有當不收斂時對應的 δmax/H在 0.05～0.1之間時,才可把不收斂時刻 t代入式(4)來確定安全系數 Fs,因為

4 算例分析

4.1 均質土坡的安全系數

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

圖 1 均質土坡的有限元網格

圖 2給出無量綱位移 δmax/H與時步 t的關系。取 δmax/H急劇增大,且其值位于 0.05～0.1之間時,對應的強度折減系數作為邊坡整體的穩定系數,可以取 t=0.3115,對應 δmax/H=0.05,安全系數 Fs=1.390。

圖 2 無量綱位移與時步t的關系

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

圖 3 對應于Fs=1.390的網格變形圖(a)節點位移矢量;(b)變形網格

4.2 非均質土坡的安全系數

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

圖 4 非均質邊坡的幾何模型

表 1 各土層的強度參數指標

圖 5給出無量綱位移 δmax/H與時步 t的關系,可以得出當 t=0.304時達到極限狀態,對應的 δmax/H=0.08,安全系數 Fs=1.377。這與推薦的裁判答案 Fs=1.39誤差 1%。此例收斂于 t=0.304,δmax/H=0.08,位于 0.05 ～ 0.1之間 ,可以直接采用計算不收斂作為判斷標準。

圖 5 無量綱位移與時步t的關系

5 結語

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

本文通過兩個算例,得出取 δmax/H在 0.05～0.1之間時,對應的 Fs作為整體安全系數這一土坡破壞標準,便于數值計算,避免采用不收斂作為破壞標準的缺點。

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

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

〔3〕鄭穎人,趙尚毅.巖土工程極限分析有限元法及其應用[J].土木工程學報,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〕欒茂田,武亞軍,年廷凱.強度折減有限元法中邊坡失穩的塑性區判據及其應用[J].防災減災工程學報,2003,23(3):1～8.

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

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