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

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

2）中國河北三河 065201 中國地震局建筑物破壞機理與防御重點實驗室

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

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

引言

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

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

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

1 顆粒流方法基本理論

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

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

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

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

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

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

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

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

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

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

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

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

表2 細觀試驗參數(shù)Table 2 Micro-scale test parameters

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

圖2 土體顆?？辜魪姸惹€Fig. 2 Shear strength curve of soil particles

3 模型建立

3.1 滑坡基本情況

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

本文作者參與了對海原特大地震誘發(fā)滑坡的再調(diào)查，實地調(diào)查1920年海原特大地震誘發(fā)黃土滑坡620個，衛(wèi)星影像補充調(diào)查滑坡605個，總計1 225個.調(diào)查內(nèi)容包括地理位置、滑動類型、地層巖性、坡形、斜坡結(jié)構(gòu)類型、滑坡基本特征以及影響因素等.調(diào)查所得滑坡分布如圖3所示.滑坡集中分布在三個區(qū)域：① 寧夏西吉與甘肅靜寧、會寧三縣交會地帶，區(qū)內(nèi)滑坡分布密集，發(fā)育特征明顯，代表滑坡為西吉縣興平鄉(xiāng)堡灣村下馬達子滑坡；② 寧夏海原九彩鄉(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

下馬達子滑坡野外編號為XJ06007，位于西吉縣興平鄉(xiāng)堡灣村喜家灣組葫蘆河支流濫泥河左岸黃土斜坡上，地理位置為（105.621°E，35.893°N），是1920年海原地震誘發(fā)的大型黃土地震滑坡.下馬達子滑坡整體形狀呈長舌狀，滑坡周界清晰，主滑方向為310°，滑坡體上陡下緩，上部平均坡度約15°，下部平均坡度小于10°，長度達到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 下馬達子滑坡現(xiàn)場調(diào)查圖Fig. 4 Field investigation photo of Xiamadazi landslide

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

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

3.2 模型建立

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

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

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

3.3 動力輸入

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

4 結(jié)果分析

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

圖8 測點分布示意圖Fig. 8 Abridged general view of measuring points distribution

4.1 滑坡過程

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

圖9 不同時間-步長滑坡位移圖（顆粒數(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個測點顆粒的速度變化，得到速度曲線示于圖10.由圖10可以看出：在整個滑坡過程中，滑床內(nèi)部顆粒產(chǎn)生的速度變化較小，趨近于零；而滑坡后緣測點（測點1，2，3）速度變化較為劇烈.說明在地震波輸入后滑坡后緣運動劇烈程度大于滑坡前緣（測點5），這與滑坡后緣失穩(wěn)后重力勢能轉(zhuǎn)化為動能有關(guān)，測點3處速度高于其它五個測點，表明其首先被破壞并持續(xù)滑動.

圖10 測點速度曲線Fig. 10 Velocity curve of measuring points

4.3 位移變化

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

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

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

4.4 滑動面情況

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

圖13 滑坡運動位移圖Fig. 13 Displacement diagram of landslide movement

5 討論與結(jié)論

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

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

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

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

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