岩石球砾与平板碰撞法向恢复系数

2024-02-28孙翰卿许龙朱成马新军黄俊光

科学技术与工程 2024年3期

孙翰卿, 许龙, 朱成, 马新军, 黄俊光

(1.华南理工大学土木与交通学院, 广州 510640; 2.广州建筑股份有限公司, 广州 510000;3.广东省有色矿山地质灾害防治中心, 广州 510080; 4.广州市建筑集团有限公司, 广州 510030;5.广州市设计院集团有限公司, 广州 510620)

落石灾害属于自然界中危险性极大的地质灾害之一,往往具有多发性、突发性、随机性等特点,对中国山区高速公路、铁路的运营及房屋安全造成了严重的危害[1-3]。了解落石灾害现象中岩块的运动轨迹、运动机理及冲击特性可为防护结构的设计提供重要理论依据[4-6]。法向恢复系数是影响落石与边坡碰撞过程和能量损耗的关键参数之一,其取值直接决定了落石的运动轨迹。因此,开展落石法向恢复系数研究对提高落石运动轨迹预测和优化防护措施布置具有重要的理论意义和工程应用价值。

落石与边坡的碰撞属于动力接触问题,其本质可以简化成岩石球砾与平板的接触碰撞[7-9]。中外学者通过碰撞试验、理论推导和数值模拟的方法探究球砾法向恢复系数的影响因素。在试验研究方面,章广成等[10]和叶四桥[11]等分别开展了落石野外碰撞试验和模型试验,研究发现坡面特征是影响落石法向恢复系数的主要影响因素,而落石形状、尺寸、碰撞速度、入射角度对法向恢复系数没有明显影响。叶阳等[12]开展了花岗岩球砾室内碰撞试验,试验结果却表明球砾的法向恢复系数存在明显的尺寸效应,且随着碰撞速率、等效弹性模量和含水率的增加而降低。由此可见,基于试验结果得到不同因素对岩石球砾恢复系数的影响仍存在争议,有待进一步深入研究。在理论研究方面,国内外学者基于Hertz弹性接触理论,推导了多种弹塑性接触碰撞模型和黏弹性接触碰撞模型[13-15],但模型理论计算结果与实际仍存在差异,部分模型参数需要试验拟合确定。在数值模拟研究方面,Guzzetti等[16]基于数值软件计算得到落石法向恢复系数依据坡面的软硬状态在0.1～0.5取值。章广成等[10]通过ANSYS/LS-DYNA软件模拟落石与坡面碰撞过程,坡面岩土体从基岩、碎屑堆积层至松散碎石土,落石法向碰撞恢复系数逐渐减小,取值范围在0.3～0.8。

综上所述,国内外学者针对坡面特征对落石法向恢复系数的影响已达成统一认识,即坡面由硬到软,法向恢复系数逐渐减小。然而,学者们对碰撞速率对落石法向恢复系数的影响尚未形成一致结论,关于岩石材料参数对法向恢复系数的影响更是分析不足。基于落石碰撞试验研究法向恢复系数难以灵活设定岩石的材料参数,且费时费力,特别是现场落石试验还存在安全性问题。相比之下,基于理论推导和数值模拟方法可更为方便、可靠地研究碰撞速率和材料参数对落石法向恢复系数的影响。

因此,现建立岩石球砾与平板碰撞理论模型,利用限元分析软件ABAQUS模拟了球砾与平板的碰撞过程,并采用上述两种方法探究了碰撞速率、弹性模量、弹性极限及泊松比对球砾法向恢复系数的影响,比较了两种方法结果的差异。相关研究成果可为落石运动轨迹分析中法向恢复系数参数取值提供理论支撑。

1 球板碰撞理论模型

基于Hertz经典弹性理论模型[17-20]与杜妍辰两颗粒弹塑性碰撞理论模型[21],建立了简化的岩石球砾与平板的接触碰撞模型。为便于推导,本文做出如下假设。

(1)将岩石球砾和平板视为各向同性的均质岩块。

(2)不考虑岩石的硬化,弹性极限和强度极限一致。

(3)不考虑碰撞过程中岩石球砾与平板之间的摩擦力。

(4)平板的尺寸无限大。

岩石球砾以一定速度v0与平板发生法向碰撞(球砾只有垂直于板的速度,切向速度为0)。由于假定平板的尺寸无限大,可将平板看作半径无穷大的球体,则有平板的质量m2远远大于球砾的质量m1,因此,整个接触碰撞过程中可近似认为平板的速度始终为0。

岩石球砾与平板开始接触瞬间作为分析的初始状态,球砾与平板的脱离瞬间作为结束状态,可以将碰撞过程分为3个阶段[14,21]。

(1)弹性变形阶段:该阶段是以球板开始接触瞬间为起始,以接触面圆心处球砾或平板发生塑性变形为结束。岩石球砾的速度由v0减小至v1,球砾和平板的变形为弹性变形,球砾的动能全部转换为球砾和平板的弹性势能。

(2)弹塑性变形阶段:该阶段是以接触面圆心处岩石球砾或平板发生塑性变形为起始,以球砾速度减小至0为结束。球砾的速度由v1减小至0,球砾和平板的变形为弹塑性变形,球砾的动能部分转换为球砾和平板的弹性势能。

(3)弹性回弹阶段:该阶段是以球砾速度减小至0为起始,以球砾与平板的脱离瞬间为结束。球砾的速度由0增大至v末,球砾和平板的变形为弹性变形,球砾和平板的弹性势能转换为球砾的动能。

1.1 弹性变形阶段

球砾与平板接触碰撞模型如图1所示。假设球砾的球心位移为δ,球砾和平板之间的总接触力为p,接触面为圆形,圆半径为R0,接触面半径为ra,接触区中心压应力为q0,其他相关参数如表1所示。

表1 球砾和平板的相关参数Table 1 Theparameters of the rock sphere and the plate

图1 球砾与平板接触碰撞模型Fig.1 The rock sphere and the plate collision model

根据Hertz弹性接触理论[19,21]有

(1)

(2)

(3)

(4)

结合式(1)和式(2),有

(5)

(6)

(7)

弹性变形阶段结束时,球砾与板接触区域中心处(r=0)即将进入塑性变形阶段,该处接触压应力σ=σp,σp为塑性段接触应力。此时将式(5)、式(6)代入(4)式得接触半径ra1为

(8)

将式(8)结果代入式(6)得球心位移δ1为

(9)

解式(7)微分方程,考虑初始条件,当球心位移δ=0时,球砾速度v=v0,有

(10)

弹性变形阶段结束时,球心位移δ=δ1,弹性阶段末球砾速度v1为

(11)

1.2 弹塑性变形阶段

(12)

结合式(4)、式(5)、式(6)和式(12),则弹塑性区总接触力为

(13)

根据牛顿第二定律,有

(14)

解式(14)微分方程,考虑初始条件δ=δ1时v=v1,得到弹塑性阶段球砾速度v的表达式为

(15)

弹塑性阶段结束时,有v=0,代入式(15),得到球砾球心近似位移δ2为

(16)

1.3 弹性回弹阶段

该阶段岩石球砾在弹性回弹力的作用下与平板逐渐分离,弹性应变能逐步转化成球砾的动能,而塑性应变不可恢复。假定球心位移δ中包含的塑性位移δp,且总接触力大小与球心位移关系与弹性变形阶段相同,则该阶段的总接触力为

(17)

由于弹塑性变形阶段末与弹性回弹阶段初总接触力大小相等,则有

(18)

求解式(18),可得

(19)

根据牛顿第二定律,有

(20)

解式(20)微分方程,考虑初始条件δ=δ2时v=0,得到弹性回弹阶段岩石球砾速度v的表达式为

(21)

该阶段末,球心位移δ=δp,因此岩石球砾与平板碰撞结束后的速度v末为

(22)

法向恢复系数e为岩石球砾与平板碰撞后的速度与碰撞前的速度比值,即

(23)

2 球板碰撞有限元模拟

利用SIMULIA公司的大型有限元软件ABAQUS模拟岩石球砾与平板的碰撞过程,并对比分析数值模拟结果与理论计算结果。数值模型如图2所示,岩石球砾的半径为50 mm,平板的半径为250 mm,厚度为100 mm。球砾底部与平板上表面中心恰好接触,接触采用动力学的接触方法、有限滑移,接触属性为无摩擦。考虑岩石球砾与平板碰撞属于轴对称问题,因此,岩石球砾与平板部件构建均采用轴对称单元,该单元可以在不牺牲计算精度的前提下大幅提高计算效率,且能在计算结束后通过ODB软件扫掠单元形成球板碰撞三维模型。在网格划分方面,球砾与平板模型均采用四边形单元(CAX4R),并在接触区域对模型网格进行加密处理,以提高计算精度;在远离接触区域网格划分较稀疏,以提高计算速度。岩石球砾轴对称面包含1 754个四边形单元,平板轴对称面包含2 419个四边形单元,材料的相关参数初始取值参考文献[12],如表2所示。边界条件设定方面,球砾和平板轴线所在边界施加对称约束以限制水平方向位移,平板底施加固定约束,以消除碰撞过程中平板弹性振动对碰撞的影响。

表2 岩石球砾和平板的材料参数Table 2 The materialparameters of the rock sphere and the plate

主要研究岩石球砾的法向恢复系数,因此,球板碰撞数值模拟中,岩石球砾碰撞前速度方向垂直于平板,碰撞速率设置为1.0～10.0 m/s。模型计算分析步采用动态显示分析步,考虑到碰撞时间和碰撞速率有关,当碰撞速率大于或等于2.5 m/s时,分析步时间设置为4×10-4;当碰撞速率小于2.5 m/s时,分析步时间设置为6×10-4。在初始分析步对球砾施加竖直方向碰撞速率v0,并在整个计算过程中实时监测球心竖直方向的速度v。模拟结束后,绘制球心竖直方向速度曲线,读取球心碰撞后的速度v末,根据球砾v末与v0的比值,计算得到球砾的法向恢复系数。

3 理论计算与数值模拟结果分析

3.1 碰撞速率对法向恢复系数的影响

岩石球砾法向恢复系数理论计算结果可通过式(11)、式(16)、式(19)和式(23)联立求解,式中参数m1、R0、v0为球砾的质量、半径和初始速度,等效弹性模量Et由球砾与平板的弹性模量和泊松比共同决定,均为已知参数,仅塑性区接触应力σp为待定参数。

岩石球砾与平板碰撞属于局部接触问题,接触面的半径远小于球砾的半径,因此接触区域受到周围约束作用,处于三向应力状态[12],σp的取值需要考虑围压对弹性极限的增大作用。Li等[22]研究发现,当接触区中心压应力至少为材料弹性极限的1.6倍时,材料才开始屈服;贾乃文[23]推导出刚性平冲头压入的近似解答,认为该增大系数取值应当在1+π/2附近。实际上,该增大系数取值和围压大小有关,围压大小又与球砾碰撞速率、材料参数相关。考虑到本模型假定球砾的初始碰撞速率在1.0～10.0 m/s之间,笔者根据数值模拟结果,读取接触面塑性区域的接触应力,并与材料的弹性极限σs比较,得到塑性区接触应力约为σs的1.95倍。因此,弹塑性变形阶段塑性区接触应力σp为

σp=min(σp1,σp2)=min(1.95σs1,1.95σs2)

(24)

为分析碰撞速率对岩石球砾法向恢复系数的影响,球砾的碰撞速率依次设置为1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0 m/s,理论计算与数值模拟选用的材料参数一致,如表2所示。理论计算与数值模拟得到的法向恢复系数e的结果如图3所示。

图3 不同碰撞速率下球砾法向恢复系数Fig.3 The normal restitution coefficient of rock sphere upon various impact velocity

数值模拟和理论计算得到的法向恢复系数十分接近,最大差值在5%以内,验证了本文数值模拟结果的可靠性。两种计算结果均表明,岩石球砾法向恢复系数随着碰撞速率的增加而逐渐降低,且下降速率逐渐减小。究其原因,随着碰撞速率的增加,球砾与平板接触区局部塑性变形的大小与范围增大,球砾的耗能增大,恢复系数降低。当球砾的碰撞速率较低时,碰撞速率是决定恢复系数主要影响因素之一,而当碰撞速率较高时,碰撞速率不再是影响恢复系数的主要因素,因此法向恢复系数下降速率随着碰撞速率的增加而逐渐减小。图3所示法向恢复系数变化曲线与何思明等[13]在Thornton理论基础上得到的法向恢复系数随碰撞速率变化曲线趋势基本一致,同样证明了本文理论计算和数值模拟结果的可靠性。

3.2 弹性模量对法向恢复系数的影响

为讨论弹性模量对法向恢复系数的影响,球砾与平板的弹性模量依次取24、26、28、30、32、34、36、38、40 GPa,理论计算与数值模拟选取的碰撞速率和材料参数如表3所示。法向恢复系数e的理论计算和数值模拟结果如图4所示。

表3 不同弹性模量下球砾和平板的相关参数Table 3 Theparameters of the rock sphere and the plate with various elastic modulus

理论计算和数值模拟结果表明,随着岩石球砾与平板的弹性模量增加,球砾的法向恢复系数降低,且两者基本呈线性关系,如图4所示。这是由于球砾与平板的弹性模量越大,材料达到弹性极限时弹性应变越小,进而导致弹性应变能减小。在弹性回弹阶段,弹性应变能再次转化为球砾的动能,因此球砾碰撞后的动能随着弹性应变能的减小而降低,球砾的恢复系数也随之降低。

3.3 弹性极限对法向恢复系数的影响

为讨论材料弹性极限对岩石球砾法向恢复系数的影响,球砾与平板的弹性极限依次取70、80、90、100、110、120、130、140和150 MPa,其余参数取值如表4所示。岩石球砾恢复系数e的理论计算和数值模拟结果如图5所示。

表4 不同弹性极限下球砾和平板的相关参数Table 4 Theparameters of the rock sphere and the plate with various elastic limit

理论计算和数值模拟表明,随着岩石球砾和平板的弹性极限增大,球砾的法向恢复系数增大,且两者基本呈线性关系,如图5所示。其原因在于弹性变形阶段和弹塑性变形阶球砾的部分动能转化为弹性应变能,且转化率随着材料弹性极限的增加而增大。在弹性回弹阶段,材料的弹性应变能再次转化为球砾的动能,因此球砾的恢复系数随着弹性极限的增加而增大。