基于颗粒流方法的黄土地震滑坡数值模拟＊

2022-08-23贾唯龙常晁瑜李佩茹张志伟徐久欢杨济源

地震学报 2022年4期

关键词：后缘黄土滑坡

贾唯龙常晁瑜李佩茹张志伟徐久欢杨济源

1）中国河北三河 065201 防灾科技学院

2）中国河北三河 065201 中国地震局建筑物破坏机理与防御重点实验室

3）中国哈尔滨 150080 中国地震局工程力学研究所

4）中国哈尔滨 150080 中国地震局地震工程与工程振动重点研究室

引言

受黄土动力易损性和南北地震带的双重影响，黄土地区发生地震滑坡的风险很高（徐张建等，2007）.历史研究表明，黄土地区一旦发生中强地震，往往诱发数量众多的滑坡，这些滑坡单体规模大、影响范围广、冲击速度快，对人民生命财产安全和城市建设的危害程度极高（常晁瑜等，2020）.实践证明，只有掌握地震滑坡发生机理，了解滑坡发生的过程和运动特征，才能采取有效的措施，避免或减轻地震滑坡的灾害（Zhouet al，2002）.

随着数值模拟方法的发展和计算机运算能力的提升，基于有限单元法、有限差分法和块体离散单元法等数值模拟方法，地震作用下斜坡的安全系数计算问题取得了极大进步（张江伟，2016；薄景山等，2019）.然而，强震作用下的滑坡是一个存在滑动、平移、转动的复杂过程，具有宏观上的不连续性和单个块体运动的随机性，当涉及到岩土体的大变形及岩土破坏问题时，上述数值模拟方法均具有一定的局限.而离散元颗粒流方法不受变形量限制（Hadjigeorgiouet al，2009；Tanget al，2009），可方便地处理非连续介质力学问题，有效地模拟岩土体的开裂、分离等非连续现象，对于探索地震滑坡的反映机理、运动过程和致灾范围具有一定优势（周健等，2000；曹文等，2017；石崇等，2018）.许多学者采用颗粒流离散元程序对边坡的滑移破坏机理开展了一系列的研究工作，例如：贺续文等（2011）讨论了节理连通率对边坡破坏形式的影响；周喻等（2016）从细观力学角度深入研究了顺层断续节理岩质边坡模型破坏过程的力学机制；李新坡和何思明（2010）对节理岩质边坡的破坏和运动过程进行了研究，分析了不同关键参数对破坏后的堆积形状和运动距离的影响；Scholtès和Donzé （2012）研究了含裂隙岩体的渐进性破坏.

通过以上分析，静力作用下滑坡运动失稳过程的研究进展较快，但针对地震滑坡失稳破坏运动过程的研究依然不成体系.鉴于PFC在模拟边坡失稳破坏方面的优势，本文拟在野外调查和室内试验的基础上，通过标定土体细观参数、模型建立、动力输入等过程，利用PFC2D程序模拟西吉县兴平乡堡湾村下马达子滑坡动力荷载下的失稳破坏运动过程，得到该地震滑坡的破坏运动机理，以期为黄土地区滑坡防治提供一定参考.

1 颗粒流方法基本理论

1971年，伦敦大学帝国学院Cundall博士在分析准静力或动力条件下岩石边坡的运动时，借鉴分子动力学理论方法，首次提出了离散单元法（discrete element method，缩写为DEM）的概念（Cundall，Strack，1979）.PFC （particle flow code）程序又称颗粒流方法（王光谦，倪晋仁，1992），是基于离散单元模型框架，由计算机引擎和图形用户界面构成的细观分析软件，主要用于模拟有限尺寸颗粒的运动和相互作用.

PFC的计算循环使用时间-步长显示算法，需要对每个颗粒反复应用牛顿第二定律，并对每个颗粒间的接触单元反复应用力-位移定律.在接触模型中，颗粒接触点的接触力和相对位移可以分解为沿法线和切线方向的两个分量，通过力-位移定律可以由法向刚度和切向刚度把接触力的法向、切向分量分别与法向、切向相对位移联系起来.第i个接触点的力Fi可以沿法向和切向进行分解，分解为法向力Fn和剪切力Fs，即

式中，ni，ti分别为接触平面的单位矢量.

两个接触实体之间的法向“重叠”量表示为Un，则法向力Fn可以表示为

式中，Kn为接触点的法向刚度，n为接触面的单位向量.

对于集合体来说，颗粒受到的剪切力与颗粒的相对运动、加载历史以及应力路径相关.因此，在颗粒流模型中，剪切力的计算以增量表示.切向力增量为

式中，ks为接触剪切刚度，µ为颗粒摩擦系数， ΔUs表示在计算时步内接触位移的切向增量，由下式得到:

式中， Δt为时间步长增量，vs为速度.

2 参数标定

利用颗粒流进行数值模拟，本质上就是从其细观力学特征出发，将材料的力学响应问题从物理域映射到数学域内进行数值求解（Wanget al，2003）.土的结构特征除土颗粒的大小、形状、表面特性及粒度级配特征外，还包括颗粒间的排列与集合关系，孔隙的大小，颗粒间联结的特点等，这些土体细观结构的变化很大程度上影响到土体宏观力学特性的变化（周健，池永，2003；Härtl，Ooi，2008；陈达等，2018）.利用PFC2D进行数值模拟之前，需要先假定土体颗粒间的本构特性，再用颗粒流方法模拟双轴压缩或直剪等土工试验过程（Park，Song，2009），对其力学性质进行数值试验，以此得到细观力学参数与宏观力学参数之间的对应关系.

通过试验，得到黄土的土工试验参数列于表1.通过细观直剪数值试验模拟宏观大型直接剪切试验，首先赋予土体颗粒微观力学参数（表2），剪切盒模型尺寸为500 m×250 m，生成颗粒总数为1 706个（图1）.颗粒间的接触模型选用平行黏结模型，在表2的细观力学参数条件下，对直剪模型分别施加100，200，300和400 kPa的竖向应力，得到剪切应力-位移关系曲线，然后拟合出土体颗粒抗剪强度曲线（图2），从而得到数值模拟所用岩土体颗粒的黏聚力c为19.25 kPa，内摩擦角φ为13.5°.经过直接剪切数值试验后，得到的土体颗粒内摩擦角和黏聚力大小与表1中宏观参数基本对应，故此土体颗粒细观参数可用于数值模拟.

表1 土工试验参数Table 1 Geotechnical test parameters

表2 细观试验参数Table 2 Micro-scale test parameters

图1 直剪数值试验模型图Fig. 1 Model diagram of direct shear numerical test

图2 土体颗粒抗剪强度曲线Fig. 2 Shear strength curve of soil particles

3 模型建立

3.1 滑坡基本情况

1920年海原MS8.5特大地震诱发了数以百计的黄土滑坡，造成了惨重的人员伤亡和经济损失，由于特殊的地理位置和气候条件，尽管经历了近百年时间，这些黄土地震滑坡的外部形态多数依然保存完好，具有极高的学术研究价值（许冲等，2018；常晁瑜等，2019）.

本文作者参与了对海原特大地震诱发滑坡的再调查，实地调查1920年海原特大地震诱发黄土滑坡620个，卫星影像补充调查滑坡605个，总计1 225个.调查内容包括地理位置、滑动类型、地层岩性、坡形、斜坡结构类型、滑坡基本特征以及影响因素等.调查所得滑坡分布如图3所示.滑坡集中分布在三个区域：① 宁夏西吉与甘肃静宁、会宁三县交会地带，区内滑坡分布密集，发育特征明显，代表滑坡为西吉县兴平乡堡湾村下马达子滑坡；② 宁夏海原九彩乡境内，区内滑坡主要沿深切河流两侧发育，代表滑坡为海原县李俊滑坡；③ 固原市区东部河川乡山地区域，区内滑坡分布较为零散，代表滑坡为固原市石碑塬滑坡.

图3 海原大地震诱发黄土滑坡分布图Fig. 3 Distribution map of loess landslide induced by Haiyuan landslide earthquake

下马达子滑坡野外编号为XJ06007，位于西吉县兴平乡堡湾村喜家湾组葫芦河支流滥泥河左岸黄土斜坡上，地理位置为（105.621°E，35.893°N），是1920年海原地震诱发的大型黄土地震滑坡.下马达子滑坡整体形状呈长舌状，滑坡周界清晰，主滑方向为310°，滑坡体上陡下缓，上部平均坡度约15°，下部平均坡度小于10°，长度达到1 280多米.由于上部土体下滑，导致上部较薄，厚度只有8—13 m，而下部滑体较厚，约15—25 m；后缘较陡，坡度约20°，该段长约190 m，滑体上发育大量落水洞.中部向下到居民区长约640 m，该段坡体平缓，滑体已被改造为梯田.下部滑坡堆积区长约450 m.滑坡上部宽约 140 m，中部较宽，最宽约270 m.滑坡总体积约360万 m3.滑坡区相对高差为220 m （图4）.

图4 下马达子滑坡现场调查图Fig. 4 Field investigation photo of Xiamadazi landslide

根据野外勘探（图5），下马达子滑坡地层岩性主要为上更新统马兰黄土、中更新统离石黄土和新近系砖红色泥岩，滑坡发育在黄土与泥岩的不整合接触上.由于黄土垂直节理极其发育，且其结构松散，当土体含水量大时，其抗拉、抗剪强度大大降低，其它条件具备时极其易发生变形.而且此处坡体为阴坡，受北部西伯利亚冷空气影响，黄土堆积厚度大，这也是影响坡体失稳的一个重要因素.

图5 下马达子滑坡剖面图Fig. 5 Profile of Xiamadazi landslide

3.2 模型建立

基于下马达子滑坡坡面图（图5），根据等体积法进行地形复原，所得滑坡发生前斜坡剖面图示于图6，而后通过颗粒流软件PFC2D进行建模，对泥岩颗粒施加重力后，颗粒做自由落体运动，循环一定次数后形成矩形，而后进行削坡操作，利用二次函数将泥岩层削出凹面，再对斜坡开挖形成基本泥岩层形状（图7），此时模型长约1 040 m，高度约200 m.对颗粒间不平衡力进行设置，多次运算后滑床颗粒静止平衡，此时锁定泥岩层颗粒速度及位移.接着采用落球法生成黄土颗粒，采用平行黏结模型，分三次下落到泥岩层指定位置，每次下落后黄土颗粒进行一次自重平衡，通过不断运算调整颗粒间不平衡力大小，直至不平衡力下降至允许范围，多次试算平衡后生成滑坡发生前基本模型，如图7所示.

图6 滑坡发生前斜坡剖面图Fig. 6 Profile before landslide occurrence

图7 滑坡开始前模型Fig. 7 Model before the beginning of landslide

3.3 动力输入

采用拟静力法对地震进行输入，拟静力法（刘红帅等，2007；Zhouet al，2015）是一种用静力学方法近似解决动力学问题的简易方法，它发展较早，迄今仍然被广泛使用.其基本思想是在静力计算的基础上，将地震作用简化为一个惯性力系附加在研究对象上，其核心是设计地震加速度的确定.由现场勘察及后期地震灾害分析报告可知，该滑坡所在地在1920年海原地震中的烈度为Ⅸ，故此次模拟采用Ⅸ度区的设计加速度即0.4g施加于模型上，使颗粒在水平方向上产生相应的惯性力.

4 结果分析

为了获得在外力输入下滑坡的运动情况，在关键位置设置6个测点（图8），以便于观察模拟中滑坡的发生机理和运动特征.

图8 测点分布示意图Fig. 8 Abridged general view of measuring points distribution

4.1 滑坡过程

滑坡模型运行60 s后停止，为了解滑坡的发生过程，截取不同时间下滑坡的位移图，示于图9.从图中可以清晰地看出：外力输入后，模型开始失稳，在20 s时，边坡坡度发生改变，上部滑体在地震波作用下产生明显位移，并挤压坡面颗粒，推动坡面颗粒向前运动，滑坡前缘颗粒在地震波作用下逐渐脱离滑床，开始运动；30 s时，坡面颗粒在上部滑体的重力和水平力的双重作用下，抗滑能力进一步降低，此时已可辨别出滑动面，滑面颗粒沿滑面不断向下运动，同时牵引后缘颗粒下滑；40 s时，随着滑坡后缘颗粒的逐渐下滑，滑坡后缘已经破坏，后缘颗粒已全部沿滑动面下滑，滑面颗粒在后缘颗粒的推移和前缘颗粒的牵引的共同作用下，速度进一步增大，位移增长迅速；50 s时，前缘颗粒在地面摩擦力作用下逐渐静止，在坡角处产生堆积，前缘颗粒位移达到40 m左右，此时坡面抗滑力逐渐增大，滑体与滑床接连处颗粒速度减小，逐渐趋于静止；60 s时，滑面颗粒运动缓慢且处于减速阶段，逐渐趋于稳定，此时滑坡后缘颗粒在滑面位移达到最大值100 m左右，边坡再次趋于稳定.

图9 不同时间-步长滑坡位移图（颗粒数22 733）Fig. 9 Landslide displacement at different time steps （The number of particles is 22 733）

整体上，滑坡的发生过程中，滑坡后缘颗粒的推挤作用和前缘颗粒的牵引作用均发挥了重要作用，滑体在推挤和牵引的双重作用下，逐渐发生失稳破坏，形成贯通滑动面，与滑床间摩擦力降低，直至发生滑动破坏.

4.2 速度变化

为了监测滑床内部与滑体的速度特征，检测了6个测点颗粒的速度变化，得到速度曲线示于图10.由图10可以看出：在整个滑坡过程中，滑床内部颗粒产生的速度变化较小，趋近于零；而滑坡后缘测点（测点1，2，3）速度变化较为剧烈.说明在地震波输入后滑坡后缘运动剧烈程度大于滑坡前缘（测点5），这与滑坡后缘失稳后重力势能转化为动能有关，测点3处速度高于其它五个测点，表明其首先被破坏并持续滑动.

图10 测点速度曲线Fig. 10 Velocity curve of measuring points

4.3 位移变化

通过对比滑坡发生前后的变化情况（图11），可以看出模型中颗粒发生了较大的位移（Behbahaniet al，2013）.为详细研究不同测点位移情况，通过对五个测点（测点1—5）的水平位移和垂直位移进行监测，得到如图12所示的位移曲线.通过对比可知，测点位置不同，水平和垂直方向的位移大小关系也不同，但大体上呈现出水平位移约为垂直位移的1—2倍.测点3处颗粒在水平和竖直方向上位移均最大，可见此处颗粒运动距离最大.测点5在运行50 s左右达到垂直位移最大值，但与其它测点不同的是，该测点处颗粒垂直位移未进一步明显发展，即该处颗粒在竖直方向不再发生明显运动，这是由于滑体在坡脚的堆积挤压作用，在竖直方向上，使得该测点颗粒克服重力，达到平衡状态，不再产生加速度，从而出现位移无明显变化的现象.在水平方向上，通过滑坡过程位移图可知，滑坡后缘处（测点1，2，3）水平位移基本较滑坡前缘处（测点5）大，这是由于滑坡前缘颗粒在破坏后期脱离滑面后产生堆积，受到地面摩擦力作用趋于稳定；而测点1处颗粒位移较小则是因为此处颗粒靠近滑坡后缘，滑坡后缘颗粒经滑面滑下后产生堆积再次趋于稳定，滑体与滑面间抗滑力逐渐恢复，此处颗粒未运动至坡脚.