王 举，杨风威，唐 旭

(1.贵州路桥集团有限公司， 贵州 贵阳 550081； 2.黄河勘测规划设计有限公司， 河南 郑州 450003；3.贵州省质安交通工程监控检测有限责任公司， 贵州 贵阳 550081)

地震荷载作用下含节理岩质边坡的破坏机理

王 举1，杨风威2，唐 旭3

(1.贵州路桥集团有限公司， 贵州 贵阳 550081； 2.黄河勘测规划设计有限公司， 河南 郑州 450003；3.贵州省质安交通工程监控检测有限责任公司， 贵州 贵阳 550081)

以汶川地震为背景,采用UDEC数值模拟，以什邡八角镇实测的“5.12”汶川地震波作为原始波形，对地震荷载作用下，含一条节理面岩质边坡滑移、拉裂破坏的过程进行了研究。结果表明：岩石边坡破坏模式为块体沿节理面的滑移破坏，并伴随着上方岩体拉裂破坏；节理刚度主要影响边坡前期相对位移大小及塑性区的产生，其后期拉裂区的扩展模式基本是一致的；节理面倾角较小时，边坡岩体仅在节理与坡顶、坡面交叉区域产生小范围的拉裂破坏，随着倾角的增大，边坡的相对位移及拉裂塑性区都显著增大。

UDEC；地震荷载；岩质边坡；岩石强度；节理刚度；节理倾角

岩石边坡中不可避免的存在节理、裂隙、夹层等不连续结构面。地震波传播到这些不连续结构面时会发生复杂的相互作用，一方面，相比于完整岩石的强度，岩体结构面的强度低很多，因此地震波传播会引起岩体结构面的张开、扩展现象，造成岩体的失稳或破坏，这为人们预测、预防岩体的失稳提供了理论基础；另一方面，各类结构面的存在会改变地震波传播的方向、幅值和频率，并使波产生反射和折射现象，阻隔波的传播，造成地震波的衰减和波形改变，这也成为了解决防灾减灾和岩土体工程稳定问题时最为关心的问题之一[1]。王帅等[2-4]人研究非贯通节理岩体在地震作用下的动力响应特征，结果表明节理密度、倾角等相关参数的变化都会对岩体在地震作用下的动力响应产生显著的影响。由于节理岩体的动力学试验，试样制备复杂，设备价格昂贵，且解析法只能求解部分不太复杂的问题，因此，数值模拟成为研究岩体动力荷载响应的一项重要工具[5]。陶连金等[6]对节理岩体边坡的动力相应进行了离散元模拟，并推倒了离散元的计算公式。刘亚群等[7]运用二维离散元程序UDEC对某以节理岩质边坡在爆破时的动力响应进行了分析，得到的斜坡底部的速度时程与观测结果吻合的较好。应用数值模拟对岩质边坡在地震荷载下动态响应的研究并不少见，边坡几何特征、荷载对响应的影响亦有分析，但对含节理的岩质边坡、节理面的刚度、节理面倾角等因素对边坡动态响应的影响则是很少涉及[8-12]。

鉴于UDEC进行应力波穿过节理面传播规律分析的良好特性，本文以汶川地震为背景，采用UDEC数值模拟，为突出节理面性质的影响，对地震波作用下，含节理面不同节理刚度及不同节理倾角的岩质边坡滑移、拉裂破坏的过程进行研究，对其他因素进行了适当的简化，不考虑水的作用。为进一步分析震裂变形边坡的稳定性及其加固治理提供铺垫。

离散元是基于块体之间的非连续性假设，它把岩体表示为离散块体的组合，把节理表示为块体之间的接触，即把节理处理为模型中的一个边界条件而不是一个特殊的单元。每一个块体与周围的块体通过可以随时间发生变化的边界接触而发生作用。与一般的解连续介质力学问题要求满足边界条件、平衡方程、变形协调方程和本构方程不同，离散单元法不需要满足变形协调方程，因为介质一开始就假定为离散块体的集合，故块与块之间没有变形协调的约束[13-15]。

1 计算模型建立及参数选取

1.1 计算模型

为了重点分析节理面性质及荷载参数对边坡动态响应的影响，将岩质边坡的坡高和坡角固定。参考肖克强[8]的工作，边坡坡高取100 m，坡角取60°。边坡底边长取600 m以减小侧边双向黏性边界对输入应力波传播精度的影响。

在边坡中自坡顶至坡面设置一条贯通的节理面，并在节理上方岩体区域，沿坡面和坡顶布设速度监测点A-J。坡面各点的高程差为20 m，坡顶各点水平距离也为20 m。在后续的影响因素分析中，坡面各点以其距离坡角的高程来描述，坡顶各点以其距边坡顶点的水平距离来描述。图1为计算模型示意图，计算中通过改变节理与坡顶交叉点的位置改变节理倾角，节理面倾角β取值范围为0°～50°。

1.2 岩石力学模型和物理力学参数

计算中岩石材料采用理想弹塑性模型，屈服准则采用Mohr-Coulomb强度准则，屈服函数如下[9]：

图1计算模型及测点布置示意图(单位:m)

(1)

ft=σ3-σt

(2)

式中：σ1、σ3分别为最大、最小主应力；φ为摩擦角；c为黏聚力；σt为岩石抗拉强度，其中

(3)

当岩体内某一点应力满足fs<0时，发生剪切破坏；当满足ft>0时，发生拉伸破坏。本次数值模拟岩石材料，其物理力学参数见表1。

表1 岩石物理力学参数

1.3 节理力学模型和力学参数

岩体节理力学模型采用库仑滑移模型，它是UDEC中的主要节理模型，其力学特性如图2所示。根据这一模型，在节理的法向方向，应力-位移关系可以看作是线性的，由节理法向刚度控制，如下式：

Δσn=-knΔun

(4)

式中：kn为节理法向刚度；Δσn为有效法向应力增量；Δun为法向位移增量。通过块体之间的叠合量的计算可以很方便的得到节理面的相对法向位移。

图2 UDEC中节理力学模型

对于节理受拉的情况，如果法向应力σn大于其抗拉强度T(即σn<-T)，则σn=0。

与此类似，对于节理受剪的情况，其位移是由节理的剪切刚度ks控制的，并且节理的剪应力τs不能大于其抗剪强度τmax，即

当|τs|≤c+σntanφ=τmax时，

(5)

当|τs|≥τmax时，

τs=sign(Δus)τmax

(6)

表2为节理面强度参数。

表2 节理面力学参数

1.4 网格尺寸、阻尼及边界条件

本文数值模拟中，输入荷载的频率最大为20 Hz，根据岩石的物理力学参数以及式Kuhlemeyer和Lysmer建议的网格单元尺寸应不大于最小波长的1/8～1/10，网格尺寸设置为4 m，可以满足精度要求。考虑岩体介质阻尼的影响，采用局部阻尼形式，局部阻尼系数取0.125。

由于输入的地震波荷载具有较长的持时，因此为了避免地震波在侧边界反射效应的影响，模型左右边界均设置x、y双向黏滞边界，同时通过试算，确定模型纵横尺寸比为1∶4，以保证输入荷载的精确性，见图1。

2 地震波选取及基线校正

本次数值模拟以什邡八角镇实测的“5.12”汶川地震波作为原始波形，选取主震段附近约12 s的振动波形，并按比例对其幅值进行缩放[10-11]。为保证数值计算结果的可靠性，不会因速度、位移的数值过大而使得计算结果难以解释甚至得出错误的结果。对截取得到的加速度波形进行基线校正。本文按比例将截取的地震加速度波峰值调整为0.6g，并采用文献[12]提出的均方加速度法对加速度时程曲线进行基线校正，得到模型输入所需的竖直向及水平向振动加速曲线见图3。

3 地震荷载下含节理岩质边坡的破坏机理

3.1 边坡的破坏模式

图4为地震波作用下，节理上方坡体内塑性区扩展过程。图中圈号代表此处网格单元发生张拉破坏，星号代表此处网格单元恰好处于屈服状态，叉号代表此处网格单元已经屈服。图5为t=1.81 s及t=2.2 s时刻边坡最小主应力分布情况。

图3输入模型振动加速度波形

从数值模拟塑性区破坏的扩展过程来看，对于含有一条节理面的岩质边坡，在地震荷载作用下，其破坏模式可以分为沿节理面滑移、滑移和上部岩体拉裂两个阶段。

在t=1.81 s之前，地震荷载的作用使得节理面上的滑动力大于抗滑力，但尚未达到岩石材料破坏的强度，因此，边坡上方岩体仅沿节理面产生滑动。随着地震荷载的增大，边坡岩体在节理面中上部附近开始出现拉裂破坏，并呈现向坡顶扩展的趋势，如图4(a)所示。此时对应的边坡最小主应力分布如图5(a)所示，边坡在节理面中间区域及节理与坡顶交叉区域出现明显的拉应力区。

另外，拉裂塑性区不仅是从节理面向坡顶扩展，也会出现由坡顶向节理面扩展，如图4(c)所示，但仍以前者为主。在t=2.2 s时刻，边坡内部已出现明显的4条竖向拉裂区，如图4(d)所示，同时，贯通塑性区的出现也引起边坡应力的重分布，节理上方岩体基本都处于受拉状态，如图5(b)所示。图6为此时刻边坡位移等值线图，由图6可见，随着坡体沿节理面的滑移，节理面上方岩体的位移值并不一致，表明拉破坏塑性区的出现致使边坡位移产生不连续。

图4 地震波作用下边坡塑性区扩展模式

图5 地震波作用下边坡最小主应力分布图

图6边坡位移等值线图(t=2.2 s)

随着地震波的持续输入，节理上方岩体中拉裂、剪切塑性区不断增多，到计算历时结束，已基本完全变为塑性区。沿节理面自上而下选取5个监测点A1—A5，其相对于下方岩体的切向相对位移-时间曲线如图7所示。

图7节理面各点切向相对位移-时间曲线

由图7可知，各点的切向位移曲线在1.9 s附近发生分离，对应上方岩体开始出现拉裂塑性区。另外，从整个历时区间可以看出，岩体沿节理面呈现逐级滑动，与输入加速度波形的几个较大峰值时刻相对应。

值得一提的是，这里仅就边坡在地震作用下发生滑移及拉裂塑性区扩展的过程进行了分析。虽然边坡发生了几级相对滑移后，相对位移继续增大并在地震波持时结束后增幅逐渐减小而趋于稳定，但此时边坡有可能失稳破坏，也可能处于稳定状态，因为节理面在经历了较大范围的相对滑动，强度会显著降低，而且由于岩体中大量塑性区的产生，其受力情况会发生改变，需要更进一步的研究加以判断。

3.2 不同节理刚度边坡的破坏模式

分别取节理刚度为0.1 GPa/m、1.0 GPa/m，其他参数保持不变，分析不同节理刚度对边坡破坏模式的影响。边坡在t=1.8 s、t=2.2 s两个时刻塑性区扩展及岩体相对位移情况见图8。由图8可知，t=1.8 s时刻，节理刚度0.1 GPa/m的边坡已出现节理面至坡面的贯通的拉裂破坏区，且坡顶也出现向节理面扩展的拉破坏单元，而此时节理刚度1.0 GPa/m的边坡仅在节理面中上部开始形成拉破坏单元。t=2.2 s时刻，两者拉破坏单元数已基本持平，节理的相对位移略有差异。因此节理刚度主要影响边坡前期相对位移和塑性区的产生，其后期拉裂区的扩展模式基本是一致的。

3.3 不同节理倾角边坡的破坏模式

取节理面刚度kn=ks=1 GPa/m，其他参数保持不变，分析不同节理倾角边坡在地震荷载下的破坏模式。t=2.0 s及t=6.0 s时边坡的塑性区分布如图9、图10所示。

图8地震荷载下不同节理刚度边坡塑性区扩展

由图9、图10可看出t=2.0 s时刻，β=20°边坡节理上方岩体仅在节理与坡顶交叉区域产生了小范围的拉裂破坏，而β=30°边坡节理上方岩体内已形成多条贯通的竖向拉裂区，β=40°及β=50°的边坡已发生明显的滑移，且节理上方岩体内出现大范围的塑性区。t=6.0 s时刻，β=20°的边坡拉裂区出现在节理与坡顶、坡面交叉处，β=30°边坡节理上方岩体内已形成大量塑性区，而β=40°及β=50°的边坡已发生大规模的滑动，最大相对位移分别为3.021 m、10.27 m。另外，节理倾角较小时，边坡在节理顶部向边坡底部会形成一条竖向塑性区，如图10(a)、图10(b)所示。这是由于SV波入射时，节理倾角较小，其反射应力波在节理下方会充分叠加从而产生破坏。

图9 不同节理倾角边坡塑性区(t=2.0 s)

图10不同节理倾角边坡塑性区(t=6.0 s)

4 结 论

本文采用UDEC数值模拟，对地震荷载作用下，含一条节理面岩质边坡滑移、拉裂破坏的过程进行了研究，得到以下主要结论：

(1) 岩石边坡破坏模式为块体沿节理面的滑移破坏，并伴随着上方岩体拉裂破坏。

(2) 节理刚度主要影响边坡前期相对位移大小及塑性区的产生，其后期拉裂区的扩展模式基本是一致的。

(3) 节理面倾角较小时，边坡岩体仅在节理与坡顶、坡面交叉区域产生小范围的拉裂破坏，随着倾角的增大，边坡的相对位移及拉裂塑性区都显著增大。

[1] Goodman R E. Methods of Geological Engineering in Discontinuous Rock[M]. New York: John Wiley and Sons, 1976:135-140.

[2] 王 帅,盛 谦,朱泽奇,等.地震荷载作用下地下洞室不利结构坍塌机制研究[J].岩土力学,2012,33(10):2897-2902,2966.

[3] 信春雷,高 波,周佳媚,等.含节理隧道洞口地震动力放大效应研究[J].木工程学报,2013,46(S1):269-274.

[4] 崔 臻,冷先伦,朱泽奇,等.非贯通随机节理岩体的结构效应对地下洞室地震响应的影响[J].岩土力学,2014,35(S2):654-652.

[5] Jing L. A review of techniques,advances and outstanding issues in numerical modelling for rock mechanics and rock engineering[J]. International Journal of Rock Mechanics and Mining Sciences, 2003,40(3):283-353.

[6] 陶连金,苏生瑞,张倬元.节理岩体边坡的动力稳定性分析[J].工程地质学报,2001,9(1):32-38.

[7] 刘亚群,李海波,李俊如,等.爆破荷载作用下黄麦岭磷矿岩质边坡动态响应的UDEC模拟研究[J].岩石力学与工程学报,2004,23(21):3659-3663.

[8] 肖克强.地震荷载作用下顺层岩体边坡变形特征及稳定性研究[D].武汉:中国科学院研究生院(武汉岩土力学研究所),2006.

[9] Itasca Consulting Group, Inc. UDEC(Universal Distinct Element Code) User's Manual Version 4.0[M]. Minneapolis:Itasca Consulting Group, Inc., 2005.

[10] 于海英,王 栋,杨永强,等.汶川8.0级地震震动特征初步分析[J].震灾防御技术,2008,3(4):322-336.

[11] 殷跃平.汶川八级地震滑坡特征分析[J].工程地质学报,2009,17(1):29-38.

[12] Yang J, Li J B, Lin G. A simple appproch to intergration of acceleration data for dynamic soil structure interaction analysis[J]. Soil Dynamics and Earthquake Engineering, 2006,26(8):725-734.

[13] Eberhardt E, Thuro K, Luginbuehl M. Slope instability mechanisms in dipping interbedded conglomerates and weathered marls—the 1999 Rufi landslide, Switzerland[J]. Engineering Geology, 2005,77(1/2):35-56.

[14] Sepulveda S A, Murphy W, Jibson R W, et al. Seismically induced rock slope failures resulting from topographic amplification of strong ground motions: the case of pacoima canyon, california[J]. Engineering Geology, 2005,80(3/4):336-348.

[15] Itasca Consulting Group, Inc. 3DEC (3 Dimensional Distinct Element Code) User’s Manual Version 3.0[M]. Minneapolis: Itasca Consulting Group, Inc., 2003.

FailureMechanismofRockSlopewithJointUndertheSeismicLoad

WANG Ju1, YANG Fengwei2, TANG Xu3

(1.GuizhouRoad&BridgeGroupCo.,Ltd.,Guiyang,Guizhou550081,China;2.YellowRiverEngineeringConsultingCo.,Ltd.,Zhengzhou,He'nan450003,China; 3.GuizhouProvinceQualityandSafetyTrafficEngineeringMonitoringandInspectionCenterCo.,Ltd.,Guiyang,Guizhou550081,China)

“5.12” Wenchuan earthquake waves measured at the site of Shifang Bajiaozhen were taken as the original waveform to study the single joint surface rock slope sliding and the process of tension damage under the seismic load by using the UDEC numerical simulation. The results shows that the rock slope failure mode is simple slip failure of rock block along joint surface which is accompanied by tension damage of the above rock; the joint stiffness mainly affects the slope relative displacement and formation of plastic zone in earlier stage, the expansion patterns of the tension fracture zone are basically consistent in later stage; when the joint surface angle is small, a small scale tension damage will occur in the cross region of joint, top and surface of slope rock, as the angle increases, the slope relative displacement and tension fracture plastic zone increase significantly.

UDEC;seismicload;rockslope;rockstrength;jointstiffness;jointangle

10.3969/j.issn.1672-1144.2017.05.031

2017-04-27

2017-06-05

贵州省交通运输厅科技项目(2015122046)

王 举(1984—)，男，贵州毕节人，工程师，主要从事公路与桥梁施工工作。E-mail: 373636753@qq.com

TU45

A

1672—1144(2017)05—0171—07