应力波入射滑移节理传播特性的数值模拟
2017-01-03吴新霞胡英国黎卫超
饶 宇，赵 根，吴新霞，李 鹏，胡英国，黎卫超
(长江科学院 水利部岩土力学与工程重点实验室，武汉 430010 )
应力波入射滑移节理传播特性的数值模拟
饶 宇，赵 根，吴新霞，李 鹏，胡英国，黎卫超
(长江科学院 水利部岩土力学与工程重点实验室，武汉 430010 )
基于弹性波理论，建立应力波入射滑移节理的数值模拟模型，并计算应力波的透射系数、反射系数和吸收系数，并据此研究应力波垂直入射滑移节理的传播特性。研究表明：纵波垂直入射滑移节理的传播特性主要受节理面的抗拉强度影响，当纵波产生的拉应力大于抗拉强度的瞬间，节理面产生分离而破坏，透射波瞬间消失，反射波增强，吸收能量增加，发生全反射；剪切波垂直入射滑移节理的传播特性主要受节理面的切向刚度与抗剪强度、等效切向刚度乘积的比值λ的影响，随着λ增大，反射波逐渐增强，透射波逐渐减弱，吸收能量先增大后减小，当λ足够大时，应力波在节理面处发生全反射而不产生透射。
应力波传播；滑移节理；反射系数；透射系数；吸收系数；能量；数值模拟
1 研究背景
天然岩体的不连续性，诸如包含的裂隙、软弱夹层以至节理、断层带等，使得天然岩体内的应力波传播更为复杂。一方面岩体不连续性导致应力波传播过程中发生复杂的透射、反射甚至能量损失；另一方面应力波传播也使得岩体内部发生调整甚至导致岩体参数的改变。因此，研究应力波在岩体节理处的传播规律对于判定岩体稳定性具有重要的工程意义。
目前，国内外学者对应力波在节理岩体内传播进行研究时，主要将节理分为黏结型和滑移型。对于黏结型的节理面，无论应力波以何种角度入射节理面，都不会引起节理面处岩体的相对滑移和能量损失，即在任何时间段内，节理面的入射波能流恒等于反射波能流加上透射波能流。Pyrak-Nolte等[1]、Li等[2]、杨风威等[3]、刘立波等[4]、刘婷婷等[5]在这方面都做了较为深入的研究。
自然状态下，岩体内节理多为滑移型节理。应力波倾斜入射在节理面处产生剪应力，当剪应力大于抗剪强度时，两侧节理面将产生相对滑动，从而消耗入射波能量。Li 等[6-7]考虑应力波在节理面处的能量损失，研究了节理岩体的等效黏弹性介质模型；宋林等[8]采用位移不连续理论，研究了节理面滑移状态下的透射系数、反射系数及能量消耗系数；王帅等[9]运用频域分析方法，研究了应力波在层状节理中的传播，并运用离散元程序UDEC进行验证；卢文波[10]以节理刚度作为判断节理黏结程度的指标，来区分完全黏结型节理和滑移型节理。
应力波入射滑移节理时的传播较完全黏结节理更为复杂，目前的研究侧重于将岩体节理视为弹性或半弹性材料来计算应力波的透射系数和反射系数，而对应力波入射滑移型节理的能量损失研究偏少。本文基于弹性波动理论，推导了应力波倾斜入射滑移节理时的透射系数、反射系数以及吸收系数的计算公式，并对节理不同抗拉强度和抗剪强度参数下，产生滑移的过程及机理采用数值模拟的方法进行验证，并探讨透射系数、反射系数以及吸收系数与节理面参数的关系，拓展了滑移型节理面的动态响应规律的研究。
2 滑移节理的非线性动力模型
2.1 透射系数、反射系数及吸收系数定义
应力波在斜入射不连续界面时将产生波型转换，即产生透射同类型波、透射转换波、反射同类型波和反射转换波。本文采用能量法来定义透射系数、反射系数及吸收系数，即取节理处单位周期内通过单位面积的透射波、反射波与入射波的能量比值来定义透射系数、反射系数，在透射系数、反射系数基础上定义吸收系数。定义如下：
(1)
(2)
(3)
(4)
(5)
(6)
式中：下标c=P,S，表示入射波类型；下标k=P，S，表示透射或反射波类型(其中：P表示纵波，S表示剪切波);ETk,ERk,EIc分别表示1个周期内透射、反射和入射波通过单位面积的能量;Tkc，Rkc分别为透射系数、反射系数；A为不连续界面处的能量吸收系数;σ=ρCv，表示应力波作用于节理界面的应力(其中：P波作用为法向应力，S波则为剪应力)；ρ为传播介质密度;C为应力波传播速度;v(vIc,vTk,vRk)为入射波、透射波和反射波的质点振动速度;T0为振动周期;t1为应力波到达节理界面后某一时刻(作为计算开始时刻)。
图1 模型示意图Fig.1 Schematic diagrams of numerical models
2.2 FLAC3D模型的建立
由于应力波(P波或SV波)倾斜入射不连续结构面时将产生波形转换，为了更好地区分不同类型的透射波、反射波，模型的尺寸需要足够大，以通过应力波的不同到达时间来对波型进行区分。如图1所示，取模型的尺寸为400 m×4 m×400 m，在模型中心位置设置一条横向贯穿节理面，其沿y轴旋转到与xy平面成一定角度，即为入射角。动力加载方式采用的是应力时程加载，应力波加载在模型的底部，沿z方向传播；其中S波入射时沿x方向加载，P波入射时沿z方向加载。固定模型的y方向，模型上、下边界设置为z向和x向的静态边界条件；模型械右边界S波入射时设置z方向的黏性边界，而P波入射时设置x方向的黏性边界。
为更好地区分透射和反射波形，在模型内部共布置4个速度监测点，监测点A 布置在模型左下角，监测点B布置在模型左上角，监测点C，D布置在模型内部，通过监测A，B，C，D的x方向和z方向的速度，以此来计算S波和P波的透射系数、反射系数。
2.3 模型参数的确定
为了更好地研究滑移节理对应力波传播的影响，考虑在岩体结构内设置一组滑移型节理面，而传播介质为完整的岩块，视为弹性材料，岩块的参数取值见表1。
表1 岩块物理力学参数
Table 1 Physical and mechanical parameters of rock mass
密度ρ/(kg·m-3)弹性模量En/GPa剪切模量Ea/GPa体积模量K/GPa265075．0530．0250．03
则P波和S波的传播速度如下：
(7)
(8)
节理面初始参数取为无滑移弹性节理的物理力学参数，如表2所示。
表2 无滑移弹性节理的物理力学参数
Table 2 Physical and mechanical parameters of non-slip elastic joint
法向刚度EV/GPa切向刚度EM/GPa内摩擦角φ/(°)黏聚力c/MPa抗拉强度T/MPa3．873．00402．5100
FLAC3D动力分析中，动荷载的施加方法有加速度曲线、速度曲线、位移曲线以及应力曲线4种方法。本文采用入射波幅值为0.1 m/s、频率为50 Hz的正弦函数波，由于采用Lysmer和 Kuhlemeyer[11]的黏性吸收边界，动荷载输入必须为应力曲线。应力输入应为2倍的黏滞边界法向黏滞力或切向黏滞力，则对于入射P波和入射S波等价施加的应力函数如下：
(9)
(10)
式中：σP，σS为入射波等价应力函数(MPa)；ρ为传播介质密度；CP，CS分别为P波、S波传播速度；ω为质点振动角频率；vP，vS为P波、S波的质点振动速度函数，vP=vS=0.1sin(ωt)。模型网格尺寸4 m，是入射P波波长的1/29，入射S波波长的1/17，符合网格尺寸要求[11]。
3 应力波垂直入射滑移节理的数值模拟
P波入射或者S波入射节理面时，在节理面法向和切线方向将分别产生拉应力和切应力，当拉应力超过抗拉强度或剪切力超过抗剪强度时，节理两侧的岩体介质将产生相对运动，此时即存在能量损失。本文选取节理面不同黏聚力以及抗拉强度参数，对P波和S波垂直入射滑移节理时进行数值模拟，获得相应的传播特性。
3.1 节理抗拉强度对应力波传播特性的影响
完整岩石的抗压性能良好，抗拉性能差，这就是动力荷载容易造成岩质边坡的拉裂破坏的原因。而岩体内由于存在不良结构面，使得其抗拉强度更差。当拉伸应力波入射到节理面时，在节理面的法向将产生拉应力。当法向拉应力小于节理面的法向抗拉强度时，接触面如同胶合模型，不会产生张开或滑移；当应力波产生的拉应力一部分超过节理面的法向抗拉强度，接触面分离产生破坏，从而影响到应力波的传播特性。滑移节理的抗拉强度按表3取值。
表3 滑移节理抗拉强度选取
Table 3 Tensile strength of slip joint
节理类型工况抗拉强度/MPa无滑移弹性节理100．012．0滑移节理20．530．0
图2 T=2 MPa时入射、反射和透射波形Fig.2 Incidence, reflection and transmission waveforms when T=2 MPa
拉伸P波入射时，在节理面处产生应力方向与P波振动方向一致，即垂直入射时，应力方向为节理面的外(内)法线方向，根据式(9)，可得其应力大小为
(11)
由式(11)可知应力峰值σmax=1.545 MPa。
图2—图4分别为抗拉强度为2，0.5，0 MPa时，测点C(即入射波形和反射波形)，D(即透射波形)速度时程曲线。
图3 T=0.5 MPa时入射、反射和透射波形Fig.3 Incidence, reflection and transmission waveforms when T=0.5 MPa
图4 T=0 MPa时入射、反射和透射波形Fig.4 Incidence, reflection and transmission waveforms when T=0 MPa
由图2—图4可知，当应力波沿底面内法向加载时，前半周期在节理面上产生压应力，此时节理面抗拉强度大小不影响应力波传播，后半个周期应力波在节理面产生拉应力；当应力波沿底面外法向加载时，前半周期在节理面上产生拉应力。
无论应力波沿内法向加载还是沿外法向加载，当抗拉强度T>应力波应力峰值σmax时，抗拉强度不影响应力波的传播，应力波在节理面处正常透射和反射；当抗拉强度T<应力峰值σmax时，在应力波产生拉应力的半周期内(内法向加载的后半周期、外法向加载的前半周期)，抗拉强度将影响应力波传播。
当应力波在节理面处产生拉应力σ
3.2 节理剪切参数对应力波传播的影响
垂直入射时，剪切波在节理面上产生切向应力。因此，节理面的抗剪强度控制着剪切波垂直入射节理面的传播。本节主要针对剪切波垂直入射滑移节理的传播特性进行探讨。定义剪切应力波入射时节理的等效切向刚度Ks为
Ks= ks/(ωzs) 。
(12)
式中：ks为节理切向刚度；zs为剪切波的波阻抗，zs=ρCS。则入射剪切波在节理面上产生的最大剪应力为τ=ksU，U为应力波入射产生的最大位移值。
定义无量纲常数λ为
(13)
式中τs为节理面的抗剪强度。
采用摩尔-库伦强度理论的抗剪强度公式，即
τs=σtanφ+c 。
(14)
式中：σ为节理面的法向应力;φ为内摩擦角;c为节理面的黏结强度。
为了更好地了解透射系数、反射系数和吸收系数与无量纲常数λ的关系，节理的初始物理力学参数可按照表2取值，考虑内摩擦角φ=0°，节理的抗剪强度全部由节理面的黏结强度c提供。采用式(10)所示的剪切应力波入射，则剪切S波入射产生的剪应力为
(15)
由式(15)可知剪应力峰值τmax=0.892 MPa。
为了研究节理面抗剪强度对应力波传播特性的影响，节理面抗剪强度取值如表4所示，图5为抗剪强度分别取表4所示值时的入射、透射和反射波形。
表4 滑移节理抗剪强度选取
Table 4 Shear strength of slip joint
节理类型工况抗剪强度/MPa无滑移弹性节理2．5011．00 滑移节理20．2030．0540．00
图5 不同抗剪强度下入射、反射和透射波形Fig.5 Incidence, reflection and transmission waveforms at different shear strength
图5表明，当入射波剪切应力峰值τmax<节理面抗剪强度τs时，应力波能很好地通过节理面，产生透射和反射，而无能量的损失，此时节理面的抗剪强度不影响应力波的透射和反射；当τmax>τs时，节理面产生滑移和能量的损失，随着抗剪强度τs减小，反射波逐渐增强，透射波逐渐减弱。当抗剪强度τs减小到0时，节理面被破坏而如同自由面，应力波在节理面处发生全反射而不产生透射。
如图1所示，通过监测点C，D的速度，采用式(1)—式(6)来计算剪切应力波入射滑移节理面时透射系数、反射系数及吸收系数。根据设置节理面的不同抗剪强度和切向刚度获得不同的无量纲常数λ，研究无量纲常数λ与剪切应力波的透射系数、反射系数和吸收系数的关系。图6分别为节理等效切向刚度Ks=0.5和Ks=5时无量纲常数λ与剪切应力波的透射系数、反射系数和吸收系数的关系。
图6 不同节理等效切向刚度下反射系数、透射系数和吸收系数Fig.6 Reflection, transmission and absorption coefficients at different normalized tangential joint stiffness
由图6可知，主要影响剪切波在节理处透反射的因素是无量纲常量λ，垂直入射剪切应力波在节理处的传播主要与剪切波的波阻抗zs、入射剪切波的角频率ω、应力波入射产生的最大位移值U以及节理面的抗剪强度τs有关。
当无量纲常量λ≤1时，节理面上不会产生滑移，从而没有能量损失，随无量纲常量λ增大，节理面的剪应力大于抗剪强度，能量吸收系数先增大后减小，反射系数逐渐增大至1，透射系数逐渐减小到零。即节理切向刚度ks→∞或节理抗剪强度τs→∞时，应力波能很好地透射节理而反射波很小，节理近似于完整岩体而对应力波传播无影响；节理切向刚度ks→0或节理抗剪强度τs→0时(节理面被破坏)，应力波几乎全部被反射而透射波很小，节理面如同自由界面而发生全反射。
4 结 论
本文建立应力波入射滑移节理的FLAC3D模型，研究表明纵波垂直入射滑移节理的传播特性主要受节理面的抗拉强度影响，剪切波垂直入射滑移节理的传播特性主要受节理面的抗剪强度和切向刚度的影响，并获得如下结论：
(1) 纵波垂直入射，当产生的峰值应力小于抗拉强度，应力波很好地通过节理面产生透射和反射，而无能量损失；当纵波产生的拉应力大于抗拉强度的瞬间，节理面产生分离而破坏，透射波瞬间消失，反射波增强，吸收能量增加，节理面如同自由界面而发生全反射。
(2) 剪切波垂直入射时，当剪切波在节理面产生的剪应力小于抗剪强度，应力波能很好地通过节理面，产生透射和反射，而无能量的损失；当剪应力大于抗剪强度时，节理面产生滑移，随着抗剪强度τs减小，反射波逐渐增强，透射波逐渐减弱，当抗剪强度为0时，应力波被全部反射而无透射波。
(3) 节理面的透射系数、反射系数及吸收系数主要受无量纲常量λ影响，随λ增大，反射系数逐渐增大至1，透射系数逐渐减小到0，吸收系数先增大后减小；当节理切向刚度或节理抗剪强度足够大时，应力波能很好地透射节理而反射波很小，节理近似于完整岩体而对应力波传播无影响；节理切向刚度或节理抗剪强度足够小时，应力波几乎全部被反射而透射波很小，节理面如同自由界面而发生全反射。
[1] PYRAK-NOLTE L J, MYER L R, COOK N G W. Transmission of Seismic Waves Across Single Natural Fractures[J]. Journal of Geophysical Research, 1990, 95(B6):8617-8638.
[2] LI J, MA G. Analysis of Blast Wave Interaction with a Rock Joint[J]. Rock Mechanics and Rock Engineering, 2010, 43(6): 777-787.
[3] 杨风威, 李海波, 李建春, 等. 斜入射线弹性节理应力波传播特征的数值模拟[J]. 岩土力学, 2013, 34(3):901-907.
[4] 刘立波, 李建春, 李海波, 等. 应力波斜入射黏弹性节理的传播规律[J]. 岩石力学与工程学报, 2012, 31(增2): 3593-3598.
[5] 刘婷婷, 李建春, 李海波, 等. 应力波通过非线性平行节理的能量分析[J]. 岩石力学与工程学报, 2013, 32(8): 1610-1617.
[6] LI J C, MA G W, ZHAO J. An Equivalent Viscoelastic Model for Rock Mass with Parallel Joints[J]. Journal of Geophysical Research: Solid Earth, 2010, 115(B3).DOI:10.1029/2008JB006241.
[7] 李建春, 李海波, MA G W, 等. 节理岩体的一维动态等效连续介质模型的研究[J]. 岩石力学与工程学报, 2010, 29(增2): 4063-4067.
[8] 宋 林, 邵珠山, 吴敏哲,等. 考虑剪切滑移效应时节理岩体中剪切S波的传播特性探析[J]. 西安建筑科技大学学报(自然科学版), 2011, 43(4):494-500.
[9] 王 帅, 盛 谦, 朱泽奇, 等. 频域分析方法在层状节理岩体波动问题中的应用[J]. 长江科学院院报, 2013, 30(4): 62-66.
[10]卢文波. 应力波与可滑移岩石界面间的相互作用研究[J]. 岩土力学, 1996, 17(3):70-75.
[11]LYSMER J, KUHLEMEYER R L. Finite Dynamic Model for Infinite Media[J]. Journal of the Engineering Mechanics Division, 1969,95(4):859-877.
(编辑：罗 娟)
Numerical Simulation of Propagation Characteristics ofStress Waves across Slip Joint
RAO Yu, ZHAO Gen, WU Xin-xia, LI Peng, HU Ying-guo, LI Wei-chao
(Key Laboratory of Geotechnical Mechanics and Engineering of Ministry of Water Resources,Yangtze River Scientific Research Institute, Wuhan 430010, China)
Based on the elastic wave theory, the numerical simulation model of stress waves across slipping joint has been built， and transmission, reflection and absorption coefficients are calculated to analyze the propagation characteristics of stress waves. Results show that propagation characteristics of normal incident longitudinal wave are affected by joint’s tensile strength, and if the tensile stress produced on the joint by normal incident longitudinal wave exceeds the tensile strength, the interface will separate, transmitted wave disappears instantaneously, and then total reflection occurs, with the reflected wave and absorbed energy increasing. The propagation characteristics of normal incident shear wave are affected by the ratio of tangential joint stiffness and the product of shear strength and normalized tangential joint stiffness. With the increasing of the ratio, the reflection coefficient increases, transmission coefficient decreases, and absorption coefficient increases at first and then decreases, then ends with total reflection.
stress waves propagation; slip joint; reflection coefficient; transmission coefficient; absorption coefficient; energy; numerical simulation
2015-09-17
国家自然科学基金青年科学基金项目(51309026)
饶 宇(1990-)，男，湖南长沙人，助理工程师，硕士，主要从事岩土动力学方面的研究工作，(电话)027-82829793(电子信箱)raoyuray@163.com。
10.11988/ckyyb.20150795
2016,33(12):94-98
TU45
A
1001-5485(2016)12-0094-05
