上海理工大学学报  2023, Vol. 45 Issue (6): 574-583   PDF    
颗粒滚动阻力本构参数标定试验研究
董朋昆, 高政国     
北京航空航天大学 交通科学与工程学院,北京 100191
摘要: 针对颗粒滚动阻力本构参数难以通过阻力测量试验来直接标定,开展了两种离散元滚动阻力模型MDEM(modified discrete element method)和HDEM(hysteresis discrete element method)参数的试验标定方法研究。利用自由滚动颗粒的摆动过程建立了颗粒间滚动刚度系数Kr与滚动阻尼系数Cr的试验标定公式,搭建了激光位移传感器测量颗粒微振动的光学试验平台,通过测量试件摆动时程曲线实现了MDEM模型中KrCr的标定;提出基于材料单轴循环拉伸试验的标定方法,利用颗粒本体材料加卸载循环作用下的滞回曲线实现了HDEM模型中的弹性滞后系数 β的标定。试验结果表明,该试验方法能够有效地标定颗粒间滚动阻力本构参数。
关键词: 滚动阻力参数     MDEM模型     HDEM模型     弹性滞后    
Experimental study for identifying constitutive parameters of rolling resistance between particles
DONG Pengkun, GAO Zhengguo     
School of Transportation Science and Engineering, Beihang University, Beijing 100191, China
Abstract: Aiming at the difficulty of directly identifying the constitutive parameters of rolling resistance between particles using physical experiment of resistance measurement, an experimental method for identifying the parameters of MDEM (modified discrete element method) model and HDEM (hysteresis discrete element method) model was proposed. Based on the process of swing back and forth of free rolling particle, two formulas for identifying rolling stiffness coefficient Kr and rolling damping coefficient Cr in MDEM model were derived, and an experimental device for measuring time-history displacement curve of swing using a laser displacement sensor was set up. Another identifying method based on material uniaxial cyclic tensile experiment was utilized to measure the hysteresis curve of particle material on cycling load and unload, and the elastic hysteresis coefficient β was identified through fitting the experimental curves. The experimental results show that these two methods can well identify the rolling resistance constitutive parameters of particles.
Key words: rolling resistance parameter     MDEM model     HDEM model     elastic hysteresis    

滚动阻力模型参数是离散元数值仿真中描述颗粒间滚动阻力矩的重要构成部分。在现有研究中,滚动阻力模型参数的取值主要依赖于经验,直接通过试验方法对颗粒间滚动阻力参数进行标定的研究成果较少。造成这一现象的主要原因是颗粒发生滚动时受到的荷载情况复杂,比如颗粒间发生相对转动时,在接触面上不仅存在滚动摩擦力,还存在滑动摩擦力等作用,并且滚动摩擦系数与滑动摩擦系数在数值上相差数个量级。除此之外,在竖向压力与颗粒形状影响下,滚动摩擦力与滑动摩擦力之间还会发生相互转换[1]。这些因素造成在滚动阻力参数标定试验中很难构造出一个仅存在滚动摩擦力的工况,使得通过滚动摩擦力测量试验识别滚动模型参数的难度较大。针对目前已有的滚动阻力参数标定试验研究,按照试验方法,大致可分为如下两类:

a. 直接标定法。基于能量守恒原理,通过试验方法直接进行标定。如张荣芳等[2]基于高速摄像技术和能量守恒方法,依据滚动摩擦力造成的能量损失在总能量中的占比,实现了滚动摩擦系数的标定。崔涛等[3]利用高速摄像机拍摄玉米种子沿导轨的滚动过程,根据种子颗粒运动过程中的能量守恒,计算得到种子与导轨材料之间的滚动摩擦系数。林榕[4]提出将铁球从倾斜铁轨滚下,与底部粘性材料发生完全非弹性碰撞。通过测量粘性材料的形变量得到铁球动能,再根据能量守恒换算得到铁球在滚动下滑过程中的滚动摩擦系数。李贝等[5]设计了小球从斜面上滚落的试验,通过记录落点与斜面终点间的水平距离,依据能量守恒原理得到滚动摩擦力所做的功,从而得到滚动摩擦系数。刘万锋等[6]也设计出了类似的试验装置,他们将倾斜的石板面换成了一个弧形导轨,小球运动到导轨末端后做平抛运动,根据水平位移和水平抛出点高度值计算小球和导轨间的滚动摩擦系数。这类试验方法简便易行,但滚动摩擦力在做功时通常伴随有滑动摩擦力的作用,使滚动摩擦系数的测量值与实际值间存在一定的偏差。

b. 间接标定法。这类方法是从颗粒体系的宏观行为出发,如通过颗粒堆积形成的休止角,反向推算颗粒的滚动摩擦系数。比如韩燕龙等[7]通过移动圆柱管的颗粒堆积试验,建立了滚动摩擦系数与椭球颗粒物料堆积角间的线性关系。王云霞等[8]使用有机玻璃板与铝质圆筒进行了玉米种子堆积仿真试验,再通过回归分析堆积角仿真数据,获得种子间的静摩擦系数与滚动摩擦系数。吴文江等[9]通过试验和仿真相结合的方法,对沙粒间的恢复系数、静摩擦系数和滚动摩擦系数进行了标定。从颗粒的宏观现象反算细观参数方法具有试验方法操作简单、易于实施的优点,但缺点是干扰影响因素较多,识别精度不稳定。

从上述研究可以看出,若要实现通过试验方法对滚动阻力模型参数的标定,需要解决两个技术难题:a. 设计寻找适用于颗粒滚动阻力参数测量的颗粒运动行为,该运动过程中颗粒主要受滚动摩擦力作用,受其他外力作用的影响较小;b. 基于该运动行为研究,设计相应的试验方法和测量装置,以满足可识别的技术条件,并应具有足够的识别精度。

目前,颗粒运动离散元方法中常用的滚动阻力模型为MDEM(modified discrete element method)模型[10],该模型包含有两个滚动阻力参数:滚动刚度系数Kr和滚动阻尼系数Cr。HDEM(hysteresis discrete element method)模型是笔者[11]基于MDEM模型并考虑材料弹性滞后[12]建立的,主要由滚动滞弹簧、滚动粘壶、非承拉节点和滚动滑阻器等元件构成。与MDEM模型相比,引入弹性滞后的HDEM模型能更为合理地表征颗粒在临近静止时的动能耗散现象,得到与试验曲线更吻合的模拟结果。针对颗粒间的滚动阻力模型参数:MDEM模型中的KrCr,HDEM模型中的弹性滞后系数β,本文分别以颗粒的自由滚动和材料单向循环拉伸为基础,提出了相应的试验标定方法。

实现KrCr试验标定的理想方法是寻找颗粒的一种特殊运动状态:纯滚动状态,即在颗粒间仅存在滚动摩擦力,而滑动摩擦力为零[10]。然而,实际中颗粒间发生相对转动时通常伴随有碰撞、分离等复杂的接触行为,在试验中难以直接构建颗粒的纯滚动状态。因此,基于纯滚动状态直接设计试验,通过测量颗粒间滚动阻力实现模型参数的试验标定难度较大。但是,作者发现在刚性平面上自由滚动的弹性颗粒达到完全静止前存在微小的往复摆动现象[13],在该过程中颗粒仅受到滚动阻力矩和惯性作用,其他外力影响可以忽略。并且颗粒的自由滚动行为在试验中易于实现,能为KrCr的试验标定提供试验基础。

本文首先通过理论推导证明了自由滚动的颗粒在摆动过程中具有与纯滚动状态类似的振动行为,利用颗粒在自由滚动中的摆动频率值建立了KrCr的识别公式。然后基于颗粒在自由滚动中的往复摆动行为,使用具有高测量精度的激光位移传感器设计了微振动试验检测平台,对圆柱形试件在刚性平面上自由滚动的时程曲线进行了测量,并依据试验曲线频率值识别了试件的KrCr。此外,β的标定需要依赖数值模拟,通过对片状料材进行单轴循环拉伸测量滞回曲线,然后使用弹性滞后本构曲线对试验数据进行拟合来实现对β的标定。

1 MDEM模型参数试验标定 1.1 参数标定公式

根据与滚动速度的相关性,现有的离散元滚动阻力模型大致可分为两类[14]。一类是滚动阻力矩与滚动速度无关,阻力矩大小正比于法向接触力,力矩方向与滚动方向相反。另一类是与滚动速度相关,表现为一种粘滞力特征。MDEM模型是一种典型的与速度相关型模型,它主要是由法向接触模型、切向接触模型以及滚动阻力接触模型3个部分构成,如图1所示。其中:下标n代表法向;s代表切向;r代表滚动。法向接触模型与切向接触模型均包含有线性弹簧、阻尼器及非承拉节点等元件,切向接触模型中还包含有表征滑动阻力临界值的滑阻器。滚动阻力模型与切向模型的基本组成元件类似,差别在于滚动阻力模型中各元件表达的是阻力矩值M。MDEM模型的本构关系为


图 1 离散元接触模型 Fig. 1 Contact model in DEM
$ {M}_{{\rm{r}}}=\mathrm{min}({K}_{{\rm{r}}}{\theta }_{{\rm{r}}}+{C}_{{\rm{r}}}{\dot{\theta }}_{{\rm{r}}}\text{,}{\mu }_{{\rm{r}}}{F}_{{\rm{n}}}) $ (1)

式中:Mr为滚动阻力矩;Fn为法向接触力;${\mu _{\rm{r}}}$为滚动摩擦系数;${\theta _{\rm{r}}}$${\dot \theta _{\rm{r}}}$分别为颗粒的滚动角位移和角速度。

1.1.1 自由滚动状态

图2给出了弹性颗粒在刚性平面上的自由滚动过程。圆柱形颗粒受初始扰动向前滚动,接触面提供的阻力矩使颗粒滚动速度逐渐减小,在静止前颗粒会在滚动的终点位置处发生左右往复摆动。这里为便于表述,将颗粒在平面上的自由滚动划分为减速滚动过程和摆动过程。在摆动过程中,颗粒与刚性平面之间的相对水平位移值为零,滑动摩擦力值为零。颗粒在该过程中仅受滚动阻力矩和惯性作用,受其他因素影响很小,如图3所示。颗粒在自重与惯性力矩作用下,在与刚性平面形成的接触面上会分布着不均匀的竖向应力,这些应力形成了滚动阻力矩值${M_{\rm{r}}}$。利用颗粒摆动角位移值能够建立摆动过程中的力矩平衡方程。假设在某一时刻颗粒圆心$O'$的摆动水平位移值为$ x $,则绕滚动终点O的摆动角位移值为${\theta _{\rm{r}}} = {x \mathord{\left/ {\vphantom {x r}} \right. } r}$。力矩以逆时针方向为正,颗粒绕点O的力矩平衡方程


图 2 自由滚动颗粒静止前的摆动阶段 Fig. 2 Swing state of free rolling particle in process

图 3 摆动过程颗粒受力示意图 Fig. 3 Diagram of forces of particle in the swing state
$ - {I_{O'}}{\ddot \theta _{\rm{r}}} - {F_x} \cdot r + {M_{\rm{r}}} = 0 $ (2)

式中:${F_x} = m\ddot x$为颗粒摆动时的惯性力,m为颗粒质量,$\ddot x$为摆动水平加速度;${\ddot \theta _{\rm{r}}}$为颗粒的滚动角加速度;r为颗粒半径。滚动阻力矩通过弹性恢复力与阻尼力可以表示为${M_{\rm{r}}} = - {C_{\rm{r}}}{\dot \theta _{\rm{r}}} - {K_{\rm{r}}}{\theta _{\rm{r}}}$,二者均由接触面提供。颗粒绕滚动终点位置O的转动惯量${I_O} = {I_{O'}} + m{r^2}$${I_{O'}}$表示颗粒绕圆心的转动惯量,其表达式为

$ {I_{O'}} = \frac{1}{2}m{r^2} $ (3)

将滚动阻力矩和转动惯量的表达式代入到式(2),式(2)可简化为

$ {I_O}{\ddot \theta _{\rm{r}}} + {C_{\rm{r}}}{\dot \theta _{\rm{r}}} + {K_{\rm{r}}}{\theta _{\rm{r}}} = 0 $ (4)

假设颗粒的摆动时程曲线为低阻尼系统振动曲线,即

$ y = {{\rm{e}}^{ - \xi \omega t}}a\sin \left( {{\omega _{\rm{r}}}t + \alpha } \right) $ (5)

式中:$\omega $${\omega _{\rm{r}}}$分别表示有阻尼圆频率和无阻尼圆频率;$\xi $表示阻尼比;$\alpha $表示初始相位角;$a$表示初始振幅。通过式(5)对式(4)进行求解,得到弹性颗粒在自由滚动状态下的滚动刚度系数Kr和滚动阻尼系数Cr的表达式为

$ {K_{\rm{r}}} = {I_O}{\omega ^2} $ (6)
$ {C_{\rm{r}}} = 2{I_O}\omega \xi $ (7)
1.1.2 纯滚动状态

图4给出了相互接触的两个半径均为r的颗粒1与颗粒2发生相对转动,产生的转动位移值可以分解为纯滑动位移与纯滚动位移,Bardet等[15]、Iwashita等[10]和Jiang等[16]已经提出了分解公式。假设在时间步长$dt$内,颗粒1与颗粒2的转动角位移值分别为$d{\theta '_1}$$d{\theta '_2}$。若$d{\theta '_1} \ne d{\theta '_2}$,则可以将转动角位移分解为正对称和反对称两个部分。其中正对称部分为相对滑动角位移$d{\theta '_{\rm{s}}}$,反对称部分为相对滚动角位移$d{\theta '_{\rm{r}}}$,即


图 4 滚动位移分解示意图 Fig. 4 Diagram of decomposition of rotational displacement
$ \left\{ \begin{gathered} d{\theta '_1} = d{\theta '_{\rm{s}}} - d{\theta '_{\rm{r}}} \\ d{\theta '_2} = d{\theta '_{\rm{s}}} + d{\theta '_{\rm{r}}} \\ \end{gathered} \right. $ (8)

$d{\theta '_{\rm{r}}}$表示的运动部分即为纯滚动部分。当无外荷载作用时,相互挤压接触的两个弹性颗粒将发生自由振动,此时接触面上不发生相对滑移,是一种纯滚动状态,如图5所示。


图 5 纯滚动过程颗粒受力示意图 Fig. 5 Diagram of forces of particle in the pure rolling state

假设某一时刻颗粒的滚动角位移值为${\theta '_{\rm{r}}}$,惯性力矩为$- {I_{O'}}{\ddot \theta '_{\rm{r}}}$,滚动阻力矩为${M'_{\rm{r}}}$,同样是由接触面的竖向应力提供。绕颗粒圆心$O'$的运动平衡方程为$- {I_{O'}}{\ddot \theta '_{\rm{r}}} + {M'_{\rm{r}}} = 0$,由于${M'_{\rm{r}}}$是由滚动阻尼力矩$- {C_{\rm{r}}}{\dot \theta '_{\rm{r}}}$和弹性恢复力矩$- {K_{\rm{r}}}{\theta '_{\rm{r}}}$组成,则运动平衡方程可表示为

$ {I_{O'}}{\ddot \theta '_{\rm{r}}} + {C_{\rm{r}}}{\dot \theta '_{\rm{r}}} + {K_{\rm{r}}}{\theta '_{\rm{r}}} = 0 $ (9)

利用式(5),并使用$\omega '$表示颗粒在纯滚动状态下的频率值,$\xi '$表示在纯滚动状态下的阻尼比值,则纯滚动状态下的KrCr可分别表示为

$ {K_{\rm{r}}} = {I_{O'}}{\omega '^2} $ (10)
$ {C_{\rm{r}}} = 2{I_{O'}}\omega '\xi ' $ (11)

由于滚动刚度与滚动阻尼是颗粒材料的固有物理属性,在自由滚动状态下与纯滚动状态下应保持一致,即

$ {K_{\rm{r}}} = {I_{O'}}{\omega '^2} = {I_O}{\omega ^2} $ (12)
$ {C_{\rm{r}}} = 2{I_{O'}}\omega '\xi ' = 2{I_O}\omega \xi $ (13)
1.1.3 两种推导过程对比

表1对比了通过弹性颗粒自由滚动中的摆动过程与纯滚动状态建立滚动阻力参数识别公式的理论推导过程。可以看出,颗粒在摆动过程中与纯滚动状态下的振动行为类似,通过纯滚动状态与摆动过程均可以从动力学角度建立KrCr的理论识别公式,区别仅在于摆动过程与纯滚动状态下颗粒的转动中心不同。在摆动过程中颗粒是绕接触面中心转动,而在纯滚动状态下是绕颗粒质心。但是,相较于纯滚动状态,弹性颗粒在刚性平面上的自由滚动状态容易获取,更利于进行试验设计。


表 1 自由滚动状态与纯滚动状态的滚动参数识别公式对比 Table 1 Comparison of the identifying formulations between free rolling and pure rolling
1.2 试验装置

为准确识别MDEM滚动模型参数,以圆柱形试件在刚性平面上的自由滚动过程为试验基础,使用图6中的微振动试验平台对试件的滚动时程曲线进行测量。该试验装置主要由试验平台、加载气囊、激光位移传感器、电源,控制器与数据采集终端等构成,其中激光位移传感器是整个试验装置的核心。图6展示了激光位移传感器发射的激光穿过亚克力板照射到圆柱试件上的轨迹。在工作状态下,激光位移传感器的激光照射到试件上会发生反射,当试件位置发生改变时,反射的激光进入传感器接收端口的角度将发生改变,传感器内的感光片将即时检测到试件位置变化信息,并对试件的位移信息进行记录。该装置能够测量的试件形状为圆柱形或球形,能够测量的试件材质为聚氨酯、亚克力、木材、铝材、铜等。各装置的详细信息与具体功能如表2所示。


图 6 微振动试验装置 Fig. 6 Micro-vibration test device

表 2 微振动试验装置元件信息 Table 2 Component information of micro-vibration test device

微振动试验平台的操作流程具体如下:

a. 首先记录待测试件的几何尺寸和质量,用于计算转动惯量${I_O}$。对试件外表面划分不同的测量区域,目的是通过对试件不同位置进行多次测量来减小颗粒表面初始缺陷的影响。

b. 将划分区域后的试件放置在调平后的大理石平台上。然后调整固定支架位置,使激光位移传感器的激光能够照射到试件中心位置。

c. 按压试验装置左侧的加载气囊,通过产生瞬时的冲击气流对试件进行加载,使试件以一定的初速度向远离加载端的方向移动。

d. 操作数据采集终端,记录试件的滚动位移时程曲线。

使用该装置进行试验时需要注意:试验加载前需首先打开位移传感器的数据采集功能,以确保试验测量精度与时程曲线的完整性。将试件放置在激光位移传感器量程范围内,并保证试件滚动距离不会过大,避免超过激光位移传感器的有效量程。

可以观察到的试验现象为:试件在加载后开始向前滚动,在滑动摩擦力作用下逐渐减速,然后在滚动终点位置处左右摆动。

1.3 试验结果

现以试验测量得到的一条摆动时程曲线为例,介绍KrCr的识别过程。

使用图6中的微振动平台对内外径尺寸分别为12 mm和77 mm、材质为聚氨酯的圆柱试件进行测量,得到图7所示的摆动时程曲线。从图中可以看出圆柱试件的摆动时程曲线具有以下特征:a. 曲线振动幅值随时间逐渐减小。这是由于材料内摩擦、环境阻尼等作用,试件动能被逐渐耗散,摆动角位移值逐渐减小;b. 曲线上存在一些不光滑的点。这是因为试件在生产时表面存在一些初始缺陷,位移传感器发射的激光在试件表面缺陷处会发生反射偏差,而具有较高测量精度的传感器能够识别出这些偏差,反映在试验曲线上即为不光滑点。从图中可以看出试件的摆动时程曲线与低阻尼系统振动曲线在形式上类似,因此可采用式(5)对试验数据进行拟合。


图 7 圆柱形试件摆动时程曲线 Fig. 7 Time history curve of cylinderical specimen

图7中的位移曲线平移至位移为零处,然后除以试件半径,得到摆动角位移曲线,如图8所示。再使用式(5)对角位移曲线进行拟合,具体过程为:首先在曲线上选取时间为${t_{A}} = 4.266\;{\text{s}}$${t_{B}} = 5.097\;{\text{s}}$和角位移${\theta _{A}} = 0.022\;6\;{\text{rad}}$${\theta _{B}} = 0.015\;2\;{\text{rad}}$AB两点,得到


图 8 试验曲线和拟合曲线 Fig. 8 Test curve and fitted curve
$ {\omega _{\rm{r}}} = \frac{{2\pi }}{{{t_B} - {t_A}}} = 7.56\;{\text{rad/s}} $ (14)
$ \xi = \frac{1}{{2\pi }}\ln \frac{{{\theta _A}}}{{{\theta _B}}} = 0.063 $ (15)

在已知${\omega _{\rm{r}}}$$\xi $的情况下,可以得到有阻尼频率$\omega $为7.58 rad/s。然后根据AB两点的时间值与角位移值,可以求得参数$a$$\alpha $分别为0.17 4 rad、0.716 rad。因此,拟合曲线的表达式为

$ {\theta _{\rm{r}}} = 0.174{{\rm{e}}^{ - 0.479t}}\sin (7.565t + 0.716) $ (16)

使用式(16)对试验曲线进行拟合,结果如图8所示,可以看出式(16)能够与试验曲线较好吻合,表明通过低阻尼系统振动曲线对颗粒的摆动时程曲线进行拟合是可行的。最后,根据圆柱试件的质量和半径,得到试件的转动惯量值${I_O} = 4.337 \times {10^{ - 4}}\;{\text{kg}} \cdot {{\text{m}}^2}$。由式(12)与式(13)得到KrCr

$ {K_{\rm{r}}} = {\omega ^2}{I_O} = 2.482 \times {10^{ - 2}}({\rm{N}}\cdot{\rm{m}}/{\rm{rad}})$ (17)
$ {C_{\rm{r}}} = 2{I_O}\omega \xi = 4.155 \times {10^{ - 4}}({\rm{N}}\cdot{\rm{m}}\cdot{\rm{s}}/{\rm{rad}}) $ (18)
2 HDEM弹性滞后系数试验标定

通过对颗粒本体材料进行单轴循环拉伸试验,测量材料在加载和卸载过程中的滞回曲线,并使用弹性滞后本构公式[11]对试验曲线进行拟合,实现对颗粒弹性滞后系数$\beta $的标定。

2.1 HDEM滚动阻力模型

使用MDEM模型模拟颗粒自由滚动时,得到的模拟时程曲线与试验曲线之间在摆动幅值衰减上存在差异,如图9所示,从图中可以看出在6~8 s时间范围内MDEM模型无法快速耗散颗粒能量。形成这一现象的主要原因是在颗粒摆动过程中,MDEM模型的能量耗散主要依赖于与速度相关的滚动阻尼器。在滚动速度趋近于零时,阻尼器耗散动能的能力会减弱;而实际中颗粒在摆动临近静止时,由于与速度无关的材料弹性滞后产生的耗能占比是逐渐增大的[11]。因此,与速度相关的MDEM滚动阻力模型无法表征由弹性滞后形成的速度无关型的滚动阻力矩。


图 9 试验曲线和数值模拟曲线 Fig. 9 Test curve and simulation curves

图10为弹性滞后曲线示意图。材料在弹性范围内加载与卸载时,由于应变滞后于应力,加载曲线与卸载曲线不能完全重合,而是形成了一个封闭回线,称为弹性滞后曲线[12]。颗粒在滚动时由弹性滞后形成滚动阻力矩的机理可以这样描述:假设在刚性平面上一个弹性颗粒在主动力矩Ma作用下发生滚动,当滚动角速度不为零时,由于弹性滞后的作用,接触面前缘的加载应力将大于后缘上的卸载应力,如图11所示。不均匀分布的竖向应力在颗粒滚动时会产生滚动阻力矩值。这个滚动阻力矩与颗粒滚动速度无关,与接触面上的竖向应力水平有关。


图 10 弹性滞后曲线 Fig. 10 Elastic hysteresis curve

图 11 接触面应力分布示意图 Fig. 11 Diagramof stress distribution on contact surface

为表征弹性滞后在颗粒摆动临近静止时的耗能特性,笔者在文献[11]中以幂函数的形式定义了弹性滞后本构关系,提出了基于弹性滞后的速度无关型HDEM滚动阻力模型,如图12所示,图9中的数值模拟结果表明该模型能够较好地模拟颗粒在临近静止前的动能耗散过程。但是表征能量耗散的弹性滞后系数仍依赖数值拟合进行确定。因此,为实现对不同颗粒材料弹性滞后系数的直接标定,本文设计了材料单轴循环拉伸试验,并使用弹性滞后公式拟合试验数据标定弹性滞后系数。


图 12 HDEM模型 Fig. 12 HDEM model
2.2 弹性滞后本构关系

颗粒材料弹性滞后的本构关系可以表示为:

加载曲线

$ {\sigma _{{\text{load}}}}(\varepsilon ,{\varepsilon _{\rm{{rev}}}}) = \left({\sigma _{\rm{e}}} - {\sigma _{\rm{{rev}}}}\right){\left(\frac{{\varepsilon - {\varepsilon _{\rm{{rev}}}}}}{{{\varepsilon _{\rm{e}}} - {\varepsilon _{\rm{{rev}}}}}}\right)^\beta } + {\sigma _{\rm{{rev}}}} $ (19)

卸载曲线

$ {\sigma _{{\text{unload}}}}(\varepsilon ,{\varepsilon _{\rm{{rev}}}}) = {\sigma _{\rm{{rev}}}}\left[ {1 - {{\left(1 - \frac{\varepsilon }{{{\varepsilon _{\rm{{rev}}}}}}\right)}^\beta }} \right] $ (20)

式中:(${\varepsilon _{\text{e}}}$, ${\sigma _{\text{e}}}$)表示滞回曲线上的弹性极限点;(${\varepsilon _{{\text{rev}}}}$, ${\sigma _{{\text{rev}}}}$)表示曲线上的转折点;$\;\beta \in $(0, 1)为材料的弹性滞后系数,$\;\beta $值越接近于0,滞回曲线越饱满,弹性滞后耗散的能量越大。

2.3 试验装置

颗粒材料弹性滞后系数的具体试验标定步骤为:首先使用拉伸试验机对与颗粒材料相同材质的片状纺锤形试件进行单轴循环拉伸试验,得到材料的应力–应变弹性滞回曲线。然后利用式(19)与式(20)对试验曲线进行拟合,进而实现对$\;\beta $的标定。

图13给出了材料的单轴循环拉伸试验装置。该装置主要由加载试验机主体、力传感器、夹具和数据采集终端构成。使用的力传感器型号为DYLY-102,量程为5 kg,精度为0.03%。传感器两端分别与试验机主体和夹具相连。加载时,力传感器测量的荷载值将与加载试验机内引伸计测量的试件位移值形成荷载–位移曲线。试验机型号为Byes-50N,数据采集速率为30次/s,使用温度为10~35 ℃。拉伸试验使用的试件形状为片状纺锤形[17],试件截面尺寸如图14所示,厚度为2 mm,拉伸断面尺寸为24 mm2。试件材质可为聚氨酯、橡胶和硅胶等,试验温度为常温25 ℃。


图 13 单轴循环拉伸试验装置 Fig. 13 Test device of uniaxial cyclic tensile

图 14 单轴循环拉伸试验试件 Fig. 14 Test specimen of uniaxial cyclic tensile

图15为单轴循环拉伸试验的加载制度。由于橡胶材料在进行单轴循环拉伸时,填充颗粒表面附近分子链发生断裂形成Mullins效应[18]。因此,为尽可能消除残余应变对试验结果的影响,对试件进行了多次循环加载,并保持最大循环荷载值Fmax不变。


图 15 循环加载制度 Fig. 15 Cyclic loading protocols

对材料进行单轴循环拉伸试验的步骤为:

a. 首先对试件施加一个预拉力,目的是消除夹具与试验加载机、夹具与试件之间连接空隙。

b. 使用位移控制对试件进行加载,加载速率为0.001 mm/s,将试件加载至最大荷载Fmax。然后在最大荷载处维持1 s,再卸载至零。重复加载和卸载至最终得到的曲线近似重合。

2.4 试验结果

图16图17为分别对橡胶材料和聚氨酯材料进行单轴循环拉伸得到的试验曲线。可以看出,加载曲线与卸载曲线形成了明显的滞回曲线,同时材料的残余应变也随着循环加载次数增加逐渐减小,最终形成了稳定的滞回曲线,这一试验现象证实了材料弹性滞后特性的存在。在图16中,试验曲线上弹性极限点坐标为(1.74, 2.0),使用式(19)和式(20)对试验曲线进行拟合,可以得到橡胶试件的$\beta $为0.75。在图17中,试验曲线上弹性极限点坐标值为(3.42, 30.0),通过拟合可以得到聚氨酯试件的$\beta $为0.94。


图 16 橡胶材料弹性滞后试验曲线 Fig. 16 Elastic hysteresis curve of rubber sample

图 17 聚氨酯材料弹性滞后试验曲线 Fig. 17 Elastic hysteresis curve of polyurethane sample
3 结 论

基于弹性颗粒的自由滚动状态,借助激光位移传感器设计了微振动试验平台,测量了圆柱形试件在自由滚动过程中的时程曲线,然后通过低阻尼系统振动曲线拟合试验数据,得到了摆动频率值,进而实现了MDEM模型参数的标定。使用材料单轴循环拉伸机建立了HDEM模型中的弹性滞后系数试验标定方法,对橡胶材料和聚氨酯材料的弹性滞后系数进行了标定。通过上述试验研究得到以下结论:

a. 弹性颗粒的自由滚动状态具有与纯滚动状态类似的振动行为,通过颗粒的自由滚动中的摆动行为能够建立MDEM模型中参数KrCr的理论标定方法,并且颗粒的自由滚动在试验中易于实现。公式标定结果与离散元模拟结果表明,使用激光位移传感器设计的微振动试验平台具有足够的测量精度,基于颗粒自由滚动状态能够实现滚动模型参数的准确识别。

b. 通过对与颗粒相同材质的片状纺锤形试件进行单轴循环拉伸试验,利用HDEM模型本构公式对试验曲线进行拟合,能够实现对颗粒材料弹性滞后系数的标定。

参考文献
[1]
孙珊珊, 苏勇, 季顺迎. 颗粒滚动–滑动转换机制及摩擦系数的试验研究[J]. 岩土力学, 2009, 30(S1): 110-115. DOI:10.3969/j.issn.1000-7598.2009.z1.023
[2]
张荣芳, 周纪磊, 刘虎, 等. 玉米颗粒粘结模型离散元仿真参数标定方法研究[J]. 农业机械学报, 2022, 53(S1): 69-77. DOI:10.6041/j.issn.1000-1298.2022.S1.008
[3]
崔涛, 刘佳, 杨丽, 等. 基于高速摄像的玉米种子滚动摩擦特性试验与仿真[J]. 农业工程学报, 2013, 29(15): 34-41. DOI:10.3969/j.issn.1002-6819.2013.15.005
[4]
林榕. 滚动摩擦系数的测量[J]. 辽宁教育行政学院学报, 1996(5): 78-79. DOI:10.13972/j.cnki.cn21-1500/g4.1996.05.030
[5]
李贝, 陈羽, 孙平, 等. 滚动摩擦系数工程测量方法与验证[J]. 工程机械, 2017, 48(4): 29-32. DOI:10.3969/j.issn.1000-1212.2017.04.006
[6]
刘万锋, 徐武彬, 李冰. 滚动摩擦系数的测定及EDEM仿真分析[J]. 机械设计与制造, 2018(9): 132-135. DOI:10.3969/j.issn.1001-3997.2018.09.036
[7]
韩燕龙, 贾富国, 唐玉荣, 等. 颗粒滚动摩擦系数对堆积特性的影响[J]. 物理学报, 2014, 63(17): 174501. DOI:10.7498/aps.63.174501
[8]
王云霞, 梁志杰, 张东兴, 等. 基于离散元的玉米种子颗粒模型种间接触参数标定[J]. 农业工程学报, 2016, 32(22): 36-42. DOI:10.11975/j.issn.1002-6819.2016.22.005
[9]
吴文江, 郭斌, 高占凤, 等. 基于DEM的沙土颗粒建模及参数标定研究[J]. 中国农机化学报, 2019, 40(8): 182-187. DOI:10.13733/j.jcam.issn.2095-5553.2019.08.32
[10]
IWASHITA K, ODA M. Rolling resistance at contacts in simulation of shear band development by DEM[J]. Journal of Engineering Mechanics, 1998, 124(3): 285-292. DOI:10.1061/(ASCE)0733-9399(1998)124:3(285)
[11]
高政国, 董朋昆, 张雅俊, 等. 一种滞弹簧耗能的新型离散元滚动阻力模型研究[J]. 力学学报, 2021, 53(9): 2384-2394. DOI:10.6052/0459-1879-21-236
[12]
GREENWOOD J A, MINSHALL H, TABOR D. Hysteresis losses in rolling and sliding friction[J]. Proceedings of the Royal Society A:Mathematical, Physical and Engineering Sciences, 1961, 259(1299): 480-507.
[13]
WANG D R, GAO Z G, ZHANG Y J, et al. Experimental method for identifying constitutive parameters of rolling resistance between particles in discrete element method[J]. Granular Matter, 2021, 23(4): 91. DOI:10.1007/s10035-021-01113-6
[14]
ZHU H P, YU A B. A theoretical analysis of the force models in discrete element method[J]. Powder Technology, 2006, 161(2): 122-129. DOI:10.1016/j.powtec.2005.09.006
[15]
BARDET J P, HUANG Q. Numerical modeling of micropolar effects in idealized granular materials[J]. Mechanics of Granular Materials and Powder Systems, 1992, 37: 85-92.
[16]
JIANG M J, YU H S, HARRIS D. A novel discrete model for granular material incorporating rolling resistance[J]. Computers and Geotechnics, 2005, 32(5): 340-357. DOI:10.1016/j.compgeo.2005.05.001
[17]
LI Z, WANG Y L, LI X, et al. Experimental investigation and constitutive modeling of uncured carbon black filled rubber at different strain rates[J]. Polymer Testing, 2019, 75: 117-126. DOI:10.1016/j.polymertesting.2019.02.005
[18]
REBOUAH M, MACHADO G, CHAGNON G, et al. Anisotropic Mullins stress softening of a deformed silicone holey plate[J]. Mechanics Research Communications, 2013, 49: 36-43. DOI:10.1016/j.mechrescom.2013.02.002