王 壯, 蘇 雷，*, 時 偉, 張擁軍, 凌賢長, 2

(1.青島理工大學(xué) 土木工程學(xué)院，青島 266525；2.哈爾濱工業(yè)大學(xué) 土木工程學(xué)院，哈爾濱 150090)

山體邊坡常常發(fā)生各種破壞現(xiàn)象，其中滑坡是最常見的破壞形式。山體滑坡具有巨大的破壞力，其發(fā)生時攜帶大量巖土傾瀉而下，對人們的生命與財產(chǎn)安全造成巨大的威脅。誘發(fā)山體滑坡的因素眾多，如降雨、地震、人工開挖和工程爆破等，其中降雨、地震等自然現(xiàn)象是無法避免的，許多滑坡的產(chǎn)生也與此相關(guān)，如2008年汶川地震發(fā)生后，誘發(fā)了數(shù)以萬計的滑坡災(zāi)害，造成巨大的人員傷亡及財產(chǎn)損失[1]；2006年的馬達嶺滑坡由連續(xù)降雨引發(fā)，該滑坡破壞了大量的農(nóng)田及道路，所幸未造成人員傷亡[2]；2017年九寨溝發(fā)生MS7.0地震，觸發(fā)大量的滑坡災(zāi)害，地震共造成25人死亡、6人失蹤、525人受傷、1萬多人受到影響，73 671間房屋受損，經(jīng)濟損失達80.43億元[3]。

滑坡的分析方法主要有極限平衡法與數(shù)值分析法，其中數(shù)值分析法考慮土體間的應(yīng)力-應(yīng)變關(guān)系，可以在短時間內(nèi)獲得較多的體系反應(yīng)量，因此眾多研究人員采用數(shù)值模擬的方式對滑坡進行研究。如楊笑男等[4]介紹FLAC在特殊性巖土滑坡災(zāi)害治理評價中的一般方法，應(yīng)用FLAC對典型黏土質(zhì)硅藻土邊坡的滑坡防治措施進行優(yōu)化設(shè)計和數(shù)值模擬分析，提出經(jīng)濟合理的綜合防治建議。張小雪等[5]采用顆粒流(PFC)模擬方法，對黏土試樣的雙軸試驗進行數(shù)值模擬，并結(jié)合摩爾-庫侖破壞準則，得到細觀參數(shù)對宏觀特性的影響規(guī)律，模擬黏性土坡在自重作用下變形破壞的全過程。王翔南等[6]克服擴展有限元法(XFEM)需要預(yù)先設(shè)定滑裂面的起始位置，以及滑裂面前端擴展方向的判斷精度較低等缺陷，使擴展有限元法可以自動判斷裂縫發(fā)生位置，提高其對滑坡破壞過程的模擬精度。王宵等[7]采用離散元軟件UDEC4.0和有限差分軟件FLAC3D，分別建立二維和三維數(shù)值模型模擬邊坡破壞過程，結(jié)果表明巖石類邊坡變形破壞主要經(jīng)歷三個時期，李守巨等[8]采用有限元強度折減方法評估含節(jié)理巖石邊坡的穩(wěn)定性，計算臨界滑移面的安全系數(shù)，通過Drucker-Prager破壞準則判斷巖石邊坡有限元模型中每個單元的滑移狀態(tài)，采用準靜態(tài)方法模擬巖石邊坡地震作用荷載，討論地震動力放大系數(shù)對邊坡抗滑穩(wěn)定性安全系數(shù)的影響。

考慮到滑坡具有大變形的特點，F(xiàn)LAC計算程序采用顯式差分法求解微分方程，在大變形問題的分析方面具有獨特的優(yōu)勢[9]，因此本文采用FLAC計算程序模擬白泥溝滑坡，計算分析其在自重、降雨及地震工況下的穩(wěn)定性。

1 工程地質(zhì)概況及滑坡特征

白泥溝滑坡位于巧家縣白鶴灘鎮(zhèn)迆博社區(qū)躍進4組唐家山斜坡處，距離巧家縣縣城南側(cè)5 km。白泥溝滑坡整體形態(tài)呈扇形，后緣清晰，兩側(cè)邊界較為明顯?；麦w寬約350 m，長約300 m，滑坡面積約7.58×104m2?；潞缶墶⑶熬墭烁吒卟?17 m，滑坡體厚度最大為45 m，平均厚度32 m，滑坡體體積約2.5×106m3，屬于巨型厚層基巖滑坡，滑坡區(qū)全貌如圖1所示。

滑坡前緣存在采石場，開挖后造成卸荷作用導(dǎo)致滑坡后緣出現(xiàn)張拉性裂縫?；聟^(qū)出露物質(zhì)為殘坡積土及灰?guī)r，其中表層殘坡積土由棕紅色粉質(zhì)黏土與碎石土組成，厚度分布不均。下覆基巖以灰?guī)r為主，由鉆孔巖芯可見，基巖具有強風(fēng)化特征，巖體受構(gòu)造擠壓極其破碎，巖體透水性強?；戮植繆A雜泥化軟弱夾層，厚度為2～5 m，分布在地下25～40 m處，夾層傾角20～25°，接近山體斜坡傾角。鉆探結(jié)果顯示，鉆孔深度內(nèi)未見地下水，但由于雨水沿滑坡后緣裂縫入滲，導(dǎo)致泥巖層含水量較高，且該泥巖層具有良好的隔水效果，在雨水浸潤下，抗剪強度明顯降低?；滤诘貐^(qū)地震活動頻繁，地震作用增大坡體的下滑力，從而促進滑體沿泥巖層產(chǎn)生滑動現(xiàn)象。

2 FLAC數(shù)值模型建立

2.1 模型的建立及參數(shù)選取

模型依據(jù)白泥溝滑坡地質(zhì)勘察剖面圖建立，模型長約389 m，模型坡頂至模型底部約215 m。模型上部覆蓋粉質(zhì)黏土，中間分布軟弱夾層，地質(zhì)勘察為黏土層，基巖為強風(fēng)化灰?guī)r，滑體由粉質(zhì)黏土及強風(fēng)化灰?guī)r兩部分組成。根據(jù)地質(zhì)勘察報告[10]，各巖土層物理力學(xué)特性參數(shù)見表1，各層材料均采用摩爾-庫侖屈服準則。網(wǎng)格劃分考慮地震波在模型中傳播的數(shù)值精度[10-12]，最大網(wǎng)格邊長為1.6 m，x，y方向分別劃分248，179個節(jié)點，共44 392個網(wǎng)格。靜力計算下，模型兩側(cè)采用x向約束，模型底部采用x，y雙向約束。FLAC數(shù)值模型如圖2所示。

圖1 白泥溝滑坡全貌[10]

圖2 FLAC數(shù)值模型

在數(shù)值模擬中，F(xiàn)LAC采用強度折減法計算坡體安全系數(shù)。自重工況下的滑坡穩(wěn)定性計算采用天然狀態(tài)下的巖土物理力學(xué)特性參數(shù)；考慮到勘察報告中指出在鉆探范圍內(nèi)未見地下水，所以本文未考慮孔隙水壓力上升及巖土體內(nèi)部滲流，在降雨工況下，將滑體材料參數(shù)改為飽和參數(shù)，用于計算坡體的穩(wěn)定性；地震作用下的計算只考慮天然狀態(tài)下的材料參數(shù)。

表1 白泥溝滑坡模型巖土物理力學(xué)特性參數(shù)

2.2 模型合理性驗證

工程中為查明滑坡的分層、分級情況及滑面深度，掌握滑坡的變形速率及滑動方向，采用測斜儀對滑坡進行了深孔位移監(jiān)測，記錄坡體內(nèi)部水平位移情況。依據(jù)實際工程中的深孔位置，在數(shù)值模型中記錄相同位置的水平位移，并在自重作用下對模型進行計算。監(jiān)測和計算的水平位移對比如圖3所示。計算與監(jiān)測的水平位移整體趨勢一致，其中滑坡表面位移量一致，滑動面處(22～25 m)位移均出現(xiàn)增大現(xiàn)象，表明數(shù)值模型能夠再現(xiàn)滑坡不同深度的位移。但是由于工程正處于施工中，受人工作業(yè)的影響，且監(jiān)測位移量在毫米量級，致使數(shù)值模擬結(jié)果與監(jiān)測結(jié)果存在一定的差異。

2.3 地震波處理

在數(shù)值模擬中，采用1940年Imperial Valley地震中的El Centro波，地震波如圖4所示。地震波最大加速度0.349g，發(fā)生在2.12 s，地震波時間間隔0.02 s，總時長53.74 s。由地震波得到的最終速度與位移均不為0，如果直接將地震波施加到模型中進行計算，將導(dǎo)致FLAC模型在地震波結(jié)束后依舊存在持續(xù)速度或殘余位移[9, 11-12]，因此需要將地震波進行基線修正處理。基線修正可以借助FLAC中的Fish語言功能實現(xiàn)，修正前、后的地震波速度、位移時程曲線如圖5、圖6所示。

2.4 地震波輸入及阻尼的選取

地震波的施加方式有兩種[9, 11]，第一種是直接在模型底部施加加速度或速度時程，該方式適用于剛性地基，即模型底部為模量較大的材料。第二種是將應(yīng)力或力時程施加到模型中，該方式需要在模型底部施加靜態(tài)邊界。第一種方式存在潛在的缺點[11]，即模型底部的運動是完全指定的，相當于在底部施加了一個位移邊界，對于一般模擬計算是不適用的。本文采用第二種方式實現(xiàn)模型在地震工況下的模擬計算，應(yīng)力時程可以由速度時程通過公式σn=factor(ρCp)vn計算得到。式中：σn為應(yīng)力；ρ為模型底部材料的密度；Cp為縱波在模型底部材料中的傳播速度；vn為速度；factor需要進行試算，將得到的應(yīng)力時程輸入到模型中，并記錄模型底部的速度響應(yīng)時程，與速度時程進行對比，得到最合適的factor值。經(jīng)過試算，本文中factor值為1.0。模型底部施加x和y向靜態(tài)邊界，使用靜態(tài)邊界可以較好地吸收地震反射波能量，防止地震波的反射；模型兩側(cè)施加自由場邊界[9, 11]。

圖4 地震波加速度時程曲線

地震工況下的模擬計算還需要設(shè)置合適的阻尼，F(xiàn)LAC提供三種阻尼類型，分別為瑞利阻尼、局部阻尼和滯后阻尼。相較另兩種阻尼，瑞利阻尼可以近似地反映巖土體具有的頻率無關(guān)性[10]，因此本次計算選用瑞利阻尼求解[13-14]，設(shè)置阻尼比為0.5%，中心頻率為1.8 Hz。

3 數(shù)值計算結(jié)果分析

3.1 自重及降雨工況下結(jié)果分析

3.1.1 水平位移

FLAC采用強度折減法計算邊坡安全系數(shù)，計算得到的自重及降雨條件下的安全系數(shù)分別為1.35，1.00，滑坡從穩(wěn)定狀態(tài)轉(zhuǎn)變?yōu)榍贩€(wěn)定狀態(tài)，說明降雨對滑坡的穩(wěn)定性產(chǎn)生了較大影響。圖7為自重工況下的水平(即x向)位移云圖。由圖可知，在自重工況下，滑坡的水平方向最大位移為0.25 m，發(fā)生在滑體表面上部的黏土層，與該位置坡度較大有關(guān)；滑體其余位置位移量在0.05～0.15 m之間。圖8為降雨工況下水平位移云圖。由圖可知，在降雨工況下，模型水平位移擴大到了0.45 m，依舊發(fā)生在滑體表面上部的黏土層；滑體表面下半部黏土層的部分區(qū)域位移量有所擴大，達到了0.2 m；降雨主要影響了滑體表面黏土層的位移。

3.1.2 剪應(yīng)力

圖9、圖10為自重及降雨工況下的剪應(yīng)力云圖。從圖中可以得出，自重與降雨工況下的剪應(yīng)力云圖非常類似，均在滑動面中部產(chǎn)生應(yīng)力增大現(xiàn)象，最大達到了2.5×105Pa，大于滑動面其他部位的5×104Pa。此處的應(yīng)力分布也比較集中，這種集中現(xiàn)象容易造成局部剪切破壞，降低坡體的穩(wěn)定性?？梢哉f，該部位受到剪應(yīng)力最大，且剪應(yīng)力較為集中，容易產(chǎn)生滑動現(xiàn)象。

圖7 自重工況下水平位移云圖(單位：m)

圖8 降雨工況下水平位移云圖(單位：m)

圖9 自重工況下剪應(yīng)力云圖(單位：Pa)

圖10 降雨工況下剪應(yīng)力云圖(單位：Pa)

3.2 地震工況下的結(jié)果分析

3.2.1 加速度響應(yīng)

圖11 監(jiān)測點地震響應(yīng)加速度時程曲線

在數(shù)值模型中，加速度監(jiān)測點分別位于滑體表面、滑體中部及滑動面(具體位置見圖1)，三點坐標依次為(115，179)，(115，164)，(115，149)。將監(jiān)測點的最大響應(yīng)加速度與輸入的最大加速度相比，得到該點的加速度放大系數(shù)。監(jiān)測點的加速度時程曲線如圖11所示。從圖中可以看出，滑體表面的加速度幅值最大，表明該處受擾動程度最激烈?；w表面、滑體中部和滑動面最大加速度分別為0.91g，0.22g和0.69g。同輸入的加速度相比，最大值分別放大了2.87，0.64和1.97倍，滑體表面及滑動面最大加速度均出現(xiàn)放大效應(yīng)，而滑體中部最大加速度卻有所減小。結(jié)合位移時程曲線(圖12)，滑體中部加速度雖然較小，但其位移量始終處于滑體表面及滑動面之間。

圖12為三個監(jiān)測點水平位移時程曲線。從圖中可以看出，三點的位移趨勢一致，隨著地震作用的施加，位移量呈波動性上升，在20 s左右達到位移最大值，隨后位移值逐漸波動減小?；w表面的位移波動十分劇烈，這應(yīng)該與滑體表面為粉質(zhì)黏土有關(guān)。從位移量上看，滑體表面黏土層的位移量始終最大，其次是滑體中部，滑動面的位移量始終最小。由圖12知，在地震結(jié)束時刻，位移隨時間仍不斷增加，呈現(xiàn)一定的增加趨勢，說明滑坡在地震作用下發(fā)生了破壞[15]。

圖13 地震結(jié)束后水平位移云圖(單位：m)

3.2.2 水平位移響應(yīng)

圖13為地震結(jié)束后的水平位移云圖。從圖中可以看出，滑坡的最大位移達到了0.75 m，比自重工況下的位移增大了0.5 m，比降雨工況下增大了0.3 m，發(fā)生位置依舊在滑體表面上部的黏土層，滑體表面下半部的黏土層位移也達到了0.50～0.75 m，滑體其他位置水平位移處于0.20～0.50 m之間，與滑床之間存在相對位移。

3.2.3 剪應(yīng)力、最大剪應(yīng)變響應(yīng)

圖14、圖15分別為地震結(jié)束后剪應(yīng)力和最大剪應(yīng)變增量云圖。最大剪應(yīng)變增量可反映巖土體的剪切變形程度，其分布范圍及貫通情況可直觀反映滑坡潛在破壞滑動面的分布位置和發(fā)展趨勢[16]。從圖14可以看出，地震工況下的模型未出現(xiàn)剪應(yīng)力增大或應(yīng)力集中現(xiàn)象。從圖15中可以看出，最大剪應(yīng)變增量最大值發(fā)生在滑動面中部，達到0.175，且最大剪應(yīng)變增量在滑動面位置呈貫通現(xiàn)象，表明滑體容易產(chǎn)生滑動現(xiàn)象。

圖14 地震結(jié)束后剪應(yīng)力云圖(單位：Pa)

圖15 地震結(jié)束后最大剪應(yīng)變增量云圖

4 結(jié)論

1) 數(shù)值模擬表明：自重和降雨工況下安全系數(shù)分別為1.35和1.00，表明降雨降低了滑坡的穩(wěn)定性。在這兩種工況下，滑動面存在剪應(yīng)力集中現(xiàn)象，容易導(dǎo)致局部產(chǎn)生剪切破壞。最大位移產(chǎn)生在滑體表面上部的黏土層，與自重工況相比，降雨使黏土層位移量明顯增大。

2) 地震工況計算表明：地震作用下滑體不同部位動力響應(yīng)差別明顯。從位移上看，滑體表面水平位移量最大。從加速度動力響應(yīng)上看，滑體表面、中部及滑動面加速度時程變化趨勢一致，滑體表面加速度放大系數(shù)最大。模型并沒有出現(xiàn)應(yīng)力集中現(xiàn)象，但在模型滑動面處最大剪應(yīng)變增量呈貫通現(xiàn)象，表明在地震工況下，滑體易產(chǎn)生滑動現(xiàn)象。

3) 綜合自重、降雨及地震工況下的模擬情況，剪應(yīng)力及最大剪應(yīng)變增量最大值均發(fā)生在滑動面中部位置，對滑坡進行防治時，應(yīng)注意該位置的位移、應(yīng)力及應(yīng)變情況。針對滑體上部黏土層位移量較大的情況，建議對其進行削坡處理。