基于顆粒流的地震液化型滑坡運動學特征分析

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

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

引言

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

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

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

1 研究方法

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

1.1 顆粒流離散元

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

1.2 計算原理

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

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

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

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

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

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

式中：μ為摩擦系數。

2 數值模型

2.1 模型概況

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

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

圖1 西吉縣下馬達子滑坡現場照片Fig.1 Photographs of Xiamadazi landslide in Xiji County

2.2 模型顆粒參數標定

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

對下馬達子滑坡后壁的原狀黃土進行土工試驗，得到的具體參數見表1。在細觀參數標定中，通過反復調試微觀力學參數，得到和土工試驗參數一致的力學行為后，認為微觀參數可以用來模擬土工試驗得到的宏觀參數。通過反復直剪數值試驗模擬，得到了軟件中賦予土體顆粒微觀力學參數，列于表2。

表1 土工試驗參數Table 1 Soil test parameters

參數標定中直剪數值試驗模型尺寸為500 mm×250 mm，生成顆粒總數為1 706個（圖2）。顆粒間的接觸模型選用平行粘結模型，對直剪模型分別施加100、200、300、400 kPa的豎向應力，得到剪切應力－位移關系曲線，通過計算得到直剪試驗中材料的宏觀力學參數。通過直剪試驗數值模擬得到的巖土體宏觀力學參數粘聚力c=19.25 kPa，內摩擦角φ=13.5°。經過直剪數值試驗后，得到的土體顆粒內摩擦角和粘聚力大小與表2中宏觀參數基本對應，故此土體顆粒細觀參數可用于數值模擬。模擬試驗的c－φ曲線如圖3所示。

圖2 直剪試驗模擬模型Fig.2 Simulation model of direct shear test

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

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

液化區的顆粒參數在原本基礎邊坡上層覆土的顆粒參數上進行調整，液化試驗組參數如表3所示。由于PFC2D中難以呈現土樣的有效應力的數據變化，故文中從土的基本強度指標之一的內聚力出發，觀察試驗組對應的c－φ關系曲線的截距，即c值的變化。通過模擬直剪試驗可以得出試驗組內聚力的值，將其視為c1，若c1值小于覆土層內聚力c0值，可以認為該組顆粒參數為發生液化后的液化區顆粒的參數。原覆土層黃土數值模擬采用顆粒參數見表2。

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

2.3 建模

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

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

圖5 下馬達子邊坡模型和監測點分布圖Fig.5 Distribution map of Xiamadazi slope model measuring points

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

2.4 動力輸入

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

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

3 計算結果分析

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

3.1 滑坡運動過程分析

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

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

3.2 各測點速度分析

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

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

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

3.3 各測點位移分析

由圖9和圖10中各測點水平位移與豎直位移變化情況可以看出，液化后相較于未液化模型測點3、測點4、測點5的水平方向位移變化量明顯高于其他測點，與速度變化量情況相吻合。測點4由于處于幾個測點中相對液化區最近的位置所受影響最強烈，其位移變化量最大與速度變化特征趨勢相吻合。值得提出的是，測點2的變化和測點6的變化較為特殊，測點2位置處于坡肩段并且坡度很大，當兩模型均運行到30 s后，測點顆粒滑動至坡中段，坡度變緩致使兩模型測點2位移曲線斜率變小。測點6位置接近坡底，液化模型運動到后20 s時，該測點已經滑出原邊坡區域，坡度變緩導致該測點水平、豎向位移斜率均變小。

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

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

4 結論

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

（1）液化模型較未液化模型更容易發生失穩破壞，速度更大、滑距更遠、破壞性更強；

（2）無論是地震液化模型還是未液化模型，坡肩位置都是地震滑坡最先失穩的部位，壩肩位置的滑體運動速度大，滑距更遠，在地震滑坡防治中應重點關注；

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

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