闞 露，張發(fā)明

(河海大學(xué) 地球科學(xué)與工程學(xué)院，南京 210098)

利用FLAC3D預(yù)測某滑坡體地震荷載作用下的穩(wěn)定性

闞 露，張發(fā)明

(河海大學(xué) 地球科學(xué)與工程學(xué)院，南京 210098)

以瀾滄江邊某滑坡工程實例為背景，選用擬靜力法對地震荷載進行模擬。利用ANSYS軟件建立數(shù)值模型，并導(dǎo)入FLAC3D軟件進行計算，預(yù)測滑坡穩(wěn)定性。對比結(jié)果表明:天然工況條件下，滑坡體塑性區(qū)零星分布，剪應(yīng)力集中帶僅位于滑坡體表面，位移場分布合理;加載地震荷載后，滑坡體塑性區(qū)上下貫通，滑動面明顯，剪應(yīng)力集中帶深入滑坡體內(nèi)部，形成完整的潛在滑動面，而且位移大幅度增加;兩種工況條件下滑坡體的破壞形式均以水平滑動為主，地震工況條件下的滑坡體會出現(xiàn)大面積失穩(wěn)。該研究為抗震方案設(shè)計提供了有益參考。

地震;滑坡;穩(wěn)定系數(shù);擬靜力法

收稿日期:2013－01－08
第一作者簡介:闞 露(1988－)，男，安徽省滁州人，碩士，研究方向:巖體結(jié)構(gòu)與工程穩(wěn)定性，E-mail:klhhuedu@126．com。

我國西南部地區(qū)溝壑縱橫、山高坡陡，是中國大陸內(nèi)強震活動頻度最高的地區(qū)。區(qū)域內(nèi)，許多重大水利水電工程正在建設(shè)中或處于論證階段。因此，研究地震對滑坡的影響不僅具有重要的理論意義，而且具有重要的實踐意義。以5.12汶川地震為例，地震誘發(fā)了1 701處滑坡，1 844處巖石崩塌，以及1 093處邊坡失穩(wěn)[1]，并對很多水庫造成了不同程度的危害。永久位移和安全系數(shù)是評價地震作用下邊坡穩(wěn)定性的兩個主要指標(biāo)。其中安全系數(shù)是國內(nèi)采用的主要指標(biāo)，計算方法有擬靜力法、動力時程分析法、動力有限元強度折減法[2]等。

擬靜力法是一種用靜力學(xué)方法近似解決動力學(xué)問題的簡易方法，具有簡便、易操作的特點，在工程實踐中得到了廣泛的應(yīng)用[3]。其基本原理是，將復(fù)雜的、不斷變化的地震荷載簡化成一個恒定的慣性力系，然后附加在研究對象上。如何確定地震的加速度是其核心問題[4]。動力分析法雖然計算精度高，但操作復(fù)雜，現(xiàn)有規(guī)范中亦沒有明確的計算結(jié)果評價標(biāo)準(zhǔn)[5]，而且計算過程中需要大量的地震相關(guān)參數(shù)，例如地震波時程曲線等，然而，大部分中小型工程中缺乏觀測資料，難以提供相關(guān)參數(shù)，甚至忽略了對地震工況的評估?；诖耍P者以瀾滄江邊某滑坡工程實例為背景，利用FLAC3D軟件，采用擬靜力法對地震荷載進行模擬，預(yù)測滑坡體的穩(wěn)定性。

1 地震荷載模擬方法與有限差分原理

1．1 地震荷載模擬方法

根據(jù)擬靜力法，利用大小和方向恒定的水平荷載來模擬地震工況下復(fù)雜多變的地震荷載。通過模擬地應(yīng)力場將水平荷載加載到計算模型的每一個單元上。根據(jù)DL 5073—2000《水工建筑物抗震設(shè)計規(guī)范》得出滑坡某質(zhì)點i的水平向地震慣性力:

Fi=0.11gGEi，

式中:GEi——集中在質(zhì)點i的重力作用標(biāo)準(zhǔn)值;

g——重力加速度。

1.2 有限差分原理

1.2.1 FLAC有限差分公式

有限差分法就是在利用數(shù)值計算方法求解偏微分方程時，用有限差分近似公式代替每一處導(dǎo)數(shù)，從而將求解偏微分方程的問題轉(zhuǎn)化為求解代數(shù)方程的問題。

由高斯散度定理有

式中:s——封閉曲面;

ni——曲面s的單位法向量;

f——標(biāo)量、向量或張量;

xi——位置向量;

ds——增量弧長;

A——表面積。

定義f在面A上的梯度平均值，

將式(2)代入式(1)后得到

對于一個三角形子單元，

式中:Δs——三角形某一條邊的邊長，求和即在三角形三條邊上加總;

＜f＞——相應(yīng)邊的平均值[6]。

1.2.2 應(yīng)力應(yīng)變計算公式

同理可以求出

由材料的本構(gòu)方程和相應(yīng)的邊界條件，就可以求得應(yīng)力增量。對于各向同性的材料，有

式中:σij——任意一點的應(yīng)力張量，用來表示該點的應(yīng)力狀態(tài);

λ、μ——拉梅常數(shù);

θ——體積應(yīng)變;

δij——系數(shù)，當(dāng)i=j時，δij為1，否則，δij為零。

這樣，通過上述各式的迭代求解，就可以得出每一迭代時步對應(yīng)的單元的應(yīng)力應(yīng)變值[7]。

2 滑坡概況及數(shù)值模型

2．1 滑坡概況

滑坡地處瀾滄江右岸，坡體主滑方向為N40°W，岸坡為陰坡，主滑方向長度為550 m，沿河長度為300 m，表面積為16.5×104m2。坡體前緣高程為750 m，后緣高程為980 m，滑坡地貌標(biāo)志明顯，平面展布呈箕形，以左右沖溝為邊界，總體坡角約為21°，自然邊坡高度285 m，岸坡植被中等發(fā)育，滑動面平均埋深約為15 m，估算體積方量為75×104m3?；麦w發(fā)育于窩拖寨逆斷層下盤，忙糯斷層上盤，下伏基巖為二疊紀(jì)中粗?；◢忛W長巖，其表層物質(zhì)以碎塊石土層為主，結(jié)構(gòu)松散(圖1)。

圖1 滑坡剖面示意Fig．1 Geological profile of landslide

2.2 幾何模型

FLAC3D軟件采用顯示差分法求解微分方程，不形成剛度矩陣，每一步計算僅需要很小的內(nèi)存。在求解過程中通過疊加每一時步的小變形獲得大變形求解，從而僅占用計算機很少的內(nèi)存即可，并且在模擬材料的塑性破壞和塑性流動方面具有明顯優(yōu)勢。而地震工況下滑坡的破壞失穩(wěn)常常伴隨著巨大的坡體變形，所以文中采用FLAC3D軟件預(yù)測滑坡在地震工況下的穩(wěn)定性。由于FLAC3D前處理功能較弱，建立復(fù)雜的三維模型非常困難，故選用ANSYS軟件進行前期建模，劃分計算網(wǎng)格，再導(dǎo)入FLAC3D程序中計算，以優(yōu)化分析過程。

模型以主滑方向的反方向為y軸正方向，豎直向上為z軸正方向，指向河流上游方向為x軸正方向。豎直方向坐標(biāo)采用實際高程坐標(biāo)，建立直角坐標(biāo)系，如圖2所示。模型分為兩個部分，共有11 518個單元，2 777個節(jié)點，其底面高程為650 m，沿平行河流方向長度為869 m，沿主滑方向長度為1 077 m，采用Solid45單元計算。計算模型邊界約束形式為，側(cè)邊界只對水平方向進行約束，底邊界在水平和豎直方向均進行約束，模型的上部邊界取為自由面。FLAC3D程序中，節(jié)點速度是主要變量，所以選取模型的邊界條件是通過約束模型邊界的節(jié)點速度實現(xiàn)的[8]，即模型底部邊界的水平、豎直方向的速度約束和四周邊界水平方向的速度約束。在計算程序中，將邊界上xvel或yvel的值設(shè)置為0，以達到約束邊界的目的。利用彈性求解法生成初始地應(yīng)力場的步驟是，將體積模量和剪切模量設(shè)置為大值并將材料的本構(gòu)模型設(shè)置為彈性模型，然后進行求解。初始地應(yīng)力場生成以后，賦予材料相關(guān)參數(shù)，設(shè)置重力，加載水平荷載，并采用摩爾－庫倫模型進行計算。

圖2 三維模型網(wǎng)格Fig．2 Three-dimensional block group model

根據(jù)實際地勘報告，計算采用的各層物理力學(xué)參數(shù)如表1所示，其中，基巖為微風(fēng)化和新鮮巖體。

表1 巖、土體物理力學(xué)指標(biāo)Table 1 Physico-mechanical parameters of rock and soil

3 結(jié)果分析

3．1 塑性區(qū)分布規(guī)律

天然工況條件下塑性區(qū)零星分布于滑坡體上游邊界的沖溝附近，如圖3a所示。塑性區(qū)分布高程(z坐標(biāo))范圍為790～830 m，平行河流方向(x坐標(biāo))范圍約為550 m。而地震工況條件下塑性區(qū)已連成大片，基本覆蓋了整個滑坡表面，滑坡體上部以張拉屈服為主，中部及下部以剪切屈服為主，如圖3b所示。通過對比得知，在地震工況條件下，滑坡體塑性區(qū)面積急劇增加。

從天然工況條件下x=560 m的塑性區(qū)剖面(圖4a)分布區(qū)域來看，滑坡體的屈服區(qū)域很小，僅在滑坡體表面零星出現(xiàn)，塑性區(qū)分布于y坐標(biāo)(主滑方向)270～480 m、z坐標(biāo)790～830 m之間，厚度為5 ～10 m。這說明滑坡體發(fā)生剪切和張拉破壞的可能性很小，即使發(fā)生，也僅是局部區(qū)域，不會出現(xiàn)整體性滑動。而從地震工況條件下同一剖面的塑性區(qū)分布圖(圖4b)來看，塑性區(qū)已經(jīng)上下貫通，距離滑坡體表面深度為10～20 m，可以明顯看出滑動面位置，很可能發(fā)生整體性滑動。

圖3 整體塑性區(qū)分布Fig．3 Plastic zones distribution on whole

圖4 x=560 m剖面塑性區(qū)Fig．4 Plastic zones distribution on X=560 m plane

3.2 剪應(yīng)力(剪應(yīng)變增量)分布規(guī)律

判斷堆積體的(潛在)滑動面(帶)，可根據(jù)其剪應(yīng)力(應(yīng)變增量)來判斷。剪應(yīng)力較為集中或剪應(yīng)變增量較大(絕對值)的部位，則為其(潛在)滑動面(帶)，變形破壞也都沿此處發(fā)生;剪應(yīng)力分布較為分散(均勻)或者剪應(yīng)變增量較小、基本上未發(fā)生變化的部位，一般不會有潛在滑動面產(chǎn)生，因此，這些部位也不會發(fā)生較大的變形和破壞。天然工況條件下剪應(yīng)力集中帶主要出現(xiàn)在滑坡體中上部、上游邊界沖溝附近，高程為880 m左右(圖5a);而地震工況條件下剪應(yīng)力集中帶的表面積約為天然工況條件下表面積的兩倍(圖5b)。

圖5 剪應(yīng)變增量立體展示Fig．5 Shear strain increment on whole

天然工況條件下x=450 m的剪應(yīng)變增量剖面見圖6a，剪應(yīng)力集中帶僅出現(xiàn)在滑坡表面以下約5 m范圍內(nèi)，坡體內(nèi)部沒有明顯的剪應(yīng)力集中。地震工況條件下同一剖面的剪應(yīng)變增量見圖6b，可以看出，潛在滑動面的位置已經(jīng)深入坡體內(nèi)部，很可能會發(fā)生整體性滑動。

3.3 位移場分布規(guī)律

天然工況條件下的滑坡體變形主要出現(xiàn)在滑坡體中部和中上部區(qū)域，如圖7、8所示。豎直方向位移總體表現(xiàn)為下沉，最大沉降區(qū)與水平方向最大位

圖6 x=450 m剪應(yīng)變增量剖面Fig．6 Shear strain increment on X=450 m plane

圖7 滑坡體豎直方向位移云圖Fig．7 Vertical displacement contour on whole

移區(qū)位置基本一致，局部沉降值最大為0.71 cm (圖7a)。出現(xiàn)局部最大沉降的區(qū)域位于滑坡體中上部、上游邊界沖溝附近，高程(z坐標(biāo))范圍為890～930 m，平行河流方向(x坐標(biāo))范圍為550～590 m，表面積小于100 m2;在滑坡體中部，高程為830～860 m出現(xiàn)隆起，最大值為0.29 cm。地震工況條件下滑坡體沉降區(qū)和隆起區(qū)的位置與天然工況條件下的基本一致，但位移量在數(shù)量級上有很大差別。局部沉降值最大為49.77 cm，在滑坡體中下部，最大隆起值為50.83 cm(圖7b)。

水平方向位移主要表現(xiàn)為傾向河流方向。天然工況的水平方向局部最大位移值為2.10 cm，出現(xiàn)區(qū)域與豎直方向局部最大沉降區(qū)一致，滑坡體其他部位的水平方向位移值大都在0～1 cm之間(圖8a)。地震工況下局部最大水平位移值為200.93 cm，位移值過大，判斷為滑坡體失穩(wěn)(圖8b)。

圖8 滑坡體水平方向位移云圖Fig．8 Horizontal displacement contour on whole

根據(jù)強度折減法得出天然工況下滑坡體穩(wěn)定系數(shù)為1.19，地震工況條件下滑坡體穩(wěn)定系數(shù)為0.88。

4 結(jié)論

文中結(jié)合工程實例，選用擬靜力分析法對一滑坡地震工況下的穩(wěn)定性進行數(shù)值分析，主要得到以下結(jié)論:

(1)天然工況條件下，滑坡體塑性區(qū)分布零星，剪應(yīng)力集中帶僅位于滑坡體表面，位移場分布合理;加載地震荷載后，滑坡體塑性區(qū)上下貫通，滑動面明顯，剪應(yīng)力集中帶深入滑坡體內(nèi)部，形成完整的潛在滑動面而且位移大幅度增加?？梢耘袛?，地震工況條件下滑坡體會出現(xiàn)大面積失穩(wěn)。

(2)天然工況和地震工況兩種條件下，滑坡體的水平位移值均約為豎直位移值的三倍，可見滑坡體的破壞形式以水平滑動為主。

[1] 徐夢珍，王兆印，漆力?。氪ǖ卣鹨l(fā)的次生災(zāi)害鏈[J]．山地學(xué)報，2012，30(4):502－512．

[2] 葉海林，黃潤秋，鄭穎芝，等．地震作用下邊坡穩(wěn)定性安全評價的研究[J]．地下空間與工程學(xué)報，2009，5(6):1248－1252．

[3] 李 亮，褚雪松，龐 峰，等．地震邊坡穩(wěn)定性分析的擬靜力方法適用性探討[J]．世界地震工程，2012，28(2):57－63．

[4] 包世華．結(jié)構(gòu)動力學(xué)[M]．武漢:武漢理工大學(xué)出版社，2004．

[5] 宮經(jīng)偉，嚴(yán)新軍，夏新利，等．?dāng)M靜力法在溢洪道閘室抗震穩(wěn)定計算中的應(yīng)用[J]．人民黃河，2010，32(1):107－108．

[6] 鄒 棟，鄭 宏．快速拉格朗日法及其在邊坡穩(wěn)定性分析中的應(yīng)用[J]．礦業(yè)研究與開發(fā)，2005，25(5):84－87．

[7] 黃潤秋，許 強．顯式拉格朗日差分分析在巖石邊坡工程中的應(yīng)用[J]．巖石力學(xué)與工程學(xué)報，1995，14(4):346－354．

[8] 陳育民，徐鼎平．FLAC/FLAC3D基礎(chǔ)與工程實例[M]．北京:中國水利水電出版社，2009．泥的鋼纖維混凝土密實度得到提高，從而使得超聲波聲速、回彈值和混凝土抗壓強度都得到提高。這與在普通混凝土中得到的結(jié)論是一致的。

(編輯 荀海鑫)

Stability prediction of landslide under seismic load using FLAC3Dsoftware

KAN Lu， ZHANG Faming
(School of Earth Sciences＆Engineering，Hohai University，Nanjing 210098，China)

This paper building on a landslide associated with a reservoir construction on Lancang River discusses the simulation of seismic load by selecting pseudo-static method，the development of numerical model using ANSYS software，the introduction of FLAC3Dsoftware for calculation，and the prediction of the landslide stability．The comparison shows that，under the natural condition，there occur the sporadic distribution of the plastic zones of the slope，the concentration of the shear stress only on the surface of the landslide，and the rational displacement of field;but under seismic load，there arise the complete connection of the plastic zones of the slope，the obvious slip plane，and the extension of the shear stress concentration area deep into the interior body of the slope，resulting in a complete potential sliding surface，marked by a significant increase in the displacement;the failure mode of the landslides under the two conditions is dominated by the horizontal slip while landslide due to the seismic condition is accompanied by large areas of instability．The study serves as a

for the anti-seismic designs．

earthquake;landslide;stability factor; pseudo-static method

10．3969/j．issn．1671－0118．2013．01．017

P642．22

1671－0118(2013)01－0077－06

A