多体系统碰撞动力学中接触力模型的研究进展1)

王庚祥 马道林 刘洋 刘才山,2)

*(西安建筑科技大学机电工程学院,西安 710055)

†(北京大学工学院,湍流与复杂系统国家重点实验室,北京 100871)

**(上海交通大学船舶海洋与建筑工程学院,上海 200240)

††(埃克塞特大学工程、数学和物理科学学院,英国埃克塞特 EX44 QF)

引言

接触碰撞行为在工程机械[1]、生物医学仪器与航空航天等领域[2]的多体系统中无处不在[3-4].典型场景如飞机起落架轮胎与地面接触、车辆碰撞测试、嫦娥五号着陆与空间飞行器之间的对接过程等.随着工程机械系统不断向轻质、高速、重载与精密方向发展[5-6],碰撞行为在短时间内引起的局部弹塑性变形与剧烈冲击力,愈发容易导致多体系统的结构损伤与工作状态失稳;对高端工程装备的安全使用造成了严重的威胁[7-8].因此,有必要研究碰撞行为引起多体系统冲击与振动的物理机制,分析碰撞与系统整体动态性能之间的耦合关系[9],为工程机械系统的安全稳定工作提供理论依据[10].然而,多体系统中碰撞行为的研究与接触力模型的发展历程密切相关[11-12],碰撞力作为衡量碰撞行为是否对多体系统造成结构损伤的重要指标高度依赖于碰撞力模型的精确性[13-14].

至目前为止,对动态接触碰撞行为计算的模型主要分为两类(如图1):(1)基于Hertz 接触模型发展的考虑能量耗散的连续接触力模型[11,15-18],该类模型又可分为阻尼因子中分母包含与非包含初始碰撞速度的接触力模型;(2)通过静态压力试验获得的准静态弹塑性接触力模型[11-12,19-20],该类模型从初始不考虑弹塑性过渡阶段的非连续接触模型发展为连续的准静态接触力模型.本文就两类模型与相关的关键碰撞参数的发展历程进行了着重介绍.

图1 接触力模型的分类Fig.1 Classification of contact force models

1 等效连续接触力模型

碰撞行为是典型的非线性与不连续的非光滑系统[21],其中由于碰撞引起系统速度的跳跃具有典型的非光滑特征[22].当前描述碰撞体接触行为的概念分为[17]:非光滑动力学方法(nonsmooth dynamics formulation)[23]和正则化方法(regularized approach)[24].其中基于几何约束的非光滑动力学方法认为接触体在碰撞过程中不发生变形,因此该方法也称之为刚体方法[25].相反地,正则化方法认为碰撞体在接触区域允许发生变形且接触力是变形量的函数,所以该方法也称之为罚方法或者连续分析方法[26].

非光滑动力学方法中最常用的接触建模方法是线性补偿技术[27](linear complementary approach),该方法的核心思想是接触刚体之间单边约束的显式表达[28].单边约束被用来处理接触问题的主要特征是利用了刚体不可穿透性的概念[29],意味着接触体在碰撞过程中其接触点不能越过碰撞体的边界[30];防止碰撞发生变形保证其接触力为正值[31].然而,该方法在处理考虑摩擦行为的接触问题时可能会出现多解和无解的可能性,或者可能出现能量不守恒的情况[17];另外,该方法不利于多体系统碰撞动力学代码的通用化.相反,正则化方法很好地避免了刚体方法遇到的诸多问题,该方法认为接触体的每个接触区域都覆盖着一些分散在其表面的弹簧阻尼元件[32],该元件的变形与刚度大小依赖于碰撞体的几何与材料属性以及相对渗透深度[21],因此该方法也称之为柔顺接触力模型(compliant contact force model).该方法的主要缺点在于难以确定刚度系数与阻尼系数以及能量指数等接触参数,以及在多体系统动力学仿真过程中引入大量高频信息影响其求解效率[33].但是该方法的数学表达形式简单直接且是渗透深度的连续函数,便于多体系统碰撞动力学的程序通用化,因此该方法被广泛应用于多体系统动力学软件[22,34-35],包括ADAMS,RecurDyn,EDEM 等.

连续接触力模型的原型为Hertz 接触力模型[36-37],考虑到Hertz 模型不能描述碰撞过程中的能量耗散;阻碍了Hertz 接触力模型在多体系统碰撞动力学中的进一步应用[38-40].为此,Kelvin 与Voigt 人为地在Hertz 接触力模型中引入了阻尼项描述碰撞过程中的能量耗散,其基本形式为F=(K为Hertz接触刚度,δ为相对接触变形量,D为阻尼系数,为相对碰撞速度),其中为弹性力项,为阻尼力项.

为了更加精确地计算多体系统中的碰撞行为,众多学者在推导阻尼系数的过程中采用了不同的近似与求解方法[21,34,41-45],相应地获得了不同性质的阻尼系数.阻尼[46]是振动过程中能量耗散的一种度量,也是理解系统动力行为的重要组成部分,对抗震结构的设计至关重要[47].阻尼一般分为两种:(1)与频率相关的阻尼称为黏性阻尼[48-50];(2)与频率无关的阻尼称为迟滞阻尼[45].

另外,考虑到Hertz 接触刚度独立于接触变形量,为了建立Hertz 刚度与相对接触变形量之间的关系,修改了接触力模型中的Hertz 刚度而忽略了能量指数与刚度系数的关系[41,51],造成了一系列的数值仿真问题.下面的内容集中讨论了接触力模型中人为阻尼的类型与推导过程,阐述了能量指数与刚度系数之间的关系;最后总结了等效连续接触力模型在计算碰撞行为时所面临的挑战.

1.1 等效迟滞阻尼

1881 年Hertz[52]在研究两个完全弹性体的接触应力时提出了著名的Hertz 接触理论,其接触力F是变形量的非线性函数

其中,δ为接触变形量;n为能量指数;K为Hertz 接触刚度,其表达式为

其中,Ei和Ej为碰撞体的杨氏模量;Ri和Rj为碰撞体的接触半径;vi和vj为碰撞体的泊松比.值得注意的是Hertz 接触刚度不同于胡克定律中弹簧的刚度,该刚度独立于接触变形量且完全依赖于接触体的几何与材料属性(如式(2))[51],且该参数的量纲为N/m1.5,而不是胡克定律中弹簧刚度系数的量纲N/m.

众所周知,碰撞过程实际上是碰撞体之间能量相互转化与耗散的过程,针对纯弹性碰撞行为的Hertz 接触理论并不能满足实际的工程需求,其原因在于式(1)并不能描述碰撞规程中的能量损耗(碰撞过程中,接触体的压缩路径与恢复路径相同)[18];于是为了描述碰撞行为的动态接触过程以及碰撞导致的能量耗散,人为地在原始Hertz 接触力模型(1)的基础上引入表示能量耗散的阻尼项,将Hertz 接触力模型扩展为可表示能量耗散的动态连续接触力模型[53].该方法不仅能描述完整的动态接触过程与能量损耗,并且建立了接触力、相对变形量和变形速度之间的耦合关系[54],以及能快速获得碰撞行为随时间的动态变化过程,避免恢复系数法不能计算接触进程的缺点[55];另外,为了达到精确计算碰撞过程中能量耗散的目的,众多国内外学者利用不同的近似推导方法从阻尼系数入手开展了一系列的研究工作[17].

为了在Hertz 接触模型中引入人为的阻尼因子描述碰撞过程中的能量耗散,1960年Kelvin 与Voigt[17]利用线性弹簧阻尼模型计算了接触过程中的能量耗散

其中式(3) 右边的第一项线弹性项(胡克定律),K1为常系数刚度(不同于Hertz 接触刚度);D1为常系数阻尼;为相对碰撞速度.然而,该模型在接触开始与结束阶段由于接触速度不为零[56],导致其接触力模型中的阻尼分量一直存在并引起接触力的不连续,甚至在恢复阶段结束时侵入深度为零而相对接触速度为负导致其接触力也为负,显然这种情况违反了接触力学的物理意义(碰撞体之间不可能存在拉力)[24,26];以上原因导致Kelvin-Voigt 模型产生的迟滞环不是一个完整的封闭区域,丧失了描述完整碰撞过程中能量耗散的能力.虽然该模型有以上缺陷,但在某些碰撞场景也有成功的应用[56-58],例如Dubowsky等[57]利用该模型计算了三维间隙关节元素之间的接触力,并提供了相关实验数据验证了该模型在一维振动冲击系统中的成功使用[58].类似地,Rogers等[59]同样将该模型应用在考虑间隙关节的多体系统动力学仿真中,并认为该模型在材料阻尼较小时能近似地表现出黏弹性的属性.Taylor等[60]利用该模型计算车辆轮胎与地面的法向接触力.以上学者均认为接触力模型的改进是精确计算多体系统碰撞动力学响应的基础.

为了克服Kelvin-Voigt 模型在压缩起始与恢复结束时刻出现的大于零与小于零的碰撞力,导致阻尼环呈现出不连续的状态,文献[61,17]基于Hertz 理论,将阻尼系数描述成迟滞阻尼因子与接触变形量之间的函数(D=χδn,χ为迟滞阻尼因子),避免了Kelvin-Voigt 模型在压缩开始阶段接触阻尼力的存在,同时也消除了恢复阶段结束时相对接触速度不为零引起的缺点[61],保证了接触力模型的连续性与迟滞环封闭的特点,其表达式为

其中m为接触参数.自此开始,在随后的近40多年时间里,国内外众多学者在追求更加精确的迟滞阻尼因子的这条路上正式拉开了序幕[11,16,22,51,60-62];实际上,连续接触力模型的发展历程本质上是阻尼因子的发展历史[34],因为,式(4)右边的第一项弹性力项(原始的Hertz 接触模型)保持不变,第二项阻尼项的发展目的是为了准确描述碰撞过程中的能量耗散.

考虑到Hunt-Crossley 模型的特殊性[61],这里有必要详细介绍该模型的推导过程,两个小球之间发生碰撞过程如图2.

图2 碰撞体接触过程Fig.2 Contact process between two contact bodies

首先,根据碰撞过程能量守恒原则,碰撞过程中耗散的能量可表示为

其中,T(-)为碰撞前系统动能;T(+)为碰撞后系统动能;m1与m2为碰撞体的质量;v1i与v2i为碰撞前接触体的初始速度;v1j与v2j为碰撞后接触体的速度.

其次,根据碰撞前后动量守恒原则有如下关系式

该模型在推导阻尼因子过程中最大的特点在于:假设衡量能量损耗的恢复系数为初始相对碰撞速度的线性函数[15]

其中 α为常数,一般取0.08 s/m 至0.32 s/m.

因此,将式(9)代入式(8),并考虑到 α 小于1,在忽略 α 的高阶项后式(8)重新写为

以上式(5)～式(10)是根据碰撞前后能量守恒的原则获得了碰撞过程的能量耗散,该过程与碰撞过程中是否发生弹塑性变形无关,即对弹性与弹塑性碰撞行为均适用.

另外,对式(4)中的阻尼力积分被视为碰撞过程中耗散的能量

为了获得阻尼项耗散的能量,式(11)中碰撞中间过程对应的相对变形速度需要表示成相对变形量的函数以确保式(11)能被积分.由此,该模型认为在压缩结束阶段碰撞体的初始动能被分为两部分:一部分被转化为弹性势能,一部分为碰撞后系统的动能;其表达式为

其中v12为压缩结束时碰撞体的共同接触速度.此时,同样根据动量守恒原则,其共同接触速度可表示为

将式(13) 代入式(12),压缩结束阶段的应变能可写为

另外,压缩结束阶段的应变能也可通过弹性力做功表示

联立式(14)与式(15),可以获得相对初始速度与最大相对接触变形量之间的关系

根据该式,式(10)改写成

然而,对于中间接触过程(相对变形量0<δ<δmax),根据能量守恒,其接触过程被认为是初始碰撞动能不断向势能转化的过程

将式(14)与式(15)代入式(18),中间接触过程对应的相对变形速度可表示为

将式(19)代入式(11),由此可以通过对阻尼项积分获得碰撞过程中的能量损耗

联立通过碰撞过程能量守恒获得的式(17)与通过对阻尼项积分获得的式(20),再利用恢复系数是初始相对碰撞速度的线性函数的假设式(9),Hunt-Crossley 模型对应的迟滞阻尼因子可表示为

迟滞阻尼因子是Hertz 接触刚度、恢复系数与初始相对碰撞速度的函数,该迟滞阻尼因子推导的关键步骤有4 点:(1) 假设恢复系数是相对初始碰撞速度的线性函数[63];(2) 在假设条件(1)的情况下,根据能量守恒原则将耗散能量表示为相对初始碰撞速度的函数[64];(3) 假设在压缩结束时刻,系统初始动能被转化为碰撞后系统的动能与势能[34];(4) 将中间碰撞过程对应的相对变形速度表示为相对变形量的函数,以便对阻尼项积分时获得碰撞过程的能量耗散[11].通过以上两种不同思路推导的碰撞过程中的能量耗散完全等价,以此获得连续接触力模型中的迟滞阻尼因子.

以Hunt-Crossley 模型[34]的推导过程为基础,其他8 种迟滞阻尼因子——包括Lankarani-Nikravesh模型[65-66]、Zhiying-Qishao 模型[67]、Flores 等模型[68](首次提出该模型中的迟滞阻尼因子的是文献[69])、Hu-Guo 模型[70]、Shen 等模型[71]、Safaeifar-Farshidianfar 模型[21]、Zhang 等模型[72-73]与Zhao 等模型[22]——在推导过程中采用了不同的假设或近似计算的方法获得的阻尼表达式不尽相同;例如Lankarani-Nikravesh 模型[65-66]认为恢复系数为常数,并未假设其为相对初始碰撞速度的函数以此获得另外的迟滞阻尼因子.Flores 等模型[68-69]则认为在压缩结束阶段系统的初始动能分为三个部分:一部分为碰撞后的系统动能,一部分为损耗的能量,另一部分能量被转化为势能;以此获得新的迟滞阻尼因子.

除此之外,其他接触力模型在推导过程中将碰撞行为等效为非受迫的弹簧阻尼系统,通过近似求解非线性振动方程的方式获得阻尼因子[26,74-75]

其中,D为阻尼系数;为相对变形加速度.由于该方程没有解析解,一般通过近似求解的方式进行求解,并在此过程中假设恢复系数为初始碰撞速度的函数,或者直接将迟滞阻尼因子假设为的形式(其中 β为恢复系数的函数).根据此方法推导迟滞阻尼因子的模型包括Herbert-McWhannell 模型[76]、Gonthier 等模型[75]、Carvalho-Martins 模型[15]、Gharib-Hurmuzlu 模型[77]与Hu 等模型[78].此外,还有一种颇具争议的推导方法,即Lee-Wang 模型[79]将碰撞过程视为线性的弹簧阻尼系统,通过线性振动方程的解析解来获得阻尼因子;然而,在其接触力模型中仍然保留了非线性的Hertz 接触刚度,前后存在明显矛盾[68,80].感兴趣的读者可以查阅相关文献[79],表1 中给出了15 种连续接触力模型的表达式[45],图3 描述了表1 中部分连续接触力模型中接触力与变形量的变化趋势.

图3 常见连续接触力模型的力与变形量之间的关系[51]Fig.3 Force-deformation relationship of the common contact force models[51]

表1 中的连续接触力模型有以下特点:(1)能量指数n与接触参数m均等于1.5[44];(2)所有迟滞阻尼因子均是接触刚度、恢复系数与相对初始碰撞速度的函数[43];(3)所有迟滞阻尼因子的分母中均含有初始相对碰撞速度[34,81].该类连续接触力模型在计算具有初始速度接触体的碰撞行为时,由于其数学结构简洁[34],并在计算碰撞问题时只需要判断接触是否发生,不需要判断接触行为处于压缩阶段还是恢复阶段[82-83],其计算流程相对简单;因此该类模型被广泛应用于多体系统动力学软件中(ADAMS 与RecurDyn 等)[84].然而,该类模型也存在不可调节的缺陷,由于该类模型中阻尼项分母中含有初始碰撞速度,限制了该模型对颗粒物质中孤立波传播特性刻画的能力[85-88].因为颗粒物质的初始状态一般处于静止状态,以至于颗粒物质间的相对碰撞速度等于零[89];该接触状态导致表1 中连续接触力模型的阻尼项无穷大,这显然违反了颗粒物质间的物理属性[8];同时也给颗粒物质的数值仿真引入了数值奇异问题[8,90-92].从表1 中连续接触力模型的迟滞阻尼因子推导过程发现,导致阻尼项分母中含有初始相对碰撞速度的原因有以下3 点:(1)利用能量守恒原则推导碰撞过程的能量耗散[93];(2)假设恢复系数是初始相对碰撞速度的函数[61];(3)直接将阻尼系数假设成分母中含有初始碰撞速度的形式[78].以上三类因素均直接或者间接与初始相对碰撞速度相关,以至于推导过程中在迟滞阻尼因子分母中引入了初始碰撞速度.

1.2 等效黏性阻尼

为了避免阻尼项分母中含有相对初始碰撞速度带来的数值奇异问题,在推导阻尼因子时绕开上述相关假设[94]与避免借助碰撞前后的能量守恒原则;将碰撞行为等效为单自由度的弹簧-阻尼模型(如图4),并直接对单自由度欠阻尼非受迫振动方程求解,就能获得分母中不含初始碰撞速度的黏性阻尼项[29,32,68].

图4 弹簧-阻尼模型Fig.4 The spring-damper model

该类模型分为线性和非线性两类系统:(1)假设系统的刚度系数为常数,其阻尼系数通过求解单自由度线性振动方程获得,该类具有代表性的接触力模型为Maxwell 模型[8];该模型与Kelvin-Voigt 模型类似,但是Maxwell 模型的阻尼系数是恢复系数的函数;(2)考虑Hertz 接触模型的非线性特征,基于非线性弹性碰撞行为特征建立的单自由度非线性振动系统的动力学模型(如式(22)),通过近似求解该动力学方程获得接触力模型中的阻尼因子;或者直接利用经验值确定其阻尼系数[74-75,95].为了清楚地诠释该类模型与表1 中模型的区别,有必要分别介绍线性与非线性系统中阻尼系数的推导过程,为该文后续的接触力模型分析奠定基础.

(1) 线性碰撞行为对应的振动方程

需要注意的是:参数KL是线性弹簧的刚度,而不是Hertz 接触刚度,所以方程(23)中弹性力项KLδ的变形量指数为1,其具体原因将在4.1 小节解释接触力模型中刚度系数与指数系数关系的重要性.

式(23)作为常见的欠阻尼非受迫振动系统的动力学方程,其解析解为

其中系统的幅值和相角可由振动系统的初始条件获得

将式(25)代入式(24),线性系统的解析解为

整个碰撞过程的持续时间为

另外,根据式(26),可以获得振动系统的变形速度

因此,当碰撞行为结束时变形量完全恢复,将式(27)代入式(28),其相对变形速度为

最后,根据Newton 恢复系数的定义

将式(29)与式(27)代入式(30),线性振动系统的阻尼系数为

Maxwell 模型的数学表达式为

(2) 非线性碰撞行为对应的振动方程

其中 α为恢复系数的函数.

为了获得相对变形速度,对上式进行无量纲化,即定义

无量纲化后,式(35)改写为

通过对上式积分可以获得相对碰撞速度,以此导出阻尼系数

式(39) 正是被广泛应用于颗粒物质仿真软件(EDEM)中的连续接触力模型[96].显然,该模型中阻尼项不涉及初始相对碰撞速度引起的数值奇异问题,有效避免了迟滞阻尼针对颗粒物质碰撞行为仿真的缺陷[97].其他接触力模型的阻尼系数推导方式均借鉴了上述线性或非线性动力学方程的求解策略,采用了不同的近似手段和假设发展了其他6 种类似的模型.表2 中给出了接触力模型的阻尼项分母中不含初始碰撞速度的接触力模型,此类模型主要应用于颗粒物质的数值仿真[86,89,98].

表2 黏性接触力模型(阻尼项分母中不含初始相对碰撞速度)Table 2 Contact force models with viscous damping factors(the denominators of the damping force do not include the initial impact velocity)

1.3 等效接触力模型中阻尼系数的讨论

显然,1.1 节中通过碰撞过程中能量与动量守恒原则获得阻尼与频率无关,但可保证碰撞前后系统的能量守恒[34];1.2 节中通过求解振动方程获得与频率相关的黏性阻尼,但无法保证碰撞前后系统的能量守恒[11,101].两类阻尼在碰撞模型中不仅推导方式与使用背景不同,而且其展现的接触动力学行为也不尽相同[102-103].为了详细说明两类阻尼的不同之处,假设两个相同的钢球在拥有不同初始速度(分别为0与8 m/s)下发生碰撞行为,其中小球的半径为0.02 m,杨氏模量为2.07×1011Pa,泊松比为0.3.图5中给出了两种不同接触力模型(Hunt-Crossley模型中的阻尼为迟滞阻尼,EDEM 软件中的接触力模型[104]对应其黏性阻尼)中阻尼力的对比情况;其中迟滞阻尼随着接触变形量的增大缓慢增加,当碰撞行为处于压缩阶段的结束阶段时,其碰撞速度趋近于0以至于其阻尼力从最大值(8000N)急剧减小至0N;在恢复起始阶段虽然相对接触速度较小但其接触变形量较大,以至于相对接触速度发生较小变化其阻尼力将急剧增加;另外,在恢复的结束阶段其接触变形量为0以至于其阻尼力也为0.因此,包含迟滞阻尼的接触力模型中描述能量损耗的阻尼环不会出现“拉力”区域(见图5),其主要原因在于恢复结束阶段其阻尼力快速趋近于0.相反,EDEM 软件中的接触力模型中黏性阻尼力虽然在整体趋势上与迟滞阻尼力保持一致(见图6),符合碰撞行为中阻尼项的表现形式;但是在细节上存在明显的区别.在压缩起始阶段与恢复结束阶段,黏性阻尼力戏剧化地增加,其原因在于黏性阻尼力中δ0.25在接触变形量较小时其值较大.因此黏性阻尼力在压缩起始阶段与恢复结束阶段明显大于迟滞阻尼,这也是为什么黏性阻尼对应的阻尼环会在恢复临近结束阶段出现“拉力”区域(见图6).黏性阻尼力在经过极速增加区域后,当接触变形量逐渐增大时,阻尼力缓慢增加;即使在靠近压缩结束阶段时其阻尼力也是缓慢减小,与迟滞阻尼力急剧减小不同.另外,关于“拉力”区域另一种直观的解释为:该区域对应的接触变形量大于0,但其碰撞力为负值;原因在于恢复结束阶段其阻尼力大于弹性力.由于该类模型主要用于颗粒物质之间的接触碰撞行为描述[105-107],由于该区域在整个阻尼环中只出现在恢复的结束阶段对大量颗粒物质之间的碰撞行为与颗粒之间孤立波的传播基本不产生影响,这也正是该模型在EDEM 软件[104]中被广泛使用的原因.图5 中两类接触力模型中力与接触变形量之间的关系在整体趋势上保持一致,但包含黏性阻尼的接触力小于其包含迟滞阻尼的接触力,其原因在于最大迟滞阻尼力大于最大黏性阻尼力(因为两类接触力模型中Hertz 弹性项相同).

图5 两类阻尼的比较Fig.5 Comparative analysis between two different damping systems

图6 拥有不同阻尼类型的接触力模型Fig.6 Contact force models with different damping systems

1.4 能量指数与接触刚度之间的关系

关于能量指数与接触刚度的关系往往被大家所忽略而造成不易察觉的错误.大家通常认为Hertz 模型中能量指数等于1时[51],其数学表达式就是弹性力学中的胡克定律;其实两者之间并没有直接的联系,并且该观点从量纲的角度出发就不正确.因为Hertz 刚度的单位为N/m1.5,当能量指数等于1 时,最后获得接触力的量纲不是N,明显违反了力学常识.这也是为什么Hertz 接触模型中能量指数等于1.5,说明接触体为金属材料且接触应力呈抛物线分布(这里需要注意的是:Hertz 接触刚度本身是一个特殊的物理量,有别于传统的刚度系数(变形量相关的系数);而Hertz 接触刚度完全取决于接触体的几何与材料属性且独立于接触变形量).因此,能量指数的大小必须与接触刚度的量纲保持一致,否则将造成不易察觉的错误,例如表2 中的Kuwabara-Kono 模型[99]与Lee-Wang 模型[79],就忽视了Hertz接触刚度量纲与能量指数的关系;很明显,Kuwabara-Kono 模型中阻尼力的量纲为N/s;Lee-Wang 模型中阻尼力的量纲为N/m0.25.以上两种模型已经在颗粒物质的动态仿真中有所应用,并取得了 “合理”的仿真结果,但依然不能避免量纲不正确的缺点.

那么在何种情况下,基于Hertz 理论扩展的连续接触力模型中的能量指数会发生变化？一般有两种情况能量指数会发生变化:(1) 接触材料不再是金属材料,而是玻璃或者高分子聚合物等材料,直接影响了接触应力的分布形式;(2)其接触形式发生变化,由点接触变化为线接触或者面接触.然而,至目前为止就表1 与表2 中的接触力模型而言,这种情况还未发生[51];也就是说,只要是基于Hertz 接触模型拓展的考虑能量耗散的接触力模型,其能量指数必定等于1.5;除非改变接触刚度的数学表达式,使得接触刚度系数的单位不再是N/m1.5.

另外,为什么要一直强调,接触刚度与能量指数的关系很容易被忽视,并且导致不易察觉的错误.因为在数值仿真当中,以考虑间隙关节的多体系统动力学仿真为例,间隙关节元素之间的碰撞行为属于典型的多重碰撞与多重压缩问题,一般采用连续接触力模型计算其碰撞力;其中间隙关节元素之间的相对接触变形量一般在[1×10-10,1×10-5] m,其能量指数的大小决定了接触力的大小处于不同的数量级;直接影响间隙关节元素之间的接触碰撞判定条件.例如文献[51,108-109]中接触模型的刚度系数单位为N/m2,但假设其能量指数仍然等于1.5,此时仍能获得“合理”的仿真结果,且该类错误很难发现.如果为了配合接触刚度的量纲,假设能量指数等于2;此时其接触力将处于另一个量级甚至约等于零,这将原本处于接触的碰撞体根据其仿真结果判定为未接触,给考虑间隙关节的多体系统动力学仿真带来了严重的数值仿真问题;其根本原因在于能量指数的大小直接改变了 δn的数量级,具体细节可以参考文献[51].

造成以上错误的根本原因在于对连续接触力模型中Hertz 接触刚度的改进,使接触力模型中的刚度与接触变形量之间产生耦合关系,从整体上改善接触力模型的精度.另外,为了改善表示能量耗散的阻尼因子的精度,甚至通过拟合能量指数的方式;然而,在此过程中忽略了接触刚度与能量指数之间的协调关系.由此可见,为了改善接触力模型的精度,最好的方式是不断改进表示能量耗散的阻尼因子的精度;而不建议在忽略接触刚度量纲与能量指数之间关系的情况下,对接触刚度进行改进;更不可基于Hertz 接触刚度的情况下对能量指数进行拟合.除非同时改变接触刚度,并保证接触刚度与能量指数之间的协调关系以确保接触力量纲的正确性.

1.5 关于等效连续接触力模型的讨论

当前所有的连续接触力模型(见表1 与表2)均是人为的在Hertz 接触力模型[16-17](F=Kδn)中引入表示能量耗散的阻尼项(D=χδn)[110],将准静态Hertz 接触力模型扩展为可表示能量耗散的动态连续接触力模型(F=Kδn+χδm)[16-17].然而,通过上述连续接触力模型描述碰撞行为存在以下四个方面的问题:(1)虽然可通过对阻尼项积分将碰撞过程中的能量耗散均匀化分布在压缩与恢复路径上[65],但是无法通过人为添加的阻尼因子解释能量耗散的物理原因[111];(2)在碰撞体几何形状与材料属性确定的情况下,Hertz 接触刚度可描述低速碰撞环境下的弹性碰撞(碰撞速度低于临界弹性碰撞速度ρ为碰撞体的密度)[66];然而当碰撞速度大于临界弹性碰撞速度时,碰撞体将发生弹塑性变形[112];此时的Hertz 接触刚度高估了弹塑性碰撞阶段的接触刚度,无法正确描述弹塑性碰撞行为[113];以至于目前连续接触力模型中的接触刚度在弹塑性接触过程中高估了压缩过程中的应变能与接触力的大小[61];(3)在碰撞过程中,一部分能量在局部接触碰撞过程中被损耗,另一部分能量在压缩过程中以应变能的形式储存在碰撞体中,并以应力波的形式影响未受冲击区域的变形状态[114],但是当前的连续接触力模型均未考虑碰撞引起的应力波在接触体中的传播效应[115];(4)为了分析碰撞引起的应力波传播效应,需要借助连续介质力学对碰撞问题进行求解,研究柔性碰撞体在吸收应变能之后的变形状态;当碰撞时间接近或者大于接触体的基础振动周期时,碰撞激发的应力波将从边界返回并达到碰撞位置[116],与还没有结束的碰撞过程产生耦合效应激发多重压缩-恢复等过程[117].然而,由于连续接触力模型关于弹塑性碰撞过程中能量损耗与接触力的计算存在偏差,影响了柔性体吸收能量的大小,以至于在接触碰撞持续时间与相对碰撞速度的预测上存在误差,使得多体系统中碰撞行为与柔性变形之间的耦合关系变得更加复杂.因此,正确表征碰撞过程中的能量耗散对于揭示能量耗散机理至关重要[118],也是揭示多体系统中碰撞与整体系统柔性特征之间耦合关系的前提条件.

更重要的是,由于压缩过程中系统能量等于初始碰撞动能减去碰撞体吸收的应变能

该等式清楚地描述了压缩行为本质上是初始动能不断向应变能转化的过程.然而,在能量的转化过程中,因为Hertz 接触刚度明显大于弹塑性变形阶段的接触刚度,根据式(40)可知,在弹塑性碰撞的压缩过程中碰撞初始动能转化为应变能的效率变快,以至于能量快速被消耗而改变了碰撞接触时间,同时也改变了接触速度的变化规律.另外,由于阻尼因子是接触刚度与恢复系数的函数,所以对阻尼因子积分计算碰撞过程中能量耗散时,会发现压缩过程中耗散能量[61]明显比实际情况多(由于接触刚度的原因),且耗散速度比实际的快,无法准确计算碰撞过程经历的时间.图3 中可以看出耗散能量多(迟滞环的面积越大)的接触力模型不仅相对侵入深度较小且相对碰撞速度变化急促.所以,当前动态连续接触力模型描述弹塑性碰撞行为能量耗散不准确的根本原因在于压缩过程中接触刚度与实际情况不符.更进一步,正因为Hertz 接触刚度高估了弹塑性碰撞阶段的接触刚度,改变碰撞体在压缩过程中吸收能量的效率,进而影响碰撞过程中的能量损耗,使多体系统中弹塑性碰撞行为与柔性变形之间能量传播与耗散机制更加复杂,无法准确获得弹塑性碰撞引起系统振动机理的真实原因.

2 准静态接触力模型的研究现状

在实际碰撞行为中,弹塑性碰撞行为相比纯弹性碰撞现象更加常见,Lankarani等[66]指出当初始碰撞速度大于或者等于(说明较小的碰撞速度就能引起接触体的局部弹塑性变形)时碰撞表面将会出现永久接触变形.魏悦广等[119]强调碰撞体在重载荷下必须考虑接触体的弹塑性变形.在当前高端工程装备与航空航天等领域,经常面临高速重载轻质的工况环境,因此,相比纯弹性碰撞现象,弹塑性接触碰撞现象更值得关注.

最早弹塑性变形现象的研究主要来源于人们对接触材料硬度与工程机械中结合面粗糙度的探索[120].其原因在于现实中根本不存在完全平整的接触面,物体的接触表面实际上是由众多几何尺寸不同的微凸体构成(如图7[121]为了便于研究一般近似假设微凸体的形状为圆柱形、球体、正弦波或者波浪形[20,52];其中球体之间的碰撞行为在实际工程应用中普遍存在),即从宏观角度上看似光滑平整的零件表面,但实际上其表面接触形貌呈现不同程度的粗糙和凹凸不平的特征[122].两个机械表面的真实接触情况是许多离散微凸体相互挤压与滑动的过程,导致其粗糙表面的真实接触面积远小于名义接触面积[44],造成了在极小的真实接触面积上承受较大的载荷,以至于产生弹塑性变形最终导致其局部表面被压溃[123],通常情况下塑性变形必定伴随着弹性变形[124].

图7 表面粗糙度形貌Fig.7 The surface roughness topography

一般而言,接触行为按照是否“贴合”分为协调接触(conforming contact)和非协调接触(nonconforming contact)[125].其中协调接触为两个固体表面在无变形时精确地或者相当接近地 “贴合” 在一起;而非协调接触为具有不相似外形且无变形的两个固体接触时,首先在一个点或沿一条线相碰,分别称为点接触和线接触.以上对协调与非协调接触的定义是从宏观接触表面出发,忽视了表面粗糙度对接触行为的影响[21,47].然而,当从微观接触表面出发时,任意接触表面均不可能出现“贴合”的接触行为,而是均为点接触或者线接触;因此,当考虑接触表面微观粗糙度时,其协调接触与非协调接触并没有明显的界限,均可视为非协调接触[123,126-127].

当忽略接触过程中应变率的变化与塑性流,并认为接触体各向同性,一般将整个接触过程分为完全弹性阶段、弹塑性阶段和完全塑性等三个阶段[128].其中在完全弹性接触阶段,Hertz 理论提供了一个封闭的力学本构关系;一旦接触行为达到材料的屈服极限,其接触行为将进入弹塑性阶段;作为完全弹性和完全塑性变形的过渡阶段,使得整个接触行为连续的弹塑性本构方程则相对难以获得[12];当达到弹塑性接触变形的极限时,接触行为将进入完全塑性阶段,该阶段的本构关系只有在简单的假设条件下才能获得其解析解[12,129].

最早对完全塑性问题的研究源自于人们对材料硬度的探索,1900年Brinell 研究了塑性变形区域的应力分布[52],并提出接触体的硬度为接触力除以接触区域的整体面积;随后,1908 年Meyer 认为接触材料的硬度等于接触力除以在压痕方向上的接触投影面积[130](值得注意的是,1933 年Abbott 与Firestone 在研究轴承表面特征的文章中并没有提及接触过程中塑性变形的工作[12]).1948 年Tabor[131]在Ishlinskii 的基础上分析了一个硬性球形压头被压入一个软性金属的表面,当外力卸载时,金属表面发生了塑性的永久变形.1972 年,Lee等[132]利用有限元与实验测量的方式获得接触材料硬度与压痕深度没有必然的联系.1985 年Sindair等[133]在忽略摩擦的情况下研究了刚性球体与弹塑性半空间接触体之间的压痕,分析了接触过程中的应力与应变之间的关系.1986 年Johnson[130]系统地研究了不同几何形状接触体在完全弹性和弹塑性以及完全塑性接触阶段的本构关系,获得了一系列适用于不同接触场合的有效弹塑性接触力模型.

直到1987 年,Chang等[134](CEB 模型)基于统计模型根据G-W 模型的假设条件研究了粗糙接触表面在大载荷下的弹塑性变形状态,首次提出了包含弹性和塑性变形阶段的准静态接触模型.然而,CEB 模型[134]在描述完全弹性阶段过渡到完全塑性阶段出现了不连续现象,即塑性阶段的本构关系并不能使接触过程从弹性阶段平滑地过渡到塑性阶段.为了弥补CEB 模型中的不连续特征,Zhao-Maletta-Chang(ZMC 模型)[83]利用多项式插值的方式描述弹塑性接触阶段的力与变形之间的关系;然而,该模型依然没有很好地解决CEB 模型的不连续问题.Thornton[135]在忽略弹塑性过渡阶段的情况下,基于Hertz 理论提出了可连续描述整个接触过程(包括完全弹性和塑性变形)的接触力模型.在此基础上,Johnson[130]推导了刚性球压入半空间的平均接触压力与下压量的关系,获得了简化的弹塑性接触力模型.Stronge等[136]基于Johnson 模型[137]通过修正压痕与接触半径之间的几何关系提出了包含弹性、弹塑性与完全塑性变形的接触力与接触变形量之间的函数关系式.为了进一步提高准静态接触力模型的精度,并考虑到有限元方法的不断发展,Jackson等[138]基于有限元法根据Mises 屈服理论提出的JG 模型,证实了ZMC 模型并没有很好解决CEB 模型的不连续问题.另外,Kogut 等(KE 模型)[139]采用有限元法验证了CEB 模型在塑性变形阶段低估了接触力的大小,而ZMC 模型在弹塑性阶段低估了接触力的大小;并提出了一种无量纲的连续弹塑性本构方程可描述整个弹塑性接触过程.Shankar等[140]利用有限元法扩展了KE 模型与JG 模型在弹塑性接触过中弹性阶段过渡到塑性阶段的临界值.Zhang等[141]利用有限元法在同时根据Hertz 理论与Mises 屈服理论情况下,忽略碰撞过程中的摩擦效应,根据不同接触阶段变形量与载荷之间的关系基于Newton-Raphson 方法获得弹塑性接触阶段的接触力模型.Etsion等[142]借助有限元法在研究两个小球之间弹塑性碰撞过程中卸载阶段的力与变形之间的关系,以此为基础,提出了可连续描述弹塑性接触行为的力学模型.Du等[143]分析了两个小球在发生弹塑性法向碰撞时的能量耗散,提出了可描述低速环境下发生的弹塑性碰撞行为,并通过有限元法验证了该模型中碰撞参数的合理性.Burgoyne等[144]在研究颗粒物质中孤立波的传播特征时,采用有限元分析方法导出了材料性能的经验函数,提出可用于弹塑性碰撞行为仿真的准静态弹塑性接触力模型,并通过Hopkinson 实验验证了该模型的正确性.Komvopoulos等[145]采用拟合有限元结果的方式将碰撞过程分为弹性、小变形的弹塑性、中等变形的弹塑性与大变形的弹塑性等四个阶段.在类似的工作中,Brake[146]做出了突出的贡献,将弹塑性接触行为分为弹性、弹塑性和塑性接触阶段,其中弹性阶段服从Hertz 理论,塑性阶段服从Johnson 理论或者利用Meyer 硬度指数描述完全塑性接触阶段;忽略接触过程中硬度的变化,利用多项式插值的方式近似弹塑性接触阶段的本构方程.马道林等[113]借助Brake 模型的思想,利用合理的边界条件改进了过渡阶段的函数关系,基于Johnson弹塑性接触理论发展了可连续从弹性变形过渡到塑性变形的准静态弹塑性分段接触力模型.

从准静态弹塑性接触力模型的发展历程可看出该类模型的发展同样经历了从不连续到连续的演化过程.建立该类模型的方法主要分为3 种:(1)以固体力学为基础,在简化假设条件的情况下推导精度较低的弹塑性解析模型;(2)借助有限元结果进行多项式拟合获得半解析解的弹塑性接触力模型;(3)在半解析弹塑性接触力模型的基础上,通过合理的假设与边界条件推导精确的解析解模型.表3 中给出了9 种常见的连续准静态弹塑性接触力模型,图8中描述了表3 中准静态接触力模型中接触力与变形量之间的关系.

表3 连续准静态弹塑性接触力模型Table 3 Continuous quasi-static elastoplastic contact force models

续表3

一般而言,弹塑性接触行为在加载时分为弹性、弹塑性与完全塑性接触阶段,而在卸载时只有弹性变形(如图8);整个弹塑性接触行为的加载与卸载阶段又称为压缩与恢复阶段.图8 中接触力模型加载与卸载曲线的差异性与研究对象、接触材料属性、边界条件以及曲线拟合与近似方式不一致有关.为了清楚描述整个弹塑性接触过程,图9 给出了表3 中一般性的接触力与弹塑性变形量之间的关系,其中弹性阶段在整个弹塑性接触行为中只会在初始接触阶段发生,并且只占整个接触过程的极短一部分,几乎可以忽略[113].当接触变形量超过临界弹性变形量 δy时,接触行为进入弹塑性接触阶段,其接触力与变形量之间的关系近似线性;当载荷继续增加,接触变形量进一步增大且超过临界塑性变形量 δp,接触行为进入完全塑性接触阶段,其接触力与变形量之间的关系基本为线性.另外,准静态弹塑性接触力模型认为在弹性阶段没有能量耗散(其加载与卸载曲线完全重合),且接触过程遵守Hertz 理论.当接触变形量大于临界弹性接触变形量时[113],由于在压缩过程中发生了弹塑性变形,以至于其加载曲线不再遵守Hertz 接触理论[139](其接触刚度也不再是Hertz 接触刚度),而卸载过程依然符合Hertz 理论,最终准静态弹塑性接触力模型通过不同加卸载曲线构成的封闭面积(如图9)描述了接触过程中的能量损耗,说明了接触过程中能量耗散是由于接触力改变了压缩阶段接触体的本构关系.以上的力学现象是表3 中所有准静态弹塑性接触力模型的共性特征(注意:表3 中所有模型设计的临界接触力与临界变形量不尽相同,具体的表达式请参考相关文献).

图8 准静态弹塑性接触力模型[146]Fig.8 Quasi-static elastoplatsic contact force model[146]

图9 一般弹塑性接触力模型中力与变形量之间的关系Fig.9 Force-deformation relationship for a general elastoplastic contact force model

此外,虽然表3 中的弹塑性接触力模型可精确捕捉弹塑性冲击过程中碰撞力与变形量之间的关系,但是,需要在仿真计算过程中区分碰撞过程是处于压缩阶段还是恢复阶段,如果处于压缩阶段,需要判定当前接触变形量是否超过临界弹性变形量或者临界塑性变形量[149-150];以此为依据,根据不同阶段的本构关系计算碰撞过程中的力学行为.

更重要的是,每一次弹塑性接触的计算中均需要保存永久变形量与最大接触变形量,为下次碰撞行为仿真做准备.如此繁琐的计算流程尤其不适合计算多重压缩和多重碰撞现象;也正因为复杂的计算流程导致该类模型虽然能准确描述弹塑性接触过程,但没有在多体系统动力学的商业软件中被广泛应用.

3 关于恢复系数的讨论

多体系统中碰撞行为导致的能量损耗不可避免[17],恢复系数作为衡量冲击过程中能量耗散的一个全局度量指标被广泛应用于工程与理论研究中,它可描述不同机制导致的能量耗散,包括阻尼损耗、弹塑性变形以及冲击波在接触体内的传播[151].常见的恢复系数有四类[152]:Newton 恢复系数,Poisson 恢复系数,Stronge 恢复系数与Hedrih 恢复系数.其中Newton 恢复系数根据碰撞过程中动量守恒,利用碰撞前后法向接触相对速度之比衡量碰撞过程中的能量损耗,其表达式为

Poisson 恢复系数将碰撞过程分为压缩与恢复两个阶段,在压缩阶段碰撞体之间的相对法向速度不断减小直至为零,在恢复阶段利用压缩阶段储存的冲量使相对法向速度为零的接触体分离,于是通过碰撞恢复阶段与压缩阶段的冲量之比衡量碰撞过程中的能量损耗,其表达式为

其中F为接触力;tc为压缩结束时间;te为整个碰撞过程结束时间;ts为碰撞行为起始时间.

Stronge 恢复系数认为压缩过程中系统的动能以应变能的形式储存在接触体中,在恢复阶段通过释放应变能将接触体分离,因此利用碰撞前后吸收的应变能来衡量碰撞过程中的能量损耗,其表达式为

其中 δe为碰撞后对应的相对变形量; δc为碰撞过程中的最大相对变形量; δs为接触开始的变形量.

Hedrih[153-154]在研究两个沿着圆轨迹滚动小球的正碰撞动力学行为时,通过碰撞前后小球的相对角速度定义了恢复系数,其表达式为

其中 ω(-)与 ω(+)为碰撞前后小球之间的相对角速度.

以上四种根据不同物理量定义的恢复系数均是为了衡量碰撞前后的能量耗散,在本质上是等效的;除非在斜碰撞行为中考虑了摩擦效应等因素.一般而言,当恢复系数等于1 时,被称之为纯弹性接触,即整个接触过程没有能量耗散;当恢复系数等于0时,被称之为完全塑性接触,即碰撞行为完成后初始动能被完全耗散.

对于表1 与表2 中的连续接触力模型而言,在利用该类模型计算碰撞行为时,必先确定恢复系数[3,155-156];而恢复系数的取值精度则决定了多体系统碰撞动力学的仿真精度.目前在利用该类模型计算碰撞行为时其恢复系数一般根据经验值或通过实验确定,而为了便于计算则一般采用经验值.然而,恢复系数在单纯表征能量耗散的情况下,未考虑碰撞引起的局部接触变形与碰撞时间,无法描述碰撞过程中接触力随时间的变化历程,在多体系统碰撞动力学的研究中受到限制[5];所以表1 与表2 中的连续接触力模型采用迟滞阻尼因子来描述碰撞过程中的能量耗散(阻尼环包围的面积如图3).然而,由于含阻尼因子的连续接触力模型的原型为Hertz 接触模型,导致该类模型在高恢复系数(大于0.85)下有较好的计算精度,一旦恢复系数小于0.8 该类模型则不能准确描述碰撞过程中的能量耗散[4,84],降低了多体系统碰撞动力学的仿真精度.所以当利用该类模型计算弹塑性碰撞行为时(EDEM 软件中的连续接触力模型(见表2)用于颗粒物质之间的弹塑性碰撞仿真),一般通过调节恢复系数的大小来平衡弹塑性碰撞过程中的能量耗散(大部分学者认为表1和表2 中的连续接触力模型只能用于弹性碰撞).然而,通过阻尼因子来描述碰撞过程中能量耗散的连续接触力模型,虽然可以通过调节恢复系数描述弹塑性碰撞过程中的能量耗散,但是Hertz 接触刚度高估了弹塑性接触阶段的刚度系数[157];即从整体上基于连续接触力模型依然不能准确获得弹塑性碰撞行为的力学特征(第四节中将详细说明该问题).

恢复系数不是关于接触材料的常数,而是描述碰撞过程中能量耗散的标志性参数[61,63,69,153,158];为了获得准确的恢复系数,众多学者基于准静态弹塑性本构关系在不同的碰撞环境中推导了一系列的恢复系数具体表达式.其原理在于准静态弹塑性模型可准确地获得加载阶段与卸载阶段接触力所做的功,所以可通过碰撞前后能量之比的平方根或者碰撞前后的速度之比来获得准确的恢复系数.其中包括Chang等[159]在CEB 弹塑性接触力模型的基础上利用接触表面粗糙度的统计学模型推导了恢复系数的另一种解析解模型.Thornton 模型[135]通过线性化处理完全塑性阶段的接触力与变形之间的本构关系,缓解了CEB 模型不连续的力学特性,并通过接触力在碰撞前后的做功比获得了封闭的恢复系数模型.Vu-Quoc等[147,160]在基于有限元法获得本构关系的基础上分析了碰撞过程的恢复系数,但并没有给出恢复系数的封闭解.另外,Wu等[150]从碰撞速度的观点基于有限元法定义了临界速度,将恢复系数分为两个部分描述碰撞过程中的能量损耗;Weir等[161]则推导了碰撞速度在低速环境下的恢复系数解析解;Jackson等[138,162-163]通过拟合有限元对恢复系数的分析结果优化了基于碰撞速度的恢复系数精度.Ma等[113]利用其本构方程获得了完全塑性阶段对应的恢复系数(表4 中仅给出了6 种具有封闭解的恢复系数模型,其他形式的恢复系数模型可参考文献[164]).从表4 中可以看出恢复系数不仅与碰撞体几何与材料属性有关,而且与碰撞速度密切相关.

表4 恢复系数模型Table 4 Coefficient of restitution(CoR) model

恢复系数作为衡量碰撞过程中能量损耗的直接参数,其中连续接触力模型在计算碰撞行为时需要依赖恢复系数的大小,而准静态弹塑性接触力模型可根据精确的本构关系确定恢复系数的大小.因此,对同一种接触行为,完全可以先通过准静态弹塑性接触力模型确定恢复系数的大小,并以此作为输入至连续接触力模型中,避免了根据经验值与实验测试方式获得恢复系数的弊端.在此基础上,要将面临一个关键的问题:连续接触力模型中阻尼因子做功代表的能量耗散是否与准静态弹塑性接触力模型的碰撞前后接触力做功之差描述的能量耗散具有内在联系,即阻尼环包围的面积(如图3)是否等于加卸载曲线包围的面积(如图9)？下面一节将以恢复系数为桥梁,详细解释两类模型之间的内在联系.

4 两类接触力模型之间的联系

考虑到表1 与表2 中连续接触力模型在计算弹塑性接触碰撞时,Hertz 接触刚度高估了弹塑性接触阶段的接触刚度;导致其连续接触力模型不能准确获得动态弹塑性碰撞行为的力学特征.为此,该部分将基于本课题组在2015 年提出的连续准静态弹塑性接触力模型[113](ML 模型)修正连续接触力模型在计算弹塑性碰撞时的接触刚度.

4.1 基于材料弹塑性本构的等效迟滞阻尼

由于当接触行为进入弹塑性阶段时其接触力与变形量之间的关系近似为线性,所以,借助ML 模型的本构关系将弹塑性阶段的刚度线性化为

基于该线性化的弹塑性接触刚度根据线性弹簧-阻尼模型,假设其考虑能量耗散的弹塑性接触力模型的形式为

其中χp为新的迟滞阻尼因子.

根据1.1 节中阻尼因子的推导方式,在压缩过程的结束时刻,最大应变能为

对阻尼项积分就可获得整个接触过程的能量耗散

阻尼项将接触过程耗散的能量均匀分布在压缩和恢复阶段,根据式(48),其整个接触过程的能量耗散

以及该时刻系统动量守恒式(13),将式(49)改写为

联立式(8)、式(49)与式(51),新的阻尼因子为[43]

因此,新的连续弹塑性接触力模型为

很明显,由于该阻尼因子的推导过程利用了系统的能量与动能守恒原则,以至于其阻尼项的分母中含有初始碰撞速度,即该阻尼项为迟滞阻尼;但并不妨碍计算碰撞过程中的能量耗散.

假设一个间隙球面关节,其接触参数如表5 所示,其初始碰撞速度为4 m/s,此时最大接触变形量(111.52 μm)已超过临界弹性变形量(1.578 8 μm),接触过程处于弹塑性接触阶段;此时的恢复系数可由ML 模型识别为0.6909.将该恢复系数作为新连续接触力模型的输入量可计算弹塑性碰撞行为的力学特征如图10.

图10 接触力-变形量(关于弹塑性阶段)[43]Fig.10 Relationship between the force and deformation(in elastoplastic phase)[43]

表5 接触参数Table 5 Contact parameters

通过相对误差分析,新的连续接触力模型中阻尼因子代表的能量耗散与ML 模型描述的能量耗散误差不超过0.25%[43],即阻尼环的面积与三角形的面积基本一致.该分析结论表明:当恢复系数保持一致时,两类模型在描述弹塑性碰撞行为的能量耗散时是等效的,说明了在弹塑性发生时,人为阻尼项表示的能量耗散就是弹塑变形做的功,与阻尼项的积分路径无关.由于两类模型能保证碰撞过程中能量耗散的一致性,所以两类模型在计算最大接触力和最大变形量方面虽然不能完全一致(毕竟两类模型的数学表达式和计算原则均不相同),但整体上可保持协调,尤其是碰撞后接触体的运动状态基本一致;主要归因于新的连续弹塑性接触力模型式(52)利用线弹塑性接触刚度代替了Hertz 接触刚度,避免了Hertz 接触刚度与弹塑性接触刚度不符的缺陷;更重要的是由于两类模型在描述碰撞过程中的能量耗散时完全等效.

因此,相比表1 与表2 中的连续接触力模型,新的弹塑性接触力模型可精确地计算弹塑性碰撞过程,弥补了Hertz 接触刚度引入的误差;另外,相比表3 中的准静态弹塑性接触力模型,新的连续接触力模型不仅在力学特征上能与准静态弹塑性接触力模型保持一致,并且简化了准静态弹塑性接触力模型在计算弹塑性碰撞过程中的复杂流程,避免了在计算过程中区分碰撞行为处于压缩或恢复阶段,更不必保存每一次碰撞行为导致的永久与最大变形量;而只需要判定接触行为是否发生.然而,该模型由于迟滞阻尼项中分母含有初始速度导致该模型在计算颗粒物质之间的碰撞接触行为时会导致数值奇异问题(因为大部分颗粒物质在初始阶段均保持静止状态,颗粒之间初始相对碰撞速度为零).

4.2 基于材料弹塑性本构的等效黏性阻尼

为了避免上述所提的连续接触力模型分母中初始相对碰撞速度导致的数值奇异问题[88],该小节依然基于线性化的弹塑性接触刚度,将弹塑性碰撞过程等效为线性的振动行为(如图4),该模型的动力学方程为单自由度非受迫阻尼自由振动方程式(23),根据该方程的解析解与初始碰撞条件式(24)～ 式(30)可推导系统的阻尼系数为[94]

所以该连续弹塑性接触力模型为

在该模型的推导过程中没有利用碰撞前后的能量与动能守恒原则,避免了获得的阻尼项分母中含初始相对碰撞速度的特征,即阻尼项为黏性阻尼;成功回避了连续接触力模型式(53)的数值奇异问题.

关于该模型的仿真精度,本文以一维球链为研究对象(如图11),以实验数据为基准,对比分析了EDEM 软件中连续接触力模型(见表2 中Tsuji等[96]模型)与新接触力模型式(55)之间的精确度.相关仿真参数与实验数据详见文献[117].

图11 一维球链的碰撞实验装置示意图[117]Fig.11 Experimental setup of the one-dimension granular chain[117]

图12 给出了新连续接触力模型式(55)在计算颗粒物质动态碰撞行为时的计算结果,孤立波在一维球链中的传播特性被很好地复现,其结果被实验数据证明是合理的.当与EDEM 软件中的连续接触力模型对比分析时,图13 中展现出两种模型基本能反映孤立波在一维球链中的传播特性,但是在孤立波波峰对应的最大接触力方面,两者与实验数据的误差百分比见表6;可以清楚地看出,新的连续接触力模型在计算最大接触力方面精度明显高于EDEM软件中被使用的连续接触力模型[94],再一次证明了Hertz 接触刚度在描述弹塑性接触行为的刚度属性时存在误差.

图12 一维球链中孤立波的传播特征Fig.12 Propagation of solitary wave in the one-dimension granular chain

图13 连续接触力模型之间的对比分析Fig.13 Comparative analysis between two continuous contact force models

表6 误差分析对比Table 6 Comparative analysis of the error percentage

5 接触力模型在多体系统动力学中的应用

多体系统中最早采用刚体模型基于冲量动量法研究碰撞问题,认为应力波在碰撞体中的传播速度无限大,几乎同时影响了碰撞体各质点速度;引起了多体系统广义速度的跳跃,造成在考虑切向摩擦的碰撞动力学方程求解中往往出现无解或者多解的情况[115];并在整个碰撞过程中,采用恢复系数描述碰撞前后速度的变化以及能量耗散.然而,该方法并不能获得碰撞力随时间的变化历程;要正确处理碰撞中的摩擦问题,与分析接触力随时间的变化规律,需要利用连续接触力模型考虑碰撞体之间的微运动规律,解决由于切向摩擦引起的数值不稳定问题,并通过阻尼因子阐述碰撞过程中的能量耗散.另外,连续接触力模型从数学的角度为碰撞问题提供了一种简单易行的方法,并且认为碰撞不再是瞬时过程,可以精细研究碰撞过程中碰撞力与侵入量以及碰撞速度之间的变化关系,这有利于研究由于碰撞而引起的结构破坏现象[166],这也是连续接触力模型在工程领域被广泛应用的原因.然而,连续接触力模型忽略了碰撞引起的应力波在碰撞体中的传播效应[115].实际多体系统中的碰撞问题不仅会在碰撞位置产生局部的弹塑性变形,而且碰撞引起的局部应力集中会以应力波的形式在碰撞体中传播;为了同时描述碰撞行为引起的局部弹塑性变形以及应力波的传播规律,必须借助连续介质力学研究多体系统中的碰撞接触行为.

以一维球链的碰撞现象为例[117](如图11),可以反映利用上述单点接触力模型在计算多体系统中碰撞问题时遇到的困难.图11 中的一维球链,采用弹簧的方式局部柔化刚体小球之间的接触碰撞问题,该球链的接触碰撞现象不同于两个小球之间的单点接触碰撞行为,球链中的碰撞现象不仅涉及多点碰撞,而且涉及接触过程中的多重压缩与多重碰撞现象.在第一对小球碰撞结束之前,第二对小球已经进入压缩阶段并发生多点碰撞现象,以此类推,第一对小球发生碰撞激发的碰撞力以孤立波的形式在球链中传递并引发多重碰撞现象;当最后一对小球产生碰撞并将碰撞力重新传递到开始的碰撞位置时将发生多重压缩现象,并且除第一次碰撞外,因为能量耗散导致其他碰撞激发的动态响应逐渐衰减[167].被撞击的一维球链中孤立波在球链中传播与反射引发多重碰撞现象,当波反射到碰撞开始位置时可能在单个接触点上发生多重压缩现象,与此同时碰撞过程中伴随着能量耗散以至于碰撞现象越来越微弱.

因此,当柔性体之间发生碰撞行为时,首先在碰撞的局部区域发生弹塑性变形,并产生高幅值且急剧变化的碰撞力[115],同时激发应力波在柔性体中传播,以及会诱发应力波与碰撞之间相互耦合作用等一系列复杂的动态行为[116].多体系统中的能量一部分由于局部接触区域发生的弹塑性变形被耗散,另一部分能量被传递到碰撞体中以瞬态应力波的形式在柔性体中传播与反射,以至于影响柔性体碰撞后的运动状态;另外,在碰撞结束前,应力波在柔性体中被反射到接触位置与碰撞行为发生耦合现象,进而诱发多重压缩与多重碰撞等现象,进一步影响接触力的大小.除此之外,在此期间需要对多次微碰撞过程中的碰撞力峰值、碰撞时间及碰撞次数进行准确预测,错综复杂的碰撞与柔性变形耦合关系增加了多体系统碰撞动力学的建模仿真难度[168];更进一步,上述复杂的动态交互行为令碰撞行为与柔性变形之间的耦合关系变得更加扑朔迷离.围绕上述复杂耦合问题国内外众多学者开展了系统的研究[41,169],按照碰撞过程中是否考虑弹塑性变形将以下研究分为两类:(1)基于动态连续接触力模型的多体系统碰撞动力学;(2)基于准静态弹塑性接触模型的多体系统碰撞动力学.

5.1 未考虑弹塑性接触变形的碰撞动力学研究

在忽略碰撞行为中弹塑性变形的情况下,众多学者对多体系统中碰撞现象与柔性变形的相互作用进行了系统的研究.文献[170-171]利用浮动坐标法建立曲柄滑块机构中柔性连杆的动力学模型,采用动态接触力模型计算刚性滑块之间的碰撞力,分析了刚柔耦合与碰撞行为之间的动态关系.文献[114,172]与胡海岩等采用绝对节点坐标法分析了多体系统中间隙关节碰撞现象与构件变形的特征,研究了柔性机器人在抓取过程中的接触碰撞问题.考虑到绝对节点坐标在描述变形体时引入了大量的弹性坐标,文献[173]采用Craig-Bampton 子结构模态综合法缩减了弹性坐标数量,分析了间隙关节元素之间碰撞对柔性系统的影响.张云清等[174]利用浮动坐标法分析了柔性构件小变形与碰撞行为的耦合现象,建立了曲柄滑块机构柔性连杆的动力学模型,讨论了关节刚度与碰撞现象对系统振动频率的影响.文献[175-176]采用绝对节点坐标法建立柔性杆与刚性孔之间的动力学模型,考虑了柔性杆与刚性孔之间的碰撞与摩擦效应对系统的动态影响;或者采用假设模态法建立柔性体的动力学模型,在弹性范围内采用动态接触力模型计算柔性杆与刚体之间的碰撞力,揭示了冲击波的传播效应.刘昊等[177]采用绝对节点坐标法建立了空间充气展开绳网捕获系统的动力学模型,基于Hertz 接触理论分析了被捕获物体与柔性绳网之间的碰撞现象.方建士等[178]利用子系统法考虑柔性梁的动力刚化效应,基于考虑能量耗散的Hertz 接触力模型研究了大范围运动柔性梁与刚性地面之间的接触碰撞现象.彭海军等[179]针对柔性多体系统中接触体之间的碰撞与摩擦问题,基于辛离散化提出了一种非光滑接触方法,避免了因互补变量过多而导致求解效率较低的困难.虞磊等[180]基于绝对节点坐标法建立了柔性梁与柔性板的动力学模型,将柔性体之间的碰撞检测转化成基本几何体的碰撞检测,并利用Hertz 接触理论计算了柔性体与刚体以及柔性体之间的碰撞力,仿真结果与RecurDyn 对比验证了该方法的正确性.

上述工作主要集中研究碰撞行为对多体系统中变形构件的动态影响[181],在弹性范围内分析了碰撞行为与柔性变形之间的耦合关系,对柔性体与碰撞行为之间的精细化建模与求解方面做出了贡献.但是较少关注多体系统中能量在经历碰撞与柔性变形之后的重新分配与耗散问题,实际上多体系统中碰撞与柔性变形的耦合关系本质就是能量相互转换与耗散的过程[182-183];同时,现有工作较少关注碰撞与应力波之间引发的多重压缩与多重碰撞问题.

5.2 考虑弹塑性接触变形的碰撞动力学研究

众多国内外学者发现多体系统在高速碰撞下,如果忽略碰撞行为中的弹塑性变形将导致其系统动态响应与实际情况存在偏差[184-186].Yigit[187]研究了柔性梁与刚性地面的碰撞动力学,利用完全弹性接触力模型计算了柔性梁与刚体之间碰撞力,仿真结果显示完全弹性接触模型并不能获得精确的结果,通过实验证明弹塑性接触力模型更加符合实际情况.蓦朋波等[116]利用动态子结构模态综合法推导了柔性构件之间碰撞激发瞬态波传播的动态子结构模型,在保持接触刚度不变的情况下,在碰撞过程中考虑了永久的弹塑性变形,基于单轴压缩(UC)模型计算了接触体之间的碰撞力,分析了弹塑性波的传播特性.文献[111,188]围绕柔性体与刚体碰撞的动力学做了系统的研究,主要采用浮动坐标法建立柔性构件的动力学模型,分别利用动态连续接触力模型与单轴压缩(UC)模型计算柔性体与刚体之间的接触力,考虑了正碰撞与斜碰撞过程中的摩擦现象,分析了接触在发生弹性变形与塑性变形时系统的动态特性,指出了动态连续接触力模型与静态弹塑性模型不仅对能量耗散的描述不同,而且通过两种模型计算的接触力大小也大相径庭.Chen等[112]采用模态综合法建立柔性体的动力学模型,提出接触区域判定方法,通过罚函数的惩罚因子计算柔性体之间的弹塑性接触力,对比分析了弹性和弹塑性碰撞后系统的动能相差大约30%,说明了考虑碰撞中弹塑性变形的必要性与能量耗散的主要来源.王检耀等[168]与洪嘉振基于有限元法同时采用罚函数法和附加约束法的接触模型计算了柔性杆之间的接触碰撞问题,发现柔性梁与碰撞之间的相互作用会引发多次微碰撞问题.

上述研究考虑了碰撞行为中发生的弹塑性变形,证实了弹塑性变形在碰撞行为中较为常见且不可忽视,主要分析了弹塑性碰撞行为引起的应力波在柔性体中的传播规律,研究了碰撞与柔性变形之间耦合作用诱发的多重碰撞问题.然而,在碰撞力的计算过程中对于弹塑性变形引起的接触刚度变化缺乏精细化的建模与分析(实际的接触刚度分为弹性、弹塑性与塑性接触刚度)[189],导致其碰撞中能量耗散的预测出现偏差,影响了碰撞的接触时间与接触力的大小,干扰了柔性体中应力波与碰撞行为之间的耦合作用.另外,利用准静态弹塑接触力模型或者惩罚因子计算碰撞体之间的接触力,其中惩罚因子难以表征碰撞体之间的真实接触性质;其次,由于准静态弹塑性接触模型对接触行为的描述完全取决于相对接触变形量的大小;因此,基于该类模型计算接触行为容易受到积分误差的干扰.相反,连续接触力模型从相对接触变形量与相对变形速度上同时控制能量耗散与接触力的大小[190-191],对积分误差不敏感,这也是连续动态弹塑性接触力模型的另一种优势.

多体系统碰撞动力学中碰撞与柔性体耦合效应容易引起多重压缩与多重碰撞现象,其中柔性体从碰撞行为中吸收的能量与耗散能量直接相关[192];而应力波的传播规律取决于柔性体吸收的能量[193],其碰撞时间也依赖于能量耗散的效率.所以,正确表征碰撞过程中的能量耗散与揭示能量耗散的本质原因是多体系统碰撞动力学的核心内容,只有保证精确计算碰撞过程中能量耗散在时间和空间上的传播规律,才能从本质上揭示碰撞与柔性变形之间的耦合关系.

6 总结与发展趋势

(1) 本文主要阐述了不考虑碰撞过程中应变率和硬度发生变化的正向接触力模型的研究进展,分析了两类不同接触力模型的发展过程,指出了连续接触力模型在弹塑性碰撞动力学计算中存在的问题,讨论了准静态弹塑性接触力模型在计算弹塑性碰撞行为时复杂的计算流程.明确了Hertz 接触刚度是导致连续接触力模型在计算弹塑性碰撞行为的主要误差来源,以此为依据,以线性的弹塑性接触刚度为基础,通过碰撞过程中的能量守恒原则推导了新的迟滞阻尼因子,证实了人为阻尼因子代表的能量损耗与弹塑性碰撞中弹塑性变形做的功等效,而与阻尼项的积分路径无关.

(2) 当前的连续接触力模型与准静态弹塑性接触力模型均是由法向碰撞与接触行为发展而来,忽略了碰撞行为中的切向摩擦效应,以至于在分析斜碰撞行为时,需要先根据接触力模型计算法向接触力,再根据摩擦模型计算接触体之间的切向接触力,整个过程忽略了法向接触力与摩擦行为之间的耦合关系;因此,未来的接触力模型的发展要从斜碰撞入手,发展考虑法向和切向耦合作用的连续接触力模型.

(3)载人返回舱着水/着陆以及航空发动机轴承的失效问题均涉及碰撞体与流体之间的流固耦合关系,为满足国家重大航天发展需求,结合Reynolds 方程将接触体之间的流体引入至连续接触力模型的阻尼因子中;因此,充分考虑流体动力学对碰撞行为的影响,建立可描述流体与法向碰撞力相耦合的连续接触力模型是未来的发展趋势.

(4) 考虑到高速重载与高精度是未来多体系统的发展趋势,间隙关节导致的接触碰撞与磨损现象将更加突出,然而当前对于间隙关节碰撞行为的研究均忽略了高速重载下导致接触体局部的弹塑性变形,给间隙关节元素之间的碰撞力计算引入了较大的误差;另一方面,在基于Archard 磨损模型预测间隙关节元素之间磨损深度时,也将引起接触应力的计算误差,最终影响间隙关节磨损表面的几何属性重构.因此,考虑间隙关节局部弹塑性变形的多体系动力学性能分析与寿命预测是未来的发展趋势.

(5) 综合考虑多体系统碰撞动力学中的多种非理想因素,包括多体系统中构件的柔性变形、间隙关节的碰撞行为、关节的润滑与摩擦效应,以及碰撞行为引起的局部弹塑性变形;建立保真度更高的多体系统碰撞动力学模型,分析局部弹塑性变形与系统柔性变形之间的耦合关系,探索应力波在柔性体中的传播特征;揭示多体系统关节中流固耦合效应与系统大范围非线性运动之间的内在联系,开发具有通用性多体系统碰撞动力学软件将是未来面临的挑战之一.