(長(zhǎng)江工程職業(yè)技術(shù)學(xué)院, 武漢 430212)

1 引 言

我國(guó)是一個(gè)地質(zhì)災(zāi)害頻發(fā)、災(zāi)害分布廣泛、災(zāi)害類型復(fù)雜的多山國(guó)家，每年因地質(zhì)災(zāi)害導(dǎo)致的人員傷亡數(shù)量眾多，經(jīng)濟(jì)損失巨大、社會(huì)影響廣泛。其中庫(kù)區(qū)滑坡災(zāi)害是地質(zhì)災(zāi)害中分布最廣、危害最嚴(yán)重的一種，因此，一直以來(lái)都是國(guó)內(nèi)外研究人員關(guān)注的熱點(diǎn)問(wèn)題。最著名的庫(kù)區(qū)滑坡災(zāi)害為1963年意大利的vajont水庫(kù)滑坡。該滑坡入水速度高，激起的波浪以150 m的高度越過(guò)壩頂，造成下游近3 000余人死亡。在我國(guó)三峽庫(kù)區(qū)，2003年7月13日，長(zhǎng)江支流青干河左岸的秭歸縣千將坪村發(fā)生了巨型基巖滑坡，滑坡災(zāi)害造成了房屋倒塌、青干河斷流，傷亡達(dá)到10余人，并造成巨大直接經(jīng)濟(jì)損失。這些滑坡事件引起了人們對(duì)于庫(kù)區(qū)滑坡災(zāi)害的極大關(guān)注。為了避免和減少滑坡給公眾生命財(cái)產(chǎn)帶來(lái)的巨大危害，揭示水庫(kù)蓄水誘發(fā)滑坡產(chǎn)生的動(dòng)力學(xué)機(jī)制、并開(kāi)展滑坡運(yùn)動(dòng)全過(guò)程的模擬研究，對(duì)滑坡災(zāi)害的預(yù)測(cè)預(yù)報(bào)具有重要的理論意義和應(yīng)用價(jià)值。

滑坡體的變形失穩(wěn)破壞是一個(gè)動(dòng)態(tài)的運(yùn)動(dòng)過(guò)程，包含了滑坡體內(nèi)部塊體的相互錯(cuò)動(dòng)、張開(kāi)和塊體轉(zhuǎn)動(dòng)等多種復(fù)雜運(yùn)動(dòng)方式。因此，傳統(tǒng)的剛體極限平衡計(jì)算方法以及基于連續(xù)介質(zhì)假定的有限元分析方法均無(wú)法模擬滑坡體的運(yùn)動(dòng)演化過(guò)程和運(yùn)動(dòng)過(guò)程中所反映的動(dòng)力學(xué)特點(diǎn)。作為一種模擬巖土體不連續(xù)變形的數(shù)值計(jì)算方法[1]。非連續(xù)變形分析(DDA)本身具備有完善的運(yùn)動(dòng)學(xué)理論，該方法在分析過(guò)程中，不僅可以模擬塊體本身位移和變形，而且還能模擬塊體間的相互錯(cuò)動(dòng)、張開(kāi)以及自身轉(zhuǎn)動(dòng)等運(yùn)動(dòng)形式。因此，非常適合于大位移、大變形滑坡體的動(dòng)力學(xué)分析。自不連續(xù)變形方法提出以來(lái)，便受到了國(guó)內(nèi)外的廣泛關(guān)注。多名學(xué)者發(fā)展了該方法并應(yīng)用到工程實(shí)踐中[1-6]。本文采用不連續(xù)變形方法對(duì)我國(guó)西南地區(qū)水庫(kù)內(nèi)古滑坡體滑坡的運(yùn)動(dòng)全過(guò)程進(jìn)行模擬，研究滑坡在整個(gè)動(dòng)態(tài)過(guò)程的不連續(xù)變形及滑坡最終堆積體幾何特征，揭示滑坡產(chǎn)生和發(fā)展的運(yùn)動(dòng)機(jī)制，為類似庫(kù)區(qū)滑坡的研究提供參考與借鑒。

2 DDA方法基本原理

DDA方法是一種基于最小勢(shì)能原理的解決離散塊體集合的數(shù)值計(jì)算方法?；趬K體理論，該方法中的每個(gè)塊體的變形和移動(dòng)變量是相互獨(dú)立的。塊體間的接觸計(jì)算則通過(guò)增加接觸彈簧來(lái)實(shí)現(xiàn)。

2.1 DDA控制方程

在二維DDA方法中，任意形狀塊體有6個(gè)自由度。其中，用于表征塊體運(yùn)動(dòng)和塊體變形的自由度分別有3個(gè)。因此，對(duì)于任意DDA塊體i，可表示為：

(1)

在DDA方法中，每一時(shí)步的迭代方程：

(2)

式中K為等效剛度矩陣

(3)

式中，Di、Fi(i=1,2，……，n)為6×1的子矩陣。Fi為作用于塊體i的荷載向量。

2.2 DDA接觸算法

相較于其他巖土類數(shù)值分析方法，DDA方法最大的特點(diǎn)在于處理離散塊體之間的接觸問(wèn)題。在二維DDA方法中有3種接觸類型：角-角接觸、角-邊接觸以及邊-邊接觸。其中，邊-邊接觸在計(jì)算中被轉(zhuǎn)換為兩個(gè)角-邊接觸。因此，在實(shí)際計(jì)算過(guò)程中，只有2種，即角角接觸及角邊接觸。

DDA方法采用罰函數(shù)方法計(jì)算塊體間的接觸力。即當(dāng)塊體之間有嵌入時(shí)，通過(guò)在接觸區(qū)域的法線和切向添加罰彈簧阻止塊體之間進(jìn)一步相互嵌入。塊體間接觸力的大小取決于塊體之間的相互嵌入程度。在每一計(jì)算時(shí)步中通過(guò)增減罰彈簧來(lái)實(shí)現(xiàn)無(wú)張拉、無(wú)嵌入的接觸條件。此外，DDA方法采用Mohr-Coulomb破壞準(zhǔn)則判別接觸面的狀態(tài)，即

kndn>0和ksds≤kndntanφ+cl→鎖定

(4)

kndn>0和ksds>kndntanφ+cl→滑動(dòng)

(5)

kndn≤0→張開(kāi)

(6)

其中，kn和ks分別法向和切向彈簧剛度。一般而言，kn一般取為塊體彈性模量E的10～1 000倍。ks一般取為kn的1/4。dn和ds則分別為法向和切向的嵌入距離。φ和c分別為為接觸面上的摩擦角和粘聚力。l為接觸區(qū)域的計(jì)算長(zhǎng)度。

DDA方法的接觸存在三種接觸狀態(tài)，包括：鎖定、滑動(dòng)以及張開(kāi)[6]。鎖定狀態(tài)時(shí)，平衡方程(2)中將同時(shí)添加法向和切向彈簧。當(dāng)接觸處于滑動(dòng)狀態(tài)時(shí)，添加法向彈簧以及切向摩擦力。當(dāng)接觸處于張開(kāi)狀態(tài)時(shí)，無(wú)需添加法向彈簧、切向彈簧或切向摩擦力。需要注意的是，當(dāng)接觸狀態(tài)從鎖定轉(zhuǎn)變?yōu)榛瑒?dòng)或張開(kāi)時(shí)，接觸面上的粘聚力c被移除。

3 計(jì)算模型和計(jì)算條件

3.1 計(jì)算模型

本次計(jì)算中斷面為滑坡的某個(gè)典型縱剖面，模型的長(zhǎng)度為1 050 m，高度約500 m。實(shí)際地質(zhì)資料表明，該滑坡為斜坡堆積物順向坡，因而在計(jì)算過(guò)程中引入了假想的不連續(xù)面剖分確定。在堆積體下部的基巖的完整性很好，因此在計(jì)算過(guò)程中塊體尺寸較大。根據(jù)滑坡可能的破壞模式，建立了如圖1所示的滑坡DDA計(jì)算模型，其高程范圍為EL1230～EL1460。計(jì)算域約束條件為：地面自由，左邊界和右邊界為水平方向約束，而底部邊界為垂直方向約束。

圖1 西南地區(qū)某庫(kù)岸滑坡DDA計(jì)算模型

3.2 材料參數(shù)的選取

滑坡體以及滑帶材料力學(xué)參數(shù)見(jiàn)表1。

表1 滑坡計(jì)算參數(shù)

3.3 荷載條件

通常來(lái)說(shuō)，作用于滑坡上的荷載包括集中荷載、分布荷載、自重、初始地應(yīng)力以及地下水壓力等。本文主要研究滑坡體在自重作用力和地下水壓力共同作用下的變形。計(jì)算過(guò)程中，在每一時(shí)步由位移、速度變化所產(chǎn)生的慣性力將傳遞到下一時(shí)步中。對(duì)于地下水壓力的模擬，則通過(guò)對(duì)所有在地下水位線以下的塊體表面施加裂隙水壓力來(lái)實(shí)現(xiàn)[3]。

4 計(jì)算結(jié)果及分析

4.1 滑坡運(yùn)動(dòng)過(guò)程模擬

作為一種隱式的離散塊體分析方法，不連續(xù)變形方法不僅可以模擬塊體的變形和運(yùn)動(dòng)，而且還能較好的反映出塊體之間的滑移、閉合以及張開(kāi)等非連續(xù)變形形式，從而可以進(jìn)行巖體大變形分析。此外由于不連續(xù)變形方法在提出之處就通過(guò)動(dòng)力學(xué)方法解決準(zhǔn)靜力學(xué)問(wèn)題，并在計(jì)算過(guò)程中考慮了時(shí)間因素，因此它也可以用于分析滑坡這種動(dòng)力學(xué)響應(yīng)問(wèn)題。

本文采用不連續(xù)變形方法模擬了西南地區(qū)某庫(kù)岸滑坡變形演化的全過(guò)程，圖2～圖5給出了在水庫(kù)蓄水和強(qiáng)降雨聯(lián)合作用下，滑坡體從開(kāi)裂、錯(cuò)動(dòng)到產(chǎn)生整體性滑動(dòng)的動(dòng)力學(xué)過(guò)程。

圖2 500時(shí)步變形(t=2.5 s)

圖3 2 000時(shí)步變形(t=10 s)

圖4 4 000時(shí)步變形(t=20 s)

圖5 5 000時(shí)步變形(t=25 s)

如圖2所示，在滑坡變形初期，滑坡體前緣表面、滑坡體底滑面(滑帶)以及滑坡體的后緣淺部處均發(fā)生了不同程度的拉裂以及錯(cuò)動(dòng)變形。這充分說(shuō)明在滑坡發(fā)生變形的前期階段，由于水庫(kù)的蓄水作用，滑坡體底滑面發(fā)生軟化。同時(shí)，在地下水的滲透作用下，滑體的前緣部位出現(xiàn)首先出現(xiàn)了局部破壞現(xiàn)象。隨著滑坡體變形進(jìn)一步增加，在第2 000時(shí)步時(shí)(如圖3所示)，滑坡體的后緣部位開(kāi)始下滑，并逐漸同基巖面分離。從圖中可以發(fā)現(xiàn)，此時(shí)滑坡體內(nèi)部各塊體之間發(fā)生了較大的錯(cuò)動(dòng)變形，并表現(xiàn)出張開(kāi)、滑動(dòng)以及傾倒變形特征。在第4 000時(shí)步時(shí)(如圖4所示)，滑坡體已呈現(xiàn)出整體滑動(dòng)的趨勢(shì)，并沖到河對(duì)岸。隨后，滑坡下滑速度逐漸減小并最終形成滑坡堆積體(如圖5所示)，整個(gè)滑坡運(yùn)動(dòng)歷時(shí)25s，滑坡體前緣塊體平均滑速為13.4 m/s。

4.2 位移時(shí)程曲線分析

為分析滑坡體在滑動(dòng)過(guò)程中的運(yùn)動(dòng)及變形特征，在計(jì)算模型中選取了3個(gè)典型塊體進(jìn)行分析，分別為445#塊體(坡腳前緣)、321#塊體(坡體中前部)以及95#塊體(滑坡體的后緣)。各典型塊體的位移時(shí)程曲線如圖6～7所示。

圖6 水平方向位移時(shí)程曲線

圖7 豎直方向位移時(shí)程曲線

由圖6～圖7可知，在滑坡的啟動(dòng)階段(0～500時(shí)步)，滑坡體前緣部位的變形較大，如第445#塊體的位移達(dá)到5.9 m。然而，位于后緣部位的變形則較小，如第95號(hào)塊體的位移為0.5 m。

這說(shuō)明庫(kù)區(qū)滑坡首先在坡腳處發(fā)生局部破壞，并在滑坡體前緣表層發(fā)生拉裂運(yùn)動(dòng)，然后牽引后緣滑體發(fā)生滑動(dòng)。由圖7中各典型塊體的豎直方向的位移時(shí)程曲線可知，在滑坡啟動(dòng)初期，塊體的豎直向位移為負(fù)，說(shuō)明滑坡體是向下滑動(dòng)的。當(dāng)滑坡變形穩(wěn)定后，位于滑坡體前緣的最大水平方向變形達(dá)到200 m，位移滑坡體后緣的最大水平變形達(dá)到128 m。

4.3 速度時(shí)程曲線分析

圖8和圖9分別為滑坡體中典型塊體(445#、321#和95#塊體)的水平及豎直方向的速度時(shí)程曲線。

圖8 水平方向速度時(shí)程曲線

圖9 豎直方向速度時(shí)程曲線

由圖8、圖9可知，滑坡體內(nèi)部各塊體在滑坡運(yùn)動(dòng)過(guò)程中的速度變化較大。此外，計(jì)算結(jié)果顯示，在滑坡體啟動(dòng)初期，其前緣部位的代表性塊體速度要較后緣部位塊體的速度大?；麦w整體滑動(dòng)后，滑坡體后緣部位塊體的速度則較前緣部位大。由于滑坡后緣的塊體處于高高程部位，具有較高的勢(shì)能，因此該部位塊體速度隨時(shí)間變化更加的劇烈。在滑坡運(yùn)動(dòng)后期，滑坡體整體滑移速度逐漸減小，直至形成滑坡堆積體。

5 主要結(jié)論

(1)作為一種隱式的離散塊體分析方法，不連續(xù)變形方法不僅可以模擬塊體的變形和運(yùn)動(dòng)，而且還能還好的反映出塊體之間的滑移、閉合以及張開(kāi)等非連續(xù)變形形式，從而可以進(jìn)行巖體大變形分析。因此在模擬滑坡運(yùn)動(dòng)時(shí)較傳統(tǒng)的有限元分析有明顯優(yōu)勢(shì)。

(2)滑坡不連續(xù)變形分析結(jié)果顯示，滑坡穩(wěn)定后，其前緣部位的水平方向位移約為200m，其后緣部位的水平方向位移為128m?；碌奈灰茣r(shí)程曲線和速度時(shí)程曲線表明，該庫(kù)區(qū)滑坡首先在坡腳處發(fā)生局部破壞，同時(shí)滑坡體前緣表層發(fā)生拉裂運(yùn)動(dòng)，從而進(jìn)一步牽引后部滑坡體發(fā)生變形，最后在強(qiáng)降雨以及水庫(kù)蓄水共同作用下引發(fā)了西南地區(qū)某庫(kù)岸滑坡發(fā)生整體滑動(dòng)。