0 引言
不同固体结构表面之间的弹塑性接触冲击现象是一个重要的力学问题，无论是在宏观还是微观尺度中均较为常见，特别是在机械系统的设计与应用中普遍存在[1-2]。接触冲击现象发生时间极短，量级在微秒与毫秒之间，并伴随有冲击过程中复杂的接触力变化和冲击后的速度变化[3]。近年来，随着固体力学与其他前沿交叉学科的快速发展，多体表面之间的接触冲击或碰撞所产生的复杂局部力学行为(如恢复系数、接触力、平均压力(即硬度)、接触区域(或接触半径)、弹塑性变形等[4])，受到了越来越多研究者和工程技术人员的广泛关注。
对于任何接触冲击事件，研究者主要关注冲击过程中和冲击后的动力学参数变化规律。一方面由于接触冲击发生时间极短，使得冲击过程中冲击物体(如刚性的球体或杆件等)的运动学特性需要经过先进的视觉技术(如高帧率高速摄像机)及数字图像处理技术(如数字图像相关法DIC、粒子图像测速法PIV等)才能获取[5-6]；另一方面接触冲击过程中的接触力与接触变形的关系、接触力和运动参数的求解以及测量与预测冲击对象的运动仍然面临着巨大的挑战[7]。
国内外许多研究者开展了大量有关接触冲击问题的研究，且已取得了较为显著的研究成果。从接触冲击变形角度，接触模型主要可分为压陷模型和扁平模型两大类。对于压陷模型，冲击物体被认为是刚性的，被冲击物体发生变形，代表性模型有Ye-Komvopoulos模型[8]、Kogut-Komvopoulos模型[9]、Brake模型[10-11]等；相反，对于扁平接触模型，被冲击物体被认为是刚性的，而冲击物体发生变形，代表性模型有Vu-Quoc模型[12]、Kogut-Etsion模型[13]、Jackson-Green模型[14]等。此外，还有一些假设冲击物体和被冲击物体均发生变形的模型，如Gheadnia模型[4,7]、Ma-Liu模型[15]等；也有学者从阻尼角度建立黏弹性/黏弹塑性接触模型来分析弹塑性接触冲击问题，如Hunt-Crossley模型[16]、Christoforou-Yigit模型[17]等。但总体来说，上述模型将接触冲击过程分为加载和卸载两个阶段，不同之处体现在对加载阶段的子阶段划分上，如Jackson-Green模型将加载阶段划分为弹性阶段和弹塑性阶段，Brake模型将加载阶段划分为弹性阶段、弹塑性阶段和塑性阶段，实际上对于塑性阶段，可认为是弹塑性阶段的一个极限；而卸载阶段即为恢复反弹阶段。其中弹性阶段和恢复阶段均被公认为满足赫兹接触理论[18]，而对于弹塑性阶段，不同接触模型的表达公式是不同的。此外，总结国内外研究现状，对两构件系统弹塑性接触冲击问题的研究工作可概括为正冲击(如球体[1,8-9,13-14,17]或杆件[3,5]垂直冲击固体平面，球体与球体正冲击[10-12,15])和倾斜冲击(如倾斜杆件冲击固体平面[2,5,7]，网球冲击倾斜球拍平面[6])。在目前的公开文献中，还没有关于垂直杆件冲击倾斜固体平面方面的研究报道。
针对垂直杆件冲击倾斜固体平面的弹塑性接触冲击问题，本文以末端圆形杆件为冲击对象、固体平面为被冲击对象，通过分析系统运动链关系来建立杆件动力学方程，并基于Brake模型建立了接触冲击各阶段瞬时接触力计算模型，旨在得到弹塑性接触冲击过程中接触力与接触变形量的关系，以精确描述其实际情况；借助商业数学软件MATLAB中ODE45函数对所建模型进行数值求解，分析了杆件动能、接触点位置和速度、杆件角度和角速度等的时变规律，恢复系数和永久变形量与初始冲击速度的关系，以及接触力与接触变形的关系；最后，对恢复系数和永久变形量两个参数进行了试验与仿真结果的对比分析，以验证所建模型的正确性。
1 动力学数学模型
图1为倾斜接触冲击类型示意图。对于杆件与固体平面的倾斜接触冲击，根据形成初始冲击角的方式可分为两种类型，分别是倾斜杆件冲击固体平面(图1a)和垂直杆件冲击倾斜固体平面(图1b)。
以固体平面为参照物，这两种冲击类型的区别是：前者在接触冲击之前，杆件相对于固体平面只有法向速度，切向速度为零；而对于后者，在接触冲击之前，杆件相对于固体平面存在初始的法向和切向速度。对于倾斜杆件冲击固体平面的冲击类型，文献[2]考虑振动响应的影响，基于假设模态法和广义动量定理建立了杆件的动力学模型，且采用数字图像处理技术，当初始冲击角θ=45°时，对恢复系数和冲击后的运动参数进行了试验与仿真对比分析；并根据Jackson-Green模型，将振动响应考虑到接触变形中，系统分析了弹塑性接触冲击过程中接触力的变化规律。在上述研究基础上，本文将对垂直杆件冲击倾斜固体平面的冲击类型进行分析研究，以期对倾斜弹塑性接触冲击问题获得更为全面的认识与了解。
(a) 倾斜杆件冲击固体平面
(b) 垂直杆件冲击倾斜固体平面
图1 倾斜接触冲击类型示意图
Fig.1 Schematic diagram of oblique contact-impact type
1.1 系统运动链
对图1b所示的系统运动链解释如下：冲击对象为末端圆形的杆件B，圆头半径为Rb，被冲击对象为固体倾斜平面S。固定绝对坐标系为
该坐标系下的正交基矢量为(i0,j0,k0)；固体倾斜平面上的局部坐标系为
其正交基矢量为(i1,j1,k1)，均满足笛卡儿坐标系；固体平面的倾斜角为θ(将其定义为初始接触冲击角)，冲击对象(杆件)的质心为C，长度为L，直径为d(d=2Rb)，接触点为E；重力为G，接触力为F，且将接触力F在
坐标系下分解为沿Y1正方向的Fn和沿X1负方向的Ft。
定义E点在下的广义坐标分别为q1、q2、q3 (均为时间t的函数qi(t)，i=1,2,3)，用来描述接触冲击过程中的瞬时位置。广义坐标q1表示E点与坐标系Y1轴的距离，q2表示E点与坐标系X1轴的距离，q3表示杆件与倾斜平面的夹角(与冲击角θ互为余角)；相应的广义速度u1、u2、u3用来描述其运动速度。
1.2 动力学方程
基矢量(i0,j0,k0)与基矢量(i1,j1,k1)的坐标变换关系为
式中,R10为转换矩阵。
接触点E的位置向量为
rE=q1i1+q2j1
(2)
质心点C的位置向量为
rC=[q1i1+(q2+Rb)j1]+
(3)
杆件的角速度向量和角加速度向量分别为
对式(2)、式(3)分别求时间t的导数，可得质心点C和接触点E的速度向量为
vE=vC+ω×(rE-rC)
(7)
质心点C和接触点E的加速度向量为
aE=aC+α×(rE-rC)+ω×[ω×(rE-rC)]
(9)
根据牛顿-欧拉(Newton-Euler)方程，得到冲击对象的运动方程:
式中,m为杆件质量；IC为杆件质心的转动惯量；μk为动摩擦系数。
接着，通过转换矩阵R10可求得下接触点E的位置向量为
2 弹塑性接触冲击模型
接触冲击过程可被细分为3个阶段：弹性阶段、弹塑性阶段和恢复阶段。采用赫兹接触理论对弹性阶段和恢复阶段进行分析；而对于弹塑性阶段，考虑碰撞两物体的材料属性，本节将基于Brake模型对瞬时接触力的弹塑性变化规律进行探讨。
2.1 接触冲击过程
图2为接触冲击过程中接触变形示意图。当接触点E开始与固体斜平面S发生接触冲击时弹性阶段开始，直到碰撞两物体中的弱材料发生屈服，屈服接触变形量和屈服接触力分别为δy、Fy，此时弹性阶段结束而弹塑性阶段开始；当接触点E的法向速度达到零时，接触变形达到最大变形量δm，对应最大接触力为Fm，此时弹塑性阶段结束且恢复反弹阶段开始；当接触力F为零时，变形区域达到永久变形量δr，接触冲击结束。
(a) 开始接触
(b) 屈服
(c) 最大变形
(d) 永久变形
图2 接触冲击过程中接触变形示意图
Fig.2 Schematic diagram of contact deflection during contact-impact
2.2 弹塑性接触力
系统的等效弹性模量E′和等效曲率半径R′分别由下式求得：
式中，Eb、Ef分别为两碰撞物体的弹性模量；μb、μf分别为两碰撞物体的泊松比；Rb、Rf分别为两碰撞物体半径，由于Rf=∞，则R′=Rb。
弹性阶段[18]接触力
式中,δ为接触变形量。
接触半径
弹性阶段开始于初始接触并直到碰撞两物体中的弱材料发生屈服现象。根据von Mises屈服准则，屈服接触变形量
式中,Fμ为应力场最大振幅;μ为碰撞两物体中弱材料的泊松比;z为接触表面之下的距离;σy为屈服应力。
对于弹塑性阶段，假设随塑性区域的增大弹性区域呈现衰减变化趋势。当δ>δy时，弹塑性阶段开始直到达到最大变形量δm，弹塑性阶段接触力
式中：nε为应变硬化指数;n为Meyer硬度指数;δp为塑性变形;p0为不考虑应变硬化的接触压强;H′为等效布氏硬度;Hb、Hf分别为两碰撞物体的布氏硬度；g为重力加速度;ap为塑性变形时的接触半径。
同理，弹塑性阶段的接触半径
对于恢复阶段(或卸载阶段)，假设没有反向屈服发生[19]，接触变形从最大变形开始直到发生永久不可恢复变形。恢复阶段接触力
式中,Rr为恢复阶段的曲率半径;δr为永久变形量。
相应地，恢复阶段的接触半径
3 数值求解与分析
采用ODE45函数对上节所建模型进行数值求解。冲击对象(末端圆形杆件)材质为20钢，长度L=150 mm，直径d=18 mm；被冲击对象(固体平面倾斜固定于大质量斜面上)材质为巴氏合金，牌号ZSnSb11Cu6。笔者前期研究结果[2-3]表明：对于上述弹塑性材料，在接触冲击过程中，巴氏合金为弱材料。表1所示为碰撞两物体的材料及几何属性。
表1 碰撞两物体的材料及几何属性
Tab.1 Material and geometric properties of two colliding bodies
属性参数20钢ZSnSb11Cu6密度(kg/m3)7 8507 737弹性模量(GPa)20648泊松比0.2800.285布氏硬度(MPa)1 528.80321.44屈服强度(MPa)24866Meyer硬度指数2.21半径(mm)9∞重力加速度g(m/s2)9.81
3.1 杆件动能
在垂直杆件与倾斜固体平面的接触冲击过程中，永久变形、摩擦等会损失一部分能量，这部分能量等于冲击前后的能量之差。杆件接触冲击过程中的动能T主要由两部分组成，分别为杆件的平动动能T1和杆件的转动动能T2，其表达式如下：
式中：
为杆件质心C在
下的速度。
在接触冲击过程中，杆件动能随时间的变化趋势见图3。分析图3可知：随着初始坠落高度Hi(即初始冲击速度)的增加，杆件与固体平面发生弹性接触冲击的时间缩短，且在此阶段，杆件动能T与其平动动能T1几乎完全一样，这表明弹性接触冲击的动能损失很小，可忽略不计；当进入弹塑性阶段后，在其接触冲击初期，杆件动能T与其平动动能T1呈现相同的变化规律，之后，接触变形逐渐出现塑性变形，杆件能量损失逐渐增大，而且随着塑性变形量的不断增大，杆件转动动能T2逐渐增大，当Hi>0.45 m时，伴随有应变硬化现象(图3c)；当进入恢复阶段瞬间，残余可恢复变形使各能量曲线的斜率均有突变发生。由此可得不同Hi下杆件动能呈现两种时变规律：当Hi≤0.45 m时，杆件动能随时间的增加而减小，见图3a、图3b；而当Hi>0.45 m时，在接触变形达到最大变形前，杆件动能将达到最小，随后出现反弹现象，见图3c。可见，随着Hi的增加，发生接触冲击的时间会缩短，且恢复反弹现象更加明显。
(a) Hi=0.25 m
(b) Hi=0.45 m
(c) Hi=0.65 m
图3 接触冲击过程中杆件动能随时间变化曲线
Fig. 3 Variation curves of kinetic energy of rod-like structure depending on time during contact-impact
3.2 接触点位置和速度
以初始坠落高度Hi=0.45 m的接触冲击事件为例，对接触冲击过程中接触点E的位置和速度变化规律进行详细分析与说明。
图4为接触点E的位置变化曲线。杆件与倾斜平面发生接触冲击的时间t=188.36 μs≈190 μs。当t≤0.104 μs时，杆件与倾斜平面发生弹性接触冲击，该阶段切向位置q1不变，法向位置q2线性增大到0.26 μm；当t≤172.19 μs时，杆件与倾斜平面发生弹塑性接触冲击，在该阶段随着时间的增加，法向位置q2呈现非线性递增，直到达到最低点(即接触变形达到最大变形量δm=274.44 μm)；当t>172.19 μs时，恢复阶段开始，直到该阶段的接触力Fr等于0，此时接触变形即为永久变形，即δr=268.27 μm。
(a) 切向位置q1
(b) 法向位置q2
图4 接触冲击过程中接触点的位置随时间变化曲线
Fig.4 Variation curve of position of contact point depending on time during contact-impact
图5为接触点E的速度变化曲线。接触点E切向速度和法向速度的曲线斜率均是在最大变形点(即t=172.19 μs)处不连续，这表明弹塑性阶段结束和恢复阶段开始时的杆件加速度存在突变现象，也间接表明塑性变形的不可恢复性；随时间t的增加，切向速度u1逐渐增大，而法向速度u2则先减小后增大；在接触冲击前后，接触点E的切向速度分别为-2.499 m/s、-3.747 m/s；同理可知，接触点E的法向速度在接触冲击前后分别为-2.499 m/s、0.598 m/s。
(a) 切向速度u1
(b) 法向速度u2
图5 接触冲击过程中接触点的速度随时间变化曲线
Fig.5 Variation curve of velocity of contact point depending on time during contact-impact
3.3 角度与角速度
图6、图7分别为接触冲击过程中杆件角度q3和角速度u3的变化曲线。接触冲击过程中，杆件角度q3的斜率在最大变形点处连续，而角速度u3的斜率在该处不连续，这表明弹塑性阶段的结束对杆件的角速度变化有显著影响；在接触冲击前后，杆件角度q3分别为0.785 4 rad、0.783 4 rad，两者相差0.002 rad；同理可知，杆件角速度u3在接触冲击前后分别为0、-31.253 rad/s。
图6 接触冲击过程中杆件角度随时间变化曲线
Fig.6 Variation curve of impact angle of rod-like structure with time during contact-impact
图7 接触冲击过程中杆件角速度随时间变化曲线
Fig.7 Variation curve of angular velocity of rod-like structure with time during contact-impact
3.4 恢复系数
恢复系数e是反映接触碰撞物体之间变形恢复能力的一个重要参数[20]，其定义为两碰撞物体接触冲击前后法向分量速度的比值，表达式如下：
其中，vEf、vEb分别为杆件接触点E冲击前后的速度向量。需要说明的是：由于被冲击物体为大质量固体平面，故其冲击前后的速度可认为是零。
图8 恢复系数随初始冲击速度变化曲线
Fig. 8 Variation curve of coefficient of restitutiondepending on initial impact velocity
图8所示为恢复系数与初始冲击速度的函数关系。可以看出，随初始冲击速度vi的增大，恢复系数e呈递减趋势，其曲线斜率随初始冲击速度的增大而减小；当vi>3.25 m/s时，e≈0.227 4，恢复系数基本趋于稳定且其平均绝对误差为0.006 5，最大相对误差为5.3%，平均相对误差为2.84%。同时，恢复系数也可间接反映出巴氏合金易变形的特性。
3.5 永久变形量
图9所示为永久变形量与初始冲击速度的关系。可以看出，随初始冲击速度vi的增大，永久变形量δr呈递增趋势；当vi≤3.25 m/s或vi≥4.25 m/s时，曲线呈近直线变化趋势；而当3.25 m/s<vi<4.25 m/s时，曲线斜率呈现先减小后增大、再减小后增大的趋势，δr的增大速率放缓，这可能是因为该区间发生了应变硬化现象。
图9 永久变形随初始冲击速度变化曲线
Fig. 9 Variation curve of permanent deformation depending on initial impact velocity
3.6 接触力与接触变形量
图10所示为接触冲击过程中接触力关于接触变形的函数变化关系。如图10a所示，随接触变形量δ的增大，接触力F逐渐增大，且当接触力达到最大接触力Fm时，接触变形量达到最大值δm，弹塑性阶段结束；之后，接触冲击进入恢复阶段，当接触力减小至0时，接触变形为永久变形。图10b为不同初始坠落高度下接触力与接触变形量的变化曲线，可以看出，随初始坠落高度Hi的增大，函数变化关系虽呈现相同的趋势，但最大接触力逐渐增大，且由最大变形量与永久变形量的差值Δδ(Δδ=δm-δr)可知，Δδ越小，永久变形量更接近于最大变形量。
4 试验验证
为验证所建模型的正确性，本文采用OLYMPUS激光共聚焦扫描显微镜(CLSM)对杆件在不同初始坠落高度下接触冲击后的永久变形量进行了测量，并将轮廓最深处的变形值与仿真结果进行了对比分析。此外，选择文献[3]的正冲击试验数据，对正冲击下恢复系数和永久变形量两个参数进行试验与仿真结果的对比分析。
(a) Hi=0.45 m
(b) 不同初始坠落高度
图10 接触力随接触变形量的变化曲线
Fig.10 Variation curves of contact force vs.contact deflection
4.1 倾斜冲击试验
每次接触冲击完成后，用记号笔对轮廓外形进行框选并标记，之后在CLSM下测量变形区的三维轮廓。由于表面粗糙度Ra不超过1.6 μm，故可忽略其影响，则横截面上的最大变形值可认为是永久变形量δr。图11所示为初始坠落高度0.20 m、初始冲击角度45°下接触冲击后变形区域的三维轮廓，可以看出，永久变形量δr为130.2 μm。
图11 倾斜冲击下永久变形量的三维轮廓
Fig.11 3D profile for permanent deformation under oblique impact
本文将图11中轮廓最深处的变形值与图9所示仿真结果进行了对比分析，如图12所示。由图12可知：永久变形量δr随初始冲击速度vi的增大而增大，且试验值略大于仿真值，其最大相对误差与平均相对误差分别为5.17%、2.22%。仿真结果与试验结果基本吻合，从而验证了所建模型的正确性。
图12 倾斜冲击下永久变形量的仿真与试验结果对比
Fig.12 Comparison between simulation and testing results for permanent deformation under oblique impact
4.2 正冲击试验
由图13可知：恢复系数e随初始坠落高度Hi的增大而减小，其最大相对误差与平均相对误差分别为11.82%、5.65%，且最大绝对误差为0.014。误差可能是杆件动力学建模过程中忽略了振动响应对运动参数的影响[2]。
图13 正冲击下恢复系数的仿真与试验结果对比
Fig.13 Comparison between simulation and testing results for coefficient of restitution under normal impact
图14 正冲击下永久变形量的仿真与试验结果对比
Fig.14 Comparison between simulation and testing results for permanent deformation under normal impact
分析图14可知：仿真结果与试验数据呈现相同的变化趋势，永久变形量δr随初始坠落高度Hi的增大而增大，且仿真值小于试验值，其最大相对误差与平均相对误差分别为7.12%、3.95%。仿真结果与试验结果基本吻合，从而验证了所建模型的正确性。
5 结论
(1) 随初始坠落高度Hi的增大，发生接触冲击的时间会缩短且恢复反弹现象更加明显；当Hi>0.45 m时，接触冲击过程中伴随有应变硬化现象。
(2) 以初始坠落高度Hi=0.45 m的接触冲击事件为例，杆件与倾斜平面发生接触冲击的时间约为190 μs。随初始冲击速度vi的增大，恢复系数e呈递减趋势，而永久变形量δr呈递增趋势；当vi>3.25 m/s时，e≈0.227 4，恢复系数基本趋于稳定；而当vi≤3.25 m/s或vi≥4.25 m/s时，δr呈近直线变化趋势，当3.25<vi<4.25 m/s时，δr的增大速率放缓。
(3) 随初始坠落高度的增大，接触力F与接触变形量δ的函数变化关系呈现相同变化趋势，但永久变形量δr更接近于最大变形量δm。
(4) 通过对恢复系数e和永久变形量δr两个参数的试验与仿真对比分析，得出仿真与试验结果基本吻合，验证了所建模型的正确性。
研究中忽略了振动响应对运动参数及接触变形的影响，在后续建立更加精确的动力学与接触力模型时应予以考虑。接触力模型中接触点实际曲率半径需通过更高帧率的高速摄像机来获取原始图像并采用数字图像处理方法求解，以进一步提高模型精度，再结合获取得冲击过程运动学参数，进而对恢复系数作出验证。
[1] DONG Xiaoyun,YIN Xiaochuan,DENG Qingming,et al.Local Contact Behavior between Elastic and Elastic-plastic Bodies[J]. International Journal of Solids & Structures,2018,150:22-39.
[2] 王尧,孟文俊,项丹,等. 基于数值图像处理方法的弹塑性倾斜接触冲击动力学试验与数值分析[J]. 机械工程学报,2019,55(1):81-90.
WANG Yao,MENG Wenjun,XIANG Dan,et al.Experimental and Numerical Analysis of the Elasto-plastic Oblique Contact-impact Dynamics Using Digital Image Processing Method[J]. Journal of Mechanical Engineering,2019,55(1):81-90.
[3] WANG Yao,LI Shujun,XIANG Dan,et al. Experimental and Theoretical Analyses of the Contact-impact Behavior of Babbitt ZChSnSb11-6[J]. Proceedings of the Institution of Mechanical Engineers,Part L:Journal of Materials Design & Applications,2019,233(7):1267-1276.
[4] GHAEDNIA H,WANG X Z,SAHA S,et al. A Review of Elastic-plastic Contact Mechanics[J]. Applied Mechanics Reviews,2017,69(6):060804.
[5] GHAEDNIA H,CERMIK O,MARGHITU D B,et al. Collision Measurements Using Digital Image Correlation Techniques[J].International Journal of Mechanical Sciences,2017,131:836-846.
[6] GHAEDNIA H,CERMIK O,MARGHITU D B. Experimental and Theoretical Study of the Oblique Impact of a Tennis Ball with a Racket[J]. Proceedings of the Institution of Mechanical Engineers,Part P:Journal of Sports Engineering & Technology,2015,229(3):149-158.
[7] GHEADNIA H,CERMIK O,MARGHITU D B. Experimental and Theoretical Analysis of the Elasto-plastic Oblique Impact of a Rod with a Flat[J].International Journal of Impact Engineering,2015,86:307-317.
[8] YE N, KOMVOPOULOS K. Indentation Analysis of Elastic-plastic Homogeneous and Layered Media:Criteria for Determining the Real Material Hardness[J].Journal of Tribology,2003,125(4):685-691.
[9] KOGUT L,KOMVOPOULOS K. Analysis of the Spherical Indentation Cycle for Elastic-perfectly Plastic Solids[J]. Journal of Materials Research,2004,19(12):3641-3653.
[10] BRAKE M R. An Analytical Elastic-perfectly Plastic Contact Model[J].International Journal of Solids and Structures,2012,49(22):3129-3141.
[11] BRAKE M R. An Analytical Elastic Plastic Contact Model with Strain Hardening and Frictional Effects for Normal and Oblique Impacts[J]. International Journal of Solids and Structures,2015,62:104-123.
[12] VU-QUOC L,ZHANG X,LESBURG L. A Normal Force-displacement Model for Contacting Spheres Accounting for Plastic Deformation:Force-driven Formulation[J]. Journal of Applied Mechanics,2000,67(2):363-371.
[13] KOGUT L,ETSION I. Elastic-plastic Contact Analysis of a Sphere and a Rigid Flat[J].Journal of Applied Mechanics,2002,69(5):657-662.
[14] JACKSON R L,GREEN I. A Statistical Model of Elasto-plastic Asperity Contact between Rough Surfaces[J].Tribology International,2006,39(9):906-914.
[15] MA Daolin,LIU Caishan. Contact Law and Coefficient of Restitution in Elasto-plastic Spheres[J]. Journal of Applied Mechanics,2015,82(12):121006.
[16] HUNT K H,CROSSLEY F R E. Coefficient of Restitution Interpreted as Damping in Vibroimpact[J].Journal of Applied Mechanics,1975,42(2):440-445.
[17] CHRISTOFOROU A P,YIGIT A S.Inelastic Impact and the Coefficient of Restitution[J]. Journal of Engineering Research,2017,4(4):194-213.
[18] HERTZ H. Üher Die Berührung Fester Elastischer Körper (on the Contact of Elastic Solids)[J].Journal Fur Die Reine und Angewandte Mathematik,1882,92:156-171.
[19] CHUN B K,JINN J T,LEE J K. Modeling the Bauschinger Effect for Sheet Metals,Part Ⅰ:Theory[J]. International Journal of Plasticity,2002,18(5/6):571-595.
[20] KHULIEF Y A. Modeling of Impact in Multibody Systems:an Overview[J]. Journal of Computational & Nonlinear Dynamics,2013,8(2):021012.
