彭雪峰，朱亞林,2，馬 馳，檀 昆

(1.合肥工業(yè)大學(xué) 土木與水利工程學(xué)院，合肥 230009； 2.土木工程結(jié)構(gòu)與材料安徽省重點實驗室，合肥 230009)

在工程中，邊坡穩(wěn)定性分析是伴隨著土壓力和地基承載力理論一起發(fā)展起來的[1-4].庫倫和朗肯建立了一系列關(guān)于土壓力的計算方法，學(xué)者們將其運用在邊坡穩(wěn)定分析方面并形成了體系，即極限平衡法[5].極限平衡法的基本特性是：在摩爾-庫倫破壞準則基礎(chǔ)上建立靜力平衡條件,針對土體破壞時力的平衡狀態(tài)來求問題的解[6].然而，在實際中，很多情況是靜不定狀態(tài).為了使原本靜不定問題變成一個靜力問題，極限平衡法會引入一些簡化：設(shè)定滑裂面形狀為折線、圓弧、對數(shù)螺旋線等；放松靜力平衡要求，求解過程中僅滿足部分力和力矩的平衡要求；對多余未知數(shù)的數(shù)值或分布形狀作假定[7-8].但由于該方法沒有考慮到土體的本構(gòu)性質(zhì)，不能求出失穩(wěn)時土體各處應(yīng)力和應(yīng)變，無法反映土坡的破壞機理.所以，采用這種方法求出的僅僅是一種綜合性近似解答.

從進入20世紀70年代之后，隨著有限元理論方法的完善和計算機水平的提高，使用嚴格的應(yīng)力應(yīng)變分析結(jié)構(gòu)成為可能[9].一般來說，數(shù)值分析都是和強度折減法一起使用.強度折減法對邊坡上的強度指標(摩擦系數(shù)φ和黏聚力c)進行折減，并且通過以力和位移作為收斂標準得到一個新的強度指標，再與原先的強度指標做比值即為安全系數(shù)[10].從強度折減法的基本概念出發(fā)，對于同一個邊坡，理論上只有一個折減后的強度指標.于是，當進行動力分析時，會導(dǎo)致無法直接通過折減強度得到不同時間步上的安全系數(shù).

為了能夠了解邊坡在地震中抵抗滑動破環(huán)的能力.文獻[11]利用矢量和法計算得到邊坡安全系數(shù)時程圖.矢量和法是基于邊坡真實應(yīng)力狀態(tài)，將已知滑動面上的滑動矢量和抗滑矢量在相同下滑方向作投影并求和，兩者的比值為安全系數(shù)[12].以往學(xué)者在利用矢量和法進行邊坡穩(wěn)定性分析時，往往只是將邊坡滑動面處理成剪切破壞.但是，由文獻[13]論述：在地震作用過程中，邊坡主要破壞形式是由坡頂向下一定深度內(nèi)的拉破壞和坡腳向上延展的剪切破壞共同組成貫通破裂面.為了更加符合在地震作用下邊坡的破壞模式，本文結(jié)合拉-剪破壞強度準則對矢量和法進行優(yōu)化，使得到的安全系數(shù)時程圖更能反映邊坡真實抗震性能.

1 基于拉-剪破壞的矢量和方法

1.1 考慮抗拉強度矢量和方法安全系數(shù)公式推導(dǎo)

1.1.1 矢量和法的討論

矢量和法的優(yōu)勢是基于邊坡真實應(yīng)力狀態(tài)，而與之相對應(yīng)的安全系數(shù)稱為矢量和法安全系數(shù)，其具體定義：在邊坡潛在滑動面上，極限抗滑力矢量在整體下滑趨勢方向上投影代數(shù)和R與總下滑力矢量在相同方向上投影代數(shù)和T的比值[14].受力分析如圖1所示，邊坡計算區(qū)域內(nèi)潛在滑動面S已知，σs、σn、στ為滑動面A點的應(yīng)力矢量、法應(yīng)力矢量和剪應(yīng)力矢量，σ′s、σ′n、σ′τ為滑動面B點的極限抗滑應(yīng)力矢量、法應(yīng)力矢量和剪應(yīng)力矢量，m為潛在滑動體的整體下滑趨勢方向.滑動面上各點的應(yīng)力矢量和極限抗滑應(yīng)力矢量在m上投影求和，并做比值就得到相應(yīng)的安全系數(shù).

具體公式如下：

(1)

式中，分母是邊坡真實應(yīng)力矢量和在整體下滑趨勢方向做投影.分子代表一個邊坡在滑動面S上極限抗滑能力，如何量化邊坡極限抗滑能力是討論重點.

先對“σ′n+σ′τ”進行分析，極限抗滑能力是由極限抗滑法應(yīng)力矢量σ′n和剪應(yīng)力矢量σ′τ共同提供.在矢量和法中，極限法應(yīng)力矢量σ′n的大小等于該點巖土法應(yīng)力矢量σn，方向與之相反；而極限剪應(yīng)力矢量σ′τ的大小滿足Mohr-Coulomb強度準則(τ=σntanφ+c)，方向為該點在滑動面上剪應(yīng)力的反方向.接著，對于整體下滑趨勢方向m而言，有3種計算方法：a)利用從潛在滑動面剪入點到剪出點的連線作為下滑方向；b)潛在滑動面剪應(yīng)力矢量和的方向；c)潛在滑動面極限抗剪應(yīng)力矢量和的反方向[15].假設(shè)：當滑動體產(chǎn)生實際的位移時，定義此時的滑動方向為整體下滑趨勢方向m.為了滿足這種情況需要將抗滑應(yīng)力達到最大值(極限抗滑強度)，所以，本文選擇方法c).

圖1 沿潛在滑動面受力分析圖

從上文中可以看出，無論是對于極限抗下滑應(yīng)力矢量σ′s還是整體下滑趨勢方向m，對邊坡潛在滑動面上的抗滑強度準則選取都有重要的作用.在靜力學(xué)邊坡穩(wěn)定分析中，潛在滑動面主要是受剪破壞[16].但經(jīng)過調(diào)查發(fā)現(xiàn)，邊坡在地震的作用下，坡頂區(qū)域會發(fā)生拉裂破壞，有些巖石甚至被拋出[13].如果僅僅認為在滑動面上極限抗滑強度全部來源于抗剪強度，無疑會造成一定的誤差.基于以上關(guān)于矢量和法的思考，下面將給出考慮抗拉強度的矢量和法的推導(dǎo)過程和計算方法.

1.1.2 拉-剪破壞下矢量和法安全系數(shù)

首先，利用數(shù)值分析手段得到邊坡應(yīng)力場和潛在滑動面，對滑動面所在單元進行應(yīng)力矢量計算.任取一個單元P，由彈性力學(xué)基本公式可得

(2)

式中：σi表示潛在滑動面上一點應(yīng)力張量，σm、σn、στ分別表示在該點應(yīng)力矢量，法應(yīng)力矢量和剪應(yīng)力矢量，n表示該點平面的單位法向量.

1.1.2.1 矢量和安全系數(shù)公式推導(dǎo)

基于矢量和法安全系數(shù)關(guān)鍵是求出下滑面極限抗滑能力.動荷載作用使得邊坡上部形成拉裂縫.如圖2所示，設(shè)滑動面區(qū)域中包含受剪區(qū)域有pi個單元，受拉區(qū)域有pu個單元.H1到H2為受剪區(qū)，按照1.1.1節(jié)所述方法計算滑裂面極限抗滑強度，其滑裂面上的極限抗滑矢量σ′mi、法向應(yīng)力矢量σ′ni、剪應(yīng)力矢量σ′τi滿足如下公式：

(3)

式中:φ、c分別為巖土內(nèi)摩擦角和黏聚力.

圖2 邊坡極限抗滑能力受力示意

H2到H3為受拉區(qū)，其極限抗滑應(yīng)力矢量σ′mu中的法向應(yīng)力矢量σ′nu大小由巖土抗拉強度確定，極限抗剪應(yīng)力矢量σ′τu大小等于該點剪應(yīng)力矢量，公式如下：

(4)

式中ftu為巖土抗拉強度.

將極限抗滑應(yīng)力矢量在整體下滑趨勢方向反方向做投影，并求和可得

(5)

由式(3)、(4)、(5)可得

(6)

對于總下滑力，由邊坡潛在滑動面上應(yīng)力在整體下滑趨勢方向上求解代數(shù)和可得

(7)

(8)

1.1.2.2 整體下滑趨勢方向確定

為了將兩個矢量和能夠在相同方向進行投影并做比較，需要判斷出滑動體在滑動面上可能的下滑趨勢方向.為了滿足“相對運動趨勢與靜摩擦力相反”的基本物理意義，當已知邊坡坡頂處是受拉破壞，原先下滑方向公式需要做相應(yīng)的調(diào)整.其中，對于受剪滑動面上，滑動體是沿著滑動平面的切線方向運動，所以，利用極限抗剪強度σ′τi作為計算整體下滑方向因子；對于受拉滑動面，由于邊坡潛在滑動體是背離邊坡方向運動，利用滑動面上極限受拉強度σ′nu作為計算下滑方向因子，即可以得

(9)

將式(3)、(4)代入式(9)可得

(10)

1.2 滑動面計算

一般邊坡滑動失穩(wěn)是沿著剪應(yīng)變增量最大部位發(fā)生，所以，通過剪應(yīng)變增量去確定剪切破壞面位置，再利用地震作用下所得受拉破壞區(qū)域來判斷受拉滑動面高度，兩者共同形成貫通復(fù)合滑動面[17].本文利用FLAC程序計算的邊坡應(yīng)力場，其應(yīng)力信息是存儲在單元格中.首先，利用剪切應(yīng)變增量云圖確定滑動面大致范圍(見圖3)，自定義單位層高h，從坡腳H1進行逐層搜索直到坡頂H3并得到每層剪切應(yīng)變增量最大值點；接著，根據(jù)單元應(yīng)力狀態(tài)找出受拉屈服單元，將這些單元所在區(qū)域最大豎直高度作為拉裂縫深度H2；最后，從H1到H3將所有最大值點通過n階多項式(如式(11))進行擬合可以得到剪切破壞模式下的滑動面曲線l1，若利用H1到H2高度范圍內(nèi)最大值點進行擬合并結(jié)合受拉裂縫豎直線段(H2～H3)便可得到拉-剪破壞模式下的滑動面曲線l2.這里有兩點需要說明：1)剪切應(yīng)變增量最大值的坐標是所在單元中心坐標，可能會出現(xiàn)擬合之后曲線所對應(yīng)坐標偏移離開原先的單元中心坐標，這時需要進行調(diào)整；2)由于FLAC軟件中應(yīng)力狀態(tài)都儲存在單元中，受拉裂縫深度的精度取決于網(wǎng)格劃分的尺寸，單元劃分越小精度越高.在實際中，巖土顆粒小于網(wǎng)格尺寸，所以，本文采用方法得到的受拉裂縫深度會稍大于實際的拉裂縫深度，提供一定冗余度.

擬合多項式表達式：

(11)

式中ai為待求系數(shù).

圖3 滑動面搜索方法示意

1.3 地震作用下邊坡穩(wěn)定性分析過程

采用數(shù)值分析程序FLAC可以精準模擬地震波在土體中傳播并且得到各個單元應(yīng)力.首先，根據(jù)邊坡物理、幾何參數(shù)建立相應(yīng)模型，在輸入地震荷載之后進行動力求解.然后，提取每個時刻各個單元上應(yīng)力大小、方向、剪切應(yīng)變增量以及單元狀態(tài).再使用MATLAB自編程序得到滑動面單元上應(yīng)力信息，最后，利用1.1中優(yōu)化后的矢量和方法求解得到每個時刻的安全系數(shù)并對邊坡的穩(wěn)定性進行評價分析，具體流程見圖4.

圖4 邊坡穩(wěn)定性分析流程

1.4 算法驗證

為了驗證本文得到的邊坡安全系數(shù)算法合理性，利用強度折減法與本文新方法進行對比分析邊坡動穩(wěn)定性.由于強度折減法局限性，采用擬靜力法計算地震作用效應(yīng)，水平地震慣性力代表值按照相應(yīng)規(guī)范標準公式，設(shè)計烈度為7度，水平向設(shè)計地震加速度代表值為0.98[18].邊坡巖土采用Mohr-Coulomb準則，詳細參數(shù)：坡高H=5 m，坡比i為1∶1.50，黏聚力c=50 kPa，內(nèi)摩擦角φ=30°，剪切模量G=10 GPa，體積模量K=30 GPa，抗拉強度T=5 kPa.

圖5為模型單元狀態(tài)分布圖.可以看出，邊坡腹部單元主要受到剪切破壞，坡頂出現(xiàn)了受拉破壞單元，這與文獻[13]分析的結(jié)果相吻合，表明在動荷載中坡頂破壞形式主要是由拉裂縫組成.

圖5 邊坡單元狀態(tài)圖

將邊坡單元剪切應(yīng)變增量數(shù)值提取出來，以自定義層高為基本單位進行逐層搜索，將搜索得到最大值點利用式(11)擬合成多項式函數(shù)，為了防止曲線的上凸，選用多項式的最高次數(shù)n為3，具體系數(shù)見表1.

表1 多項式系數(shù)

圖6為兩種不同破壞模式下擬合曲線.其中，剪切破壞模式下的滑裂面曲線擬合高度范圍為從坡腳一直到坡頂，即豎向范圍：h1=0.2 m～h3=5 m；對于拉-剪破壞模式下的滑裂面曲線擬合從坡腳到受拉裂縫深度處，即豎向范圍為h1=0.2 m～h2=4.3 m，剩余部分按照受拉豎直裂縫處理，即范圍為h2=4.3 m～h3=5 m.

由擬合曲線逐層提取相應(yīng)單元應(yīng)力場，并帶入公式可以得到下滑趨勢角度和矢量和法安全系數(shù)，再利用FLAC自帶的強度折減法計算安全系數(shù)作為參考對比，具體數(shù)據(jù)見表2.可以看出，本文改進后拉-剪破壞矢量和法的安全系數(shù)和強度折減法只差了0.01，表明了本文方法的正確性，而剪切破壞的安全系數(shù)相對來說，比另外兩種方法數(shù)值偏大.

圖6 擬合曲線示意

表2 邊坡安全系數(shù)與下滑趨勢角

2 算例分析

2.1 邊坡模型概況

邊坡模型為某一巖質(zhì)邊坡，邊坡高度為41 m，邊坡按照1∶1.00放坡，詳見圖7.為了滿足精度要求，模型左邊界離坡腳距離20 m，坡頂?shù)接疫吔缇嚯x為35 m，下邊界到坡腳的距離為20 m.

圖7 邊坡二維計算模型

巖土材料采用Mohr-Coulomb強度準則，第1巖土層(1#)為全風化花崗片麻巖，第2巖土層(2#)為強風化花崗片麻巖，具體材料參數(shù)如表3所示.其中模型四周設(shè)有自由邊界場，阻尼采用瑞利阻尼，其中通過計算得到阻尼比為0.05，最小中心頻率為4.065 Hz，計算過程參見文獻[19-20].

表3 土體力學(xué)參數(shù)

通過對某地區(qū)100年超越頻率為2%的基巖期望反應(yīng)譜曲線進行擬合而得到人工波(如圖8)，其加速度峰值為0.28g.在計算過程中，在模型底部以水平方向輸入加速度峰值為0.28g，豎直方向輸入加速度峰值為水平方向的2/3，地震總共歷時20 s.

圖8 輸入地震時程曲線

2.2 邊坡滑動界面分析

圖9為邊坡在不同地震時刻單元狀態(tài)圖.可以看出，當t=1 s時，單元狀態(tài)主要是在坡腳處出現(xiàn)剪切破壞單元并且向坡內(nèi)方向延伸，只有在坡頂處出現(xiàn)個別的拉破壞單元.當t=2.2 s時，地震波震蕩明顯加劇，坡頂拉破壞單元不斷擴大并且與剪切破壞單元連同形成貫通界面.在圖9(c)中，t=6 s時地震波達到峰值，邊坡的腹部大部分已經(jīng)發(fā)生了剪切破壞，并且坡頂?shù)睦茐膮^(qū)域明顯增加.圖9(d)為地震t=12 s時刻的單元狀態(tài)圖，整個邊坡大部分單元都發(fā)生過了剪切破壞，坡頂區(qū)域也從地震最開始的2.2 s出現(xiàn)的個別拉破壞單元一直到地震持續(xù)時間末端被拉破壞單元的完全覆蓋.

通過前文對邊坡單元狀態(tài)的分析，可以根據(jù)拉破壞單元得到受拉裂縫深度.利用受拉裂縫深度選取受拉區(qū)域高度，結(jié)合相應(yīng)的單元剪切應(yīng)變增量最大值點可以分別擬合出兩種破壞模式下的潛在滑動面曲線，擬合多項式最高次數(shù)為3次.由于篇幅問題，僅僅展示了1,5,10,15,20 s的擬合曲線系數(shù)，見表4,5.

表4 擬合曲線系數(shù)(剪切破壞)

表5 擬合曲線系數(shù)(拉-剪破壞)

圖9 邊坡單元狀態(tài)動態(tài)圖

圖10,11為5個時間節(jié)點兩種破壞模式下的邊坡滑動面擬合曲線.通過對小部分曲線群放大可以看出，隨著時間的推移，地震滑弧有一定的變化.這是由于地震持續(xù)過程中，滑動體的位移是一個逐漸增大的過程，使得剪切應(yīng)力增量整體向坡內(nèi)平移，進而求得的擬合曲線也有向坡內(nèi)走移的趨勢.

圖10 剪切破環(huán)擬合曲線示意

圖11 拉-剪破壞擬合曲線示意

根據(jù)不同的擬合曲線，分別提取相應(yīng)單元并求出邊坡下滑趨勢角度(如圖12).可以很容易看出，整體上根據(jù)拉-剪滑裂面求出的下滑矢量角要大于剪切滑裂面的下滑矢量角；在地震過程中，最大差值達4.5°，出現(xiàn)在19.7 s時刻；最小差值為0.13°.這是由于坡頂拉裂縫的出現(xiàn)改變了原先滑動面的軌跡使得下滑角度增大，不利于邊坡穩(wěn)定性.

圖12 邊坡整體下滑趨勢角度時程

2.3 邊坡穩(wěn)定安全系數(shù)時程分析

圖13,14分別為剪切破壞和拉-剪破壞的滑動合力和抗滑合力.從整體趨勢看，兩種不同方法的合力變化趨勢和地震波波動趨勢一致，在4～6 s、和8～10 s波動最大；在數(shù)值上，拉-剪破壞下滑動合力和抗滑合力相較剪切破壞都要偏小.一方面是由于投影矢量角度變大，另一方面是由于在拉裂破壞面引入了抗拉強度條件.

圖13 邊坡的動態(tài)響應(yīng)規(guī)律(剪切破壞)

圖14 邊坡的動態(tài)響應(yīng)規(guī)律(拉-剪破壞)

利用兩種不同的方法分別對邊坡的安全系數(shù)進行計算，圖15為兩種方法的安全系數(shù)時程圖.整體上看出，兩種方法所得到安全系數(shù)趨勢一致，表明新方法在時程分析上的可靠性；其中拉-剪破壞的安全系數(shù)大小波動范圍為1.047～1.574，剪切破壞的安全系數(shù)大小波動范圍為1.077～1.622，前者比后者小了2.8%～3%.比較兩者大小可以得到，增加了拉裂面的安全系數(shù)比原先方法得到安全系數(shù)要小，表明以往用矢量和方法進行邊坡動穩(wěn)定性分析時，會高估邊坡的抗震性能，不利于邊坡的抗震設(shè)計.

圖15 邊坡安全系數(shù)時程圖

當?shù)卣饡r間為11.4 s時，兩種方法所得安全系數(shù)為最小值，拉-剪破壞的安全系數(shù)為1.05，剪切破壞的安全系數(shù)為1.07；當?shù)卣饡r間為9.4 s時，所得安全系數(shù)為最大值，其中前者為1.57，后者為1.62.圖16為t=11.4和9.4 s時刻邊坡的加速度矢量圖.在地震時間為11.4 s時刻，邊坡的加速度矢量均朝向邊坡外側(cè)，這是對邊坡動穩(wěn)定性最不利的狀態(tài)，觀察地震波加速方向，在下一秒加速度方向出現(xiàn)了明顯的改變，安全系數(shù)得到一定的增加，所以，此時計算得到的安全系數(shù)最小，是邊坡最不利的狀態(tài).同樣，當t=9.4 s時，加速度矢量的方向大致是朝向邊坡內(nèi)部，這種狀態(tài)有利于邊坡的穩(wěn)定性，并且在下一秒地震波的加速度發(fā)生了改變，安全系數(shù)又得以恢復(fù)，所以，t=9.4 s時計算的安全系數(shù)達到最大值.

2.4 邊坡動力穩(wěn)定性綜合評價

一般來說，邊坡動力穩(wěn)定性評價的具體指標在國內(nèi)外還沒有形成統(tǒng)一的標準，其中比較主流兩種指標為永久位移和安全系數(shù).當前，基于安全系數(shù)時程圖進行邊坡穩(wěn)定性評價方法大體分為平均安全系數(shù)[21]、最小動安全系數(shù)[22]和基于可靠度理論的動力安全系數(shù)[21].平均安全系數(shù)是指利用安全系數(shù)時程圖各個點的平均安全系數(shù)作為評價的指標；最小動安全系數(shù)就是以安全系數(shù)時程圖中的最小安全系數(shù)作為評價的指標；而基于可靠度理論的動力安全系數(shù)Kf的表達式如下：

(12)

式中：N為安全系數(shù)時程的時間點總數(shù)，K(ti)為ti時間點上對應(yīng)的安全系數(shù)，β為邊坡失效概率的可靠度指標.

本文將失效概率設(shè)置為1%，在其失效概率情況下對應(yīng)的可靠性指數(shù)為β=2.33，再結(jié)合各個時間節(jié)點的安全系數(shù)帶入式(12)后得到相應(yīng)的邊坡動力安全系數(shù)，并將得到結(jié)果編制成表6.可以清晰看出，3種對邊坡動力穩(wěn)定性評價指標均大于1，說明在當前地震荷載情況下，邊坡處于安全范圍內(nèi)不會發(fā)生失穩(wěn)破壞；對比拉-剪破壞的安全系數(shù)和剪切破壞的安全系數(shù)，所有指標均是前者較小，其中動力安全系數(shù)前者為1.07，后者為1.11，兩者相差3%.

表6 邊坡整體動安全系數(shù)計算結(jié)果

3 結(jié) 論

1)通過算例將強度折減法、剪切破壞矢量和法和拉-剪破壞矢量和法3種方法對比分析，結(jié)果表明，在動荷載情況下拉-剪破壞矢量和法得到安全系數(shù)更加精確.

2)在地震作用下，邊坡潛在滑動面會發(fā)生變化，并且隨著地震持續(xù)時間增加，滑動面向著邊坡內(nèi)側(cè)走移.

3)拉-剪破壞模式下邊坡整體下滑趨勢角度相對剪切破壞的矢量角度有所增大，表明坡頂?shù)睦芽p改變了原先的下滑趨勢角度.

4)拉-剪破壞的安全系數(shù)要低于剪切破壞的安全系數(shù)，說明原先的矢量和方法高估在動荷載作用下邊坡的抗震性能，這會不利于邊坡的抗震設(shè)計.