《1、 引言》

1、 引言

地震是诱发滑坡的主要原因之一[12]。地震诱发的山体滑坡,如岩石边坡垮塌和泥石流(图1)[36]造成的人员伤亡和经济损失甚至比地震直接造成的损失还要大。目前遥感和地理信息系统(GIS)已广泛应用于地震诱发滑坡的地表监测和分析[710]。深入理解地震诱发滑坡的动力机制是预测滑坡、防治滑坡灾害的基础。为了更加系统地了解地震诱发滑坡的动力机制,研究人员采用多种研究方法评估了地震诱发的永久位移,分析了边坡的受力状态[11]。此外,研究人员还研究了地震波荷载作用下的裂纹起裂、扩展和贯通过程及破坏模型。将动态剪切裂纹模型应用于断层地震过程的研究。Dalguer等[13]模拟了地震中三维(3D)动态剪切破裂扩展引起的拉伸裂纹起裂、扩展过程。

《图1》

图1 地震诱发的滑坡案例。(a)罐滩滑坡[3];(b)青川县董家村滑坡[4];(c)头寨滑坡[5];(d)大光包滑坡[6]。

上述研究主要集中在分析地震滑坡的现象、滑坡分布,或者破坏模型。然而,这些研究并没有明确地震波荷载作用下岩体内部裂纹起裂和扩展的机理。裂纹的起裂、扩展和贯通以及局部变形破坏造成的强度退化是节理岩质边坡破坏的主要原因[14],同时岩体内部的不连续结构(如天然裂缝、层理、断层、节理)也会显著影响岩石破裂过程。在各种外部地质力的作用下(如沉积作用、河流侵蚀、风化作用[1819]),张开型裂隙广泛存在于岩体中(图2)[1517]。当地震发生时,裂纹通常在这些张开裂隙周围起裂、扩展,导致岩体强度降低。主震发生后通常还会发生一定数量的余震,在余震过程中这些起裂岩体容易受到进一步破坏而导致滑坡[2021]。

《图2》

图2 岩体中的张开型裂隙。(a)岩体中的类圆孔张开裂隙[15];(b)岩体中张开裂隙[16];(c)岩石边坡中张开裂隙[17]。

本文采用颗粒元黏结模型(bonded particle model, BPM),研究了地震波荷载对含单条预制裂纹的模型中的裂纹起裂和扩展的影响。黏结模型不仅可以模拟岩石裂纹起裂扩展和累积损伤,还可以模拟地震波荷载作用下岩石内部的波速的变化[22]。对于在实验室或现场研究中无法直接测量或观察的现象,该模型可以较好地模拟重现如裂纹起裂、扩展、贯通、应力波传播、地震波传播和速度变化等过程[2324]。

地震波荷载作用下的裂纹起裂和扩展机制与准静态荷载作用下的具有较大不同[2526]。研究人员利用黏结模型研究了不同类型岩土边坡在地震波荷载作用下的地震反应。 例如,有学者采用二维黏结模型模拟了Chi-Chi地震引发的Tsaoling滑坡及其发生机理[2728]。有学者采用三维黏结模型研究了干燥状态下颗粒土中的波传播过程[29]。这些研究主要集中在主震引起的地震滑坡,结果表明在地震作用下边坡容易发生破坏。然而,在外力或地震作用下,边坡中如果仅有少量离散的裂纹并不会立即导致边坡的破坏。而如果反复对已经存在少量离散裂纹的边坡施加力或地震波,这些裂纹将会进一步地扩展、贯通,最终导致边坡破坏。为了深入研究含初始裂纹的岩石在循环地震波作用下的破坏机制,本文采用二维黏结模型研究了地震波在同一方向和两个正交方向上反复施加地震波载荷条件下裂纹的起裂、扩展过程及破坏机制。

《2、 模拟方法》

2、 模拟方法

《2.1 黏结模型》

2.1 黏结模型

既有研究已对黏结模型的原理做过详细介绍[3031],因此,本文仅对其主要特点进行简要介绍。根据接触特点,可以将黏结模型分为两类:接触黏结模型和平行黏结模型。平行黏结模型中颗粒之间既可以传递力也可以传递力矩,多用于岩石材料的模拟。可将平行黏结模型的黏结关系视为相邻颗粒间有一组弹性弹簧,弹簧位于接触面上并以接触点为中心,弹簧均匀分布在一个具有恒定法向刚度和剪切刚度的矩形截面上。这组弹簧类似于梁,可以抵抗张拉、剪切或由颗粒旋转引起的力矩。当作用在黏结上的力超过黏结的法向或剪切强度时,黏结断裂,形成微裂纹(图3)[32]。该模型可以呈现不同的破坏模式,如张拉或剪切破坏。Zhang和Wong [3334]采用平行黏结模型研究了不同颗粒尺寸及加载速率(0.005~0.600 m·s-1)对岩石材料损伤演化过程的影响。研究结果表明,平行黏结模型能够真实模拟岩石和类岩石材料在准静态压缩荷载下的裂纹起裂、扩展和贯通过程。

《图3》

图3 平行黏结破坏准则。(a)法向力与法向位移;(b)剪切力与剪切位移;(c)强度包络线[32]。

《2.2 黏结模型中的波传播》

2.2 黏结模型中的波传播

地震波是指在地震发生时,积聚的能量以弹性波的形式向四面八方释放出来。由于地球内部的非均质性,波的传播路径不是线性的,经过多次反射和透射后呈现复杂的传播路径。为了真实地模拟地震波引起的裂纹扩展,首先要阐明试样边界的反射和透射现象。在本文中,除非另有说明,仅使用最简单的入射角形式,即法向入射。Joyner和Chen [35]采用了一种方法,可以精确满足波从下面垂直入射。该方法类似于Lysmer和Kuhlemeyer [36]所使用的方法,即假定在弹性介质中存在垂直入射的平面横波。基于这一假设,可以推导出入射波的质点速度和边界处剪切应力的表达式。质点速度(vB)和剪应力(τB)公式如下:

vB=vR+vI(1)

τB=ρvs(2vI-vB)(2)

式中,vIvR分别为入射波和反射波的质点速度;vs为介质中的剪切波速;ρ为介质密度。当弹性波通过介质Ⅰ到达介质Ⅱ的边界时,入射波将变为两种类型的波:反射波和透射波[37]。反射波振幅(R)和透射波振幅(T)如下:

R=ZI-ZIIZI+ZII(3)

T=2ZIZI+ZII(4)

式中,Z为声阻抗,Z = ρPw = (ρEe)1/2。其中ρ为介质密度,Pw为波速,Ee为弹性模量。如果两种介质的声阻抗不同,则透射波和反射波的振幅和速度也不同。例如,如果介质II的声阻抗ZII趋于0,透射波的振幅是入射波的两倍。反射波振幅和速度与入射波相同。当介质II的声阻抗趋于无穷大时,透射将不再发生,入射波将全部反射回介质I中。

为了确保模拟结果的准确性,首先在黏结模型中模拟上述两种情况下弹性波的反射和透射情况,以及它们的振幅和速度之间的关系。图4为由50个颗粒组成的一维杆体(bar)。为了吸收入射波,将杆体的左端设置为静止边界,右端设置为自由边界,即模拟介质II的声阻抗为零的情况。公式(5)描述了输入速度脉冲与左侧边界上颗粒速度的关系,即模拟外力作用在左侧球体上。外力F与颗粒速度vB的关系为:

F=πRp2Pwρ(2U.(t)-vB)(5)

式中,Rp为颗粒半径;Pw为波速;U.(t)为给定的速度脉冲。

《图4》

图4 50个颗粒组成的一维杆体。

速度脉冲公式为:

U.(t)=A2(1-cos (2πft))(6)

式中,A是脉冲的振幅;f为频率。颗粒半径为0.5 m,波速为100 m·s-1,给定速度脉冲时间t为0.25 s,频率为4 s-1。杆体左侧、中间和右侧颗粒的速度如图5所示,其中黑色线段表示左侧颗粒的速度曲线,红色线段表示中间颗粒的速度曲线,绿色线段表示右侧颗粒的速度曲线。第一个曲线和第二个曲线分别为左侧颗粒和中间颗粒的速度,第三个曲线为右侧颗粒的速度。由图可知,右侧颗粒速度是左侧颗粒和中间颗粒速度的两倍。第四个曲线为中间颗粒的反射波速度,第五个曲线为左侧颗粒的反射波速度。由于左侧颗粒设置为静止边界,反射波传播至左侧边界后能量被完全吸收,弹性波没有再次发生反射[38]。图6为右侧颗粒声阻抗为无穷大时的左侧颗粒、中间颗粒和右侧颗粒的速度曲线。右侧颗粒声阻抗无限大,弹性波传播至右侧颗粒时,未发生透射;而中间颗粒和左侧颗粒的反射波速度与入射波速度大小一致,但方向相反。

《图5》

图5 自由边界条件下颗粒速度曲线。

《图6》

图6 固定边界条件下颗粒速度曲线。

上述例子分别模拟了两种特殊情况,即阻抗为零的自由边界和阻抗为无限大的边界情况。结果表明,模拟结果与Jaeger等[37]理论结果一致,说明黏结模型可以较好地模拟弹性波的传播。

《2.3 数值模型》

2.3 数值模型

本文中模拟试样尺寸为76 mm × 152 mm,该尺寸的岩石试验在实验室中已被广泛应用于岩石材料中裂纹起裂和扩展的研究[3941]。数值模型大约由38 000个颗粒组成,最小颗粒半径为0.21 mm,最大颗粒半径为0.35 mm,试样颗粒尺寸在最小和最大尺寸之间均匀分布。在黏结模型中,颗粒半径对模型的宏观力学性质有显著影响;然而,在黏结模型中颗粒的大小无法与岩石内部的矿物颗粒尺寸对应。为了更加真实地模拟岩石损伤演化和破坏过程,数值模型中要求颗粒数量不能太少,但目前无具体的数量要求[20]。黏结模型中的微观参数见表1。为了更加真实地模拟地震波在黏结模型中的传播过程,首先将数值模型的力学性能与实验室试验得到的力学性能进行对比、标定(表2)。其中,数值模型的单轴抗压强度(UCS)与平行黏结的抗拉和剪切强度成正比;弹性模量与颗粒的刚度和平行黏结强度成正比;泊松比μnp/sp(颗粒的法向刚度与剪切刚度的比值)和nb/sb(平行黏结法向与剪切刚度的比值)成正比;数值模型的破坏模式与σn/σs(法向强度与剪切强度之比)有关。然而,在模拟过程中,无法直接校准σn/σs的比值。因此,在参数校准过程中,需在保持其他参数不变的情况下,首先以不同的σn/σs比值进行一系列单轴压缩试验,将模拟结果与室内试验结果进行比较,得到数值模型的UCS、弹性模量、泊松比和破坏模式。随后,改变σn/σs比值,进行一系列的数值试验,将裂纹起裂位置、起裂角度和裂纹形态与室内试验结果进行比较,最终确定数值模型中的微观参数,本文中数值模型使用的微观参数如表1所示。有关黏结模型参数校准更具体的过程和细节,请参考文献[31]。完整试验标定完成后,在数值模型样品中设置一个长12.6 mm、宽1.3 mm的单条裂隙(图7),裂隙位于试验的几何中心,裂隙倾角β = 30°。由于不能再进一步划分黏结模型中的颗粒,因此裂隙表面局部较粗糙(图7)。

《表1》

表1 黏结模型中的微观参数

Particle parametersParallel bond parameters
VariableMeaningValueVariableMeaningValue
EcYoung’s modulus of particle (GPa)3.95EpYoung’s modulus of parallel bond (GPa)3.95
np/spStiffness ratio1.00nb/sbStiffness ratio1.00
ϕFriction coefficient0.10σnNormal strength (MPa)24.5 ± 6.5
Rmax/RminRadius ratio1.66σsShear strength (MPa)34.5 ± 6.5
RminMinimum radius (mm)0.21λRadius multiplier1.00

《表2》

表2 数值模型与岩石室内试验力学性质对比[]

ParameterPhysical experimentPFC simulation
Density, ρ (g·cm-3)1.541.54
Young’s modulus, E (GPa)5.966.02
Poisson’s ratio, μ0.150.16
UCS (MPa)33.8533.30
Tensile strength, σt (MPa)3.204.33

《图7》

图7 数值试样(76 mm × 152 mm)。裂隙长12.6 mm,宽1.3 mm,裂隙倾角为30°。

《2.4 边界条件》

2.4 边界条件

当地震波从一种介质传播到另一种介质时,在两种介质交界处通常发生透射和反射。Lysmer和Kuhlemeyer [36]指出为了防止透射和反射的发生,在边界处通常建立黏性边界用来吸收地震波能量。然而,如果数值模型的尺寸与P波波长相比非常小时,则在数值模型中观察不到反射和透射现象。P波速度vP定义为:

vP=E(1-μ)ρ(1+μ) (1-2μ)(7)

式中,E是弹性模量。本文数值模型中弹性模量E = 6.02 GPa,泊松比μ = 0.16,材料密度ρ = 1540 kg·m-3,P波波速vp = 2040 m·s-1。地震波的平均周期约为0.6 s,因此,P波波长约为1224 m。数值模型长度为152 mm,占P波波长的0.012%。因此,数值模型中无需设置黏性边界。当地震波垂直加载时,将数值模型的上下边界设为固定边界,其中边界颗粒的速度为定值,在加载过程中不变。将模型侧向边界设置为自由边界,其颗粒速度可以变化[图8(a)]。当侧向施加地震波荷载时,侧向边界为固定边界,模型上下边界为自由边界[图8(b)]。

《图8》

图8 数值模型边界条件。(a)垂直加载地震波,上下边界为固定边界,侧向为自由边界;(b)侧向加载地震波,上下边界为自由边界,侧向为固定边界。

《2.5 地震波加载》

2.5 地震波加载

地震波首先从数值模型顶部加载,之后再从模型右侧进行加载,加载过程中需监测模型中颗粒速度及模型中裂纹起裂和扩展过程。本文选用的地震波为1995年1月17日日本神户地震的地震波波形,神户地震震级7.3级,持续时间约为21 s。图9图10给出了地震波中的速度(最大-90 cm·s-1)和加速度(最大801.63 cm·s-2)。颗粒元软件说明书[38]建议在进行压缩试验时,试样加载速率为0.02 m·s-1,同时还指出了加载速率需足够慢,以确保模型保持准静态平衡。Zhang和Wong [33]建议在进行单轴压缩试验和巴西试验时,加载速率分别取0.02 m·s-1和0.01 m·s-1;当加载速率超过0.08 m·s-1时,多余的能量会转化为动能,如果加载速率再大,模型将处于动态加载状态。因此,为保证模型保持准静态平衡,所输入地震波的加载速率比神户地震的地震波的加载速率小一个数量级,输入波形的速度大部分小于0.08 m·s-1(地震波最大速度为0.09 m·s-1)。颗粒元黏结模型计算过程是基于牛顿第二定律。为了确保模型在计算过程中保持稳定,每个计算周期中的时间步长(∆t)不超过临近阈值,临近阈值取决于颗粒的刚度、密度和几何形状。在本文采用的数值模型中,每个计算周期中的时间步长较小,约为10-8 s,即数值模型中每一个时间步长对应物理时间10-8 s。0.08 m·s-1的加载速率在数值模型中为每个时间步长产生8 × 10-7 mm的位移,即每移动1 mm需要超过100万步。

《图9》

图9 地震波速度曲线。

《图10》

图10 地震波加速度曲线。

《3、 结果分析》

3、 结果分析

《3.1 裂纹起裂与扩展》

3.1 裂纹起裂与扩展

图11(a)给出了模型在轴向地震波加载作用下的裂纹起裂情况,黄色线段表示微张拉裂纹。在时间步长为1͏ 644͏ 627时,少量微裂纹从预制裂隙尖端开始起裂,裂纹的扩展方向与地震波加载方向垂直。随着地震波荷载的继续,裂纹没有进一步扩展。第一次地震波加载完成后,一共产生了60条微裂纹。第一次加载完成后,在轴向再次施加地震波载荷[图11(b)]。在第二次地震波加载期间,微裂纹数量没有发生变化,即第二次地震波加载过程中未产生新的微裂纹。图11(c)为两次地震波加载的时间-速度和时间-位移曲线,时间-位移曲线由时间-速度曲线积分变换获得。由图11可知,在第二次地震波加载过程中,裂纹没有进一步扩展,这种效应被称为Kaiser效应。Kaiser效应是指在首次加载过程中监测到岩石破裂产生的声发射信号,如果在随后的加载过程中,当加载应力低于先前的峰值应力时,岩石将不会发生破坏,不再产生声发射信号;当后续加载应力超过先前的峰值应力时,岩石发生破坏,产生声发射信号[42]。

《图11》

图11 地震波循环加载。(a)第一次地震波轴向加载;(b)第二次地震波轴向加载;(c)地震波时间-速度曲线和时间-位移曲线。

作为对比,在第一次地震波加载后[图12(a)],第二次加载时改变加载方向,即从模型右侧加载。在第二次加载过程中,时间步长为7 772 490,在预制裂隙的左下角和右上角开始出现新的微裂纹,裂纹扩展方向垂直于侧向加载方向[图12(b)]。之后再次在轴向施加地震波载荷,裂纹沿水平方向开始扩展,最终试样被完全破坏[图12(c)]。Lavrov等[43]在巴西劈裂试验中观察到了Kaiser效应,并指出Kaiser效应对加载方向较敏感,当加载方向旋转角度为15°或更大时,Kaiser效应将消失。此外,一些声发射实验也证明了Kaiser效应对加载方向有较强的依赖性[4445]。在本次加载过程中,在第二次侧向加载地震波荷载时,Kaiser效应消失,预制裂隙尖端形成少量微裂纹。在第二次加载过程中出现的新裂纹改变了原有预制裂隙周围的应力分布,因此,在第三次地震波加载过程中,裂纹沿起裂方向扩展,最终试样被完全破坏。该现象说明地震波荷载引起的裂纹起裂及扩展与地震波荷载方向有关。在同一方向反复施加地震波荷载时,已经起裂的裂纹,无法进一步扩展,而将地震波荷载的方向由轴向变为侧向进行多次加载后,裂纹将进一步扩展。

《图12》

图12 地震波循环加载。(a)第一次地震波轴向加载;(b)第二次地震波侧向加载;(c)第三次地震波轴向加载;(d)地震波时间-速度曲线和时间-位移曲线。

在三次地震波加载过程中,裂纹的起裂和扩展发生在张拉位移峰值处[图12(d)中的点A、B、C]。关于裂纹为何会在最大张拉位移出起裂,而不是在最大压缩位移或峰值速度处起裂,将会在3.2节中详细讨论。由图12可知,地震波加载过程中,模型中产生的裂纹均为张拉裂纹(模型中黄色线段)。如果将地震波(加载的地震波波形相同但振幅大)施加到完整数值模型试样上,模型破坏产生的裂纹也是张拉裂纹(图13),表明地震波在地震中主要通过张拉裂纹破坏岩石。

《图13》

图13 完整试样在地震波荷载下的破坏形态及裂纹分布。

《3.2 裂纹起裂时机》

3.2 裂纹起裂时机

为了准确确定裂纹的起裂时机,在裂纹起裂位置附近选取一个颗粒[在图12(a)中标注为绿色;x = -7.4 mm, y = -2.89 mm ],监测其位移和速度。图14为监测颗粒在第一次地震波加载时的位移和速度曲线图,裂纹起裂时间(图14点A)与最大速度和最大压缩位移(-0.15 mm)并不对应,而是与其最大张拉位移(0.11 mm)对应。本文张拉位移为正,压缩位移为负。在黏结模型中,最大拉应力或剪切应力与颗粒的位移密切相关。简而言之,在最大张拉位移处,颗粒所受张拉应力和剪切应力也会达到峰值。同时,岩石材料的抗压强度一般是其抗拉强度的4~25倍。因此,裂纹在最大张拉位移处起裂,而不在最大压缩位移处起裂。监测结果表明,裂纹起裂后,颗粒的速度有较大波动,这是由裂纹周围应力集中突然释放引起的,经过几个时间步长后,应力重新分布并达到新的平衡,颗粒速度恢复正常(图14)。

《图14》

图14 第一次地震波加载过程中监测颗粒的位移和速度曲线。

《3.3 模型中速度分布》

3.3 模型中速度分布

图15图16给出了不同位置颗粒的速度-时间步长曲线:模型顶部边界颗粒(x = -0.03 mm, y = 74.41 mm)、模型1/4处颗粒(x = -0.02 mm, y = 37.15 mm)、模型1/2处颗粒(x = -0.02 mm, y = 2.73 mm)和模型3/4处颗粒(x = -0.02 mm, y = -37.45 mm)。如图16(b)所示,为了更加直观地对比和观察不同位置的速度,图中未绘制剧烈波动时期的粒子速度。如2.3节所述,地震波加载过程中模型中未发生反射或透射。4个颗粒的速度曲线较类似,但速度大小具有较大不同,呈线性分布vAvBvCvD = 4∶3∶2∶1。在物理试验中,应变速率通常比加载速率更常用,应变速率(ε¯)指应变的变化率[33],具体表达式为:

ε¯=dεdt=ddtL-L0L0=1L0dLdt=vlL0(8)

式中,ε为应变;t为时间;L为模型试样长度;L0为模型试样原始长度;vl为加载速率。应变速率在整个试样中是均匀的。例如,在边界附近的颗粒(点A)的速度是位于中间区域颗粒(点C)的两倍。由上可知,在离散元模型中,速度沿加载方向呈线性分布,模型中的颗粒速度逐渐变化,类似于连续材料中的速度分布。

《图15》

图15 4个颗粒的位置分布。

《图16》

图16 4个颗粒速度曲线,其中图(b)为图(a)的放大图。

《4、 讨论》

4、 讨论

地震发生时,作用在岩体上的地震波荷载均可分为水平分量和垂直分量。岩体受到张拉和剪切的共同作用。地震波荷载作用下的岩体破坏崩塌是一个复杂的过程,涉及地震波传播、能量衰减、裂纹起裂、扩展和贯通。为了更好地理解地震波荷载方向与裂纹起裂方向之间的关系,基于岩石断裂力学和地震能量角度已经开展了一些研究[4647]。研究结果表明,随着地震波输入方向从0°逐渐增加至90°时,裂纹起裂角从0°逐渐增加至70.5°(图17)[47],这一现象与本文研究结果较一致。

《图17》

图17 地震波加载方向与裂纹起裂角度关系[47]。

本文研究表明,裂纹的破坏形式主要是张拉破坏,这一结果与滑坡[48]和实验室测试[49]的观察结果基本一致(图18)。地震过程中张拉应力使得滑坡岩体断裂面呈锯齿状或较为粗糙,这与由剪切应力引起的重力滑坡的光滑弧形断裂面不同。地震波荷载在地震滑坡的早期阶段起着重要作用,它产生的拉应力导致边坡产生张拉裂纹(图19)[6]。拉应力在滑坡变形中起主导作用,同时地震波荷载导致岩体破裂,岩体黏聚强度降低,随着裂纹的不断扩展和岩体强度的降低而诱发滑坡。

《图18》

图18 滑坡过程中观察到的张拉裂纹[4849]。

《图19》

图19 滑坡破坏机制示意图[6]。

上述结果解释了为何大量的山体滑坡发生在余震期间。由于野外大型边坡和岩体内部发生反射和透射,岩石反复受到地震波载荷作用,同时地震波的作用方向不同,导致岩石边坡初始裂纹进一步扩展、贯通,最终发生破坏。需要说明的是,本文利用一维模型研究了地震波的传播和岩石破裂过程,但地震波荷载作用下岩石损伤演化和破坏过程还受到预制裂隙的大小、位置、倾角、岩体含水率、围压、岩石类型等因素的影响。因此后续还需通过实验室试验、数值模拟和现场调查等多种手段,考虑地形、地质条件、地形放大和地应力状态的影响,研究地震波在大尺度边坡的传播过程、岩体破坏过程及滑坡诱发机制。

《5、 结论》

5、 结论

本文基于平行黏结模型,从两个正交方向模拟了循环地震波荷载作用下的裂纹起裂、扩展过程,讨论了地震波荷载作用下裂纹起裂扩展的方向和裂纹破坏机理。主要结论总结如下:

(1)在地震波加载过程中,张拉破坏是震动荷载条件下的主要破裂机制。

(2)在地震波加载序列中,裂纹起裂和扩展不是发生在加速度最大或速度最大处,而是发生在张拉位移最大值处。

(3)如果地震波荷载持续反复作用在一个方向上,第一次加载导致部分离散裂纹的产生,后续重复加载并不能使裂纹继续扩展,这与岩石的Kaiser效应是一致的。如果改变加载方向,进行不同方向地震波的反复交替加载,在地震波波形不变的情况下,也会导致初始起裂裂纹继续扩展、贯通,从而导致岩石破坏。