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

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

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

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

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

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

1 理論及方法

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

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

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

(1)

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

(2)

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

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

(3)

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

(4)

式中:A為加速度;t為振動(dòng)時(shí)間。

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

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

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

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

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

1.2 計(jì)算方法

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

(5)

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

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

(6)

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

(7)

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

(8)

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

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

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

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

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

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

2.1 計(jì)算條件

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

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

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

表2 地震動(dòng)時(shí)程特征信息Tab. 2 Ground motion time history characteristic information

圖6 五個(gè)地震加速度時(shí)程Fig.6 Five seismic acceleration time histories

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

2.2 計(jì)算結(jié)果

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

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

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

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

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

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

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

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

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

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

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

Gmax=κpa(σm/pa)n

(10)

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

c=c0(σm/pa)n

(11)

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

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

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

4 結(jié) 論

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

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

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

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