杨 博, 田文通,2,3, 孙军杰, 刘 琨,2, 徐舜华,2

(1. 中国地震局兰州地震研究所 中国地震局(甘肃省)黄土地震工程重点实验室， 甘肃 兰州 730000；2. 甘肃省 岩土防灾工程技术研究中心, 甘肃 兰州 730000;3. 兰州理工大学 土木工程学院， 甘肃 兰州 730050；4. 温州大学 建筑与土木工程学院， 浙江 温州 325035)

0 引言

黄土是一种多孔隙、弱胶结的第四纪沉积物，广泛分布于我国西北地区，由于其微观结构所具有的独特动力性质，使其表现出很高的地震易损性[1]。大量的震害分析表明[2-3]，地震诱发滑坡是主要的地震地质灾害类型之一。西部黄土地区地震频发，加上国家“一带一路”新兴战略的不断推进，催生了大量公路、铁路等生命线工程毗邻的人工高路堑黄土边坡，严重影响了生命线工程的安全运营， 因此，开展地震作用下高路堑黄土边坡稳定性研究具有重要意义。

针对黄土高路堑边坡的失稳，业界普遍认为是由于开挖卸荷引起应力重分布所致。研究手段大多为室内试验及数值模拟，并以宏观力学研究为主。学者通过边坡开挖变形的离心机模型试验研究，认为路堑边坡破坏从边坡顶部开始，属于蠕滑-压张拉裂[4-6]；王浩等[7]基于应力增量差值的变化解释了路堑边坡的开挖卸荷机理；张哲[8]以蒙内铁路为例，分析了膨胀土路堑边坡的失稳机理；胡晋川等[9]通过离心机试验和ANSYS有限元模型分析得出，黄土公路阶梯状高路堑边坡的变形存在先垂直、后水平的特点；张帅等[10]通过详细野外地质灾害调查，研究了我国黄土梁峁区地震诱发边坡——尉家湾巨型边坡；王鼐等[11]通过数值计算，分析了地震作用下黄土斜坡的动力稳定性和位移变形问题，推导出黄土滑坡运动的滑速和滑距公式。张灿灿[12]以徐州沙虎山边坡为例，运用FLAC3D结合强度折减法理论，分析了不同工况下边坡的稳定性。以上研究主要基于宏观力学方法，而针对土体内部细观力学机理研究仍显薄弱。

离散元方法能以颗粒为单位，应用牛顿第二定律，可模拟颗粒接触位置应力传递和位移的变化过程，对于研究地震作用下的大变形提供了基础。王家鼎等[13]认为地震力作用在每个土颗粒上，由于各土粒质量、排列方式、起始应力的不同，使各点所受作用力也存在明显差异，从而在土颗粒的接触点引起新的应力；石崇等[14]采用PFC2D模拟分析了地震作用下江坪河水电站陡岩的崩塌灾害，研究边坡破坏扩展过程；蒋明镜等[15]和廖优斌等[16]采用一种能够考虑水软化与化学风化作用的复合微观接触模型，将其植入PFC2D中，进行节理岩质边坡运动渐进直至失稳全过程的模拟，得出了边坡关键部位的位移、孔隙率、应变等时空演化规律。以上研究对象多基于岩质边坡，黄土属于散粒体，颗粒结构较为均匀，若能应用离散元方法进行细观力学分析，对于研究地震作用下高路堑黄土边坡失稳机理大有裨益。

本文以甘肃省内某在建高速公路典型深挖路堑边坡为对象，运用FLAC3D结合强度折减理论，对开挖后该边坡的稳定性进行判定后，通过PFC2D建立了边坡的颗粒模型，通过施加地震波荷载，研究了边坡地震动力响应及稳定性，揭示了地震过程中边坡土体颗粒间应力应变随时空的变化规律，对于认识地震作用下高路堑黄土边坡的失稳机理具有重要意义。

1 边坡区域工程背景

该高速路边坡工程位于甘肃省境内，工程区属黄土丘陵地貌，地形起伏较大，自然坡度35°～45°，斜坡自然稳定性较好；本工程土层大部分为风积黄土，分布在塬梁区坡顶，披盖式覆盖在基岩顶面，土质较均匀，孔隙发育，具有柱状节理，基岩主要是砂岩，砂质结构，抗风化能力差，遇水易软化，为软岩-极软岩。地质调查结果表明：该段区域地质构造稳定，无构造带通过；挖方段落内勘探深度内亦未见地下水出露；主要特殊性岩土为湿陷性黄土；该段工程场地土分类属于软土-中硬土，覆盖层厚度>5 m，工程场地类别划分为Ⅱ类场地。

2 黄土高路堑边坡初始状态

由于边坡开挖卸荷使得坡体内应力重分布，影响边坡整体稳定性，故进行开挖后边坡稳定性判定。

2.1 边坡尺寸以及物理力学参数

本研究依托甘肃省内高速公路高路堑边坡，边坡最高处高度约为54 m，长度为152 m，坡角约为40°，现场工程地质横断面如图1所示。室内试验得到本次数值模拟各岩土体的物理力学参数：黄土物理力学参数统计列于表1，岩土体物理力学参数指标列于表2。

图1 区域工程地质剖面图Fig.1 Regional engineering geological profile

表1 坡体黄土物理力学参数统计表

表2 坡体岩土体物理力学参数

2.2 边坡开挖后稳定性评价

当前边坡稳定性评价以数值模拟和极限平衡分析法为主。FLAC3D数值模拟由于无需假定滑移面位置和形状，较适用于边坡稳定性分析。结合强度折减理论，可分析处于极限状态下的土体应力变形特征。本文数值模拟法是基于强度折减理论，在未假定滑移面条件下，考虑岩土体的弹塑性，依靠折减岩土体的物理力学参数，而且能够对处于极限状态下的应力变形特征进行分析。为简化分析，本文选取左半边坡进行建模分析，建立的模型最高处约为55 m，水平长度约为112 m，边界条件设置为左右两侧约束水平位移，底边固定，上部自由，本构模型采用Mohr-Coulomb模型，初始地应力场是自重应力场。通过强度折减法计算出该边坡稳定性结果为1.12，如图2，该边坡在开挖后基本处于稳定状态。

从开挖后边坡位移云图(图3)可以看出，在自重条件下，总体位移向临空方向移动。沿着X方向的位移可以看出坡中发生明显的滑移，根据安全系数可知，该处边坡在开挖后处于稳定状态，但边坡开挖后要及时采取支挡防护措施，加强坡脚必要的加固措施，并对岩质坡面采取抗风化措施。数值计算与勘察资料结果一致，说明数值模拟具有一定的可靠性。

图4呈现了滑坡过程中沿着不同方向的应力分布图。通过应力分布图得出，滑坡过程中的最大剪应力集中在坡底，临空面附近也是应力集中带。在坡顶呈现出较小的水平应力，出现张应力。应力从坡底向坡体内部变化的过程中，应力状态由两相变为三相。

图2 安全系数Fig.2 Factor of safety

图3 开挖后边坡位移云图Fig.3 Nephogram of slope displacement after excavation

图4 开挖后边坡应力分布Fig.4 Stress distribution of slope after excavation

3 高路堑边坡动力响应分析

在分析该路堑边坡在开挖后初始状态的基础上，为进一步研究地震作用下细观力学过程，明晰动力响应过程中应力、应变、颗粒配位数、孔隙率等变化规律，采用颗粒流分析软件PFC2D建立数值计算模型，并对模拟结果进行了分析。

3.1 模型建立

从现场取样进行粒度分析和电镜SEM分析，粒度分析结果如图5所示，可以看出原状黄土粒径具有多样性，因粒径范围较大，故粒度采用对数坐标，粒度总体呈正态分布。室内电镜试验如图6，为不同放大倍数下黄土微观结构照片。由图可知，黄土颗粒粒径具有多样性，颗粒间形成点接触形式，孔隙较大，排列疏松，致使黄土骨架结构不稳定，在外力作用下容易发生破坏。

图5 粒度分析结果Fig.5 Particle size analysis result

根据场地工程勘察报告中地形地貌内容，该黄土边坡的坡角分布在35°～45°，因此本文建立了坡角为40°的边坡模型，为了便于观察边坡过程中颗粒位移的变化，对整个模型赋予了不同颜色，如图7所示。在边坡模型的的坡肩、坡体和坡角分别布设5个测量圆(图8)，以便对滑坡过程中颗粒位移进行监测。

图6 黄土SEM照片Fig.6 SEM photos of loess

图7 边坡离散元计算模型Fig.7 Discrete element calculation model of slope

图8 测量圆布设Fig.8 Layout of measuring circle

3.2 细观力学参数

模型分为黄土和砂岩两层，对于颗粒天然状态下的参数，根据数值直剪试验和物理直剪试验的参数标定拟合，由拟合试验结合经验参数，得出本文所用的颗粒细观参数取值如表3所列。

3.3 输入地震波荷载

根据国家地震局颁布的《中国地震动峰值加速度区划图》《中国地震动反应谱特征周期区划图》(G18306-2001图A和图B)，工程区地震动峰值加速度为0.10g，地震反应谱特征周期为0.45 s，相当于地震基本烈度Ⅶ度区。采用白银地区地震波，原始地震波荷载如图9所示。为了方便加载将其积分为速度，选取模型的底部边界作为地震施加源，传播方向自下而上。对于岩体边坡，岩体阻尼参数对结果往往有重要影响，颗粒流中一般采用局部阻尼考虑，其值按照0.7考虑。图10为地震波施加初期的位移和速度矢量图，由于地震荷载是通过模型底部施加，速度矢量图呈现出底部速度较大。而位移矢量图表明坡肩部分具有较大的位移，因此在坡肩部分发生滑坡的可能性较大，布设测量圆观察坡肩位移在滑坡过程中的变化趋势。

表3 数值试样细观力学参数

图9 地震波荷载输入Fig.9 Seismic wave load input

3.4 模拟结果分析

通过布设测量圆，在滑坡过程中监测了颗粒位移、应力、应变率及颗粒之间接触力随着时间变化的关系曲线。在地震荷载下，模型破坏以后位移、接触力、速度图，如图11～15所示。从图11可以看出，随着时间步长的增大，边坡变形逐渐增大，边坡破坏始于后缘的拉张裂隙，在重力作用下裂隙扩大推动滑动体滑移。边坡变形经历了应力重分布、后缘拉裂，蠕滑和加剧变形几个阶段。图12、图13反映了滑坡前后颗粒接触力的变化，由于滑坡过程中坡体中部土体被挤密，坡体中部较坡肩和坡底呈现出更大的接触力。图14是滑坡结束后的位移变化，可以明显看出滑坡前后，坡肩颗粒位移较大的部分滑移到坡底。图15呈现出模型破坏后速度变化，整体模型速度呈现出坡底速度较坡肩和坡中更大，这是由于颗粒在滑坡过程中滑移到坡底引起的。

图10 地震波荷载施加瞬间的矢量图Fig.10 Vector diagram of the moment of seismic wave loading

图11 地震荷载下破坏后模型Fig.11 Post failure model under earthquake load

图12 模型破坏前接触力Fig.12 Contact force before model failure

图13 模型破坏后接触力Fig.13 Contact force after model failure

颗粒滑动过程中配位数、孔隙比、应力和应变率的变化如图16～19所示。通过在边坡体的底部、中部、顶部和后缘布设测量圆，监测了这四个部位颗粒配位数、孔隙比、应力和应变率的变化。图16呈现了颗粒配位数的变化，可以看出整体配位数的变化都表现出先减小后不变的趋势，坡底随着时间步长的增大配位数中间一段呈现出增大的趋势，这是由于滑坡结束后坡肩和后缘部分的颗粒滑移至底部，造成底部颗粒配位数增加。同时边坡后缘呈现出急剧下降的趋势，是由于滑坡经历的最后一个阶段加剧变形导致颗粒迅速下滑。孔隙率变化(图17)中坡中部位增加更快，伴有颗粒的弹性变形、挤密的产生与迁移，颗粒变化复杂，说明边坡在坡中的变形更加剧烈，相比坡中，坡底、坡肩和后缘部位的颗粒孔隙率较小，边坡在这三个部位的变形较小。应力变化(图18)在初期阶段基本保持不变，随着时间步长的增大，呈现出先减小后增大的趋势，这是由于在滑坡蠕变阶段边坡体中颗粒重分布产生较多的自由颗粒导致的。应变的变化(图19)呈现出先减小后不变的趋势，是因为滑坡过程中发生大变形造成的这种结果。

图15 模型破坏后速度变化Fig.15 Velocity change after model failure

图16 颗粒配位数变化Fig.16 Change of particle coordination number

图18 颗粒间应力变化Fig.18 Stress change between particles

图19 颗粒应变率变化Fig.19 Change of particle strain rate

4 结论

通过有限差分法和离散元方法，进行了黄土地区高路堑边坡的动力失稳机制分析，主要结论如下:

(1) 通过FLAC3D强度折减法分析了开挖后边坡的初始状态，认为坡体最大剪应力集中在坡底，临空面也存在应力集中带，坡顶有较小的张应力。

(2) 地震作用下，边坡位移存在明显的顶部和临空面放大效应；地震波荷载施加瞬间，坡肩附近发生较大位移，主要是由于施加地震荷载时坡面表层加速度放大效应引起较大位移；根据模拟结果，可将滑坡过程分为应力重分布、后缘拉裂、蠕滑和加剧变形几个阶段。

(3) 根据接触力、颗粒配位数、孔隙率、应力等的监测结果表明，在应力重分布阶段，配位数等均保持稳定，颗粒几乎没有位移；在后缘拉裂阶段，坡肩处接触力增大，颗粒产生位移，斜坡模型速度变化不大；在蠕滑阶段，坡肩和坡体后缘下滑，坡体中部土体由于被持续挤密而变形剧烈，颗粒重分布产生自由颗粒，使得坡底颗粒配位数出现增加趋势；在急剧变形阶段，滑动面完全贯通，坡体整体失稳，边坡后缘处颗粒迅速下滑导致该位置配位数急剧下降。