李健¹, 高微¹, 张亚双², 刘欲诺³
1. 东北大学 理学院, 辽宁 沈阳 110819;
2. 中国航发沈阳黎明航空发动机有限责任公司, 辽宁 沈阳 110043;
3. 东北大学 机械工程与自动化学院, 辽宁 沈阳 110819
收稿日期：2017-04-20
基金项目：国家自然科学基金资助项目(11672072, 11502050)。
作者简介：李健(1972-), 男, 辽宁沈阳人, 东北大学教授。

摘要：提出一种颗粒弹塑性接触的通用建模方法来研究颗粒系统的振动及耗能特性.构造颗粒法向塑性接触本构方程基本形式, 并采用有限元法(FEM)获得无量纲本构关系; 提出颗粒弹塑性细观接触加-卸载多路径模型, 给出了颗粒塑性接触能量耗散的计算公式.采用离散单元法(DEM)对简谐激励下颗粒系统的振动和耗能特性进行了数值仿真分析.结果表明, 颗粒使系统产生非线性振动, 颗粒介质在系统进入共振区的前后均呈现出类似混沌的复杂非线性运动状态.颗粒间的法向塑性接触没有改变颗粒阻尼效应的主要耗能方式, 但对颗粒动态接触行为及颗粒系统动力学特性产生重要影响.
关键词：颗粒介质塑性接触有限元法离散单元法无量纲本构方程多路径模型
A Modeling Approach of Elastoplastic Contact for Analyzing the Vibration and Energy Dissipation Characteristics of Particles
LI Jian¹, GAO Wei¹, ZHANG Ya-shuang², LIU Yu-nuo³
1. School of Sciences, Northeastern University, Shenyang 110819, China;
2. Aero Engine Group of China Shenyang Liming Aero-Engine Corporation Ltd., Shenyang 110043, China;
3. School of Mechanical Engineering & Automation, Northeastern University, Shenyang 110819, China
Corresponding author: LI Jian, E-mail: jianli@mail.neu.edu.cn
Abstract: A general contact modeling method of elastoplastic contact between particles is proposed to study the vibration and energy dissipation characteristics of the particle system. The basic form of constitutive equation concerning normal plasticity contact between particles is established and a dimensionless constitutive relation can be obtained based on the finite element method accordingly. A loading-unloading multipath model of elastoplastic contact between particles is presented. A theoretical formula of normal plastic energy dissipation is deduced. Based on the discrete element method, numerical simulations on the characteristics of vibration and energy dissipation for the particle system under harmonic excitation are carried out. The results show that the linear system exhibits the characteristics of nonlinear vibration caused by particles. Granular media appears to be a complicated nonlinear state of motion similar to the chaos near the resonance region. Although the way of energy dissipation on particle damping effect remains unchanged, normal plastic contact alters the dynamic contact behavior between particles and has great influence on the dynamic characteristics of the vibrating particle system.
Key words: granular mediaplastic contactFEMDEMdimensionless constitutive equationmultipath model
振动系统中的颗粒介质会呈现复杂运动形式并消耗大量系统能量^[1], 研究人员将其应用于系统减振并称之为“颗粒阻尼”^[2-3].颗粒耗散振动系统能量的同时也会改变系统的动力学响应, 前期研究中已发现竖直振动颗粒介质对边界的冲击力出现分岔和混沌现象^[4]，并在理论上进行了解释^[5].同时颗粒介质本身的耗能特性、运动行为^[6]以及系统因颗粒运动产生的非线性动力响应之间也应该存在密切关联；目前, 这方面的研究鲜有报道，有待进一步展开, 基于离散单元法(DEM)的颗粒物质体系细观力学行为研究有助于理解这种关联并进一步揭示颗粒介质的振动耗能机理及其对振动系统动力学响应的影响规律.由于颗粒阻尼效应的一个重要实现方式为稀疏快速流颗粒间的弹塑性碰撞, 故需引入颗粒法向塑性接触行为来完善颗粒物质细观接触模型.
本文目的是给出一种适用于任意材料的颗粒动态黏-弹-塑性细观接触建模方法, 并将该模型应用于DEM来研究振动颗粒系统的动力学行为和能量耗散.颗粒细观接触中的黏性和弹性模型已较为成熟, 本文主要工作为法向塑性接触的引入.主要思路为:考虑到有限元(FEM)技术已能方便、高效地研究任意材料的单颗粒复杂弹塑性接触行为, 基于FEM结果, 首先构造颗粒法向塑性接触无量纲本构方程；由于考虑塑性后颗粒接触行为与载荷路径相关, 故进一步提出颗粒弹塑性细观接触加-卸载多路径模型, 并推导塑性接触耗能的理论公式.将建立的颗粒黏-弹-塑性细观接触模型引入DEM并开发数值计算程序, 采用完全非弹性蹦球分析法、相平面法和谱分析方法研究简谐激励下颗粒群运动形式、耗能特性以及振动结构的非线性动力学响应.
1 颗粒塑性接触本构方程颗粒塑性接触本构关系近年开始受到部分学者关注.由于不同材料的塑性行为非常复杂, 为简化计算, 这些研究通常引入理想弹塑性等假设并局限于特定材料^[7].FEM作为一种成熟高效的数值算法可用来研究任意材料颗粒间的弹塑性细观接触行为.本文基于ABAQUS软件构建颗粒间法向塑性接触本构关系, 采用轴对称半球模型并在局部加密.颗粒选用铝材料(弹性模量E=70 GPa, 泊松比μ=0.33)并最初只考虑完全弹性情况, 两颗粒半径均为R=5 mm.计算得到的法向接触力和接触半宽的关系与赫兹理论解吻合良好, 说明了颗粒接触FEM模型和计算方法的有效性.
以颗粒完全弹性接触FEM模型为基础, 在5 mm和10 mm两球接触模型中引入铝材料真实塑性参数, 获得两铝球颗粒发生弹塑性接触时, 其接触力-接触半宽本构关系曲线如图 1所示, 图中还同时给出了两球完全弹性接触的对应曲线.由图中可知, 考虑铝材料塑性参数后, 铝球颗粒在逐渐增大的法向接触力作用下很快进入塑性, 其接触区抵抗变形的能力也随之有较大程度的下降.
图 1(Fig. 1)

图 1 5 mm和10 mm铝球颗粒接触力-接触半宽本构关系曲线Fig.1 Force-half width curves for 5 mm and 10 mm aluminum balls

5种工况对应的颗粒塑性接触本构关系曲线如图 2所示.其中, 塑性接触力通过在法向总接触力中去除由Mises准则确定的弹性力临界值获得.图 2中5条曲线形状非常接近, 反映了不同工况下的塑性接触都具有相同的塑性材料参数.同时, 曲线倾斜程度与等效半径变化趋势保持一致, 说明等效半径是颗粒塑性接触本构关系中的重要变量, 这一点与颗粒弹性接触经典理论相同^[8].FEM结果还表明, 塑性接触带来颗粒半径的不可逆改变, 这与文献[7]的研究相吻合, 这种改变将对等效半径产生影响.考虑到塑性接触带来的粒径改变很小并发生在范围很小的局部，且颗粒频繁碰撞, 故本文假设颗粒间每次碰撞的接触点不重合, 颗粒间重叠量的计算不考虑颗粒变形.
图 2(Fig. 2)

图 2 颗粒塑性接触本构关系曲线Fig.2 Constitutive relation curves of plastic contact between particles 图中5-10(3.5)表示等效半径3.5 mm的5 mm和10 mm两个铝球接触, 5-板(5)表示等效半径5 mm的5 mm铝球与平板接触.其余工况类似.

为了获得更具一般性的颗粒塑性接触本构关系, 将ABAQUS得到的接触半宽与法向总接触力均去除弹性部分后, 引入弹塑性临界值(即弹性接触阶段的最大值)对数据进行无量纲化处理, 根据图 2的分析结果并参考赫兹接触理论, 构造颗粒塑性接触无量纲本构方程基本形式:

(1)

式中:为无量纲法向塑性接触力; 为无量纲塑性接触半宽; 为无量纲等效半径(引入典型等效半径R₀=2.5 mm进行无量纲化).函数G(, m_i)表示颗粒塑性接触力与等效半径和材料塑性力学性能无量纲参数m_i有关, 这其中也包含了材料塑性力学性能对颗粒等效半径的影响.的指数项待定常数H由材料塑性力学性能决定.具体应用时, 函数G(, m_i)可展开成如式(2)所示的多项式形式并根据FEM结果进行参数识别.函数G(, m_i)选取的多项式阶次通常与材料塑性参数的复杂程度成正比.

(2)

将图 2中不同工况下FEM结果代入式(1), 取G(, m_i)为二阶多项式, 基于最小二乘法进行参数识别, 获得m_i(i=0, 1, 2)和H后可得铝材料颗粒塑性接触无量纲本构方程:

(3)

图 3以两种工况为例, 给出了式(3)的理论计算值与ABAQUS数值仿真结果的比较, 二者吻合良好, 说明G(, m_i)函数的二阶多项式形式已经可以满足铝颗粒的塑性材料参数要求.
图 3(Fig. 3)

图 3 不同工况下理论计算与ABAQUS数值结果比较Fig.3 Comparison of results obtained from the theoretical formula and ABAQUS for different conditions (a)—5 mm和10 mm铝球；(b)—10 mm铝球和平面铝板.

将式(1)与有限元数值计算相结合, 可以获得任意材料颗粒的塑性接触本构关系, 这种方法为离散元研究塑性接触下颗粒介质的动力学行为提供了有效途径.
2 弹塑性接触加-卸载模型及塑性耗能引入塑性后颗粒的接触行为与载荷路径相关, 任意时刻的法向接触本构关系由加-卸载历史和接触状态决定.本文在给出颗粒塑性接触本构关系的基础上, 提出了图 4所示的颗粒法向弹塑性细观接触加-卸载多路径模型.纵坐标f_N表示法向弹塑性细观接触力；横坐标δ表示两颗粒接触法向重叠量, δ_y，δ_max和δ_R分别代表弹性区临界值、本卸载路径的最大变形及残余变形.颗粒弹塑性细观接触加-卸载路径分为4段:1(OA)→2(AB)→3(BC)→4(BD).其中路径1和3为加-卸载段, 2和4仅为加载段.
图 4(Fig. 4)

图 4 颗粒法向弹塑性细观接触加-卸载多路径模型Fig.4 Loading-unloading multipath model of elastoplastic contact between particles in the normal direction

路径1为弹性接触阶段, 其法向接触力f_e采用赫兹经典理论计算^[8]:

(4)

式中:E和R分别为两接触颗粒的等效弹性模量和等效半径, 下标1和2为颗粒号, υ为泊松比.该路径对应的弹性接触力最大值f_y可通过由Mises准则确定的式(5)得到^[9]：

(5)

式中σ_y为材料屈服强度.由式(4)、式(5)可得到弹性区内颗粒法向接触重叠量临界值δ_y:

(6)

同样, 由几何关系a²=Rδ可得δ_y对应的弹性临界接触半宽a_y:

(7)

当δ>δ_y时, 颗粒接触进入路径2和4的塑性阶段, 此时法向总接触力中除弹性力外还包括塑性部分.为便于计算和表达, 根据式(8)恢复铝球颗粒法向塑性接触无量纲本构方程(3)的量纲后, 可以得到颗粒塑性接触力f_p的表达式:

(8)

(9)

其中a为总接触半宽.此时颗粒间的法向弹塑性细观接触力f_N为

(10)

路径3代表颗粒进入塑性接触后的卸载及再加载, 本文假定该过程中颗粒接触区域出现硬化, 其接触本构关系符合赫兹经典弹性解.
颗粒进入塑性接触后, 因其塑性变形不可恢复, 该过程中弹塑性细观接触力将耗散振动系统能量, 能耗U可通过式(11)的精确积分计算(即图 4中OABC所围成的面积):

(11)

3 颗粒黏-弹-塑性细观接触模型颗粒接触模型的精细程度决定数值模拟方法能否对颗粒真实接触状态和耗能特性准确描述^[10].本文在颗粒黏-弹性模型^[3]基础上引入弹塑性细观接触加-卸载多路径模型，并由此提出颗粒黏-弹-塑性细观接触模型如图 5所示.
图 5(Fig. 5)

图 5 球形颗粒黏-弹-塑性细观接触模型Fig.5 Viscosity-elastic-plastic contact models of spherical particle at the meso-level (a)—法向接触；(b)—切向接触；(c)—滚动接触.

图 5中连接器l用来判断是否发生接触, 且不引入任何附加力.图 5a为法向接触细观模型, 其中弹性接触力和黏滞力分别采用经典赫兹理论^[8]和黏性阻尼模型计算, k_n和c_n代表法向弹性接触刚度和黏滞阻尼系数.根据图 5, 当法向弹性接触力大于临界值f_y时, 颗粒法向接触进入弹塑性阶段，此时, 塑性接触选择器s_p闭合, 法向接触合力中包含塑性接触刚度k_p的贡献.图 5b表示基于摩尔-库伦摩擦定律及明德林增量理论^[11]的颗粒切向接触模型, μ_t和s_t分别代表滑动摩擦系数和选择器, k_t和c_t表示切向接触刚度和黏滞阻尼系数^[3].图 5c为颗粒滚动方向的细观接触模型^[3], μ_r和s_r表示滚动摩阻系数和选择器, 滚动力偶矩采用扭簧-旋转黏滞阻尼器力学模型计算^[12], k_r和c_r分别表示滚动方向的抗扭刚度和黏滞阻尼系数.颗粒与边界的接触力计算可采用上述相同方法, 任意时刻某颗粒位置和运动参数(速度和位置矢量等)可由牛顿第二定律及其数值积分获得.由上可知, 图 5中的参数可由理论公式计算得到^{[3, 8, 10-12]}，并依赖于粒径、密度、弹性模量、泊松比、恢复系数、摩擦系数等颗粒材质及表面属性参数.这些基本参数可能对颗粒介质的动力学行为产生重要影响并有待展开进一步研究.
由颗粒黏-弹-塑性细观接触模型可知, 弹性接触力不产生能耗, 法向塑性力和黏滞力、切向和滚动方向的黏滞力、摩擦力总计6种接触力损耗系统能量.
4 数值算例本文基于提出的颗粒黏-弹-塑性细观接触模型, 在Fortran平台开发了DEM数值仿真程序, 以简谐激励下含颗粒介质振动系统为例(如图 6所示), 从宏观和细观不同尺度研究了颗粒及振动结构的动力学行为及耗能特性.图 6中箱体长42 mm, 宽42 mm, 高24 mm, 与颗粒材质相同；质量0.3 kg并做刚体运动.只考虑系统在z方向的振动, 弹簧刚度系数k为9 720 N/m, 不考虑颗粒影响的系统固有频率为28.7 Hz.
图 6(Fig. 6)

图 6 含颗粒介质振动系统示意图Fig.6 Schematic of the vibration system with particles

颗粒为铝材料, 弹性模量70 GPa, 泊松比0.33, 密度2 700 kg/m³.颗粒恢复系数0.5, 摩擦系数0.5, 填充率0.2, 颗粒数量总计2 000个, 平均粒径D为2 mm, 颗粒直径在0.99D~1.01D之间呈正态分布.颗粒细观接触模型中各参数(详见图 5)取值见表 1.颗粒总质量0.222 kg, 若将颗粒质量全部视为系统的附加质量, 则系统固有频率为21.7 Hz.沿z向施加幅值9.72 N的简谐激励, 在10~50 Hz范围的简谐激励下逐点扫频.考虑到颗粒介质运动的复杂性, 引入颗粒群体质心和箱体质心的运动参数来研究系统的振动特性, 即将此时的颗粒介质视为一个完全非弹性蹦球来研究其运动.
表 1(Table 1)

表 1 颗粒细观接触模型中的参数值Table 1 Parameter values in mesoscale contact model kn/(N·m-3/2) 1.66×1010 cn/(N·s·m-1) -5.27×102 ct/(N·s·m-1) -5.27×102 kt/(N·m-3/2) 1.76×1010 kr/(N·m1/2) 7.04×105 cr/(N·m·s) 8.66×10－6 μt 0.3 μr/m 4.44×10－4 sp , st, sr 选择器, 不引入附加力

表 1 颗粒细观接触模型中的参数值 Table 1 Parameter values in mesoscale contact model

图 7、图 8分别为激振频率为22.3 Hz时, 箱体和颗粒质心的频域和相平面图.图中参数均经无量纲化处理.其中, z, f, v分别表示位移、频率和速度, A₀，f_i分别表示系统静变形和激振力频率(以下同).由图可知, 箱体除激励频率下的主振动外开始出现其他非线性频率, 而颗粒质心运动呈现了较明显的三倍周期现象, 表明此频率激振下, 颗粒运动开始进入混沌状态.
图 7(Fig. 7)

图 7 激振频率为22.3 Hz时, 箱体与颗粒群质心的频域曲线Fig.7 Frequency-amplitude curves of the box and particles for excitation frequency 22.3 Hz 图中数字为峰值的横坐标值 (a)—箱体；(b)—颗粒群质心.

图 8(Fig. 8)

图 8 激振频率22.3 Hz时, 箱体与颗粒群质心的相平面图Fig.8 Phase planes of the box and particles for excitation frequency 22.3 Hz

图 9为10~50 Hz范围内, 几个典型频率简谐激励下箱体与颗粒群质心的相平面图.由图可知, 激振频率为25.5 Hz时, 箱体与颗粒均表现为激励频率下的单频线性振动, 且幅值基本保持一致, 但存在恒定相位差.颗粒质心相图近似矩形, 表明颗粒在一个运动周期内与箱体存在4次明显碰撞, 但这些碰撞对箱体的运动形式影响不大.颗粒系统与箱体位移响应幅值在激励频率为30.2 Hz时达到极值, 该数值与系统固有频率28.7 Hz较接近, 颗粒对系统的附加质量特征几乎消失, 此时系统呈现共振状态, 箱体与颗粒均为单频主振动, 颗粒运动幅值超过箱体并在相图上更接近于矩形.
图 9(Fig. 9)

图 9 不同激励频率下箱体与颗粒群质心的相平面图比较Fig.9 Phase planes of the box and particles for different excitation frequencies (a)—12.7 Hz；(b)—22.3 Hz； (c)—25.5 Hz；(d)—30.2 Hz； (e)—35.0 Hz；(f)—47.7 Hz.

另外, 箱体和颗粒在低频激励下振动时均表现出了较明显的非线性, 随着频率增大, 二者振动形式不再保持一致, 颗粒相对箱体开始出现更大的运动幅值和较低的运动速度.箱体和颗粒在高频激励下振动时主要表现为单频线性主振动, 颗粒振动的位移幅值和速度都明显小于箱体, 颗粒群体在与箱体的高频碰撞中聚集并悬浮于空中做微幅振动.值得注意的是, 系统在进入共振区的前后, 颗粒总要经历一段类似混沌的运动状态(如22.3 Hz及35 Hz), 这与不考虑颗粒间塑性接触的情况有很大不同^[3].
图 10为10~50 Hz激振频率范围内, 各分项耗能占总耗能的比例.由图可知, 相比不考虑塑性接触^[3], 法向黏滞耗能在整个激励频带下仍为颗粒阻尼的最主要耗能方式, 只是在低频激励下比例略低.这些现象说明, 虽然考虑颗粒间的塑性接触后, 塑性耗能占比很小, 没有改变颗粒阻尼的主要耗能方式, 但颗粒间的法向塑性接触影响了颗粒的动态接触行为, 进而对颗粒群体和结构的动力学特性产生影响.
图 10(Fig. 10)

图 10 不同激振频率下各分项耗能占总耗能之比Fig.10 Ratio of subitem energy dissipation to total energy dissipation for different excitation frequencies

5 结论1) 材料复杂的弹塑性接触力学特性及其与加载路径的相关性，使颗粒塑性接触理论模型的建立非常困难.基于有限元数值模拟建立颗粒弹塑性细观接触模型的方法能为几乎任意材料的颗粒系统离散元多尺度仿真分析提供一条有效途径.
2) 颗粒介质是受激振动系统产生非线性的来源.在共振区, 颗粒运动幅值超过系统并与系统保持一致的线性振动特性; 而在共振区附近, 颗粒呈现出类似混沌的高度非线性运动状态.
3) 虽然考虑法向塑性接触时, 法向黏滞耗能仍然为系统颗粒阻尼效应的最主要耗能方式, 但颗粒间的塑性接触对颗粒动态接触行为以及颗粒介质和结构的动力学特性均产生重要影响.
参考文献

[1]	Shah B M, Pillet D, Bai X, et al. Construction and characterization of a particle-based thrust damping system[J].Journal of Sound and Vibration, 2009, 326(3/4/5): 489–502.
[2]	Panossian H V. Structure damping enhancement via non-obstructive particle damping technique[J].Journal of Vibration and Acoustics, 1992, 114(1): 101–105.DOI:10.1115/1.2930221
[3]	李健, 刘璐, 严颖, 等. 基于离散单元法的颗粒阻尼耗能减振特性研究[J].计算力学学报, 2013, 30(5): 664–670. ( Li Jian, Liu Lu, Yan Ying, et al. Study on energy dissipation and vibration reduction characteristics particle damping based on DEM[J].Chinese Journal of Computational Mechanics, 2013, 30(5): 664–670.)
[4]	姜泽辉, 郑瑞华, 赵海发, 等. 完全非弹性蹦球的动力学行为[J].物理学报, 2007, 56(7): 3727–3732. ( Jiang Ze-hui, Zheng Rui-hua, Zhao Hai-fa, et al. Dynamical behavior of a completely inelastic ball bouncing on a vibrating plate[J].Acta Physica Sinica, 2007, 56(7): 3727–3732.DOI:10.3321/j.issn:1000-3290.2007.07.017)
[5]	姜泽辉, 李斌, 赵海发, 等. 竖直振动颗粒物厚层中冲击力分岔现象[J].物理学报, 2005, 54(3): 1273–1278. ( Jiang Ze-hui, Li Bin, Zhao Hai-fa, et al. Phenomena of impact bifurcations in vertically vibrated granular beds[J].Acta Physica Sinica, 2005, 54(3): 1273–1278.DOI:10.3321/j.issn:1000-3290.2005.03.047)
[6]	周宏伟, 陈前, 林莎. 垂直简谐激励下阻尼颗粒动态特性研究[J].振动与冲击, 2007, 26(9): 124–127. ( Zhou Hong-wei, Chen Qian, Lin Sha. An experimental investigation on NOPD coupled system in vertical direction[J].Journal of Vibration and Shock, 2007, 26(9): 124–127.DOI:10.3969/j.issn.1000-3835.2007.09.030)
[7]	Vu-Quoc L, Zhang X, Lesburg L. Normal and tangential force-displacement relations for frictional elasto-plastic contact of spheres[J].International Journal of Solids and Structures, 2001, 38(36/37): 6455–6489.
[8]	Hertz H. The contact of elastic solids[J].Journal für die Reineund Angewandte Mathematik, 1882, 92(4): 156–171.
[9]	Johnson K L. Contact mechanics[M]. Cambridge: Cambridge University Press, 1985: 96-119.
[10]	Saeki M. Impact damping with granular materials in a horizontally vibrating system[J].Journal of Sound and Vibration, 2002, 251(1): 153–161.DOI:10.1006/jsvi.2001.3985
[11]	Mindlin R D, Deresiewicz H. Elastic spheres in contact under varying oblique forces[J].Journal of Applied Mechanics, 1953, 20(3): 327–344.
[12]	Ji S Y, Daniel M H, Hayley H S. Comparisons of physical experiment and discrete element simulations of sheared granular materials in an annular shear cell[J].Mechanics of Materials, 2009, 41(6): 764–776.DOI:10.1016/j.mechmat.2009.01.029