賈唯龍 常晁瑜 李佩茹張志偉 徐久歡 楊濟(jì)源
1)中國(guó)河北三河 065201 防災(zāi)科技學(xué)院
2)中國(guó)河北三河 065201 中國(guó)地震局建筑物破壞機(jī)理與防御重點(diǎn)實(shí)驗(yàn)室
3)中國(guó)哈爾濱 150080 中國(guó)地震局工程力學(xué)研究所
4)中國(guó)哈爾濱 150080 中國(guó)地震局地震工程與工程振動(dòng)重點(diǎn)研究室
受黃土動(dòng)力易損性和南北地震帶的雙重影響,黃土地區(qū)發(fā)生地震滑坡的風(fēng)險(xiǎn)很高(徐張建等,2007).歷史研究表明,黃土地區(qū)一旦發(fā)生中強(qiáng)地震,往往誘發(fā)數(shù)量眾多的滑坡,這些滑坡單體規(guī)模大、影響范圍廣、沖擊速度快,對(duì)人民生命財(cái)產(chǎn)安全和城市建設(shè)的危害程度極高(常晁瑜等,2020).實(shí)踐證明,只有掌握地震滑坡發(fā)生機(jī)理,了解滑坡發(fā)生的過(guò)程和運(yùn)動(dòng)特征,才能采取有效的措施,避免或減輕地震滑坡的災(zāi)害(Zhouet al,2002).
隨著數(shù)值模擬方法的發(fā)展和計(jì)算機(jī)運(yùn)算能力的提升,基于有限單元法、有限差分法和塊體離散單元法等數(shù)值模擬方法,地震作用下斜坡的安全系數(shù)計(jì)算問(wèn)題取得了極大進(jìn)步(張江偉,2016;薄景山等,2019).然而,強(qiáng)震作用下的滑坡是一個(gè)存在滑動(dòng)、平移、轉(zhuǎn)動(dòng)的復(fù)雜過(guò)程,具有宏觀上的不連續(xù)性和單個(gè)塊體運(yùn)動(dòng)的隨機(jī)性,當(dāng)涉及到巖土體的大變形及巖土破壞問(wèn)題時(shí),上述數(shù)值模擬方法均具有一定的局限.而離散元顆粒流方法不受變形量限制(Hadjigeorgiouet al,2009;Tanget al,2009),可方便地處理非連續(xù)介質(zhì)力學(xué)問(wèn)題,有效地模擬巖土體的開(kāi)裂、分離等非連續(xù)現(xiàn)象,對(duì)于探索地震滑坡的反映機(jī)理、運(yùn)動(dòng)過(guò)程和致災(zāi)范圍具有一定優(yōu)勢(shì)(周健等,2000;曹文等,2017;石崇等,2018).許多學(xué)者采用顆粒流離散元程序?qū)吰碌幕破茐臋C(jī)理開(kāi)展了一系列的研究工作,例如:賀續(xù)文等(2011)討論了節(jié)理連通率對(duì)邊坡破壞形式的影響;周喻等(2016)從細(xì)觀力學(xué)角度深入研究了順層斷續(xù)節(jié)理巖質(zhì)邊坡模型破壞過(guò)程的力學(xué)機(jī)制;李新坡和何思明(2010)對(duì)節(jié)理巖質(zhì)邊坡的破壞和運(yùn)動(dòng)過(guò)程進(jìn)行了研究,分析了不同關(guān)鍵參數(shù)對(duì)破壞后的堆積形狀和運(yùn)動(dòng)距離的影響;Scholtès和Donzé (2012)研究了含裂隙巖體的漸進(jìn)性破壞.
通過(guò)以上分析,靜力作用下滑坡運(yùn)動(dòng)失穩(wěn)過(guò)程的研究進(jìn)展較快,但針對(duì)地震滑坡失穩(wěn)破壞運(yùn)動(dòng)過(guò)程的研究依然不成體系.鑒于PFC在模擬邊坡失穩(wěn)破壞方面的優(yōu)勢(shì),本文擬在野外調(diào)查和室內(nèi)試驗(yàn)的基礎(chǔ)上,通過(guò)標(biāo)定土體細(xì)觀參數(shù)、模型建立、動(dòng)力輸入等過(guò)程,利用PFC2D程序模擬西吉縣興平鄉(xiāng)堡灣村下馬達(dá)子滑坡動(dòng)力荷載下的失穩(wěn)破壞運(yùn)動(dòng)過(guò)程,得到該地震滑坡的破壞運(yùn)動(dòng)機(jī)理,以期為黃土地區(qū)滑坡防治提供一定參考.
1971年,倫敦大學(xué)帝國(guó)學(xué)院Cundall博士在分析準(zhǔn)靜力或動(dòng)力條件下巖石邊坡的運(yùn)動(dòng)時(shí),借鑒分子動(dòng)力學(xué)理論方法,首次提出了離散單元法(discrete element method,縮寫為DEM)的概念(Cundall,Strack,1979).PFC (particle flow code)程序又稱顆粒流方法(王光謙,倪晉仁,1992),是基于離散單元模型框架,由計(jì)算機(jī)引擎和圖形用戶界面構(gòu)成的細(xì)觀分析軟件,主要用于模擬有限尺寸顆粒的運(yùn)動(dòng)和相互作用.
PFC的計(jì)算循環(huán)使用時(shí)間-步長(zhǎng)顯示算法,需要對(duì)每個(gè)顆粒反復(fù)應(yīng)用牛頓第二定律,并對(duì)每個(gè)顆粒間的接觸單元反復(fù)應(yīng)用力-位移定律.在接觸模型中,顆粒接觸點(diǎn)的接觸力和相對(duì)位移可以分解為沿法線和切線方向的兩個(gè)分量,通過(guò)力-位移定律可以由法向剛度和切向剛度把接觸力的法向、切向分量分別與法向、切向相對(duì)位移聯(lián)系起來(lái).第i個(gè)接觸點(diǎn)的力Fi可以沿法向和切向進(jìn)行分解,分解為法向力Fn和剪切力Fs,即
式中,ni,ti分別為接觸平面的單位矢量.
兩個(gè)接觸實(shí)體之間的法向“重疊”量表示為Un,則法向力Fn可以表示為
式中,Kn為 接觸點(diǎn)的法向剛度,n為接觸面的單位向量.
對(duì)于集合體來(lái)說(shuō),顆粒受到的剪切力與顆粒的相對(duì)運(yùn)動(dòng)、加載歷史以及應(yīng)力路徑相關(guān).因此,在顆粒流模型中,剪切力的計(jì)算以增量表示.切向力增量為
式中,ks為接觸剪切剛度,μ為 顆粒摩擦系數(shù), ΔUs表示在計(jì)算時(shí)步內(nèi)接觸位移的切向增量,由下式得到:
式中, Δt為時(shí)間步長(zhǎng)增量,vs為速度.
利用顆粒流進(jìn)行數(shù)值模擬,本質(zhì)上就是從其細(xì)觀力學(xué)特征出發(fā),將材料的力學(xué)響應(yīng)問(wèn)題從物理域映射到數(shù)學(xué)域內(nèi)進(jìn)行數(shù)值求解(Wanget al,2003).土的結(jié)構(gòu)特征除土顆粒的大小、形狀、表面特性及粒度級(jí)配特征外,還包括顆粒間的排列與集合關(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)過(guò)程(Park,Song,2009),對(duì)其力學(xué)性質(zhì)進(jìn)行數(shù)值試驗(yàn),以此得到細(xì)觀力學(xué)參數(shù)與宏觀力學(xué)參數(shù)之間的對(duì)應(yīng)關(guān)系.
通過(guò)試驗(yàn),得到黃土的土工試驗(yàn)參數(shù)列于表1.通過(guò)細(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ù)條件下,對(duì)直剪模型分別施加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)過(guò)直接剪切數(shù)值試驗(yàn)后,得到的土體顆粒內(nèi)摩擦角和黏聚力大小與表1中宏觀參數(shù)基本對(duì)應(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)P蛨DFig. 1 Model diagram of direct shear numerical test
圖2 土體顆??辜魪?qiáng)度曲線Fig. 2 Shear strength curve of soil particles
1920年海原MS8.5特大地震誘發(fā)了數(shù)以百計(jì)的黃土滑坡,造成了慘重的人員傷亡和經(jīng)濟(jì)損失,由于特殊的地理位置和氣候條件,盡管經(jīng)歷了近百年時(shí)間,這些黃土地震滑坡的外部形態(tài)多數(shù)依然保存完好,具有極高的學(xué)術(shù)研究?jī)r(jià)值(許沖等,2018;常晁瑜等,2019).
本文作者參與了對(duì)海原特大地震誘發(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)容包括地理位置、滑動(dòng)類型、地層巖性、坡形、斜坡結(jié)構(gòu)類型、滑坡基本特征以及影響因素等.調(diào)查所得滑坡分布如圖3所示.滑坡集中分布在三個(gè)區(qū)域:① 寧夏西吉與甘肅靜寧、會(huì)寧三縣交會(huì)地帶,區(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á)子滑坡野外編號(hào)為XJ06007,位于西吉縣興平鄉(xiāng)堡灣村喜家灣組葫蘆河支流濫泥河左岸黃土斜坡上,地理位置為(105.621°E,35.893°N),是1920年海原地震誘發(fā)的大型黃土地震滑坡.下馬達(dá)子滑坡整體形狀呈長(zhǎng)舌狀,滑坡周界清晰,主滑方向?yàn)?10°,滑坡體上陡下緩,上部平均坡度約15°,下部平均坡度小于10°,長(zhǎng)度達(dá)到1 280多米.由于上部土體下滑,導(dǎo)致上部較薄,厚度只有8—13 m,而下部滑體較厚,約15—25 m;后緣較陡,坡度約20°,該段長(zhǎng)約190 m,滑體上發(fā)育大量落水洞.中部向下到居民區(qū)長(zhǎng)約640 m,該段坡體平緩,滑體已被改造為梯田.下部滑坡堆積區(qū)長(zhǎng)約450 m.滑坡上部寬約 140 m,中部較寬,最寬約270 m.滑坡總體積約360萬(wàn) m3.滑坡區(qū)相對(duì)高差為220 m (圖4).
圖4 下馬達(dá)子滑坡現(xiàn)場(chǎng)調(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
基于下馬達(dá)子滑坡坡面圖(圖5),根據(jù)等體積法進(jìn)行地形復(fù)原,所得滑坡發(fā)生前斜坡剖面圖示于圖6,而后通過(guò)顆粒流軟件PFC2D進(jìn)行建模,對(duì)泥巖顆粒施加重力后,顆粒做自由落體運(yùn)動(dòng),循環(huán)一定次數(shù)后形成矩形,而后進(jìn)行削坡操作,利用二次函數(shù)將泥巖層削出凹面,再對(duì)斜坡開(kāi)挖形成基本泥巖層形狀(圖7),此時(shí)模型長(zhǎng)約1 040 m,高度約200 m.對(duì)顆粒間不平衡力進(jìn)行設(shè)置,多次運(yùn)算后滑床顆粒靜止平衡,此時(shí)鎖定泥巖層顆粒速度及位移.接著采用落球法生成黃土顆粒,采用平行黏結(jié)模型,分三次下落到泥巖層指定位置,每次下落后黃土顆粒進(jìn)行一次自重平衡,通過(guò)不斷運(yùn)算調(diào)整顆粒間不平衡力大小,直至不平衡力下降至允許范圍,多次試算平衡后生成滑坡發(fā)生前基本模型,如圖7所示.
圖6 滑坡發(fā)生前斜坡剖面圖Fig. 6 Profile before landslide occurrence
圖7 滑坡開(kāi)始前模型Fig. 7 Model before the beginning of landslide
采用擬靜力法對(duì)地震進(jìn)行輸入,擬靜力法(劉紅帥等,2007;Zhouet al,2015)是一種用靜力學(xué)方法近似解決動(dòng)力學(xué)問(wèn)題的簡(jiǎn)易方法,它發(fā)展較早,迄今仍然被廣泛使用.其基本思想是在靜力計(jì)算的基礎(chǔ)上,將地震作用簡(jiǎn)化為一個(gè)慣性力系附加在研究對(duì)象上,其核心是設(shè)計(jì)地震加速度的確定.由現(xiàn)場(chǎng)勘察及后期地震災(zāi)害分析報(bào)告可知,該滑坡所在地在1920年海原地震中的烈度為Ⅸ,故此次模擬采用Ⅸ度區(qū)的設(shè)計(jì)加速度即0.4g施加于模型上,使顆粒在水平方向上產(chǎn)生相應(yīng)的慣性力.
為了獲得在外力輸入下滑坡的運(yùn)動(dòng)情況,在關(guān)鍵位置設(shè)置6個(gè)測(cè)點(diǎn)(圖8),以便于觀察模擬中滑坡的發(fā)生機(jī)理和運(yùn)動(dòng)特征.
圖8 測(cè)點(diǎn)分布示意圖Fig. 8 Abridged general view of measuring points distribution
滑坡模型運(yùn)行60 s后停止,為了解滑坡的發(fā)生過(guò)程,截取不同時(shí)間下滑坡的位移圖,示于圖9.從圖中可以清晰地看出:外力輸入后,模型開(kāi)始失穩(wěn),在20 s時(shí),邊坡坡度發(fā)生改變,上部滑體在地震波作用下產(chǎn)生明顯位移,并擠壓坡面顆粒,推動(dòng)坡面顆粒向前運(yùn)動(dòng),滑坡前緣顆粒在地震波作用下逐漸脫離滑床,開(kāi)始運(yùn)動(dòng);30 s時(shí),坡面顆粒在上部滑體的重力和水平力的雙重作用下,抗滑能力進(jìn)一步降低,此時(shí)已可辨別出滑動(dòng)面,滑面顆粒沿滑面不斷向下運(yùn)動(dòng),同時(shí)牽引后緣顆粒下滑;40 s時(shí),隨著滑坡后緣顆粒的逐漸下滑,滑坡后緣已經(jīng)破壞,后緣顆粒已全部沿滑動(dòng)面下滑,滑面顆粒在后緣顆粒的推移和前緣顆粒的牽引的共同作用下,速度進(jìn)一步增大,位移增長(zhǎng)迅速;50 s時(shí),前緣顆粒在地面摩擦力作用下逐漸靜止,在坡角處產(chǎn)生堆積,前緣顆粒位移達(dá)到40 m左右,此時(shí)坡面抗滑力逐漸增大,滑體與滑床接連處顆粒速度減小,逐漸趨于靜止;60 s時(shí),滑面顆粒運(yùn)動(dòng)緩慢且處于減速階段,逐漸趨于穩(wěn)定,此時(shí)滑坡后緣顆粒在滑面位移達(dá)到最大值100 m左右,邊坡再次趨于穩(wěn)定.
圖9 不同時(shí)間-步長(zhǎng)滑坡位移圖(顆粒數(shù)22 733)Fig. 9 Landslide displacement at different time steps (The number of particles is 22 733)
整體上,滑坡的發(fā)生過(guò)程中,滑坡后緣顆粒的推擠作用和前緣顆粒的牽引作用均發(fā)揮了重要作用,滑體在推擠和牽引的雙重作用下,逐漸發(fā)生失穩(wěn)破壞,形成貫通滑動(dòng)面,與滑床間摩擦力降低,直至發(fā)生滑動(dòng)破壞.
為了監(jiān)測(cè)滑床內(nèi)部與滑體的速度特征,檢測(cè)了6個(gè)測(cè)點(diǎn)顆粒的速度變化,得到速度曲線示于圖10.由圖10可以看出:在整個(gè)滑坡過(guò)程中,滑床內(nèi)部顆粒產(chǎn)生的速度變化較小,趨近于零;而滑坡后緣測(cè)點(diǎn)(測(cè)點(diǎn)1,2,3)速度變化較為劇烈.說(shuō)明在地震波輸入后滑坡后緣運(yùn)動(dòng)劇烈程度大于滑坡前緣(測(cè)點(diǎn)5),這與滑坡后緣失穩(wěn)后重力勢(shì)能轉(zhuǎn)化為動(dòng)能有關(guān),測(cè)點(diǎn)3處速度高于其它五個(gè)測(cè)點(diǎn),表明其首先被破壞并持續(xù)滑動(dòng).
圖10 測(cè)點(diǎn)速度曲線Fig. 10 Velocity curve of measuring points
通過(guò)對(duì)比滑坡發(fā)生前后的變化情況(圖11),可以看出模型中顆粒發(fā)生了較大的位移(Behbahaniet al,2013).為詳細(xì)研究不同測(cè)點(diǎn)位移情況,通過(guò)對(duì)五個(gè)測(cè)點(diǎn)(測(cè)點(diǎn)1—5)的水平位移和垂直位移進(jìn)行監(jiān)測(cè),得到如圖12所示的位移曲線.通過(guò)對(duì)比可知,測(cè)點(diǎn)位置不同,水平和垂直方向的位移大小關(guān)系也不同,但大體上呈現(xiàn)出水平位移約為垂直位移的1—2倍.測(cè)點(diǎn)3處顆粒在水平和豎直方向上位移均最大,可見(jiàn)此處顆粒運(yùn)動(dòng)距離最大.測(cè)點(diǎn)5在運(yùn)行50 s左右達(dá)到垂直位移最大值,但與其它測(cè)點(diǎn)不同的是,該測(cè)點(diǎn)處顆粒垂直位移未進(jìn)一步明顯發(fā)展,即該處顆粒在豎直方向不再發(fā)生明顯運(yùn)動(dòng),這是由于滑體在坡腳的堆積擠壓作用,在豎直方向上,使得該測(cè)點(diǎn)顆??朔亓Γ_(dá)到平衡狀態(tài),不再產(chǎn)生加速度,從而出現(xiàn)位移無(wú)明顯變化的現(xiàn)象.在水平方向上,通過(guò)滑坡過(guò)程位移圖可知,滑坡后緣處(測(cè)點(diǎn)1,2,3)水平位移基本較滑坡前緣處(測(cè)點(diǎn)5)大,這是由于滑坡前緣顆粒在破壞后期脫離滑面后產(chǎn)生堆積,受到地面摩擦力作用趨于穩(wěn)定;而測(cè)點(diǎn)1處顆粒位移較小則是因?yàn)榇颂庮w粒靠近滑坡后緣,滑坡后緣顆粒經(jīng)滑面滑下后產(chǎn)生堆積再次趨于穩(wěn)定,滑體與滑面間抗滑力逐漸恢復(fù),此處顆粒未運(yùn)動(dòng)至坡腳.
圖11 滑坡前后對(duì)比(虛線部分為滑前模型)Fig. 11 Comparison before and after landslide (The dotted line represnets the pre-slip model)
圖12 測(cè)點(diǎn)x方向(a)和y方向(b)的位移圖Fig. 12 Displacement diagram of measuring points in x (a)and y (b)direction
結(jié)合運(yùn)動(dòng)終止時(shí)刻的顆粒位移圖像,可以確定下馬達(dá)子滑坡的滑動(dòng)面如圖13所示,這與現(xiàn)場(chǎng)勘察結(jié)果顯示滑坡體上陡下緩(圖5)基本吻合,表明顆粒流方法可以用于地震滑坡滑動(dòng)面的確定.
圖13 滑坡運(yùn)動(dòng)位移圖Fig. 13 Displacement diagram of landslide movement
本文通過(guò)標(biāo)定土體細(xì)觀參數(shù)、模型建立、動(dòng)力輸入等過(guò)程,利用PFC2D程序模擬了西吉縣興平鄉(xiāng)堡灣村下馬達(dá)子滑坡的失穩(wěn)破壞運(yùn)動(dòng)過(guò)程,得到如下結(jié)論:
1)下馬達(dá)子滑坡的失穩(wěn)機(jī)制是地震作用下的斜坡前緣牽引、后緣推擠、坡肩受拉破壞,由于坡肩高度高且具有臨空面,失穩(wěn)后坡肩位置速度和位移均較大,是地震滑坡破壞力強(qiáng)、致災(zāi)范圍大的主要原因;
2)由于黃土特有的垂直節(jié)理,黃土高原區(qū)常常存在高陡的斜坡,且重力滑坡通常具有高陡的后壁,而地震作用下滑坡的后壁相對(duì)平緩,這是地震滑坡區(qū)別于重力滑坡的重要特征之一;
3)通過(guò)顆粒流模擬得到滑坡前后相對(duì)高差為250 m左右,滑坡后總長(zhǎng)度達(dá)到1 080 m,滑坡后坡角約13.0°;與實(shí)際情況相對(duì)高差220 m,滑坡總長(zhǎng)度1 040 m,滑坡后坡角約11.94°相比較為符合,顆粒流方法可以用于地震滑坡滑距的預(yù)測(cè).
值得指出,擬靜力法是將地震期間最大慣性力施加在土體上,認(rèn)為土體中各點(diǎn)的最大加速度同時(shí)出現(xiàn),且沒(méi)有考慮土體隨時(shí)間變化的非線性,定量分析還不夠精準(zhǔn).今后的研究中,需要利用動(dòng)力時(shí)程分析法考慮地震動(dòng)的時(shí)間過(guò)程對(duì)滑坡失穩(wěn)破壞過(guò)程的影響,計(jì)算滑坡的演化過(guò)程.