基于顆粒流的地震液化型滑坡運(yùn)動學(xué)特征分析

常晁瑜，徐久歡，薄景山，楊濟(jì)源，孫雪晨，顧駿

（1.防災(zāi)科技學(xué)院，河北 三河 065201；2.中國地震局建筑物破壞機(jī)理與防御重點實驗室，河北 三河 065201；3.中國地震局工程力學(xué)研究所地震工程與工程振動重點實驗室，黑龍江 哈爾濱 150080）

引言

地震液化型滑坡是指斜坡土體在地震荷載的作用下，斜坡土體亞穩(wěn)定結(jié)構(gòu)破壞，產(chǎn)生超孔隙水壓力，降低抗剪強(qiáng)度，導(dǎo)致斜坡失穩(wěn)產(chǎn)生流態(tài)化運(yùn)動的過程。地震液化型滑坡發(fā)生坡度低，滑移距離遠(yuǎn)，破壞能力強(qiáng)，具有發(fā)生的隱蔽性和致災(zāi)的嚴(yán)重性雙重屬性，是對人民的生命財產(chǎn)安全造成危險的極大隱患。例如：1303年山西省大地震中洪洞縣郇堡即洪洞縣東北部坡度為2°左右的山前洪積扇上的郇堡地區(qū)超過20 km2的范圍內(nèi)發(fā)生了由于下伏地層液化而導(dǎo)致的大規(guī)模地滑，造成上部黃土塬支離破碎［1］。山西臨汾大地震1695年山西臨汾東堡村東堡村一帶發(fā)生液化滑移，致使上部地層“逆層滑移”［2］。1989年蘇聯(lián)塔吉克斯坦地震造成四處低角度液化滑移；致使220名村民死亡或失蹤［3］。1920年海原特大地震石碑塬［4］、黨家岔、夏家大路等地發(fā)生液化；坡度為2°～5°；滑距最遠(yuǎn)至6 km；致使多個村莊被毀。2008年汶川地震導(dǎo)致甘肅省清水縣田川村發(fā)生了黃土液化滑移，造成了較為嚴(yán)重的震害［5］。2013年甘肅岷縣漳縣永光村，永光村滑坡分為東西兩處，發(fā)育于同一坡體之上，物源基本一致。東側(cè)滑坡掩埋6戶人家，無人員傷亡；西側(cè)泥流狀滑坡掩埋8戶人家，造成12人被埋死亡［6］。2018年印度尼西亞東部蘇拉威西省首府帕盧市附近地震引發(fā)海嘯，震中位于帕盧市以北約78 km處，震源深度為20 km。地震導(dǎo)致了大規(guī)模液化滑移，對蘇拉威西島的中西部造成了嚴(yán)重破壞，2 256人死亡，1309人失蹤，近70 000座房屋損毀，超過22萬人無家可歸，總經(jīng)濟(jì)損失超過9.2億美元［7］。數(shù)次發(fā)生在黃土高原區(qū)的中強(qiáng)震均誘發(fā)過大規(guī)模液化型滑坡，造成了極大的人員傷亡和經(jīng)濟(jì)損失。因此，開展地震液化型滑坡的研究，提出針對性的防治措施，具有重要的科學(xué)價值和現(xiàn)實意義。

自發(fā)現(xiàn)地震液化型滑坡的存在后，眾多學(xué)者對這一類型滑坡的機(jī)理展開研究。白銘學(xué)等［8］提出了由黃土含砂層液化上涌推擠土層造成的黃土低角度滑移。王家鼎等［9］提出蠕動液化理論，并指出對于滑坡滑面平緩而又遠(yuǎn)程滑動最為普遍的解釋是震動液化。何開明等［10］總結(jié)了黃土液化研究的現(xiàn)狀與設(shè)想，同時建議對其開展5個方面的研究。孫萍等［11］通過對西吉黨家岔滑坡進(jìn)行環(huán)剪試驗，表明“滑動面液化”現(xiàn)象是高速遠(yuǎn)程滑坡形成的重要因素。張茂省等［12］提出對于地下水位淺或黃土層下部有泥巖隔水層的滑坡，地震液化是其破壞和運(yùn)動的主要機(jī)理。鄧龍勝等［13］通過研究隨機(jī)地震荷載下黃土中孔隙水壓力增長基本特性、影響效應(yīng)和增長規(guī)律，得到黃土在不同變形情況下孔隙水壓力的發(fā)展水平及液化程度。許領(lǐng)等［14］對涇陽南塬黃土滑坡的運(yùn)動特征進(jìn)行統(tǒng)計分析，發(fā)現(xiàn)該類滑坡的運(yùn)動特征受液化程度影響較大。Sassa等［15］研發(fā)了一系列高速不排水環(huán)剪儀，通過精準(zhǔn)監(jiān)測滑坡運(yùn)動過程中超孔隙水壓力的行為，從而定量驗證了滑坡液化減阻機(jī)制的正確性。Wang等［16］通過調(diào)查黨家岔滑坡，并對黃土試樣進(jìn)行了不排水三軸壓縮和環(huán)剪試驗，得出發(fā)生滑坡因素是黃土層液化而不太可能是黃土中的空氣這一結(jié)論。金艷麗等［17］通過ICU、ACU試驗得到了其穩(wěn)態(tài)性和潛在不穩(wěn)定區(qū)，認(rèn)為液化發(fā)生必要條件是土體狀態(tài)處于或被動處于潛在性不穩(wěn)定區(qū)域內(nèi)；蒲丹等［18］利用數(shù)值模擬對3種工況穩(wěn)定性系數(shù)進(jìn)行計算分析，認(rèn)為黃土低角度液化型滑坡運(yùn)動過程可分為震裂及液化、震陷形成滑帶和滑移3個階段；張曉超等［19］對西吉縣黨家岔滑坡的形成和運(yùn)動過程進(jìn)行了研究，確認(rèn)在強(qiáng)震過程中滑帶土產(chǎn)生了液化，并在地震力的作用下高速運(yùn)動。馬星宇等［20］歸納綜合物探、顆粒分析及動三軸試驗結(jié)果，總結(jié)了海原地震中石碑塬黃土地層液化滑移的形成機(jī)制。

綜合以上分析，以往的地震液化型滑坡研究主要針對滑坡機(jī)理展開，對地震液化型滑坡運(yùn)動特征研究不足，評價方法較少，文中嘗試離散元數(shù)值模擬進(jìn)行黃土液化型地震滑坡模擬，對地震液化型滑坡的發(fā)生過程進(jìn)行重現(xiàn)，重點研究地震液化型滑坡的運(yùn)動特征，通過與未液化滑坡的對比，提出地震液化型滑坡的運(yùn)動特點，為研究液化型滑坡的防治提供參考依據(jù)。

1 研究方法

地震誘發(fā)滑坡的發(fā)生過程存在滑動、平移、轉(zhuǎn)動等運(yùn)動，具有宏觀上的不連續(xù)性和運(yùn)動的隨機(jī)性，涉及到巖土體的大變形及巖土破壞問題，因此，文中利用顆粒流離散元方法進(jìn)行滑坡的運(yùn)動模擬。

1.1 顆粒流離散元

1971年，倫敦大學(xué)帝國學(xué)院Peter Cundall博士借鑒了分子動力學(xué)理論方法分析準(zhǔn)靜力或動力條件下巖石邊坡的運(yùn)動時，首次提出離散單元法（discrete element method，DEM）這一概念。而顆粒流（particle follow code，簡稱PFC）程序就是典型離散元法之一。顆粒流是從細(xì)觀結(jié)構(gòu)、不連續(xù)的角度出發(fā)，選取有代表性的數(shù)百個至上萬個顆粒單元，由平面內(nèi)的平動和轉(zhuǎn)動運(yùn)動方程來確定每一時刻顆粒的位置和速度，通過圓形顆粒介質(zhì)的運(yùn)動及其相互作用來模擬顆粒材料的力學(xué)特性。PFC2D不但能用來模擬固體力學(xué)靜態(tài)、動態(tài)問題，還能夠用于參數(shù)預(yù)測、在原始資料詳細(xì)情況下的實際模擬［21］。

1.2 計算原理

PFC采用顯示積分算法，對每個顆粒重復(fù)運(yùn)用牛頓第二定律，同時在接觸模型中運(yùn)用力－位移定律，將顆粒點的相對位移和顆?；ハ嗟慕佑|力分為法向分量和切向分量。利用力－位移定律通過法向剛度和切向剛度將法向接觸力、切向接觸力、法向相對位移和切向位移建立關(guān)系。若第i個顆粒接觸點的力Fi于法向切向分量為Fn、FS，則有

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

法向接觸力Fn可以按式（2）計算：

式中：kn為法向剛度；Un為顆粒重合長度。

切向應(yīng)力按增量方法計算，將接觸力與位移增量建立關(guān)系。接觸初始形成時切向應(yīng)力賦值為零。切向接觸力增量可按式（3）計算：

式中：ks為接觸面切向剛度；ΔUS為切向位移增量。

式中：μ為摩擦系數(shù)。

2 數(shù)值模型

2.1 模型概況

文中選取1920年海原8.5級特大地震誘發(fā)的大型滑坡——西吉縣興平鄉(xiāng)堡灣村下馬達(dá)子滑坡作為模擬對象。坐標(biāo)位置為N35°53′，E105°37′。下馬達(dá)子滑坡位于一條近南北走向的黃土梁西側(cè)，地形起伏較大，滑坡形態(tài)呈典型的上稍窄中寬長舌狀，滑坡周界清晰，主滑方向為310°?；驴傞L度約為1 280 m，上部寬度約為140 m，中部最寬達(dá)到270 m，下部寬度約為120 m，滑坡體總體積約為360萬m3?；麦w上部平均坡度約15°，厚度較薄僅約8～13 m；下部平均坡度小于10°，滑體較厚，據(jù)勘察資料，滑體厚度約15～25 m。

滑坡區(qū)整體與滑坡體上陡下緩類似，后緣較陡，平均坡度約為20°，后緣段長度約190 m且該部分上的滑坡體落水洞發(fā)育。滑坡區(qū)中部坡體平緩，中部向下至居民區(qū)長度約為640 m，下部滑坡體堆積區(qū)長度約為450 m。左右兩側(cè)均為滑坡下滑形成的陡坎，高5～20 m，陡坎中部較高，上部和下部相對較低，滑坡左側(cè)陡坎較高，右側(cè)稍低。兩側(cè)陡坎下部坡腳都有大量沿陡坎滑落或崩落堆積物?；w在下滑時遇到右前方山坡阻擋，轉(zhuǎn)向左側(cè)沿溝底滑出，因此滑坡右側(cè)下半部分以右前方斜坡坡腳為邊界?；轮芙绫容^清晰，滑坡后側(cè)為滑坡滑動時牽引形成的陡峭斜坡，最上方為高2～5 m的錯落陡坎，呈半封閉圈椅狀，其滑坡現(xiàn)場情況如圖1所示［22］。

圖1 西吉縣下馬達(dá)子滑坡現(xiàn)場照片F(xiàn)ig.1 Photographs of Xiamadazi landslide in Xiji County

2.2 模型顆粒參數(shù)標(biāo)定

由于顆粒流離散元軟件PFC2D中的參數(shù)不能由土工試驗直接得到，在運(yùn)用PFC2D進(jìn)行數(shù)值模擬時，需要標(biāo)定顆粒流模型的細(xì)觀參數(shù)。對于非動力學(xué)問題的分析，需要定義顆粒流的摩擦系數(shù)、顆粒密度、法向剛度、切向剛度等其它參數(shù)［23］。文中采用模擬直剪試驗進(jìn)行顆粒參數(shù)標(biāo)定。

對下馬達(dá)子滑坡后壁的原狀黃土進(jìn)行土工試驗，得到的具體參數(shù)見表1。在細(xì)觀參數(shù)標(biāo)定中，通過反復(fù)調(diào)試微觀力學(xué)參數(shù)，得到和土工試驗參數(shù)一致的力學(xué)行為后，認(rèn)為微觀參數(shù)可以用來模擬土工試驗得到的宏觀參數(shù)。通過反復(fù)直剪數(shù)值試驗?zāi)M，得到了軟件中賦予土體顆粒微觀力學(xué)參數(shù)，列于表2。

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

參數(shù)標(biāo)定中直剪數(shù)值試驗?zāi)Ｐ统叽鐬?00 mm×250 mm，生成顆?？倲?shù)為1 706個（圖2）。顆粒間的接觸模型選用平行粘結(jié)模型，對直剪模型分別施加100、200、300、400 kPa的豎向應(yīng)力，得到剪切應(yīng)力－位移關(guān)系曲線，通過計算得到直剪試驗中材料的宏觀力學(xué)參數(shù)。通過直剪試驗數(shù)值模擬得到的巖土體宏觀力學(xué)參數(shù)粘聚力c=19.25 kPa，內(nèi)摩擦角φ=13.5°。經(jīng)過直剪數(shù)值試驗后，得到的土體顆粒內(nèi)摩擦角和粘聚力大小與表2中宏觀參數(shù)基本對應(yīng)，故此土體顆粒細(xì)觀參數(shù)可用于數(shù)值模擬。模擬試驗的c－φ曲線如圖3所示。

圖2 直剪試驗?zāi)M模型Fig.2 Simulation model of direct shear test

圖3 模擬所得c－φ曲線Fig.3 Simulation of c－φ curve

表2 試驗組的模擬顆粒參數(shù)Table 2 Simulation particle parameters of test group

液化區(qū)的顆粒參數(shù)在原本基礎(chǔ)邊坡上層覆土的顆粒參數(shù)上進(jìn)行調(diào)整，液化試驗組參數(shù)如表3所示。由于PFC2D中難以呈現(xiàn)土樣的有效應(yīng)力的數(shù)據(jù)變化，故文中從土的基本強(qiáng)度指標(biāo)之一的內(nèi)聚力出發(fā)，觀察試驗組對應(yīng)的c－φ關(guān)系曲線的截距，即c值的變化。通過模擬直剪試驗可以得出試驗組內(nèi)聚力的值，將其視為c1，若c1值小于覆土層內(nèi)聚力c0值，可以認(rèn)為該組顆粒參數(shù)為發(fā)生液化后的液化區(qū)顆粒的參數(shù)。原覆土層黃土數(shù)值模擬采用顆粒參數(shù)見表2。

表3 液化試驗組的模擬顆粒參數(shù)Table 3 Simulation particle parameters of liquefaction test group

2.3 建模

根據(jù)下馬達(dá)子滑坡的地質(zhì)剖面圖（圖4），在PFC2D中進(jìn)行等體積法仿真建模，首先生成矩形區(qū)域顆粒球，并賦予顆粒參數(shù)使其對應(yīng)于室內(nèi)土工試驗土顆粒參數(shù)，運(yùn)用削坡法削坡并對顆粒施加重力，當(dāng)顆粒球的不平衡力小于10－5時，可以認(rèn)為模型自平衡，平衡后形成泥巖層斜坡。上層覆蓋黃土層采用落球法生成，進(jìn)行平衡循環(huán)后形成泥巖層上方的覆土層。模型總顆粒數(shù)24 812個，在模型表面顆粒上設(shè)置6個測點（圖5），記錄滑坡過程中運(yùn)動位移和速度等信息。根據(jù)高九龍的研究結(jié)果［24］，設(shè)置液化區(qū)如圖6所示，對液化區(qū)內(nèi)顆粒賦值液化參數(shù)，并再次平衡。

圖4 滑坡前下馬達(dá)子邊坡剖面圖（單位：m）Fig.4 Profile of Xiamadazi slope before landslide（Unit：m）

圖5 下馬達(dá)子邊坡模型和監(jiān)測點分布圖Fig.5 Distribution map of Xiamadazi slope model measuring points

圖6 含液化區(qū)的邊坡模型Fig.6 Slope model with liquefaction zone

2.4 動力輸入

采用擬靜力法模擬地震輸入［25］，擬靜力法實質(zhì)就是將地震動作用簡化成施加于計算對象重心的水平方向、豎直方向的恒定加速度作用，大小常用地震系數(shù)kh、kv表示，作用方向取最不利于坡體穩(wěn)定的方向。研究對象所承受水平和豎向地震作用力等于水平、豎向地震系數(shù)乘以其重量。

由現(xiàn)場勘察及后期地震災(zāi)害分析報告可知，該滑坡所在地在1920年海原地震中烈度為Ⅸ，故此次模擬地震橫波輸入模擬烈度Ⅸ度區(qū)，采用0.4 g施加于水平方向。

3 計算結(jié)果分析

文中采用顆粒流數(shù)值模擬是通過模擬滑坡模型與現(xiàn)場勘察資料中的滑坡剖面坡角、坡高、坡長等數(shù)據(jù)進(jìn)行比對驗證模型可行性，文中選用此方法時主要考慮該模型的滑坡前后相對高差、滑坡后總體坡長和滑坡后坡角。通過顆粒流模擬得到滑坡前后相對高差為250 m左右，滑坡后總長度達(dá)到1 080 m，滑坡后坡角約13.0°；與現(xiàn)場勘察資料中相對高差220 m，滑坡總長度1 040 m，滑坡后坡角約11.94°相比，較為符合真實情況［26］。

3.1 滑坡運(yùn)動過程分析

為探究地震液化滑坡模型的滑坡體運(yùn)動特征，運(yùn)用位移云圖對比分析地震液化型滑坡和未液化滑坡的運(yùn)動過程，如圖7所示。不難發(fā)現(xiàn)，地震液化型滑坡模型和未液化模型的位移有較大差別，模型運(yùn)行至10 s時，地震液化型滑坡模型出現(xiàn)較大面積的失穩(wěn)，而未液化模型僅在坡肩位置出現(xiàn)小面積失穩(wěn)運(yùn)動。模型運(yùn)行至20 s時，2個模型均存在較大面積區(qū)域土體顆粒發(fā)生位移，位移主要發(fā)生在坡肩位置，相較于未液化模型，液化模型中發(fā)生位移的顆粒數(shù)量更多，位移更大，表明滑坡在地震作用下液化后滑體的體積更大，滑動距離更遠(yuǎn)。模型運(yùn)行30～40 s時，所有坡段皆出現(xiàn)超出100 m大位移土顆粒，主要集中于坡肩段和坡中段。宏觀坡型變化上，地震液化模型較為明顯，其發(fā)生大位移的滑坡體占比也明顯高于未發(fā)生液化的模型。模型運(yùn)行50～60 s，滑坡過程結(jié)束，兩模型均發(fā)生較大改變。3個坡段滑坡體中大位移占比均達(dá)到90%，由圖7（e）、（f）中易見地震液化模型中大位移滑坡體占比要略高于未液化模型，同時也可看出液化模型的坡底段由于滑坡沖出的土顆粒數(shù)量明顯大于未液化模型，延伸部分也較長，可以說明地震動導(dǎo)致部分黃土液化后滑坡體更易滑出。

圖7 液化模型與未液化模型宏觀滑坡體位移云圖Fig.7 Macroscopic landslide displacement nephogram of liquefaction model and non－liquefaction model

3.2 各測點速度分析

據(jù)圖8可知，兩組曲線中峰值速度最大均為測點2，由圖5邊坡監(jiān)測點分布圖知該測點處于坡肩段位置，并且兩組曲線中均可看出測點2的速度大部分時間高于其他測點，這一現(xiàn)象與宏觀位移云圖中最先出現(xiàn)大位移區(qū)域為坡肩段相吻合。從圖8中可以看出各測點的速度液化后的運(yùn)動過程中速度均明顯提升。此外，未液化模型中除測點2外，測點1的速度最大，但在液化模型中，測點5的速度峰值明顯高于測點1，測點3、測點4的速度多數(shù)時步下也高于測點1速度。

圖8 兩類模型運(yùn)動過程各測點速度曲線Fig.8 Velocity curve of each measuring point in motion process of two models

將測點位置與其所處坡段位置結(jié)合來看，可以發(fā)現(xiàn)測點3、測點4、測點5均處于坡中段，這3個測點正下方為液化區(qū)設(shè)置區(qū)域，即液化后對滑坡體影響最大的區(qū)域就是與其距離最近的坡段。而相較之下測點1、測點6距離液化區(qū)較遠(yuǎn)，其運(yùn)動速度受到的影響相對較小。

3.3 各測點位移分析

由圖9和圖10中各測點水平位移與豎直位移變化情況可以看出，液化后相較于未液化模型測點3、測點4、測點5的水平方向位移變化量明顯高于其他測點，與速度變化量情況相吻合。測點4由于處于幾個測點中相對液化區(qū)最近的位置所受影響最強(qiáng)烈，其位移變化量最大與速度變化特征趨勢相吻合。值得提出的是，測點2的變化和測點6的變化較為特殊，測點2位置處于坡肩段并且坡度很大，當(dāng)兩模型均運(yùn)行到30 s后，測點顆?；瑒又疗轮卸危露茸兙徶率箖赡Ｐ蜏y點2位移曲線斜率變小。測點6位置接近坡底，液化模型運(yùn)動到后20 s時，該測點已經(jīng)滑出原邊坡區(qū)域，坡度變緩導(dǎo)致該測點水平、豎向位移斜率均變小。

圖9 兩類模型運(yùn)動過程各測點水平位移對比Fig.9 Comparison of horizontal displacement of each measuring point in the motion process of two types of models

圖10 兩類模型運(yùn)動過程各測點豎向位移對比Fig.10 Comparison of vertical displacement of each measuring point in the motion process of two types of models

4 結(jié)論

文中基于顆粒流數(shù)值模擬軟件對黃土地震滑坡和地震液化型滑坡進(jìn)行對比研究，得到以下結(jié)論：

（1）液化模型較未液化模型更容易發(fā)生失穩(wěn)破壞，速度更大、滑距更遠(yuǎn)、破壞性更強(qiáng)；

（2）無論是地震液化模型還是未液化模型，坡肩位置都是地震滑坡最先失穩(wěn)的部位，壩肩位置的滑體運(yùn)動速度大，滑距更遠(yuǎn)，在地震滑坡防治中應(yīng)重點關(guān)注；

（3）地震作用下發(fā)生液化后會改變滑坡各位置運(yùn)動過程中的位移和速度，液化區(qū)附近區(qū)域位移和速度變化量明顯增加，距離液化區(qū)較遠(yuǎn)部位位移和速度變化量較小。在液化區(qū)附近的區(qū)域會推動或者沖擊下部滑坡體，造成斜坡更大面積的失穩(wěn)，這也是地震液化型滑坡具有低角度，強(qiáng)破壞力特征原因之一。

值得指出的是，文中所得結(jié)論適用于黃土高原地區(qū)，砂質(zhì)黃土含量較高，或者黃土斜坡含水率較高、地下水分布較淺、覆蓋層厚度較小的砂質(zhì)黃土層。但是由于文中模型是對下馬達(dá)子地震誘發(fā)液化滑移的理論模擬，模擬動力輸入使用的擬靜力法是將地震期間最大慣性力施加在土體上，認(rèn)為土體中各點的最大加速度同時出現(xiàn)，且沒有考慮土體隨時間變化的非線性，定量分析還不夠精準(zhǔn)，呈現(xiàn)的現(xiàn)象具有一定的局限性。今后的研究中，將會逐步完善采用動力時程動力輸入法進(jìn)行更準(zhǔn)確地研究。