基于自由落体的落石冲击土层的理论及LS-DYNA模拟研究

基于自由落体的落石冲击土层的理论及LS-DYNA模拟研究王星; 周天跃; 师江涛; 夏永旭
【期刊名称】《《北京交通大学学报》》
【年(卷),期】2019(043)004
【总页数】9页(P9-17)
【关键词】防护结构; 落石冲击; 理论研究; 数值模拟
【作者】王星; 周天跃; 师江涛; 夏永旭
【作者单位】长安大学公路学院西安710064; 陕西交通建设集团西安710086【正文语种】中文
【中图分类】TU45
滚石灾害已经成为目前中国的3大地质灾害问题之一，尤其是雨季时隧道洞口处
的滚石灾害情况更为严峻.各种滚石防护结构便应需而设，而防护结构设计的核心
即是要准确解析并计算出滚石的冲击力.
国内外学者针对滚石冲击力的计算做出了大量的研究：其中传统的Hertz接触理
论及Thornton弹塑性理论被广泛引用[1-2].文献[3]深入分析了物体间在静态接触、滑动、滚动、冲击等作用下的力学响应规律，并给出了相应的计算公式.文献[4]通
过采用室内试验的方法，得出了滚石冲击力的半经验半理论算法.文献[5]以室内试
验及动量定理为基础，给出了一种滚石冲击力的半理论半经验算法.文献[6]从现场
测试的角度出发，获得了滚石冲击力及冲击深度的量纲为1的计算表达式.文献[7]
研究了滚石直接冲击混凝土板时的力学响应情况.文献[8]对滚石冲击力进行了研究，并得出了滚石冲击力呈现出脉冲式变化规律的结论.文献[9]同样证明了滚石冲击过
程属于一种脉冲式的碰撞过程.文献[10]以理论计算为手段，分别推导出了滚石冲
击力的空腔膨胀算法以及能量守恒算法.文献[11]采用实验及动量定理的方法，归
纳出了滚石冲击力的计算方法.文献[12]同样在动量定理的基础上总结出了滚石冲
击力计算方法.文献[13]研究了滚石冲击荷载下垫层材料的动力响应情况.文献[14]
分析了滚石冲击柔性防护网的力学响应规律.文献[15]结合理论算法推导出了滚石
冲击防护结构垫层的冲击力.
针对滚石冲击力，国内外学者研究成果较多，意义深远.然而目前国内的相关设计
规范及其他滚石冲击力的算法，其所得结果多为滚石的平均冲击力，同时日本、瑞士学者所给出的经典算法均是基于试验所得出的半理论半经验算法，且算式中各字母变量的意义尚不十分明确.
为此，本文作者首先基于滚石冲击力脉冲式的变化特点，推导出滚石最大冲击力的脉冲式算法；其次再考虑滚石冲击垫层土体的弹塑性特点，得出滚石最大冲击力的弹塑性算法；再次通过采用LS-DYNA软件模拟分析了滚石的冲击力及侵彻深度随时间的变化规律，并归纳出滚动冲击力的LS-DYNA算法.最终将所获得的3种算
法与已有算法进行对比分析，以期为滚石防护设计提供参考.
1 滚石冲击力的脉冲算法
文献[11]提出了一种滚石冲击力的计算方法，定义为关宝树算法，其冲击力P及冲击时间tc的表达式分别为
(1)
(2)
式中：m为滚石的质量；a为滚石冲击加速度；Q为滚石的重量；g为重力加速度；h为垫层材料的厚度；He为滚石的下落高度.
从而有滚石冲击过程中的平均冲击力Pa为
(3)
由于Pa产生的冲量与实际冲击力P(t)产生的总冲量It相等，则有
(4)
式中：2t0为滚石冲击时间.
滚石冲击过程属于一种脉冲式行为，且冲击力随时间基本呈现出近似正弦变化的趋势[8-10,16]，故考虑P(t)近似满足图1中的正弦变化规律.
图1 滚石冲击力的变化曲线Fig.1 Variation curve of rockfall impact force
则有滚石实际冲击力P(t)与滚石最大冲击力Pmax之间的关系为
P(t)=Pmaxsin(ω0t)
(5)
根据正弦曲线的性质可得
(6)
式中：T为周期；ω0为角速度.
联立式(4)～(6)可得
(7)
联立式(3)及式(7)，便可求出滚石的最大冲击力Pmax为
(8)
考虑到滚石在冲击垫层土体后会产生一定的回弹，假定滚石的初始速度为v初，反弹的速度为v弹，v初与v弹的方向相反，则有恢复系数κ为
κ=v弹/v初
(9)
根据动量定理，则有考虑恢复系数影响后，对滚石最大冲击力Pmax进行修正的滚石冲击力P修的表达式为
P修t=mv初-(-mv弹)
(10)
联立式(9)～(10)，可以得出
P修t=(1+κ)mv初=(1+κ)Pmaxt
(11)
从而有
P修=(1+κ)Pmax
(12)
考虑恢复系数κ的影响，联立式(8)及式(12)，则可得出滚石冲击力的脉冲算法的表达式，定义为脉冲算法
P修
(13)
参数κ可参考文献[17]的实验结果进行取值.
进而求取滚石侵彻深度的理论计算值.假定v(t)为随冲击时间变化的动态滚石冲击速度，则结合牛顿第二定律有
(14)
联立式(5)～(6)、(14)，则有
(15)
对式(15)进行积分后，则有
(16)
(17)
式中：C为常数.
此处存在边界条件t=2t0,v=0，则有
(18)
从而可得出滚石的冲击深度为
(19)
对式(19)进行积分后，则有
(20)
联立式(2)、(13)、(20)，可得
(21)
2 滚石冲击力的弹塑性算法
滚石在冲击下部垫层土体后，最先与滚石接触的正下方土体区域将经历：接触-弹
性变形-弹性变形的临界状态-塑性变形，而距离滚石较远的区域则一直保持着弹性变形的状态.在冲击行为发生后，滚石与垫层土体之间将会形成球冠接触面，为了对滚石冲击力进行计算，需通过转换为滚石与垫层土体之间的接触面积大小进行考虑并计算，最终形成图2所示的接触面积的区域划分.图2中，滚石与垫层土体之间形成接触半径为fp的塑性区，而弹性接触面积外缘的半径为fr.
图2 垫层的弹塑性区域分布Fig. 2Regional distribution of elastoplasticity of cushion
滚石球体冲击垫层土体时的弹性接触压力Pe与变形量Se之间的关系式为
(22)
式中：Ra为等效半径；Ea为等效弹性模量.
Ra与Ea计算方式为
(23)
式中：υ1、υ2分别为滚石、垫层土体的泊松比； E1、E2分别为滚石、垫层土体的弹性模量；R1、R2分别为滚石、垫层土体的半径.
滚石冲击垫层土体后会形成一个球冠接触面，其接触面积Ac为
Ac=πf2=πR1S
(24)
式中： f为滚石的接触半径；S为对应于接触半径f下的侵彻深度.
从而有
f2=R1S
(25)
土体在弹性变形阶段的极限屈服应力py及对应的接触半径fy之间满足关系式[18]
(26)
py与屈服应力系数λ、材料屈服强度χM之间的关系，以及λ与υ2之间的关系为[19]
(27)
联立式(25)～(27)，便可以求出滚石的弹性极限冲击深度Sel为
(28)
联立式(22)及式(28)，即可算出滚石的极限弹性冲击力Pel为
(29)
结合图2可得出相应的力学平衡关系式为
(30)
式中：fp为土体垫层的塑性区半径.
从而弹塑性冲击条件下的力学平衡关系式为
(31)
式中：fmax为滚石冲击过程中的最大接触半径；fr为弹性接触半径；Pep为弹塑性冲击力.
联立式(21)、(25)及式(31)便可得出滚石冲击力的弹塑性算法的表达式，定义为弹塑性算法
Pep=Pel+πpyR1
(32)
3 LS-DYNA动态数值模拟
采用LS-DYNA软件模拟滚石冲击垫层土体的动态过程，并选取典型计算工况进行呈现：垫层的材料为砂土，尺寸为5 m(长)×5 m(宽)×1 m(高)，滚石半径为0.5 m，冲击速度为20 m/s.滚石与垫层土体间设置为“面-面接触”，垫层土体中心1.5 m×1.5 m的区域进行了网格加密.垫层土体的下底面设置了3个方向的全约束，滚石的冲击时间设置为0.06 s，计算步设置为200步，此次计算时间约为46 min.材料参数具体见表1.
表1 垫层土体材料参数Tab.1 Material parameters of cushion soil土体材料重度/(kN/m3)内摩擦角φ/(°)弹性模量E/MPa泊松比υ屈服强度χM/kPa碎石土20.436290.15239砂土15.228160.20152黏土12.72080.2495滚石25.0—320.12—
LS-DYNA的数值计算模型及滚石冲击垫层土体后0.015 s、0.045 s及0.06 s时
的Von mises应力云图如图3所示，根据图3中的应力等级色条便可得出土体垫层各部分的应力分布情况，同时也可判别相应的弹塑性区域.
图3 LS-DYNA数值模拟结果Fig.3 LS-DYNA numerical simulation results
滚石的冲击力随时间的变化曲线如图4所示.由图4可见，在滚石冲击0.004 s后，滚石的冲击力即达到峰值，而在此后的0.004～0.048 s期间，滚石的冲击力在逐
步的震荡过程中减小至0，滚石冲击过程完成.
图4 冲击力的变化曲线Fig. 4Variation curve of impact force
滚石冲击过程中侵彻深度随时间的变化曲线如图5所示.可见在0～0.034 s期间，滚石的侵彻深度将会随着冲击时间的增加而呈现出近抛物线的增加趋势，且在
0.034 s时形成最大冲击坑深0.301 m，而在0.034～0.048 s期间，产生了约2.6 cm的回弹，并最终形成了0.275 m的冲击坑深.
图5 侵彻深度的变化曲线Fig. 5Variation curve of penetration depth
滚石冲击结束后，垫层土体结构中轴线上的各节点的最终位移如图6所示.
图6 冲击坑的位移曲线Fig.6 Displacement curve of impact crater
可见冲击坑横剖面的左右两侧形状对称，冲击坑中心所形成的最大位移变形为0.275 m，并在坑的周围形成了约8 cm的隆起，垫层土体发生变形的主要区域是坑中心处半径为0.5 m圆的范围.
4 计算结果对比分析
为对各算法下滚石冲击力的计算结果进行对比分析，选取滚石半径为0.5 m，滚石质量为1 308.25 kg，垫层厚度取为1 m，其余的计算参数见表1及表2，计算工况中的各计算参数参考文献[20-21]进行选取，各种算法的计算表达如下[22].
表2 滚石冲击力计算工况表Tab.2 Working condition table for calculating rockfall impact force工况冲击速度/(m/s)下落高度/m冲击能量/J140.8210 466261.8423 548.5383.2741 8644105.1065 412.55127.3594 19461613.06167 45672020.41261 65082429.39376 776
1)Hertz算法：
(33)
2)瑞士算法：
(34)
3)日本算法：
(35)
4)隧道手册算法：
(36)
(37)
(38)
式中：ME为垫层土体弹性模量；ε为拉梅常数；v0为滚石初始冲击速度；vc为滚石冲击垫层的压缩波速.
碎石土、砂土、黏土垫层情况下垫层各算法所得出的滚石冲击力如图7～图9所示.通过对图7～图9进行分析，可得出如下结论：
图7 碎石土垫层的“冲击速度-冲击力”曲线Fig.7 “Impact speed - impact force” curves of gravel soil cushion
图8 砂土垫层的“冲击速度-冲击力”曲线Fig. 8 “Impact speed - impact force” curves of sandy soil cushion
图9 黏土垫层的“冲击速度-冲击力”曲线Fig. 9 “Impact speed - impact force” curves of clay cushion
1)Hertz算法所得出的滚石冲击力值无疑是偏大的，这主要是由于该算法是基于完全弹性理论所得出.瑞士算法结果略微偏大，因为该算法考虑了垫层材料的弹性模量的影响，并取弹性模量值的2/3次方作为因子，同时瑞士算法明显满足P碎石土>P砂土>P黏土，而当垫层土体为黏土材料时，瑞士算法冲击力已小于日本算法，由此可见，瑞士算法结果整体较大，同时垫层土体的弹性模量对瑞士算法结果影响较大.
2)日本算法所得结果整体较为偏大，且冲击力值整体较为稳定，是因该算法考虑将滚石重力值的2/3次方及下落高度的3/5次方作为滚石冲击力的计算因子，而拉梅常数值一般取为106.
3)关宝树算法所得出的计算结果整体偏小，这是由于其算法所得出的冲击力为整个冲击过程中的平均冲击力值，且实验过程中选用的是5 N及8 N的试验锤，或有可能存在一定的局限性，当试验中的滚石质量达到13 082.5 N时，关宝树算法可能会存在一定的偏差，从而导致计算结果偏低.
4)隧道手册算法所计算出的冲击力值普遍偏小，是因为在该算法中滚石冲击力时间的计算是通过2倍的垫层厚度与压缩波速之比所得出的，垫层厚度直接影响到冲击时间，故而冲击时间势必存在一定的偏差，最终致使冲击力与实际存在偏差. 5)LS-DYNA数值模拟结果与日本算法、瑞士算法、弹塑性算法及脉冲算法之间的吻合度较好，实现了各算法之间相互验证.脉冲算法整体与日本算法、瑞士算法及数值模拟算法保持一致，但仍略有偏小.弹塑性算法与日本算法、瑞士算法、数值模拟算法及脉冲算法整体一致性较好，从而证明了该算法的可靠性，这或可能是由于该算法充分考虑了滚石冲击过程中垫层土体的弹塑性变化情况，从而得出了较为准确的计算公式.
5 滚石冲击力的LS-DYNA算法
由于LS-DYNA数值模拟计算结果为滚石的最大冲击力，而隧道手册算法所得出的计算结果为滚石的平均冲击力，且隧道手册算法充分考虑了垫层土体的厚度、弹性模量、密度、泊松比及滚石的重量及下落高度等因素.故而尝试定义将滚石的最大冲击力与平均冲击力之比作为冲击力的放大系数，且给出放大系数的建议值，并通过计算获取碎石土、砂土及黏土垫层的滚石冲击力放大系数ζ与滚石冲击速度v0之间的关系，建立v0-ζ散点图，并进而拟合出相应的滚石冲击力放大系数曲线，具体见图10.
图10v0-ζ曲线Fig.10v0-ζ curves
对于碎石土、砂土及黏土垫层，其滚石冲击力的放大系数曲线表达式为
(39)
式中：ζ碎石、ζ砂土、ζ黏土分别为碎石土垫层、砂土垫层、黏土垫层时的冲击力放大系数.
鉴于式(39)中3条拟合曲线之间存在较小差异，故在综合考虑滚石冲击过程中的多方因素后，得出最终统一的冲击力放大系数ζe的表达式
ζe=0.655+0.241ln(v0+2)
(40)
联立式(36)～(38)、(40)，并考虑滚石恢复系数的影响，从而得出滚石冲击力为
(41)
在考虑滚石冲击角度及滚石自重因素的影响后，得出滚石最大冲击力的LS-DYNA 算法的最终表达式，定义为LS-DYNA算法
(42)
式中：θ为滚石冲击速度与垫层土体表面的夹角；β为滚石重力的方向与垫层表面的夹角.
为验证式(42)的准确性，将其计算结果与Hertz算法、日本算法等8种算法的计算结果进行对比，并分别考虑碎石土、砂土及黏土垫层的情况，对比结果如图11～13所示.
图11 碎石土垫层的验证结果Fig. 11Verification results of gravel soil cushion 图12 砂土垫层的验证结果Fig.12 Verification results of sandy soil cushion
图13 黏土垫层的验证结果Fig. 13Verification results of clay cushion
可见LS-DYNA算法结果与日本算法、瑞士算法、数值模拟结果、弹塑性算法及脉冲算法结果具备较好的一致性，脉冲算法及隧道手册算法结果仍略微偏小，而Hertz算法仍明显偏大，而关宝树算法无疑偏小，与之前的对比验证结果保持一致，其较好地验证了式(42)的准确性.
在得出滚石最大冲击力值后，此时便可结合文献[12]计算出滚石冲击力作用在防护结构顶板上的压力q，见图14.
图14 防护结构顶板上的冲击压力Fig. 14Impact pressure on roof of protective structure
由图14可知q为
(43)
假定滚石冲击力作用于防护结构顶板上部的垫层的面积为滚石球体的等效截面，即落石等效直径为2R1的区域.其中在隧道手册算法中中假定θ为30°，而在路基规范中假定θ为35°，此处由文献[22]可得
θ=45°-φ/2
(44)
式中：φ为垫层土体的内摩擦角.
6 结论
1)基于关宝树算法以及滚石冲击力的脉冲式变化规律，同时考虑了滚石冲击恢复系数的影响，得出了滚石冲击力以及侵彻深度的脉冲算法.结合滚石冲击过程中垫层
土体的弹塑性变化过程，建立了滚石冲击力的力学平衡关系式，并基于侵彻深度与接触半径之间的关系，推导出了滚石冲击力的弹塑性算法.
2)通过采用LS-DYNA软件分析研究了滚石冲击垫层土体的微观力学响应机理，得
出了滚石冲击力及侵彻深度随冲击时间的变化规律，并给出了最终塑性冲击坑的形状.
3)结合具体的滚石冲击计算工况，对各算法进行了对比分析.研究结果表明：完全弹性的Hertz算法结果明显偏大，日本算法、瑞士算法、数值模拟算法、脉冲算法、弹塑性算法结果吻合性较好，但相互之间存在少量差异，脉冲算法及隧道手册算法略微偏小，关宝树算法则明显偏小.
4)以数值模拟结果及隧道手册算法为基础，建立了 v0- ζ之间的拟合曲线，并在考虑κ、θ、滚石自重影响的基础上，得出了滚石最大冲击力的LS-DYNA算法，最终结合8种算法结果对LS-DYNA算法的适用性进行了验证.
参考文献
【相关文献】
[1] THORNTON C. Coefficient of restitution for collinear collisions of elastic-perfectly plastic spheres[J]. Journal of Applied Mechanics, 1997, 64(2): 383-386.
[2] THORNTON C, NING Z M. A theoretical model for the stick/bounce behaviour of adhesive, elastic-plastic spheres[J]. Powder Technology, 1998, 99(2): 154-162.
[3] JOHNSON H L. Contact mechanics [M].Cambridge:Cambridge University Press, 1985.
[4] KAWAHARA S, MURO T. Effects of dry density and thickness of sandy soil on impact response due to rockfall[J]. Journal of Terramechanics, 2006, 43(3): 329-340.
[5] LABIOUSE V, DESCOEUDRES F, MONTANI S. Experimental study of rock sheds impacted by rock blocks[J]. Structural Engineering International, 1996, 6(3): 171-176. [6] PICHLER B, HELLMICH C, MANG H A. Impact of rocks onto gravel design and evaluation of experiments[J]. International Journal of Impact Engineering, 2005, 31(5): 559-578.
[7] MOUGIN J P, PERROTIN P, MOMMESSIN M, et al. Rock fall impact on reinforced concrete slab: an experimental approach[J]. International Journal of Impact Engineering, 2005, 31(2): 169-183.
[8] YU B, YI W, ZHAO H B. Experimental study on the maximum impact force by rock fall[J].
Landslides, 2018, 15(2): 233-242.
[9] AZZONI A, BARBERA G L, ZANINETTI A. Analysis and prediction of rockfalls using a mathematical model[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1995, 32(7): 709-724.
[10] WANG X, XIA Y X, ZHOU T Y. Theoretical analysis of rockfall impacts on the soil cushion layer of protective structures[J]. Advances in Civil Engineering, 2018(2): 1-18. [11] 杨其新, 关宝树. 落石冲击力计算方法的试验研究[J]. 铁道学报, 1996, 18(1): 101-106. YANG Qixin, GUAN Baoshu. Test and research on calculating method of falling stone impulsive force[J]. Journal of the China Railway Society, 1996, 18(1): 101-106.(in Chinese) [12] 廖朝华, 郭小红. 公路隧道设计手册[M]. 北京: 人民交通出版社, 2012.
LIAO Chaohua, GUO Xiaohong. Manual for design of highway tunnel[M]. Beijing: China Communications Press, 2012.(in Chinese)
[13] 向波, 何思明, 庄卫林, 等. 滚石冲击荷载下垫层材料的动力响应研究[J]. 公路交通科技, 2015, 32(1): 45-49, 56.
XIANG Bo, HE Siming, ZHUANG Weilin, et al. Research on dynamic response of cushion under rockfall impact[J]. Journal of Highway and Transportation Research and Development, 2015, 32(1): 45-49, 56.(in Chinese)
[14] 周晓宇, 陈艾荣, 马如进. 滚石柔性防护网耗能规律数值模拟[J]. 长安大学学报(自然科学版), 2012, 32(6): 59-66.
ZHOU Xiaoyu, CHEN Airong, MA Rujin. Numerical simulation of energy dissipation mechanism on falling rocks protection nets[J]. Journal of Chang’an University(Natural Science Edition), 2012, 32(6): 59-66.(in Chinese)
[15] 王永东, 周天跃, 柴伦磊, 等. 基于能量分析的落石对棚洞垫层冲击力学研究[J]. 地下空间与工程学报, 2018, 14(增1): 148-153.
WANG Yongdong, ZHOU Tianyue, CHAI Lunlei, et al. Impact mechanics of rock-fall to shed cushion based on energy analysis[J]. Chinese Journal of Underground Space and Engineering, 2018, 14(S1): 148-153.(in Chinese)
[16] 何思明, 沈均, 罗渝, 等. 滚石坡面法向冲击动力响应特性研究[J]. 工程力学, 2011, 28(6): 118-124.
HE Siming, SHEN Jun, LUO Yu,et al. Study on the characteristics of normal impact of post-earthquake rock-fall on slope[J]. Engineering Mechanics, 2011, 28(6): 118-124.(in Chinese) [17] 袁进科, 黄润秋, 裴向军. 滚石冲击力测试研究[J]. 岩土力学, 2014, 35(1): 48-54.
YUAN Jinke, HUANG Runqiu, PEI Xiangjun. Tes tresearch on rockfall impact force[J]. Rock and Soil Mechanics, 2014, 35(1): 48-54.(in Chinese)
[18] 刘茂. 基于弹塑性修正Hertz接触理论的落石冲击力计算方法[J]. 中国地质灾害与防治学报, 2012, 23(3): 21-27.
LIU Mao. Calculation of impact force of falling rock based on elastoplastic modified Hertz
contact theory[J]. The Chinese Journal of Geological Hazard and Control, 2012, 23(3): 21-27.(in Chinese)
[19] BRIZMER V, KLIGERMAN Y, ETSION I. The effect of contact conditions and material properties on the elasticity terminus of a spherical contact[J]. International Journal of Solids and Structures, 2006, 43(18/19): 5736-5749.
[20] 叶四桥, 陈洪凯, 唐红梅. 落石冲击力计算方法[J]. 中国铁道科学, 2010, 31(6): 56-62.
YE Siqiao, CHEN Hongkai, TANG Hongmei. The calculation method for the impact force of the rockfall[J]. China Railway Science, 2010, 31(6): 56-62.(in Chinese)
[21] 叶四桥. 隧道洞口段落石灾害研究与防治[D]. 成都: 西南交通大学, 2008.
YE Siqiao. Research and mitigation of rockfall hazards at tunnel entrance and exit[D]. Chengdu: Southwest Jiaotong University, 2008.(in Chinese)
[22] 叶四桥, 陈洪凯, 唐红梅. 落石冲击力计算方法的比较研究[J]. 水文地质工程地质, 2010, 37(2): 59-64.
YE Siqiao, CHEN Hongkai, TANG Hongmei. Comparative research on impact force calculation methods for rockfalls[J]. Hydrogeology & Engineering Geology, 2010, 37(2): 59-64.(in Chinese)。