Dynamic shear diffusion behavior at rock interfaces
-
摘要: 动载荷下剪切失稳控制的扩散行为是岩石局部大变形发展和宏观力学性能劣化的诱因。基于广义变分原理建立剪切载荷作用下界面动态失稳的力学模型,得到了关于界面失稳的判别式和扩散方程。基于判别方程得到了剪切力和动力效应对失稳界面角度的影响,结果表明:随着外部剪切作用力的增大,剪切变形带角度有一定程度的增大;随着局部动力系数的增大,即局部惯性作用力的增大,剪切带角度明显减小。结合本征位移求解扩散方程,初步得到其位移解析表达式,位移随加载时间的增加逐渐增大。为了验证理论模型的可靠性,并进一步研究界面失稳的变形行为和对波传播的影响,建立了数值分析模型。分析结果表明:界面失稳为局部剪切破坏滑移的先导条件;界面厚度和剪切力越大,局部位移越大;界面剪切扩散行为极大降低了透射波的幅值,同时也改变了透射波的频率。研究结果可为岩石局部化变形、岩石动态强度等研究提供理论参考。Abstract: The diffusion behavior of shear instability control under dynamic load is the inducement for the development of local large deformation and the deterioration of macroscopic mechanical properties of rock. Firstly, the energy function of the unstable interface was established; and then based on the generalized variational principle, the interface disturbance analysis was carried out, while the first and second-order variances of the function were taken into consideration. Thus the governing equation of dynamic instability of the interface under shear load was established. Based on the discriminant equation, the influence of shear force and dynamic effect on the angle of the unstable interface is obtained. The results show that the angle of the shear deformation zone increases to a certain extent with the increase of external shear force. With the increase of the local dynamic coefficient, that is, the increase of the local inertial force, the shear band angle decreases obviously. By solving the diffusion equation with the edge displacement, the analytical expression of displacement was obtained, showing that the displacement increases gradually with the increase of loading time. To verify the reliability of the theoretical model and further study the deformation behavior of interface instability and its influence on wave propagation, the evolution of the fine and microscopic morphology of the contact surface during dynamic shear was described by combining with the SHPB rod-beam experiment technique, and an evaluation method for the influence of the evolution of surface contact parameters on mechanical parameters during interface instability was proposed. The numerical analysis model was established, and its result shows that interface instability is the leading condition of local shear failure slip. With the increase of interface thickness and shear force, the local displacement increases. The interfacial shear diffusion behavior greatly reduces the amplitude and changes the frequency of the transmitted wave. This study provides a good theoretical reference for the study of localization deformation and dynamic strength of rocks.
-
Key words:
- rock interface /
- finite deformation /
- shear localization /
- inertia effect /
- dynamic instability
-
构造地震很少由新剪切断裂的产生而传播,而是沿着已有断层边界面突然滑动扩散,地震是板块粘滑摩擦失稳的结果。此过程中脆性断裂只起次要作用,起主要作用的是断裂面和摩擦磨损的交替发展[1-3]。地震过程由表征断层短暂运动状态的“滑”和表征断层长期准静状态的“粘”组成,后者是地震事件间弹性应变积累的时期。地震过程中岩石界面动摩擦的变形机理主要表现为[4]在临界温度之上的晶体塑性变形和碎裂流动。虽然都显示出较强的塑性特征,但两种变形有不同的流变学机理。从动力学角度,岩石界面动摩擦过程大致可以分为波致摩擦、惯性摩擦以及滑移三个阶段[5-7]。本文中进行岩石界面动摩擦过程中剪切失稳行为的研究,主要关注与局部大变形相关的地震事件的“粘”这一过程的形成和发展,贯穿于波致摩擦和惯性摩擦两个阶段,但仍属于界面动摩擦的初期。相关研究可增强对孕震过程的认识。
对岩石界面动摩擦性能的实验研究较多。在低速冲击下,Di Toro等[8]采用一种伺服控制压-扭装置进行了一系列的快速滑动实验;Rice[9]基于旋转圆盘研究了动摩擦因数与滑动速度之间的关系,发现在滑动速度为4 m/s时界面出现了零摩擦强度的现象;张磊等[10-11]基于SHPB(split Hopkinson pressure bar)杆束技术研究了含斜节理面岩石的动摩擦行为;Zhou等[12]基于SHPB斜截面垫块技术研究了岩石的压剪动力学特性,发现在此过程中,岩石界面滑移速度可以控制在1 m/s量级,相对滑移位移可以控制在毫米量级。Okada等[13]采用平板斜撞击动摩擦实验,得到了压剪联合高速冲击下界面滑移时的临界滑动速度;Xu等[14]利用轻气炮压剪联合冲击加载装置对具有斜节理的花岗岩进行了界面动滑移特性研究,此过程的岩石界面滑移速度可以控制在10 m/s量级,相对滑移位移可以控制在0.1 mm量级。当岩石界面发生动态滑移时,岩石界面细微观结构必然发生不可逆变化,如局部损伤、局部断裂等,会产生一定频率的声发射信息和红外辐射信息,对这些信息的追踪分析可得到局部结构发生变化的位置和规模。因此,研究这些物理信息可以推测伴随岩石界面动摩擦过程发展的物理现象,主要包括:(1) 波致摩擦过程中界面Ⅱ型裂纹表面的亚瑞利波、超剪切波的传播[15-18],粗糙界面上跨声速的波、亚瑞利波,以及一种慢波的传播[19-20],这两组波有较大的差异,但在各自物理现象的描述中均有较好的效果,其机制目前尚不明确;(2) 波致摩擦过程中的红外温升和声发射信息[6,21-23],前者反映岩石界面滑移破裂时的规模,后者可以对破裂区域进行定位。岩石界面动摩擦是一个典型的失稳过程,界面作用非常复杂。Huang等[24]采用多尺度PFC离散元方法跟踪了岩石界面动摩擦因数的演化;Gao等[25]采用有限元与离散元结合的方法分析了岩石界面的破碎特性。以上研究关注的都是界面破碎带来动摩擦因数的降低及相关机理,但是此过程是如何形成的,则需要探讨岩石界面动态变形发展与界面失稳的联系,尤其是与剪切失稳扩散的关系。
对岩石失稳现象的研究已有较多[26]。基于线弹性和弹塑性小应变的理论一般不会产生本构失稳现象[27]。产生失稳现象的情况主要基于下述三种理论:(1)非关联流动理论[28],此理论认为体积应变与偏量应变具有不同的流动规律,假定界面失稳,应用弹塑性理论进行分析,由此研究平面应变条件下失稳特性与岩石扩容角和内摩擦角的关系,此理论可以较好地应用于脆性剪切带分析;(2)应变梯度塑性理论[29],此理论认为材料力学行为不仅与应变、应变历史有关,而且与应变梯度有关,由此采用屈服函数的应变梯度研究岩石的局部化变形失稳问题;(3)有限变形理论[30-33],此理论采用有限变形表述,导出一个简洁的四阶多项式判别方程,以确定产生局部化变形失稳的临界应力状态,结合有限群表征,可以得到平面压缩[34]和剪切[35]条件下岩石失稳发展过程,结合扩散型流函数分析,进一步则可以分析失稳时速率扰动的局部分布形状[36],主要适用于中等或较低强度岩石在围压作用下剪切带局部化变形失稳分析。
上述失稳分析主要基于准静态加载情况,产生剪切变形带等特征失稳变形的边界条件仍采用Hill等[30]的零应力边界,这对于动力加载情况不再适合。在岩石界面动态剪切变形的过程中,局部存在复杂的变形特征[11]。因此,本文中在有限变形理论的基础上,引入惯性力和剪切边界条件的影响,研究岩石动摩擦过程中的剪切扩散行为,结合实验观测和有限元分析探讨岩石界面动态剪切失稳的变形特征,以及其对波动传播的影响。
1. 岩石动态剪切失稳理论分析
1.1 基于有限变形的动态剪切失稳控制方程
在体积为V(t)、边界为∂V的岩石界面系统中,总能量为
ˉE 的位移表达式为[37-38]:ˉE(ui)=Eint+Ekin+Eext (1) 式中:Eint为系统内能,Ekin为系统动能,Eext为系统外力势能,其表达式分别为:
Eint=∫W(εij)dV (2) Ekin=∫ρv2idV (3) Eext=−∫VρbiuidV−∫∂VfiuidΓ (4) 式中:
W(εij) 为应变能密度,vi 为系统变形过程中的速度,bi 为系统的体力分量,ui为系统变形过程中的位移,fi 为作用于系统边界的面力分量,i=1, 2。一般情况下,忽略体力分量bi 。图1所示为剪切加载下,岩石局部剪切变形带发展示意图,图1(b)中n为剪切失稳面的外法线方向,α为x1轴与n方向的夹角,Q代表局部失稳位置。图1(b)中ρan 为法向惯性力,ρaτ 为切向惯性力。由于岩石微结构的影响,局部变形不会沿着加载路径,而是会产生分叉变形,形成剪切变形带。应变能密度W的表达式为:
W=12Lijklεijεkl (5) 式中:
Lijkl 为刚度张量,相应的本构关系满足σij=Lijklεkl 。εij 为拉格朗日应变的有限变形表达,具体形式为:εij=12(∂ui∂xj+∂uj∂xi+∂uk∂xi∂uk∂xj) (6) 对式(1)取一阶变分可以得到系统的运动方程,即:
δˉE=δ∫VWεdV=∫VσijδεijdV−∫∂VfiδuidΓ+∫Vρ∂2ui∂t2δuidV (7) 在稳定系统中,总势能不但满足式(7),而且它的二阶变分应大于零[39],即:
δ2ˉE=[ˉE,uu(u(Λ))δu]δu>0 式中:
Λ 为载荷参数。当系统达到失稳临界状态时,载荷参数
Λ 达到临界值oΛ ,位移场u(oΛ) 和对应的应力场σ(oΛ) 均达到临界状态,此时:δ2ˉE=[ˉE,uu(u(oΛ))δu]δu=0 (8) 将位移变分
δu 记为给定特征位移模式的Δu ,考虑到应力场σ(0)ij=Lijklεkl (上角标(0)表示远场),则由系统失稳的临界条件可得到系统的稳定性函数为:S(Λ)=∫V[LijklΔεklδεij+σ(0)ijΔδεij]dV−∫∂VfiδΔuidΓ−∫Vρ∂2ui∂t2δΔuidV (9) 式中:应变变分的表达式为:
δεij=12δ(∂u0i∂xj+∂u0j∂xi+∂u0k∂xi∂u0k∂xj)=(δui,j+u0k,iδuk,j)Q (10) 式中:上角标0表示临界位移,Q代表局部失稳位置。
系统局部剪切带状失稳分析如图1所示。远场应力
σ(0)ij 作用下产生局部剪切变形区,Hill等[30]采用零应力条件来处理此区的边界,但是,在动力学过程中,在局部剪切变形区的边界显然存在惯性力。因此,本文中的分析将考虑边界惯性力的影响,如图1(c)所示。基于上述分析,剪切失稳系统的运动方程为:σ(0)ij,j=ρ∂2ui∂t2i,j=1,2 (11) 根据高斯散度定理,局部运动方程可进一步表达为:
[LijklΔuk,l+σ(0)pjΔui,p],j=(ρ∂2ui∂t2)Q (12) 下面将基于式(12)和剪切边界条件,即:x2=0时,
σ(0)21=τ0 ,σ(0)1i=0 (i=1, 2)进行分析。1.2 动态剪切局部化失稳的扩散分析
由式(12),表征界面失稳临界状态的运动方程可进一步表达为:
L1212Δu1,22+L1111Δu1,11+σ(0)21Δu1,21+L1122Δu2,21+L1221Δu2,12=ρ∂2Δu1∂t2 (13) L2211Δu1,12+L2112Δu1,21+L2222Δu2,22+L2121Δu2,11+σ(0)21Δu2,21=ρ∂2Δu2∂t2 (14) 式中:本构关系
σij=Lijklεkl 采用Jaumann率下塑性可膨胀本构关系[40],则刚度张量Lijkl 表达为:Lijkl=E1+μ[12(gikgjl+gilgjk)+gijgkl(μ−E/3Eptm1−2μ+E/Eptm)−32σ2eEEpteSijSkl23(1+μ)+EEpte] (15) 式中:等效塑性切线模量
Epte=dσe/dεpe ,无损情况下,1Epte=1Et−1E ;单轴加载下的切线模量
Et=dσ1/dε1 ,dσe=32σeSkldσkl ;体积模量
Eptm=dσkk/dεpkk=dσm/dεpm ;等效应力
σe=√32SijSji ;等效塑性应变εpe=√32epijepji 。基于实验对岩石应力应变全过程[34-35]进行拟合,塑性阶段:
E/Epte =0~3,E/Eptm =0~120。结合文献[26,34-35]中的实验过程分析和局部化分叉的现象选取材料参数,如表1所示,其中λ 和μ 为拉梅系数。取塑性应变εpe =0.001,参数E/Eptm 和E/Epte 采用文献[40]推荐的方法进行试算,反复迭代确定其内部应力状态而得到。表 1 岩石材料参数Table 1. Rock material parametersρ/(kg·m−3) E/GPa ν λ/GPa μ/GPa E/Eptm E/Epte 2600 30 0.23 10.4 12.2 50 2 求解剪切局部化失稳的扩散行为。Hill等[30]给出了直角坐标系下速率变分的表达式为
δVi=Vi(n0) ,n0=n0ixi ,i=1,2 。其梯度为(图1):δVi,j=dVidnnj=cinj ,n1=cosα ,n2=sinα 。进一步,采用流函数ψ=v(x2)cos(c1x1) 来求解速度扰动。李国琛等[40]对平面应变问题采用的速率扰动函数表达式为:δVα=sin(πxαL)φΔ2 ,δVα=sin(πxαL)φΔ4 ,其中:Δ2 、Δ4 为待求的广义速率参量,φ=1+(1−ϖ)cos2πx3h−ϖcos4πx3h ,ϖ 为获取最低临界值指标的调节变量。Miles等[41]、Chau[33]则采用一种对称形位移函数χ(r,z)=X(r)cos(ηz) ,其中:η=kπ/L,k=1, 2, ···。对于剪切问题,Chau[42]采用一种反对称势函数形式:φ=C1ψ1(r)sin(ηz)cos(nθ) ,ψ=C2ψ2(r)cos(ηz)sin(nθ) 。其中η=mπ/mπLL , C1、C2为常数, m, n=0, 1, 2, ···。对于剪切局部化失稳的扩散行为,参照Chau[42]的反对称势函数形式,可设本征位移为:
Δu1=U1(x2)sin(ωx1)sinh(−k1t) (16) Δu2=U2(x2)cos(ωx1)sinh(−k2t) (17) 式中:Ui为关于x2的位移函数,k1、k2为参数。
将式(16)~(17)代入(13)~(14),为简化计算近似认为两个方向的动力影响相同,可以得到描述局部化失稳的扩散行为的控制方程组为:
L1212U″1+σ(0)21ωtan(ωx1)U′1−(ω2L1111+ρk21)U1−ω(L1122+L1221)U′2=0 (18) L2222U″2−σ(0)21ωtan(ωx1)U′2−(ω2L2121+ρk22)U2+ω(L2211+L2112)U′1=0 (19) 进一步整理,可得到:
a1U(4)1+b1U‴1+c1U″1+d1U′1+e1=0 (20) a2U(4)2+b2U‴2+c2U″2+d2U′2+e2=0 (21) 式中:
a1=L2222L1212ω(L1122+L1221),a2=−L2222L1212ω(L2211+L2112) b1=σ(0)21/tan(ωx1)L2222L1122+L1221−σ(0)21tan(ωx1)L1212L1122+L1221,b2=−σ(0)21/tan(ωx1)L2222L2211+L2112+σ(0)21tan(ωx1)L1212L2211+L2112 c1=−L2222(ω2L1111+ρk21)ω(L1122+L1221)−(σ(0)21)2ωL1122+L1221−(ω2L2121+ρk22)L1212ω(L1122+L1221)+ω(L2211+L2112) c2=L2222(ω2L1111+ρk21)ω(L2211+L2112)+(σ(0)21)2ωL2211+L2112+(ω2L2121+ρk22)L1212ω(L2211+L2112)−ω(L1122+L1221) d1=σ(0)21tan(ωx1)(ω2L1111+ρk21)L1122+L1221−σ(0)21/tan(ωx1)(ω2L2121+ρk22)L1122+L1221 d2=−σ(0)21tan(ωx1)(ω2L1111+ρk21)L2211+L2112+σ(0)21/tan(ωx1)(ω2L2121+ρk22)L2211+L2112 e1=(ω2L1111+ρk21)(ω2L2121+ρk22)ω(L1122+L1221),e2=−(ω2L1111+ρk21)(ω2L2121+ρk22)ω(L2211+L2112) 此控制方程组即为动力学过程中局部失稳存在的判断方程。若
tan(ωx1) 存在非平凡解,则系统中存在剪切局部化变形带,剪切带内的变形分布将由tan(ωx1) 决定。对其直接求解比较困难,下面将进行简化分析。2. 岩石动态剪切失稳求解
2.1 基于剪切变形带的求解
岩石动摩擦滑移主要发生在界面两侧1~5 mm范围内[9],剪切变形带作为界面破坏滑移的先导过程,一般小于此范围,岩石界面失稳主要发生在非常狭窄的区域。因此,直角坐标系下速率[30,40]变分的平面应力分量可采用Hill等[30]的表达式,即为:
δVi,j=Vi(n)n=nixi,i=1,2 (22) 式中:Vi为关于xi的速度函数,ni为参数。
其梯度为:
δVi,j=dVidnnj=cinj,n1=cosα,n2=sinα (23) 类似地,对于剪切局部化失稳的扩散行为,引入动力学项影响,可设本征位移为:
Δui=Xi(n,k)n=nix,k=kt,i=1,2 (24) 式中:Xi为关于xi的位移函数,k为局部动力参数。
将式(24)代入式(12),可得到描述剪切局部化失稳的扩散行为的控制方程组:
L1212X″1n22+L1111X″1n21+σ(0)21X″1n1n2+(L1122+L1221)X″2n1n2=ρX″1k2 (25) (L2211+L2112)X″1n1n2+L2222X″2n22+L2121X″2n21+σ(0)21X″2n1n2=ρX″2k2 (26) 设
r=X″1X″2 ,s=n2n1 ,z=k2n1n2=k2(s+1s) 代入式(25)~(26),可得:L1212L2222s4+(σ(0)21−ρz)(L1212+L2222)s3+(L1111L2222+L2121L1212+(σ(0)21−ρz)2−(L1122+L1221)(L2211+L2112))s2+(σ(0)21−ρz)(L1111+L2121)s+L1111L2121=0 (27) 此方程即为冲击过程中剪切局部化失稳存在的判断方程。其中,局部动力参数k体现了时间效应的影响。此方程中,若
s=tanα (α为剪切带角度)存在非平凡解,则系统中存在剪切带。为研究冲击过程惯性力和载荷对局部化剪切变形带的影响特征,对式(27)无量纲化,取
σ∗=σ(0)21/G ,a∗=ρz/σ(0)21 ,其中G为剪切模量。则式(27)可以改写为:(L1212L2222/G2)s4+σ∗(1−a∗)[(L1212+L2222)/G]s3+[L1111L2222/G2+L2121L1212/G2+(σ∗(1−a∗))2−(L1122+L1221)(L2211+L2112)/G2]s2+σ∗(1−a∗)[(L1111+L2121)/G]s+L1111L2121/G2=0 (28) 式(28)即为剪切局部化失稳存在的判别式。此式存在两组对称解(也称为共轭解),取对称部分换算为角度进行分析。图2所示为剪切带角度
α 随无量纲剪切力σ*和局部动力参数k的变化特征。剪切变形带角度α 为5°~9°,与岩石直剪实验结果[35]相近。由此可见:局部动力系数k的影响比外部剪切作用力的影响更为显著;随着外部剪切作用力的增大,剪切变形带角度α 有一定程度的增大,但增幅不是很大;随着局部动力系数k的增大,即局部惯性作用力的增大,剪切变形带角度α 明显减小。较强的动态作用过程,惯性效应明显,局部化程度加大,此时的剪切失效也集中于局部区域。2.2 控制方程的近似求解
对含惯性力作用的式(18)~(21)进行
tan(ωx1) 的直接求解比较困难,因此,可寻求其近似解,即:先对式(20)~(21)求解其静力学解,然后在静力学解的基础上进行动力学修正。基于此,对式(13)~(14),令加速度项等于零,即:L1212Δu1,22+L1111Δu1,11+σ(0)21Δu1,21+(L1122+L1221)Δu2,21=0 (29) (L2211+L2112)Δu1,21+L2222Δu2,22+L2121Δu2,11+σ(0)21Δu2,21=0 (30) 将式(24)代入式(29)~(30),可得:
L1212X″1n22+L1111X″1n21+σ(0)21X″1n1n2+(L1122+L1221)X″2n1n2=0 (31) (L2211+L2112)X″1n1n2+L2222X″2n22+L2121X″2n21+σ(0)21X″2n1n2=0 (32) 设
r=X″1X″2 ,s=n2n1 ,代入(31)~(32),可以得到:L1212L2222s4+σ(0)21(L1212+L2222)s3+[L1111L2222+L2121L1212+(σ(0)21)2−(L1122+L1221)(L2211+L2112)]s2+σ(0)21(L1111+L2121)s+L1111L2121=0 (33) 同理,取
σ∗=σ(0)21/σ(0)21GG 对方程无量纲化,可得到:(L1212L2222/G2)s4+σ∗[(L1212+L2222)/G]s3+[L1111L2222/G2+L2121L1212/G2+(σ∗)2−(L1122+L1221)(L2211+L2112)/G2]s2+σ∗[(L1111+L2121)/G]s+L1111L2121/G2=0 (34) 上式为关于s的多项式,一般有4个根(包括实根和虚根)。一旦得到实根s0,则表明确实存在剪切带失稳,此时得到静力学的解为:
(tanα)static =s0。然后,令s=s0+s1,将其回代到式(28),并忽略高阶小量,可得到:s1=σ∗a∗[(L1212+L2222)/G]s30−(σ∗)2[(a∗)2−2(a∗)]s20+σ∗a∗[(L1111+L2121)/G]s04(L1212L2222/G2)s30+3σ∗[(L1212+L2222)/G](1−a∗)s20+2[L1111L2222/G2+L2121L1212/G2−(L1122+L1221)(L2211+L2112)/G2+(σ∗)2(1−a∗)2]s0+σ∗(1−a∗)[(L1111+L2121)/G] (35) 此即为一阶近似解。再次令
s=s0+s1+s2 ,将其回代到式(28),并忽略高阶小量,可得到s2 的表达式,即为二阶近似解。以此类推,反复迭代,可以得到s在动力条件下的n阶近似表达式,即:s=s0+n∑i=1si 。其二阶近似解计算过程如图3所示。由此可见,惯性效应的存在使得剪切带的角度减小,反映出冲击过程变形高度局部化的趋势。3. 岩石动态剪切失稳变形的扩散分析
3.1 失稳变形的扩散分析
若方程组存在剪切局部化变形解,对于剪切变形区域,基于式(20)~(21)的扩散型位移解中的U1和U2具有如下形式:
U1=A1e−p1x2+A2e−p2x2+A3e−p3x2+A4ep4x2 (36) U2=B1e−p1x2+B2e−p2x2+B3e−p3x2+B4ep4x2 (37) 式中:A1、A2、A3、A4为待定系数;
p1=Z3+b14a1+Z2,p2=Z3+b14a1−Z2,p3=b14a1−Z3+Z1,p4=Z3−b14a1+Z1; Z1=√−9Z2/36√Z5−12Z8√Z5−Z29√Z5−3√6Z10√27Z210−72Z8Z9+2Z39+3√3Z7−12Z9Z1/36√Z5Z4, Z2=√3√6Z10√27Z210−72Z8Z9+2Z39+3√3Z7−12Z8√Z5−Z29√Z5−9Z2/36√Z5−12Z9Z1/36√Z5Z4, Z3=√Z56Z1/66,Z4=6Z1/66Z1/45,Z5=Z29+9Z2/36−6Z9Z1/36+12e1a1−9b4164a41+3b21c14a31−3b1d1a21, Z6=Z2102−4Z8Z93+Z3927+√3Z718,Z7=√27Z410−256Z38−16Z8Z49+4Z39Z210+128Z28Z29−144Z8Z9Z210, Z8=e1a1−3b41256a41+b21c116a31−b1d14a21,Z9=c1a1−3b218a21,Z10=d1a1+b318a31−b1c12a21; B1=[L1212ω(L1122+L1221)p1+σ(0)21/sL1122+L1221−ω(L1111+ρk21)L1122+L12211p1]A1,B2=[L1212ω(L1122+L1221)p2+σ(0)21/sL1122+L1221−ω(L1111+ρk21)L1122+L12211p2]A2, B3=[L1212ω(L1122+L1221)p3+σ(0)21/sL1122+L1221−ω(L1111+ρk21)L1122+L12211p3]A3,B4=[L1212ω(L1122+L1221)p4+σ(0)21/sL1122+L1221−ω(L1111+ρk21)L1122+L12211p4]A4 剪切作用下,即:x2=0时,
σ(0)21=τ0 ,σ(0)1i=0 (i=1, 2)。此问题的解以剪切变形为主,因此,为简化问题起见,在剪切变形带的边界,忽略法向惯性力,即令ρ∂2Δu1∂t2 =0。同时,考虑其切向惯性力,取ρ∂2Δu2∂t2 =τ 。则有:L2112Δu1,2+L2121Δu2,1+σ(0)21Δu2,2=τ (38) L1111Δu1,1+L1122Δu2,2+σ(0)21Δu1,2=0 (39) 此过程中,
τ/σ(0)21=β 。其中,β为无量纲的切向惯性力。为简化讨论过程,β取0.1。由于岩石动摩擦滑移主要发生在界面两侧1~5 mm范围内,此过程为破坏滑移的先导过程,失稳界面变形的位移小于滑移位移。剪切变形带通常处于一个狭窄区域,靠近基体侧的位移更小。因此,可增加位移边界条件:
x2=h,Δui=0i=1,2 (40) 式中:h为岩石动摩擦影响厚度。式(36)的待定系数A1、A2、A3、A4可由上述边界条件确定。
另外,参数k1、k2由下式所列的初始条件确定:
t=0,Δui=0,∂Δui∂t=0i=1,2 (41) 定义式(33)表达的等效位移u表征岩石界面失稳的位移。
u=√(Δu1)2+(Δu2)2 (42) 基于此,由式(36)~(42),结合表1所列岩石参数,可以进行岩石界面的动态剪切失稳变形的扩散分析。图4所示为局部化失稳变形扩散过程的初步计算结果,失稳界面的厚度h取为0.3 mm, 此时,
tanα =0.1659 ,σ∗ =0.002,β=0.1。由此可见:随着深度的增加,剪切扩散作用下的失稳位移越来越小;随着加载时间的增加,局部位移越来越大,直到局部产生破坏形成滑移系。计算结果与实验观察现象[10-11]基本一致。3.2 剪切力的影响
图5(a)所示为失稳表面的位移特征随着剪应力增大而变化的规律。在所加载的范围内,随着剪切应力的增大,剪切变形带的角度变化不大,与剪切失稳对应的5个峰的峰位比较接近;但是最大局部位移随剪切应力的增大增加较快。图5(b)所示为将失稳界面的最大位移进行投影后的图,图中剪切力较大时用蓝色和淡蓝色线表示。由此可见:与剪切失稳对应的5个峰的峰位比较接近;界面失稳位移增加幅度随着界面上剪应力的增大而增大。这也表明:当岩石所受剪切作用力过大时,材料会在局部产生失稳的基础上逐渐发展到材料破坏,例如产生岩石类材料常见的局部的剪切破裂。
3.3 界面厚度的影响
图6(a)为失稳表面位移随界面厚度变化的情况,界面厚度分别设置为0.3、0.5和1.0 mm,此时,令
σ∗ =0.002。随着界面厚度的增大,可变性位移增大,为了研究界面厚度对变形的影响,应用下式定义的相对变形来分析:εh=umaxh (43) 等效位移和应变随界面厚度变化的关系如图6(b)所示,随着界面厚度的增大,可变形区域增大;由于变形区的累积和几何改变,界面表面褶皱加深,即界面越厚,失稳位移越大。
4. 岩石界面剪切失稳数值分析
4.1 有限元模型和材料参数
应用ABAQUS软件建立如图7所示的双层结构模型进行有限元模拟计算。模型尺寸为12.5 mm×40 mm×100 mm。模型中岩石失稳界面等效成一个厚度较小的弹性层,其上层材料与下层基底完全粘结,右端施加完全约束条件,左端施加压应力,底部与侧面限制离面位移。动态载荷在ABAQUS软件载荷界面使用表(tabular)添加,如图7所示,上表面局部区域施加恒值剪力载荷作用,以模拟岩石受压过程中的局部剪切作用力。有限元计算网格类型采用C3D8R,同时对受局部剪力部分的网格进行加密。由此数值分析局部剪应力处界面失稳特征。
基体和界面层的材料参数选取主要基于张磊等[10-11]对含斜节理岩石试样进行SHPB杆束冲击的实验结果。图8(a)所示为冲击前后岩石表面粗糙系数的变化规律,图8(b)所示为岩石界面摩擦应力与表面粗糙系数对比度的关系。冲击后岩石表面的粗糙系数有一定程度的升高,冲击后粗糙表面的高程有0.2~0.4 mm的降低,均反映出一定程度的局部破碎的加剧。由表面的粗糙系数的变化幅度估计摩擦界面区材料性能的折减系数约0.2~0.4,为讨论界面影响,有限元计算参数选取如下:基体密度为
2600 kg/m3,弹性模量为30 GPa,泊松比为0.23;界面区按照基体材料属性进行折减,折减系数为0.1、0.2、0.3、0.4和0.5。4.2 有限元计算结果
图9所示为局部剪切力为20 MPa、折减系数为0.1时的位移形貌。计算结果表明:当试样受到外部载荷时,若局部薄弱处存在集中剪力,则会产生褶皱行为;随着剪应力的增大,局部变形增加;试样将会沿着失稳界面处产生滑移等大变形。利用数值模拟得到了类似现象,印证了上述理论的合理性。在有限变形的条件下,随着加载进一步施加,即会产生大变形,进而产生破坏。局部应力过大导致岩石材料分叉,为材料破坏的先导阶段。
4.3 界面参数的影响
对于小尺度的界面薄层,确定其具体的力学参数是十分困难的,类似于混凝土界面过渡区(interfacial transition zone,ITZ),大部分是根据模量进行折减,并经过有限元反复试算得到[43]。近年来,随着纳米压痕技术的发展,对中微尺度薄弱界面的力学性能参数有一定程度的评估,但获得其准确参数仍然是个挑战。下面将基于杆束实验[10-11]前后粗糙度的变化程度,讨论参数对界面起皱形貌的影响,为详细讨论,按照折减系数0.1、0.2、0.3、0.4和0.5给予界面层参数。图10所示为最大位移(umax)与折减系数的关系,此时保持界面厚度为0.3 mm, 局部剪切力为20 MPa。
从图10可以看出,界面参数对数值计算结果影响较大,因此,在数值计算中应选择合适的计算参数来讨论其影响。当界面参数与基体参数差异较小时不易引起起皱,差异较大时,在局部受载时容易形成失稳扰动,影响表面位移形貌。为了方便讨论此现象,以下分析中的折减系数取为0.1。
4.4 剪切力的影响
图11所示为界面厚度为0.3 mm,局部剪切力分别为10、20、30 MPa条件下的有限元计算结果。图11结果表明,当试样受到外部载荷时,试样局部薄弱位置产生分叉变形,表现为起皱现象;随着剪应力的增大,形貌变化,局部变形加剧;随着加载时间的延长,局部变形逐渐由小变形向大变形转变,导致试样将沿着失稳截面处产生滑移,从而产生损伤破坏。
4.5 界面厚度的影响
图12所示为剪切力20 MPa,界面厚度分别为0.3、0.5、1 mm条件下的有限元计算结果。由图12可以看出,随着失稳界面的增厚,局部力所产生的响应将被失稳界面层的变形进一步消耗掉,峰值整体位移表现为增大。
岩石SHPB杆束实验[10-11]中,由于岩石的不均匀性,发生分叉时各处界面厚度变化不一,因此裂纹萌生发展并不仅仅产生在一处,而是会形成多处裂纹路径,致密性较好的岩石可能呈规律的片状破坏(沿着分叉处滑移破坏),非均匀程度较大的岩石或呈无规律性的剪切破坏(分叉处较多整体表现较为无序)[44-45]。
4.6 局部化变形对波传播的影响
图13上半部分所示为剪切力为20 MPa,界面厚度分别为0.3、0.5、1.0 mm条件下,界面中心位置上层和下层应力变化的计算结果。结果表明,局部剪切力的大小随着界面层的起皱,载荷输入能量会被消耗,体现出衰减现象;弱界面越厚,相对衰减幅度越大。图13下部分所示为对失稳界面上层和下层两条波进行傅里叶变换得到的频域图,发现相关频域的幅值发生了改变,即界面前后不仅波的幅值改变了,波的频率也改变了。
岩石SHPB杆束实验[10-11]中,由于岩石的不均匀性,发生分叉时各处界面厚度变化不一,这种变形会消耗加载波的能量使之产生衰减,若局部化失稳持续产生,则会演化成为断裂行为;若仅局部化到一定程度继续加载,宏观应力-应变曲线上表现为:在塑性阶段,随着应变的增大,应力先微小下降再上升。此现象为断裂破坏的先导过程,表现出剪切变形带耗能的特点。局部化剪切变形带的形成和发展,为断裂滑移的先导,研究此现象及其动力学特性,可为评估局部化失稳对岩石局部化断裂、地震波传播过程中滑移界面等的影响提供参考。
5. 结 论
利用有限变形理论,对岩石界面进行了失稳分析,建立了岩石材料界面剪切局部化失稳的理论分析模型;通过实验和数值分析,说明了模型的可靠性;基于此分析了岩石界面失稳的影响因素,发现变形局部化对波传播具有衰减作用;得到的主要结论如下。
(1)基于能量法和扰动分析,提出了界面受剪力影响下局部化失稳的判别式和局部变形扩散方程。
(2)较强的动态作用过程中,惯性效应的影响明显,岩石界面局部化程度较深。剪切变形带角度随着外部剪切作用力的增大而增大,随着局部动力系数的增大明显减小。
(3)剪切力越大,局部位移越大,岩石越容易发生破坏;界面厚度越小,局部位移越大,更容易沿界面破坏。界面剪切扩散行为极大降低了透射波的幅值,同时也改变了透射波的频率。
-
表 1 岩石材料参数
Table 1. Rock material parameters
ρ/(kg·m−3) E/GPa ν λ/GPa μ/GPa E/Eptm E/Epte 2600 30 0.23 10.4 12.2 50 2 -
[1] SCHOLZ C H. Wear and gouge formation in brittle faulting [J]. Geology, 1987, 15(6): 493–495. DOI: 10.1130/0091-7613(1987)15<493:WAGFIB>2.0.CO;2. [2] SCHOLZ C H. Earthquakes and friction laws [J]. Nature, 1998, 391(6662): 37–42. DOI: 10.1038/34097. [3] COWIE P A, SCHOLZ C H. Growth of faults by accumulation of seismic slip [J]. Journal of Geophysical Research: Solid Earth, 1992, 97(B7): 11085–11095. DOI: 10.1029/92jb00586. [4] BEELER N M, TULLIS T E, WEEKS J D. The roles of time and displacement in the evolution effect in rock friction [J]. Geophysical Research Letters, 1994, 21(18): 1987–1990. DOI: 10.1029/94GL01599. [5] 徐松林, 单俊芳, 王鹏飞. 脆性材料高应变率压缩失效机制综述与研究进展 [J]. 现代应用物理, 2020, 11(3): 030101. DOI: 10.12061/j.issn.2095-6223.2020.030101.XU S L, SHAN J F, WANG P F. Review and research progress of dynamic failure mechanism for brittle materials under high strain rate [J]. Modern Applied Physics, 2020, 11(3): 030101. DOI: 10.12061/j.issn.2095-6223.2020.030101. [6] 单俊芳, 徐松林, 张磊, 等. 岩石节理动摩擦过程中的声发射和产热特性研究 [J]. 实验力学, 2020, 35(1): 41–57. DOI: 10.7520/1001-4888-19-121.SHAN J F, XU S L, ZHANG L, et al. Investigation on acoustic emission and heat production characteristics on joint surfaces due to dynamic friction [J]. Journal of Experimental Mechanics, 2020, 35(1): 41–57. DOI: 10.7520/1001-4888-19-121. [7] 徐松林, 郑文, 刘永贵, 等. 冲击下花岗岩界面动态摩擦特性实验研究 [J]. 高压物理学报, 2011, 25(3): 207–212. DOI: 10.11858/gywlxb.2011.03.003.XU S L, ZHENG W, LIU Y G, et al. Experimental investigation on interface dynamic friction of granite under combined pressure and shear impact loading [J]. Chinese Journal of High Pressure Physics, 2011, 25(3): 207–212. DOI: 10.11858/gywlxb.2011.03.003. [8] DI TORO G, PENNACCHIONI G. Superheated friction-induced melts in zoned pseudotachylytes within the Adamello tonalites (Italian Southern Alps) [J]. Journal of Structural Geology, 2004, 26(10): 1783–1801. DOI: 10.1016/j.jsg.2004.03.001. [9] RICE J R. Heating and weakening of faults during earthquake slip [J]. Journal of Geophysical Research: Solid Earth, 2006, 111(B5): 311–340. DOI: 10.1029/2005JB004006. [10] 张磊, 徐松林, 施春英. 应用杆束系统研究水泥砂浆节理面的压剪动特性 [J]. 实验力学, 2016, 31(2): 175–185. DOI: 10.7520/1001-4888-15-220.ZHANG L, XU S L, SHI C Y. On the dynamic compression-shear characteristics of cement mortar joint surface based on a bunched bar system [J]. Journal of Experimental Mechanics, 2016, 31(2): 175–185. DOI: 10.7520/1001-4888-15-220. [11] 张磊, 王文帅, 苗春贺, 等. 花岗岩粗糙表面动摩擦形态演化 [J]. 高压物理学报, 2021, 35(3): 031201. DOI: 10.11858/gywlxb.20200640.ZHANG L, WANG W S, MIAO C H, et al. Rough surface morphology of granite subjected to dynamic friction [J]. Chinese Journal of High Pressure Physics, 2021, 35(3): 031201. DOI: 10.11858/gywlxb.20200640. [12] ZHOU L J, XU S L, SHAN J F, et al. Heterogeneity in deformation of granite under dynamic combined compression/shear loading [J]. Mechanics of Materials, 2018, 123: 1–18. DOI: 10.1016/j.mechmat.2018.04.013. [13] OKADA M, LIOU N S, PRAKASH V, et al. Tribology of high-speed metal-on-metal sliding at near-melt and fully-melt interfacial temperatures [J]. Wear, 2001, 249(8): 672–686. DOI: 10.1016/S0043-1648(01)00698-6. [14] XU S, HUANG J, WANG P, et al. Investigation of rock material under combined compression and shear dynamic loading: an experimental technique [J]. International Journal of Impact Engineering, 2015, 86(1): 206–222. DOI: 10.1016/j.ijimpeng.2015.07.014. [15] BEN-DAVID O, FINEBERG J. Static friction coefficient is not a material constant [J]. Physical Review Letters, 2011, 106(25): 254301. DOI: 10.1103/PhysRevLett.106.254301. [16] DI TORO G, HAN R, HIROSE T, et al. Fault lubrication during earthquakes [J]. Nature, 2011, 471(7339): 494–498. DOI: 10.1038/nature09838. [17] RUBINO V, ROSAKIS A J, LAPUSTA N. Understanding dynamic friction through spontaneously evolving laboratory earthquakes [J]. Nature Communications, 2017, 8: 15991. DOI: 10.1038/ncomms15991. [18] PASSELÈGUE F X, SCHUBNEL A, NIELSEN S, et al. From sub-Rayleigh to supershear ruptures during stick-slip experiments on crustal rocks [J]. Science, 2013, 340(6137): 1208–1211. DOI: 10.1126/science.1235637. [19] RUBINSTEIN S M, COHEN G, FINEBERG J. Detachment fronts and the onset of dynamic friction [J]. Nature, 2004, 430(7003): 1005–1009. DOI: 10.1038/nature02830. [20] RUBINSTEIN S M, COHEN G, FINEBERG J. Dynamics of precursors to frictional sliding [J]. Physical Review Letters, 2007, 98(22): 226103. DOI: 10.1103/PhysRevLett.98.226103. [21] GEOBEL T H W, SCHORLEMMER D, BECKER T W, et al. Acoustic emissions document stress changes over many seismic cycles in stick-slip experiments [J]. Geophysical Research Letters, 2013, 40(10): 2049–2054. DOI: 10.1002/grl.50507. [22] JOHNSON P A, FERDOWSI B, KAPROTH B M, et al. Acoustic emission and microslip precursors to stick-slip failure in sheared granular material [J]. Geophysical Research Letters, 2013, 40(21): 5627–5631. DOI: 10.1002/2013GL057848. [23] MENDES R S, MALACARNE L C, SANTOS R P B, et al. Earthquake-like patterns of acoustic emission in crumpled plastic sheets [J]. Europhysics Letters, 2010, 92(2): 29001. DOI: 10.1209/0295-5075/92/29001. [24] HUANG J Y, XU S L, HU S S. Numerical investigations of the dynamic shear behavior of rough rock joints [J]. Rock Mechanics and Rock Engineering, 2014, 47(5): 1727–1743. DOI: 10.1007/s00603-013-0502-8. [25] GAO K, GUYER R, ROUGIER E, et al. From stress chains to acoustic emission [J]. Physical Review Letters, 2019, 123(5): 048003. DOI: 10.1103/PhysRevLett.123.048003. [26] 徐松林, 吴文. 岩土材料局部化变形分岔分析 [J]. 岩石力学与工程学报, 2004, 23(20): 3430–3438. DOI: 10.3321/j.issn:1000-6915.2004.20.007.XU S L, WU W. Bifurcation analysis on deformation localization of geomaterial [J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(20): 3430–3438. DOI: 10.3321/j.issn:1000-6915.2004.20.007. [27] 钱伟长. 非线性力学的新进展-稳定性、分叉、灾变、浑沌 [M]. 武汉: 华中理工大学出版社, 1988.QIAN W C. New progress in nonlinear mechanics stability, bifurcation, catastrophe and chaos [M]. Wuhan: Huazhong University of Technology Press, 1988. [28] OTTOSEN N S, RUNESSON K. Properties of discontinuous bifurcation solutions in elasto-plasticity [J]. International Journal of Solids and Structures, 1991, 27(4): 401–421. DOI: 10.1016/0020-7683(91)90131-X. [29] FLECK N A, HUTCHINSON J W. A phenomenological theory for strain gradient effects in plasticity [J]. Journal of the Mechanics and Physics of Solids, 1993, 41(12): 1825–1857. DOI: 10.1016/0022-5096(93)90072-N. [30] HILL R, HUTCHINSON J W. Bifurcation phenomena in the plane tension test [J]. Journal of the Mechanics and Physics of Solids, 1975, 23(4/5): 239–264. DOI: 10.1016/0022-5096(75)90027-7. [31] RUDNICKI J W, RICE J R. Conditions for the localization of deformation in pressure-sensitive dilatant materials [J]. Journal of the Mechanics and Physics of Solids, 1975, 23(6): 371–394. DOI: 10.1016/0022-5096(75)90001-0. [32] YOUNG N J B. Bifurcation phenomena in the plane compression test [J]. Journal of the Mechanics and Physics of Solids, 1976, 24(1): 77–91. DOI: 10.1016/0022-5096(76)90019-3. [33] CHAU K T. Non-normality and bifurcation in a compressible pressure-sensitive circular cylinder under axisymmetric tesion and compression [J]. International Journal of Solids and Structures, 1992, 29(7): 801–824. DOI: 10.1016/0020-7683(92)90017-N. [34] 徐松林, 吴文, 李廷, 等. 三轴压缩大理岩局部化变形的试验研究及其分岔行为 [J]. 岩土工程学报, 2001, 23(3): 296–301. DOI: 10.3321/j.issn:1000-4548.2001.03.008.XU S L, WU W, LI T, et al. Experimental studies on localization and bifurcation behaviors of a marble under triaxial compression [J]. Chinese Journal of Geotechnical Engineering, 2001, 23(3): 296–301. DOI: 10.3321/j.issn:1000-4548.2001.03.008. [35] 徐松林, 吴文, 张华, 等. 直剪条件下大理岩局部化变形研究 [J]. 岩石力学与工程学报, 2002, 21(6): 766–771. DOI: 10.3321/j.issn:1000-6915.2002.06.002.XU S L, WU W, ZHANG H, et al. Testing study on localization of marble under direct shear [J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 21(6): 766–771. DOI: 10.3321/j.issn:1000-6915.2002.06.002. [36] 徐松林, 吴文, 张奇华, 等. 大理岩有限变形分岔分析 [J]. 岩土工程学报, 2002, 24(1): 42–46. DOI: 10.3321/j.issn:1000-4548.2002.01.009.XU S L, WU W, ZHANG Q H, et al. Bifurcation analyses of finite/large deformation for a marble [J]. Chinese Journal of Geotechnical Engineering, 2002, 24(1): 42–46. DOI: 10.3321/j.issn:1000-4548.2002.01.009. [37] IKEDA K, MUROTA K. Imperfect bifurcation in structures and materials [M]. 3rd ed. Cham: Springer, 2019: 201–291. DOI: 10.1007/978-3-030-21473-9. [38] IKEDA K, MURAKAMI S, SAIKI I, et al. Image simulation of uniform materials subjected to recursive bifurcation [J]. International Journal of Engineering Science, 2001, 39(17): 1963–1999. DOI: 10.1016/s0020-7225(01)00038-6. [39] LEE D, TRIANTAFYLLIDIS N, BARBER J R, et al. Surface instability of an elastic half space with material properties varying with depth [J]. Journal of the Mechanics and Physics of Solids, 2008, 56(3): 858–868. DOI: 10.1016/j.jmps.2007.06.010. [40] 李国琛, 耶纳 M. 塑性大应变微结构力学 [M]. 北京: 科学出版社, 1993: 165–282.LI G C, JENNA M. Plastic large strain microstructure mechanics [M]. Beijing: Science Press, 1993: 165–282. [41] MILES J P, NUWAYHID U A. Bifurcation in compressible elastic/plastic cylinders under uniaxial tension [J]. Applied Scientific Research, 1985, 42(1): 33–54. DOI: 10.1007/BF00382529. [42] CHAU K T. Antisymmetric bifurcations in a compressible pressure-sensitive circular cylinder under axisymmetric tension and compression [J]. Journal of Applied Mechanics, 1993, 60(2): 282–289. DOI: 10.1115/1.2900791. [43] KIM S M, AL-RUB R K A. Meso-scale computational modeling of the plastic-damage response of cementitious composites [J]. Cement and Concrete Research, 2011, 41(3): 339–358. DOI: 10.1016/j.cemconres.2010.12.002. [44] 袁良柱, 苗春贺, 单俊芳, 等. 冲击下混凝土试样应变率效应和惯性效应探讨 [J]. 爆炸与冲击, 2022, 42(1): 013101. DOI: 10.11883/bzycj-2021-0114.YUAN L Z, MIAO C H, SHAN J F, et al. On strain-rate and inertia effects of concrete samples under impact [J]. Explosion and Shock Waves, 2022, 42(1): 013101. DOI: 10.11883/bzycj-2021-0114. [45] CHEN M D, XU S L, YUAN L Z, et al. Influence of stress state on dynamic behaviors of concrete under true triaxial confinements [J]. International Journal of Mechanical Sciences, 2023, 253: 108399. DOI: 10.1016/j.ijmecsci.2023.108399. 期刊类型引用(1)
1. 姚梦圆,董春亮,郝建平,谢雅. 基于Poyting-Thomson的岩石黏弹塑性时效损伤蠕变研究. 河南城建学院学报. 2024(06): 25-34 . 百度学术
其他类型引用(0)
-