地震作用下邊坡及滑坡穩(wěn)定性及破壞過程分析

李兆宇, 馬宗源, 魏睿真, 焦 凱

(1.中國水利水電第三工程局有限公司, 陜西 西安 710024; 2.西安理工大學(xué) 巖土工程研究所, 陜西 西安 710048)

地震是影響邊坡穩(wěn)定性的主要外部因素之一,而強震引發(fā)的山體滑坡是地震誘發(fā)的主要震害類型之一,對地震后期的災(zāi)后重建工作造成巨大影響。通過總結(jié)分析2008年汶川地震中形成的滑坡特征及分布情況,發(fā)現(xiàn)地震誘發(fā)滑坡的破壞力及危害性極大,有些甚至超過了地震本身所造成的損失[1]。

滑坡穩(wěn)定性分析方法主要由土力學(xué)中的邊坡穩(wěn)定分析方法改進(jìn)發(fā)展而來。上世紀(jì)初提出的圓弧滑動面的極限平衡方法,用來分析靜力情況下邊坡穩(wěn)定問題,后來基于土條間力平衡假定提出了條分方法[2-3]。我國學(xué)者基于二維情況下Spencer條分方法提出了一種新的三維極限平衡邊坡穩(wěn)定分析方法且具有一定實用性[4]。上世紀(jì)70年代研究人員提出利用強度折減的方式結(jié)合有限元方法分析邊坡穩(wěn)定,后來稱為強度折減法[5-6]。與傳統(tǒng)條分方法相比,有限元強度折減法不需要假定滑動面及土條間的作用力,可自動搜索邊坡的滑動面并計算安全系數(shù)[7]。

目前擬靜力方法是邊坡地震穩(wěn)定性問題的主要分析方法。該方法主要基于靜力條件下邊坡穩(wěn)定性分析法基本思路,即通過極限平衡方法或有限元強度折減法確定邊坡的地震安全系數(shù)[8-9]。該方法雖然簡化了計算,但是卻忽略了地震產(chǎn)生的慣性力影響以及邊坡的動力反應(yīng)過程,所以無法反映地震過程中邊坡的動力響應(yīng)及邊坡失穩(wěn)后的破壞過程。此外,擬靜力方法使用恒定的體力荷載考慮地震影響,使得邊坡安全系數(shù)的計算結(jié)果過小,往往偏于保守。Newmark基于剛塑性滑塊假定提出了地震條件下邊坡永久位移的計算公式并廣泛應(yīng)用于邊坡地震穩(wěn)定性問題的計算分析[10-11]。該方法能夠合理考慮地震效應(yīng)的影響,但仍然基于小變形及極限平衡假定進(jìn)行計算,無法考慮土體的動力特性及分析地震過程中邊坡的失穩(wěn)破壞過程。此外研究人員采用分型方法考慮巖土體力學(xué)特性、邊坡形態(tài)特征及地震烈度等因素來綜合評價邊坡的地震穩(wěn)定性[12]。

本文基于顯式有限元及動力學(xué)大變形方法提出一種適合邊坡地震穩(wěn)定性及其破壞過程的非線性分析方法[13-14],分析了地震情況下滑坡的穩(wěn)定性及滑坡失穩(wěn)后的破壞過程。

1 理論及方法

1.1 土動力學(xué)基本理論

采用土動力學(xué)模型分析地震情況下邊坡穩(wěn)定性問題。選取Hardin-Drnevich骨干加載曲線模型預(yù)測土的動模量和阻尼比隨動應(yīng)變的衰減關(guān)系。動三軸試驗數(shù)據(jù)說明Hardin-Drnevich模型更適于預(yù)測飽和黃土的動應(yīng)力應(yīng)變關(guān)系[10]。根據(jù)Hardin-Drnevich模型,土的動剪切模量和阻尼比可寫為動剪應(yīng)變的函數(shù)[15-17]:

G/Gmax=1/(1+γ/γref)

(1)

D/Dmax=(γ/γref)/(1+γ/γref)

(2)

式中:G為土體的動剪切模量;D為土體的阻尼比;Gmax為土體的初始剪切模量;Dmax為土體的最大阻尼比;γ為土體的動剪應(yīng)變幅值;γref為參考剪應(yīng)變(γref=τmax/Gmax,即G/Gmax=0.5時動剪切模量對應(yīng)的動剪應(yīng)變幅值γ)。γref越大,土的動剪切模量隨動應(yīng)變衰減越快,土的阻尼越大。

三維情況下一般將γ取為廣義剪應(yīng)變形式,即動力過程中土體的剪切模量和阻尼比與廣義剪應(yīng)變相關(guān)。但是廣義剪應(yīng)變是恒為正的,不能反映循環(huán)加載時剪應(yīng)變的正反變化,即不能準(zhǔn)確反映反向加載時的應(yīng)力路徑。本文采用將循環(huán)加載過程中卸載到反向加載階段的廣義剪應(yīng)變增量△γ乘以一個折減系數(shù)a,以抵消反向加載過程中廣義剪應(yīng)變不能為負(fù)的影響。根據(jù)黏彈性理論構(gòu)建阻尼和彈簧并聯(lián)型的黏彈性本構(gòu)模型:

(3)

使用Visual Fortran語言編制黃土動本構(gòu)關(guān)系計算程序,由VUMAT子程序接口導(dǎo)入ABAQUS有限元計算軟件。采用一個四節(jié)點減縮積分單元計算單元的剪應(yīng)力及剪應(yīng)變關(guān)系,并與黃土動三軸試驗測試數(shù)據(jù)進(jìn)行對比。對該單元的頂部約束水平自由度,單元底部水平方向上加載正弦加速度時程曲線,正弦時程函數(shù)見式(4),峰值加速度為0.04 g。

(4)

式中:A為加速度;t為振動時間。

圖1為不同類型土的動剪切模量、阻尼比與動應(yīng)變的關(guān)系,其中陜北黃土數(shù)據(jù)取自文獻(xiàn)[16],砂土和黏土數(shù)據(jù)取自文獻(xiàn)[18]和[19]。黃土動力學(xué)參數(shù)為: 最大剪切剛度Gmax=29.2MPa, 最大阻尼比Dmax=0.159,參考剪應(yīng)變γref=0.03?？梢钥闯龈鶕?jù)Hardin-Drnevich骨干加載曲線建立黃土動力本構(gòu)模型能夠很好模擬黃土加卸載過程中的滯回曲線形狀的應(yīng)力應(yīng)變關(guān)系。

圖1 土的動剪切模量、阻尼比與動剪應(yīng)變的關(guān)系Fig.1 Relationship between dynamic shear modulus, damping ratio and dynamic shear strain of soil

圖2顯示了輸入正弦幅值不同系數(shù)a取值情況下廣義剪應(yīng)變的時程變化曲線。如a=1.0情況下動剪切模量和阻尼比完全由廣義剪應(yīng)變確定,a=0情況下動剪切模量和阻尼比只與廣義剪應(yīng)變正向加載幅值有關(guān),與卸載及反向加載過程無關(guān)。圖3黃土的動應(yīng)力與動應(yīng)變關(guān)系可以看出a=1.0情況下在反向加載過程中動應(yīng)力應(yīng)變曲線切線斜率明顯偏大,而a=0情況下的動應(yīng)力應(yīng)變曲線呈橢圓形狀,a=0.5時動應(yīng)力應(yīng)變曲線與試驗結(jié)果相近。

圖2 廣義剪應(yīng)變時程曲線Fig.2 Generalized shear strain time history curve

圖3 黃土動應(yīng)力應(yīng)變關(guān)系理論預(yù)測與試驗數(shù)據(jù)的對比Fig.3 Comparison of prediction and experimental data on dynamic stress-strain relationship of loess

1.2 計算方法

使用強度折減法計算邊坡的安全系數(shù)(FOS),邊坡抗剪強度參數(shù)折減方式可寫為[7,20]:

(5)

式中:c和φ分別為邊坡土體原始的黏聚力及內(nèi)摩擦角;cf和φf分別為經(jīng)過折減之后的邊坡土體的黏聚力及內(nèi)摩擦角。邊坡位移突然增大(折減系數(shù)與邊坡最大位移關(guān)系曲線拐點)時刻的折減系數(shù)確定為邊坡安全系數(shù)。

顯式有限元方法和隱式有限元方法最大的區(qū)別在于求解非線性問題是否迭代,是否所有待求物理量在同一時刻獲得解答。隱式方法處理非線性問題時一般無法直接求解,需要采用迭代方法求解雅可比矩陣,對于高度非線性問題迭代求解可能會不收斂。顯式方法依靠時間積分求解控制方程,無需迭代直接進(jìn)行求解,求解高度非線性問題具有一定優(yōu)勢。顯式有限元計算過程存在波動引起求解誤差,需要控制時間步長保證求解穩(wěn)定性。根據(jù)節(jié)點力的平衡方程,加速度可寫為:

(6)

式中:M為節(jié)點質(zhì)量矩陣;P為外力;I為單元內(nèi)力。位移表示的平衡方程的顯式積分形式可寫為:

(7)

式中:t為時間;Δt為時間增量。保證求解穩(wěn)定性的時間步長由系統(tǒng)最高頻率及系統(tǒng)阻尼確定,穩(wěn)定步長計算如下式:

(8)

式中:ωmax為時間;ξ為系統(tǒng)阻尼。

選用ABAQUS有限元軟件顯式分析模塊及大變形模式,采用二維平面應(yīng)變問題假定進(jìn)行邊坡穩(wěn)定性問題的計算分析。采用C3D8R減縮積分單元對邊坡模型進(jìn)行網(wǎng)格劃分。以剛性基底的邊坡算例計算靜力情況下邊坡的安全系數(shù)對本文方法進(jìn)行對比驗證。邊坡坡高30 m,坡比1∶2。分別使用顯式及隱式有限元方法計算邊坡的安全系數(shù),其中顯式有限元方法采用大變形,隱式有限元方法采用小變形模式進(jìn)行計算。邊坡土體計算參數(shù)為:彈性模量E=100.0MPa,Poisson比υ=0.25,重度γ=20 kN/m3,黏聚力為c=30 kPa,內(nèi)摩擦角φ=20°,剪脹角ψ=0(不考慮剪脹性的影響),重力加速度g=9.81 m/s2。圖4(a)為邊坡計算網(wǎng)格劃分圖,圖4(b)～(c)為不同方法確定的邊坡極限狀態(tài)時網(wǎng)格變形計算結(jié)果對比。圖5為不同方法確定的邊坡最大位移與強度折減系數(shù)SRF的關(guān)系,給出了本文使用顯式有限元方法(EFEM)及使用隱式有限元(IFEM)和顯式有限差分方法(EFDM)的計算結(jié)果對比,其中隱式有限元方法采用SLOPE64邊坡穩(wěn)定分析程序進(jìn)行計算[21],顯式有限差分方法使用FLAC軟件進(jìn)行計算[22]。

圖4 顯式有限元、隱式有限元及顯式有限差分法計算出的邊坡安全系數(shù)及網(wǎng)格變形結(jié)果對比Fig.4 Comparison of safety factor and mesh deformation calculated by different methods

圖5 邊坡最大位移與強度折減系數(shù)SRF的關(guān)系Fig.5 Relationship between maximum displacement of slope and strength reduction factor SRF

從圖5可以看出邊坡最大位移隨強度折減系數(shù)SRF的增加逐漸增大,當(dāng)位移突然增大, 曲線出現(xiàn)明顯拐點時確定為邊坡的臨界破壞狀態(tài)。隱式有限元方法確定的邊坡安全系數(shù)FOS=1.36,顯式有限差分方法確定的邊坡安全系數(shù)FOS=1.43。對比兩種方法可看出,使用顯式有限元及強度折減方法同樣可以確定邊坡安全系數(shù),但由于顯式有限元方法采用大變形模式,因此邊坡位移比隱式方法要大。相比隱式方法,顯式有限元方法采用顯式時間積分方案求解動力學(xué)方程,材料及幾何非線性問題不需要進(jìn)行迭代,因此更適用于動力學(xué)及大變形問題的計算分析[23]。

2 邊坡地震穩(wěn)定性及破壞過程分析

2.1 計算條件

計算設(shè)置同上節(jié),土的計算參數(shù)取值詳見表1。在折減強度參數(shù)的同時將邊坡土體的動力學(xué)參數(shù)也進(jìn)行折減,如下式所示[17]:

(9) 表1 計算參數(shù)取值表Tab.1 Computation list of the parameter value

地震動輸入選取天然地震加速度時程,分別來自日本Kobe、美國Imperial Valley、Loma Prieta、Northridge以及土耳其Kocaeli地震[24]。將邊坡模型底面設(shè)置為加速度邊界,為了比較地震動特性對邊坡穩(wěn)定性的影響,將地震輸入的水平及豎向分量選取同一種加速度時程,其中豎向分量為水平向分量的2/3。邊坡兩側(cè)設(shè)置零加速度邊界(在該處加速度強制為零),顯式計算中約束加速度同樣也會約束位移自由度,因此邊坡兩側(cè)邊界盡量設(shè)置遠(yuǎn)離邊坡坡體區(qū)域。地震動具體信息詳見表2,地震加速度時程見圖6。圖7為5%阻尼情況下五種地震加速度時程的反應(yīng)譜曲線,可以看出五個地震波中美國的Northridge地震波峰值反應(yīng)加速度最大,土耳其Kocaeli地震波特征周期Tg最大。

表2 地震動時程特征信息Tab. 2 Ground motion time history characteristic information

圖6 五個地震加速度時程Fig.6 Five seismic acceleration time histories

圖7 五個地震時程的反應(yīng)譜曲線Fig.7 Response spectrum curves of five seismic time histories

2.2 計算結(jié)果

分別采用五種地震加速度時程進(jìn)行計算,得到地震之后邊坡最大位移與折減系數(shù)的關(guān)系以及最終確定的安全系數(shù)量值(見圖8)。從圖8結(jié)果可以看出加載Kocaeli地震加速度時程確定的邊坡的安全系數(shù)的數(shù)值最小(FOS=1.20),其他四種地震情況確定的邊坡安全系數(shù)結(jié)果相近。從五個地震加速度時程對應(yīng)的反應(yīng)譜曲線可以看出,Kocaeli地震的峰值反應(yīng)加速度最小,但其特征周期Tg最大,即地震峰值能量持時較其他四個地震時間長,因此在其他條件不變的情況下Kocaeli地震確定的邊坡安全系數(shù)是最小的。

圖8 五個地震得到的邊坡最大位移與折減系數(shù)的關(guān)系Fig.8 Relationship between maximum displacement and reduction coefficient obtained from five earthquakes

圖9為不同地震過程中邊坡坡頂水平向位移隨時間的變化,當(dāng)折減系數(shù)SRF=1.21,地震時間到20 s時邊坡坡體出現(xiàn)持續(xù)滑動破壞。

圖9 地震過程中邊坡坡頂水平向位移時程變化曲線Fig.9 The time-history curve of horizontal displacement of slope top during earthquake

圖10為Kocaeli地震過程中邊坡網(wǎng)格變形圖,可看出當(dāng)?shù)卣饡r間為26 s時刻邊坡坡體出現(xiàn)滑裂面,坡腳出現(xiàn)破壞之后出現(xiàn)整體滑動破壞。

圖10 Kocaeli地震過程中邊坡的破壞過程Fig.10 Progressive failure of slope during Kocaeli earthquake

3 滑坡地震穩(wěn)定性及破壞過程分析

滑坡地處甘肅省甘南州舟曲縣,滑坡地層結(jié)構(gòu)自上而下主要為第四系上更新統(tǒng)馬蘭黃土,中、上更新統(tǒng)滑坡堆積碎石土及殘坡積碎石土,中上志留統(tǒng)千枚巖、板巖及中泥盆統(tǒng)灰?guī)r為主。

圖11所示為北山上部滑坡體地層剖面及計算模型。將該滑坡視為二維平面應(yīng)變問題,分析滑坡在地震條件下的穩(wěn)定性及其破壞過程。

圖11 滑坡計算材料分區(qū)及單元網(wǎng)格劃分圖Fig.11 Slide calculation material partition and unit grid division diagram

邊坡一般為均質(zhì)土層構(gòu)成且高度較低,而滑坡的地層結(jié)構(gòu)更為復(fù)雜且高度較大,因此土體的剛度和強度需考慮初始應(yīng)力的影響,土體最大剪切剛度由下式確定[25]:

Gmax=κpa(σm/pa)n

(10)

式中:κ為初始剛度系數(shù);σm=(σ1+σ2+σ3)/3;pa為標(biāo)準(zhǔn)大氣壓;n為冪次系數(shù)。對于強風(fēng)化的殘坡積層土體強度參數(shù)c采用下式計算:

c=c0(σm/pa)n

(11)

式中c0為地表土體的黏聚力,隨著深度增加土體強度逐漸增大直至達(dá)到基巖的強度。地震荷載分別選取土耳其Kocaeli和美國Northridge地震加速度時程作為為橫向及豎向輸入分量,模型兩側(cè)為零加速度邊界條件。假定基巖在地震過程中不會出現(xiàn)塑性變形為完全彈性體?；鶐r上部風(fēng)化殘積層強度參數(shù)為初始黏聚力c0=40kPa,內(nèi)摩擦角φ=30°。土動力學(xué)參數(shù)為初始剛度系數(shù)κ=800,n=0.3,最大阻尼比Dmax=0.2,參考剪應(yīng)變γref=0.167。此外假定滑動面的摩擦系數(shù)μ=0.3。

圖12和圖13分別為折減系數(shù)SRF達(dá)到1.05后,地震過程中滑坡的廣義剪應(yīng)變及位移云圖。從計算結(jié)果可以看出,滑坡從形成滑動面到滑動是一個完整的破壞過程。該滑坡屬于牽引式滑動破壞,即前緣首先滑動破壞卸荷之后帶動后緣產(chǎn)生滑動。滑坡體并未完全滑動到坡腳溝谷底部,而是一半停留在山坡中部,只有一部分滑入谷底(最大滑動距離為1.46 km)。通過分析可知,顯式有限元及動力學(xué)大變形方法的結(jié)果可分析地震情況下誘發(fā)滑坡的破壞過程。

圖13 地震過程中滑坡體的位移云圖Fig.13 Contourof landslide displacement during earthquake

4 結(jié) 論

本文提出了一種地震條件下邊坡及滑坡穩(wěn)定性分析方法,并對邊坡及滑坡的地震穩(wěn)定性及破壞過程進(jìn)行了計算分析。

1) 地震情況下邊坡的失穩(wěn)破壞是一個從坡體變形到出現(xiàn)局部破壞直至整體失穩(wěn)滑動的發(fā)展過程,并非是達(dá)到極限狀態(tài)后突然發(fā)生的,須基于土動力學(xué)理論使用動力學(xué)大變形計算方法分析該問題。

2) 通過本文研究證明,顯式有限元及大變形計算方法可直接確定邊坡的地震安全系數(shù),同時還可分析地震過程中邊坡的動力響應(yīng)及破壞過程。相對于傳統(tǒng)隱式分析方法,本文所提出的方法在靜力及小變形條件下分析邊坡穩(wěn)定性具有較大的誤差,但在研究地震作用下邊坡的穩(wěn)定性及地震誘發(fā)滑坡災(zāi)害問題具有顯著優(yōu)勢。

3) 通過一個滑坡工程實例的分析,說明本文提出的方法同樣適用于地震條件下滑坡穩(wěn)定性問題的計算分析,并可進(jìn)一步確定滑坡的滑動距離及破壞范圍。本文未考慮降雨、地下水及三維地形等因素對邊坡及滑坡穩(wěn)定性的影響,下一步將建立三維地形模型考慮降雨及地下水位的影響對地震誘發(fā)滑坡災(zāi)害問題進(jìn)行計算分析工作,以期為防災(zāi)減災(zāi)工作提供參考依據(jù)。