賈唯龍 常晁瑜 李佩茹張志偉 徐久歡 楊濟(jì)源

1）中國河北三河 065201 防災(zāi)科技學(xué)院

2）中國河北三河 065201 中國地震局建筑物破壞機(jī)理與防御重點(diǎn)實(shí)驗(yàn)室

3）中國哈爾濱 150080 中國地震局工程力學(xué)研究所

4）中國哈爾濱 150080 中國地震局地震工程與工程振動重點(diǎn)研究室

引言

受黃土動力易損性和南北地震帶的雙重影響，黃土地區(qū)發(fā)生地震滑坡的風(fēng)險(xiǎn)很高（徐張建等，2007）.歷史研究表明，黃土地區(qū)一旦發(fā)生中強(qiáng)地震，往往誘發(fā)數(shù)量眾多的滑坡，這些滑坡單體規(guī)模大、影響范圍廣、沖擊速度快，對人民生命財(cái)產(chǎn)安全和城市建設(shè)的危害程度極高（常晁瑜等，2020）.實(shí)踐證明，只有掌握地震滑坡發(fā)生機(jī)理，了解滑坡發(fā)生的過程和運(yùn)動特征，才能采取有效的措施，避免或減輕地震滑坡的災(zāi)害（Zhouet al，2002）.

隨著數(shù)值模擬方法的發(fā)展和計(jì)算機(jī)運(yùn)算能力的提升，基于有限單元法、有限差分法和塊體離散單元法等數(shù)值模擬方法，地震作用下斜坡的安全系數(shù)計(jì)算問題取得了極大進(jìn)步（張江偉，2016；薄景山等，2019）.然而，強(qiáng)震作用下的滑坡是一個(gè)存在滑動、平移、轉(zhuǎn)動的復(fù)雜過程，具有宏觀上的不連續(xù)性和單個(gè)塊體運(yùn)動的隨機(jī)性，當(dāng)涉及到巖土體的大變形及巖土破壞問題時(shí)，上述數(shù)值模擬方法均具有一定的局限.而離散元顆粒流方法不受變形量限制（Hadjigeorgiouet al，2009；Tanget al，2009），可方便地處理非連續(xù)介質(zhì)力學(xué)問題，有效地模擬巖土體的開裂、分離等非連續(xù)現(xiàn)象，對于探索地震滑坡的反映機(jī)理、運(yùn)動過程和致災(zāi)范圍具有一定優(yōu)勢（周健等，2000；曹文等，2017；石崇等，2018）.許多學(xué)者采用顆粒流離散元程序?qū)吰碌幕破茐臋C(jī)理開展了一系列的研究工作，例如：賀續(xù)文等（2011）討論了節(jié)理連通率對邊坡破壞形式的影響；周喻等（2016）從細(xì)觀力學(xué)角度深入研究了順層斷續(xù)節(jié)理巖質(zhì)邊坡模型破壞過程的力學(xué)機(jī)制；李新坡和何思明（2010）對節(jié)理巖質(zhì)邊坡的破壞和運(yùn)動過程進(jìn)行了研究，分析了不同關(guān)鍵參數(shù)對破壞后的堆積形狀和運(yùn)動距離的影響；Scholtès和Donzé （2012）研究了含裂隙巖體的漸進(jìn)性破壞.

通過以上分析，靜力作用下滑坡運(yùn)動失穩(wěn)過程的研究進(jìn)展較快，但針對地震滑坡失穩(wěn)破壞運(yùn)動過程的研究依然不成體系.鑒于PFC在模擬邊坡失穩(wěn)破壞方面的優(yōu)勢，本文擬在野外調(diào)查和室內(nèi)試驗(yàn)的基礎(chǔ)上，通過標(biāo)定土體細(xì)觀參數(shù)、模型建立、動力輸入等過程，利用PFC2D程序模擬西吉縣興平鄉(xiāng)堡灣村下馬達(dá)子滑坡動力荷載下的失穩(wěn)破壞運(yùn)動過程，得到該地震滑坡的破壞運(yùn)動機(jī)理，以期為黃土地區(qū)滑坡防治提供一定參考.

1 顆粒流方法基本理論

1971年，倫敦大學(xué)帝國學(xué)院Cundall博士在分析準(zhǔn)靜力或動力條件下巖石邊坡的運(yùn)動時(shí)，借鑒分子動力學(xué)理論方法，首次提出了離散單元法（discrete element method，縮寫為DEM）的概念（Cundall，Strack，1979）.PFC （particle flow code）程序又稱顆粒流方法（王光謙，倪晉仁，1992），是基于離散單元模型框架，由計(jì)算機(jī)引擎和圖形用戶界面構(gòu)成的細(xì)觀分析軟件，主要用于模擬有限尺寸顆粒的運(yùn)動和相互作用.

PFC的計(jì)算循環(huán)使用時(shí)間-步長顯示算法，需要對每個(gè)顆粒反復(fù)應(yīng)用牛頓第二定律，并對每個(gè)顆粒間的接觸單元反復(fù)應(yīng)用力-位移定律.在接觸模型中，顆粒接觸點(diǎn)的接觸力和相對位移可以分解為沿法線和切線方向的兩個(gè)分量，通過力-位移定律可以由法向剛度和切向剛度把接觸力的法向、切向分量分別與法向、切向相對位移聯(lián)系起來.第i個(gè)接觸點(diǎn)的力Fi可以沿法向和切向進(jìn)行分解，分解為法向力Fn和剪切力Fs，即

式中，ni，ti分別為接觸平面的單位矢量.

兩個(gè)接觸實(shí)體之間的法向“重疊”量表示為Un，則法向力Fn可以表示為

式中，Kn為 接觸點(diǎn)的法向剛度，n為接觸面的單位向量.

對于集合體來說，顆粒受到的剪切力與顆粒的相對運(yùn)動、加載歷史以及應(yīng)力路徑相關(guān).因此，在顆粒流模型中，剪切力的計(jì)算以增量表示.切向力增量為

式中，ks為接觸剪切剛度，μ為 顆粒摩擦系數(shù)， ΔUs表示在計(jì)算時(shí)步內(nèi)接觸位移的切向增量，由下式得到:

式中， Δt為時(shí)間步長增量，vs為速度.

2 參數(shù)標(biāo)定

利用顆粒流進(jìn)行數(shù)值模擬，本質(zhì)上就是從其細(xì)觀力學(xué)特征出發(fā)，將材料的力學(xué)響應(yīng)問題從物理域映射到數(shù)學(xué)域內(nèi)進(jìn)行數(shù)值求解（Wanget al，2003）.土的結(jié)構(gòu)特征除土顆粒的大小、形狀、表面特性及粒度級配特征外，還包括顆粒間的排列與集合關(guān)系，孔隙的大小，顆粒間聯(lián)結(jié)的特點(diǎn)等，這些土體細(xì)觀結(jié)構(gòu)的變化很大程度上影響到土體宏觀力學(xué)特性的變化（周健，池永，2003；H?rtl，Ooi，2008；陳達(dá)等，2018）.利用PFC2D進(jìn)行數(shù)值模擬之前，需要先假定土體顆粒間的本構(gòu)特性，再用顆粒流方法模擬雙軸壓縮或直剪等土工試驗(yàn)過程（Park，Song，2009），對其力學(xué)性質(zhì)進(jìn)行數(shù)值試驗(yàn)，以此得到細(xì)觀力學(xué)參數(shù)與宏觀力學(xué)參數(shù)之間的對應(yīng)關(guān)系.

通過試驗(yàn)，得到黃土的土工試驗(yàn)參數(shù)列于表1.通過細(xì)觀直剪數(shù)值試驗(yàn)?zāi)M宏觀大型直接剪切試驗(yàn)，首先賦予土體顆粒微觀力學(xué)參數(shù)（表2），剪切盒模型尺寸為500 m×250 m，生成顆粒總數(shù)為1 706個(gè)（圖1）.顆粒間的接觸模型選用平行黏結(jié)模型，在表2的細(xì)觀力學(xué)參數(shù)條件下，對直剪模型分別施加100，200，300和400 kPa的豎向應(yīng)力，得到剪切應(yīng)力-位移關(guān)系曲線，然后擬合出土體顆粒抗剪強(qiáng)度曲線（圖2），從而得到數(shù)值模擬所用巖土體顆粒的黏聚力c為19.25 kPa，內(nèi)摩擦角φ為13.5°.經(jīng)過直接剪切數(shù)值試驗(yàn)后，得到的土體顆粒內(nèi)摩擦角和黏聚力大小與表1中宏觀參數(shù)基本對應(yīng)，故此土體顆粒細(xì)觀參數(shù)可用于數(shù)值模擬.

表1 土工試驗(yàn)參數(shù)Table 1 Geotechnical test parameters

表2 細(xì)觀試驗(yàn)參數(shù)Table 2 Micro-scale test parameters

圖1 直剪數(shù)值試驗(yàn)?zāi)Ｐ蛨DFig. 1 Model diagram of direct shear numerical test

圖2 土體顆粒抗剪強(qiáng)度曲線Fig. 2 Shear strength curve of soil particles

3 模型建立

3.1 滑坡基本情況

1920年海原MS8.5特大地震誘發(fā)了數(shù)以百計(jì)的黃土滑坡，造成了慘重的人員傷亡和經(jīng)濟(jì)損失，由于特殊的地理位置和氣候條件，盡管經(jīng)歷了近百年時(shí)間，這些黃土地震滑坡的外部形態(tài)多數(shù)依然保存完好，具有極高的學(xué)術(shù)研究價(jià)值（許沖等，2018；常晁瑜等，2019）.

本文作者參與了對海原特大地震誘發(fā)滑坡的再調(diào)查，實(shí)地調(diào)查1920年海原特大地震誘發(fā)黃土滑坡620個(gè)，衛(wèi)星影像補(bǔ)充調(diào)查滑坡605個(gè)，總計(jì)1 225個(gè).調(diào)查內(nèi)容包括地理位置、滑動類型、地層巖性、坡形、斜坡結(jié)構(gòu)類型、滑坡基本特征以及影響因素等.調(diào)查所得滑坡分布如圖3所示.滑坡集中分布在三個(gè)區(qū)域：① 寧夏西吉與甘肅靜寧、會寧三縣交會地帶，區(qū)內(nèi)滑坡分布密集，發(fā)育特征明顯，代表滑坡為西吉縣興平鄉(xiāng)堡灣村下馬達(dá)子滑坡；② 寧夏海原九彩鄉(xiāng)境內(nèi)，區(qū)內(nèi)滑坡主要沿深切河流兩側(cè)發(fā)育，代表滑坡為海原縣李俊滑坡；③ 固原市區(qū)東部河川鄉(xiāng)山地區(qū)域，區(qū)內(nèi)滑坡分布較為零散，代表滑坡為固原市石碑塬滑坡.

圖3 海原大地震誘發(fā)黃土滑坡分布圖Fig. 3 Distribution map of loess landslide induced by Haiyuan landslide earthquake

下馬達(dá)子滑坡野外編號為XJ06007，位于西吉縣興平鄉(xiāng)堡灣村喜家灣組葫蘆河支流濫泥河左岸黃土斜坡上，地理位置為（105.621°E，35.893°N），是1920年海原地震誘發(fā)的大型黃土地震滑坡.下馬達(dá)子滑坡整體形狀呈長舌狀，滑坡周界清晰，主滑方向?yàn)?10°，滑坡體上陡下緩，上部平均坡度約15°，下部平均坡度小于10°，長度達(dá)到1 280多米.由于上部土體下滑，導(dǎo)致上部較薄，厚度只有8—13 m，而下部滑體較厚，約15—25 m；后緣較陡，坡度約20°，該段長約190 m，滑體上發(fā)育大量落水洞.中部向下到居民區(qū)長約640 m，該段坡體平緩，滑體已被改造為梯田.下部滑坡堆積區(qū)長約450 m.滑坡上部寬約 140 m，中部較寬，最寬約270 m.滑坡總體積約360萬 m3.滑坡區(qū)相對高差為220 m （圖4）.

圖4 下馬達(dá)子滑坡現(xiàn)場調(diào)查圖Fig. 4 Field investigation photo of Xiamadazi landslide

根據(jù)野外勘探（圖5），下馬達(dá)子滑坡地層巖性主要為上更新統(tǒng)馬蘭黃土、中更新統(tǒng)離石黃土和新近系磚紅色泥巖，滑坡發(fā)育在黃土與泥巖的不整合接觸上.由于黃土垂直節(jié)理極其發(fā)育，且其結(jié)構(gòu)松散，當(dāng)土體含水量大時(shí)，其抗拉、抗剪強(qiáng)度大大降低，其它條件具備時(shí)極其易發(fā)生變形.而且此處坡體為陰坡，受北部西伯利亞冷空氣影響，黃土堆積厚度大，這也是影響坡體失穩(wěn)的一個(gè)重要因素.

圖5 下馬達(dá)子滑坡剖面圖Fig. 5 Profile of Xiamadazi landslide

3.2 模型建立

基于下馬達(dá)子滑坡坡面圖（圖5），根據(jù)等體積法進(jìn)行地形復(fù)原，所得滑坡發(fā)生前斜坡剖面圖示于圖6，而后通過顆粒流軟件PFC2D進(jìn)行建模，對泥巖顆粒施加重力后，顆粒做自由落體運(yùn)動，循環(huán)一定次數(shù)后形成矩形，而后進(jìn)行削坡操作，利用二次函數(shù)將泥巖層削出凹面，再對斜坡開挖形成基本泥巖層形狀（圖7），此時(shí)模型長約1 040 m，高度約200 m.對顆粒間不平衡力進(jìn)行設(shè)置，多次運(yùn)算后滑床顆粒靜止平衡，此時(shí)鎖定泥巖層顆粒速度及位移.接著采用落球法生成黃土顆粒，采用平行黏結(jié)模型，分三次下落到泥巖層指定位置，每次下落后黃土顆粒進(jìn)行一次自重平衡，通過不斷運(yùn)算調(diào)整顆粒間不平衡力大小，直至不平衡力下降至允許范圍，多次試算平衡后生成滑坡發(fā)生前基本模型，如圖7所示.

圖6 滑坡發(fā)生前斜坡剖面圖Fig. 6 Profile before landslide occurrence

圖7 滑坡開始前模型Fig. 7 Model before the beginning of landslide

3.3 動力輸入

采用擬靜力法對地震進(jìn)行輸入，擬靜力法（劉紅帥等，2007；Zhouet al，2015）是一種用靜力學(xué)方法近似解決動力學(xué)問題的簡易方法，它發(fā)展較早，迄今仍然被廣泛使用.其基本思想是在靜力計(jì)算的基礎(chǔ)上，將地震作用簡化為一個(gè)慣性力系附加在研究對象上，其核心是設(shè)計(jì)地震加速度的確定.由現(xiàn)場勘察及后期地震災(zāi)害分析報(bào)告可知，該滑坡所在地在1920年海原地震中的烈度為Ⅸ，故此次模擬采用Ⅸ度區(qū)的設(shè)計(jì)加速度即0.4g施加于模型上，使顆粒在水平方向上產(chǎn)生相應(yīng)的慣性力.

4 結(jié)果分析

為了獲得在外力輸入下滑坡的運(yùn)動情況，在關(guān)鍵位置設(shè)置6個(gè)測點(diǎn)（圖8），以便于觀察模擬中滑坡的發(fā)生機(jī)理和運(yùn)動特征.

圖8 測點(diǎn)分布示意圖Fig. 8 Abridged general view of measuring points distribution

4.1 滑坡過程

滑坡模型運(yùn)行60 s后停止，為了解滑坡的發(fā)生過程，截取不同時(shí)間下滑坡的位移圖，示于圖9.從圖中可以清晰地看出：外力輸入后，模型開始失穩(wěn)，在20 s時(shí)，邊坡坡度發(fā)生改變，上部滑體在地震波作用下產(chǎn)生明顯位移，并擠壓坡面顆粒，推動坡面顆粒向前運(yùn)動，滑坡前緣顆粒在地震波作用下逐漸脫離滑床，開始運(yùn)動；30 s時(shí)，坡面顆粒在上部滑體的重力和水平力的雙重作用下，抗滑能力進(jìn)一步降低，此時(shí)已可辨別出滑動面，滑面顆粒沿滑面不斷向下運(yùn)動，同時(shí)牽引后緣顆粒下滑；40 s時(shí)，隨著滑坡后緣顆粒的逐漸下滑，滑坡后緣已經(jīng)破壞，后緣顆粒已全部沿滑動面下滑，滑面顆粒在后緣顆粒的推移和前緣顆粒的牽引的共同作用下，速度進(jìn)一步增大，位移增長迅速；50 s時(shí)，前緣顆粒在地面摩擦力作用下逐漸靜止，在坡角處產(chǎn)生堆積，前緣顆粒位移達(dá)到40 m左右，此時(shí)坡面抗滑力逐漸增大，滑體與滑床接連處顆粒速度減小，逐漸趨于靜止；60 s時(shí)，滑面顆粒運(yùn)動緩慢且處于減速階段，逐漸趨于穩(wěn)定，此時(shí)滑坡后緣顆粒在滑面位移達(dá)到最大值100 m左右，邊坡再次趨于穩(wěn)定.

圖9 不同時(shí)間-步長滑坡位移圖（顆粒數(shù)22 733）Fig. 9 Landslide displacement at different time steps （The number of particles is 22 733）

整體上，滑坡的發(fā)生過程中，滑坡后緣顆粒的推擠作用和前緣顆粒的牽引作用均發(fā)揮了重要作用，滑體在推擠和牽引的雙重作用下，逐漸發(fā)生失穩(wěn)破壞，形成貫通滑動面，與滑床間摩擦力降低，直至發(fā)生滑動破壞.

4.2 速度變化

為了監(jiān)測滑床內(nèi)部與滑體的速度特征，檢測了6個(gè)測點(diǎn)顆粒的速度變化，得到速度曲線示于圖10.由圖10可以看出：在整個(gè)滑坡過程中，滑床內(nèi)部顆粒產(chǎn)生的速度變化較小，趨近于零；而滑坡后緣測點(diǎn)（測點(diǎn)1，2，3）速度變化較為劇烈.說明在地震波輸入后滑坡后緣運(yùn)動劇烈程度大于滑坡前緣（測點(diǎn)5），這與滑坡后緣失穩(wěn)后重力勢能轉(zhuǎn)化為動能有關(guān)，測點(diǎn)3處速度高于其它五個(gè)測點(diǎn)，表明其首先被破壞并持續(xù)滑動.

圖10 測點(diǎn)速度曲線Fig. 10 Velocity curve of measuring points

4.3 位移變化

通過對比滑坡發(fā)生前后的變化情況（圖11），可以看出模型中顆粒發(fā)生了較大的位移（Behbahaniet al，2013）.為詳細(xì)研究不同測點(diǎn)位移情況，通過對五個(gè)測點(diǎn)（測點(diǎn)1—5）的水平位移和垂直位移進(jìn)行監(jiān)測，得到如圖12所示的位移曲線.通過對比可知，測點(diǎn)位置不同，水平和垂直方向的位移大小關(guān)系也不同，但大體上呈現(xiàn)出水平位移約為垂直位移的1—2倍.測點(diǎn)3處顆粒在水平和豎直方向上位移均最大，可見此處顆粒運(yùn)動距離最大.測點(diǎn)5在運(yùn)行50 s左右達(dá)到垂直位移最大值，但與其它測點(diǎn)不同的是，該測點(diǎn)處顆粒垂直位移未進(jìn)一步明顯發(fā)展，即該處顆粒在豎直方向不再發(fā)生明顯運(yùn)動，這是由于滑體在坡腳的堆積擠壓作用，在豎直方向上，使得該測點(diǎn)顆粒克服重力，達(dá)到平衡狀態(tài)，不再產(chǎn)生加速度，從而出現(xiàn)位移無明顯變化的現(xiàn)象.在水平方向上，通過滑坡過程位移圖可知，滑坡后緣處（測點(diǎn)1，2，3）水平位移基本較滑坡前緣處（測點(diǎn)5）大，這是由于滑坡前緣顆粒在破壞后期脫離滑面后產(chǎn)生堆積，受到地面摩擦力作用趨于穩(wěn)定；而測點(diǎn)1處顆粒位移較小則是因?yàn)榇颂庮w?？拷潞缶墸潞缶夘w粒經(jīng)滑面滑下后產(chǎn)生堆積再次趨于穩(wěn)定，滑體與滑面間抗滑力逐漸恢復(fù)，此處顆粒未運(yùn)動至坡腳.

圖11 滑坡前后對比（虛線部分為滑前模型）Fig. 11 Comparison before and after landslide （The dotted line represnets the pre-slip model）

圖12 測點(diǎn)x方向（a）和y方向（b）的位移圖Fig. 12 Displacement diagram of measuring points in x （a）and y （b）direction

4.4 滑動面情況

結(jié)合運(yùn)動終止時(shí)刻的顆粒位移圖像，可以確定下馬達(dá)子滑坡的滑動面如圖13所示，這與現(xiàn)場勘察結(jié)果顯示滑坡體上陡下緩（圖5）基本吻合，表明顆粒流方法可以用于地震滑坡滑動面的確定.

圖13 滑坡運(yùn)動位移圖Fig. 13 Displacement diagram of landslide movement

5 討論與結(jié)論

本文通過標(biāo)定土體細(xì)觀參數(shù)、模型建立、動力輸入等過程，利用PFC2D程序模擬了西吉縣興平鄉(xiāng)堡灣村下馬達(dá)子滑坡的失穩(wěn)破壞運(yùn)動過程，得到如下結(jié)論：

1）下馬達(dá)子滑坡的失穩(wěn)機(jī)制是地震作用下的斜坡前緣牽引、后緣推擠、坡肩受拉破壞，由于坡肩高度高且具有臨空面，失穩(wěn)后坡肩位置速度和位移均較大，是地震滑坡破壞力強(qiáng)、致災(zāi)范圍大的主要原因；

2）由于黃土特有的垂直節(jié)理，黃土高原區(qū)常常存在高陡的斜坡，且重力滑坡通常具有高陡的后壁，而地震作用下滑坡的后壁相對平緩，這是地震滑坡區(qū)別于重力滑坡的重要特征之一；

3）通過顆粒流模擬得到滑坡前后相對高差為250 m左右，滑坡后總長度達(dá)到1 080 m，滑坡后坡角約13.0°；與實(shí)際情況相對高差220 m，滑坡總長度1 040 m，滑坡后坡角約11.94°相比較為符合，顆粒流方法可以用于地震滑坡滑距的預(yù)測.

值得指出，擬靜力法是將地震期間最大慣性力施加在土體上，認(rèn)為土體中各點(diǎn)的最大加速度同時(shí)出現(xiàn)，且沒有考慮土體隨時(shí)間變化的非線性，定量分析還不夠精準(zhǔn).今后的研究中，需要利用動力時(shí)程分析法考慮地震動的時(shí)間過程對滑坡失穩(wěn)破壞過程的影響，計(jì)算滑坡的演化過程.