A NOVEL DISCRETE ELEMENT ROLLING RESISTANCE MODEL BASED ON HYSTERESIS SPRING ENERGY DISSIPATION
-
摘要: 颗粒间滚动阻力对颗粒体系的稳定性起着重要作用. 在传统的离散元法中, 滚动阻力模型通常由转动弹簧、转动黏壶和摩擦元件表达, 颗粒滚动动能由黏滞力(矩)和摩擦力做功耗散. 由于黏滞力(矩)与滚动速度相关, 临近静止状态的颗粒滚动速度变小, 动能耗散减弱, 传统的离散元模拟得到颗粒由滚动到静止耗费的时间比试验观测的结果要长. 为解决这一问题, 基于摩擦学理论分析了滚动阻力产生的材料滞弹性机理, 将其引入离散元滚动阻力模型, 提出了一种速度无关型动能耗散的滞弹簧, 给出了滞弹簧的弹性恢复力计算公式, 建立了一种新型的离散元滞弹性滚动阻力模型(HDEM). 为验证新型滚动阻力模型的正确性, 通过一个光学物理试验对单个圆形颗粒试件的自由滚动过程进行了测量, 将测量数据与新型的滞弹型离散元模型和传统离散元模型计算结果进行了对比. 结果显示, 基于滞弹性滚动阻力模型HDEM计算结果与试验数据吻合程度更高, 而且模拟得到的颗粒摆动频率更符合试验现象.Abstract: The rolling resistance between particles plays an important role in the stability of the particulate systems. In a conventional discrete element method, the rolling resistance model between particles is usually made of springs, dashpots, and sliders in the rotational direction. The particles rolling kinetic energy is dissipated by the viscous (moment) and friction forces. With this model, the viscous force (moment) is directly related to the rolling velocity. Consequently, the dynamic dissipation capacity of particles close to the static state becomes weaker with the rolling velocity decreasing. It is known that the time required to simulate a particle rolling with a velocity close to zero by using the traditional discrete element method is longer than the experimental results. To solve this problem, the mechanism of rolling resistance caused by material hysteresis is analyzed based on tribological principle, and a new discrete element model of hysteresis rolling resistance (HDEM) is established. A hysteresis spring with velocity-independent kinetic energy dissipation is proposed, and its constitutive law’s formula is derived. To verify the new rolling resistance model, the free-rolling of a single round particle specimen on a flat surface is measured through a physical experiment. The measured data are compared with the results simulated by the new rolling resistance model HDEM and the conventional rolling resistance model. The results show that the results based on HDEM are more consistent with the experimental data, and the particle oscillation frequency is in better agreement with the experimental phenomenon observed.
-
Keywords:
- rolling resistance /
- particle system /
- DEM model /
- elastic hysteresis /
- hysteresis spring
-
引 言
滚动阻力问题涉及车辆工程[1], 土木工程[2]和农业[3-4]等许多领域. 与滑动阻力相比滚动阻力通常较小[5], 如硬质圆形或球型颗粒相对滚动时滚动摩擦系数一般为10−5 ~ 10−3, 滑动摩擦系数通常为0.1 ~ 1. 从数值来看, 滚动摩擦系数比滑动摩擦系数要小得多, 但是滚动摩擦对颗粒的宏观力学特性有着非常重要的影响. Bardet等[6]在离散元模型中最早考虑了滚动约束的作用, 他们发现不考虑滚动约束作用, 离散元模拟结果得到的力学参数在理论值范围之外. 他们认为接触力偶矩可能起着重要作用, 他们指出当约束颗粒滚动自由度后颗粒体系的内摩擦角增大. Morgan等[7]直接将转动阻尼引入离散元模型进行断层泥的模拟, 他们得到了与实验室评估接近的结果. Sakaguchi等[8]将“滚动摩擦”的概念引入离散元模型, 进行谷仓清空过程颗粒流阻塞试验与数值模拟研究. 他们在离散元计算程序中引入了一个滚动摩擦力矩, 发现常规离散元数值模拟中圆盘颗粒成拱不稳定, 易于破坏, 考虑滚动摩擦后能有效的模拟出物理试验中得到的阻塞成拱现象.
目前在滚动阻力理论及工程应用方面许多学者开展了相关研究[9-16]. 然而, 颗粒体系稳定过程中滚动阻力作用机制仍不十分清楚, 与材料[17-18]及滚动阻力有关的许多力学模型和方法仍有待研究.
滚动阻力通常被学者们称为“滚动摩擦”. 根据摩擦学理论, 滚动阻力的主要来源是接触表面上的微小滑移、塑性变形、材料的黏滞性、表面附着力和形状效应[19].
1875年, Reynolds[20]发现, 当金属圆柱体在橡胶表面上滚动时, 在竖向压力作用下接触面的切向位移会发生微小差异, 称为微滑移. 类似地, Tabor等[21]研究了弹性范围内, 硬质球体和硬质圆柱体在橡胶软基上滚动时的摩擦作用. 他们发现接触面上的切向力始终小于在滚动过程中存在润滑剂的情况下产生的界面滑移值. 因此, 他们认为微滑移并非弹性范围内滚动摩擦的主要原因. 为了解释这种情况下的滚动特性, 他们提出了弹性滞后是滚动摩擦的主要原因.
Flom和Bueche[22]采用具有黏滞特性的Voigt模型对弹性范围内材料的弹性滞后进行了研究. 根据该模型, 当坚硬球体或圆柱体在软基上以一定速度滚动时, 接触面的后缘将与软基脱离形成接触面压力的不均匀分布. 而滚动阻力力矩正是源自这种不对称分布的接触压力. 他们认为, 相对坚硬的球体在较软的基材上的滚动摩擦力会随滚动速度的变化而改变. 对于黏弹性材料的滚动颗粒, 滚动阻力主要是由接触面上体积变形产生的能量耗散引起的, 而与表面附着力的关系较小[23]. 在该研究中, 弹性滞后模型被当作具有黏弹特性的力学模型, 是速度相关型的. 当滚动速度非常小时, 滚动阻力接近于零. 同样的, 对于在硬质基体材料上滚动的软黏弹性球, 滚动摩擦阻力取决于滚动速度, 当滚动速度为零时, 滚动阻力也将为零[24]. 但是, 当采用这种速度相关型黏弹性模型研究颗粒堆积稳定时, 若颗粒处于稳定状态, 滚动阻力将会消失. 显然, 这种稳定机制是有问题的. 当滚动速度为零时, 由于阻力的消失颗粒体系稳定性会降低, 较颗粒滚动速度不为零时更容易崩塌是不合理的.
Greenwood等[23]对橡胶薄壁管进行了扭转和拉伸共同作用试验, 来研究材料的弹性滞后现象. 试验结果显示, 橡胶具有与加载速度无关的弹性滞后. 当在弹性范围内进行加载和卸载时, 由于应变变化落后于应力, 使应力应变曲线的加载和卸载过程不一致, 形成了闭环, 进而产生了能量耗散. 因此, 对于弹性材料, 弹性滞后引起的滚动阻力包括速度相关型的和速度无关型的两部分.
1979年, Cundall和Strack[25]提出了离散元方法(DEM)并迅速引起关注[26-27], 它可以方便地再现物理试验过程细节, 并且能有效降低试验成本, 目前已成为一种被广泛接受的数值计算方法. 离散元模型中滚动阻力力学模型的建立在颗粒体系力学行为模拟中至关重要.
为了研究在剪切带试验中观察到的颗粒间的巨大空隙和高旋转梯度, Iwashita和Oda[28]在DEM模型中考虑了滚动阻力, 建立了改进的离散元模型(MDEM). 该模型滚动阻力由转动方向的弹簧、黏壶、非承拉节点和摩擦型滑阻器元件表达, 其中由弹簧和黏壶组成的Voigt模型, 是一种速度相关型的阻力元件. 而速度无关型的滚动阻力通常由摩擦型滑阻器元件表达. 研究结果表明该模型可以有效预测剪切带的剪胀行为. MDEM是目前应用比较成功的离散元模型, 在此基础上, 许多学者开展了相关的计算应用与模型改进工作[9,29-30]. 在本文称MDEM模型为常规DEM模型.
目前对于滚动阻力与滚动速度之间的相关性仍没有统一的认识, 在离散元模拟中通常采用两种典型的滚动阻力公式[31-33]: 第一种滚动力矩与滚动速度无关, 力矩大小正比于法向接触力, 与滚动方向相反. 当颗粒间接触力恒定时, 滚动力矩是一个恒定值; 第二种滚动力矩与相对角速度成正比, 表现为一种黏滞力特征. 离散元模拟结果显示按照速度相关型黏滞滚动阻力公式不能很好地模拟颗粒的稳定堆积. 而单独按照速度无关型恒定的滚动阻力公式虽然模拟颗粒堆积稳定性有所提高, 但颗粒临近稳定静止时会在平衡位置往复振动, 此时滚动阻力大小不变, 方向正负变化, 微观上无法达到静止状态[19]. 因此, 在离散元算法中, 摩擦型滑阻器元件表达的速度无关型的滚动阻力当颗粒临近静止时要退出工作, 只有速度相关型黏滞力元件工作. 所以, 在颗粒堆积问题离散元模拟中, 弹性滞后引起的速度无关型滚动阻力是建立颗粒临近静止状态颗粒接触力学模型的一个重要思路.
直接采用滚动阻力公式表达的力学元件建立离散元模型应用方便, 但由于机理上缺乏深入认识, 一些滚动阻力参数确定往往通过经验或试算得到. 滚动阻力通常较小, 通过试验方法直接识别滚动阻力参数的难度较大[34].
笔者提出了通过圆形颗粒在刚性平面上滚动停止前的往复摆动现象测量滚动阻力参数的方法,研发了一个颗粒微动力光学试验检测系统[35], 可通过测量颗粒的往复摆动曲线识别Voigt模型的转动刚度与黏滞阻尼参数. 试验研究发现通过识别参数的常规DEM模型计算得到的颗粒摆动位移曲线与试验曲线在临近静止时刻吻合程度剧烈下降, 计算得到的颗粒稳定时间较试验长, 且试验曲线在临近静止时刻出现摆动频率的改变, 常规DEM模型不能从机理上解释这一现象.
为此, 本文基于弹性滞后理论研究建立了一种滞弹簧力学元件, 将与速度无关的材料弹性滞后特性引入, 提出一种滞弹性滚动阻力模型, 以此建立对试验中颗粒临近静止状态滞弹性耗能机理的解释. 与传统DEM模型相比, 改进后的滞弹性滚动阻力模型计算结果与试验结果更为符合, 验证了滞弹簧滚动阻力模型的有效性.
1. 弹性滞后
弹性范围内材料加载卸载时, 应变往往落后于应力, 使得应力−应变加载线与卸载线不重合围成一封闭回线, 形成弹性滞后现象. 与速度相关的弹性滞后现象通常表达为弹性材料的黏性行为, 而与速度无关的弹性滞后源于材料加载与卸载过程应力应变不能同步, 这一现象的主要原因在于材料分子间的相互作用和弛豫时间, 与加载速度无关. 通过薄壁橡胶管的纯拉伸试验, Tabor[21]获得应力−应变滞回曲线. 当外力没有达到极限应变卸载时, 应力−应变曲线上将形成一个转折点(εrev1, σrev1), 如图1所示. 当卸载不完全, 再次加载时, 应变将在最后一次卸载的终点(εrev2, σrev2)继续加载. 定义弹性滞后上升曲线为加载曲线, 下降曲线定义为卸载曲线, 并以(εe, σe)表示弹性极限点. 因此, 可以通过如下幂指数表达式定义应力和应变之间的关系
加载
$$ {\sigma }_{\rm{load}}(\varepsilon,{\varepsilon }_{\rm{rev}})=({\sigma }_{\rm{e}}-{\sigma }_{\rm{rev}})\frac{\varepsilon -{\varepsilon }_{\rm{rev}}}{{\varepsilon }_{\rm{e}}-{\varepsilon }_{\rm{rev}}}+{\sigma }_{\rm{rev}}$$ (1) 卸载
$${\sigma _{{\rm{unload}}}}(\varepsilon,{\varepsilon _{{\rm{rev}}}}) = {\sigma _{{\rm{rev}}}}\left[ {1 - {{\left(1 - \frac{\varepsilon }{{{\varepsilon _{{\rm{rev}}}}}}\right)}^\beta }} \right]$$ (2) 若加载起始点是(0, 0), 加载公式可简化为
$${\sigma _{{\rm{load}}}}(\varepsilon,{\varepsilon _{{\rm{rev}}}}) = {\sigma _{\rm{e}}}{\left(\frac{\varepsilon }{{{\varepsilon _{\rm{e}}}}}\right)^\beta }$$ (3) 当应变ε ∈ (0, εe), 表示材料处在弹性范围. 定义参数β ∈ (0, 1), 表示应变滞后于应力的程度. β值越接近0, 加载曲线和卸载曲线所包围的面积越大, 弹性滞后所引起的能量耗散也就越大.
为定义在复杂应力下的弹性滞后的能量耗散过程, 可根据加载和卸载构造相应的能量耗散过程. 以单个完整加载或单个完整卸载定义为一个子过程, 将整个荷载过程分为多个子过程组合, 可表达为
$$\varepsilon {\rm{ = }}\Delta {\varepsilon _{\rm{1}}}{\rm{ + }}\Delta {\varepsilon _{\rm{2}}}{\rm{ + }} \cdots {\rm{ + }}\Delta {\varepsilon _n}{\rm{ + }}\Delta \varepsilon {\rm{ = }}\sum\limits_{k = 1}^n {\Delta {\varepsilon _{{k}}} + \Delta \varepsilon } $$ (4) 其中
$$\Delta {\varepsilon _{{k}}} = {\varepsilon _{{{{\rm{rev}}(k)}}}} - {\varepsilon _{{{{\rm{rev}}(k - 1)}}}}$$ (5) Δε是不足一个完整子过程的多余应变. 子过程的能量密度表达式可写为
加载
$${e_{{k}}} = \int_{{\varepsilon _{{\rm{rev}}({{k - 1}})}}}^{{\varepsilon _{{\rm{rev}}({{k - 1}})}} + \Delta {\varepsilon _{{k}}}} {{\sigma _{{\rm{load}}}}(} \varepsilon ,{\varepsilon _{{{{\rm{rev}}}}({{k}})}}){\rm{d}}\varepsilon $$ (6) 卸载
$${e_{{{k + 1}}}} = \int_{{\varepsilon _{{{{\rm{rev}}(k)}}}}}^{{\varepsilon _{{{{\rm{rev}}(k)}}}} + \Delta {\varepsilon _{{{(k + 1)}}}}} {{\sigma _{{\rm{unload}}}}(\varepsilon,{\varepsilon _{{{{\rm{rev}}(k + 1)}}}}){\rm{d}}\varepsilon } $$ (7) 而整个过程的能量密度可通过每个子过程累加得到
$$ \begin{split} & {e_{{\rm{sum}}}} = \int_0^{\Delta {\varepsilon _1}} {{\sigma _{{\rm{load}}}}{\rm{d}}\varepsilon + \int_{{\varepsilon _{{\rm{rev(1)}}}} + \Delta {\varepsilon _2}}^{\Delta {\varepsilon _2}} {{\sigma _{{\rm{unload}}}}{\rm{d}}} } \varepsilon + \cdots +\\ & \int_{{\varepsilon _{{{{\rm{rev}}(n - 2)}}}}}^{{\varepsilon _{{{{\rm{rev}}(n - 2)}}}} + \Delta {\varepsilon _{{{n - 1}}}}} {{\sigma _{{\rm{load}}}}{\rm{d}}\varepsilon + \int_{{\varepsilon _{{{{\rm{rev}}(n - 1)}}}}}^{{\varepsilon _{{{{\rm{rev}}(n - 1)}}}} + \Delta {\varepsilon _2}} {{\sigma _{{\rm{unload}}}}{\rm{d}}} } \varepsilon+ \int_{{\varepsilon _{{{{\rm{rev}}(n)}}}}}^{\Delta \varepsilon } {\sigma {\rm{d}}} \varepsilon \\ \end{split} $$ (8) 进一步可简化为
$$ \begin{split} {e_{{\rm{sum}}}} = &\sum\limits_{k = 1}^n {\left[ {\int_{{\varepsilon _{{{{\rm{rev}}(k - 1)}}}}}^{{\varepsilon _{{{{\rm{rev}}(k - 1)}}}} + \Delta {\varepsilon _{{k}}}} {{\sigma _{{\rm{load}}}}(\varepsilon,{\varepsilon _{{{{\rm{rev}}}}}}){\rm{d}}\varepsilon } } \right.}+\\ & \left. {{\rm{ }}\int_{{\varepsilon _{{{{\rm{rev}}(k)}}}}}^{{\varepsilon _{{{{\rm{rev}}(k)}}}} + \Delta {\varepsilon _{{{k + 1}}}}} {{\sigma _{{\rm{unload}}}}(\varepsilon,{\varepsilon _{{{{\rm{rev}}}}}}){\rm{d}}\varepsilon } } \right] + \int_{{\varepsilon _{{{{\rm{rev}}(n)}}}}}^{\Delta \varepsilon } {\sigma {\rm{d}}\varepsilon } \end{split} $$ (9) $e$ 表示每个荷载过程的能量密度,$\displaystyle\int_{{\varepsilon _{{{{\rm{rev}}(n)}}}}}^{\Delta \varepsilon } {\sigma {\rm{d}}\varepsilon }$ 表示不能构成一个子过程的多余应变的能量密度.2. 滞弹簧与HDEM模型
根据接触方向, 常规DEM接触模型[28]可分为法向接触模型, 切向接触模型和滚动阻力模型, 如图2所示. 滚动方向阻力模型由滚动弹簧、滚动黏壶、摩擦器和非承拉节点组成. 滚动模型所提供的滚动阻力可表示为
$$ M=\mathrm{min}({K}_{\rm{r}}\theta +{c}_{\rm{r}}\dot{\theta },\;\; {\mu }_{\rm{r}}{F}_{\rm{n}})$$ (10) 式中Kr是弹簧刚度系数, cr是滚动阻尼系数, μr是滚动摩擦系数.
从上式可看出, 颗粒发生持续滚动时, 滚动力矩总是等于最大摩擦力矩μrFn, 当颗粒滚动不能持续时, 滚动阻力小于最大摩擦力矩μrFn, 常规DEM模型所提供的阻力值与滚动角速度有关. 滚动角速度越大, 滚动模型所提供的阻力值越大. 滚动弹簧是线弹性的, 不能耗散能量, 动能的耗散仅依靠黏壶黏滞力做功实现. 为表征滚动摩擦中速度无关的摩擦力, 根据上述建立弹性滞后的应力−应变关系, 提出图3所示的滞弹簧元件.
滞弹簧的角位移与弹性力不同于常规DEM滚动模型中滚动弹簧的线弹性关系. 参照弹性滞后应力应变表达, 滞弹簧的加载和卸载遵循以下关系
加载
$$ {F}_{\rm{load}}(\varDelta,{\varDelta }_{\rm{rev}})=({F}_{\rm{e}}-{F}_{\rm{rev}})\left(\frac{\varDelta -{\varDelta }_{\rm{rev}}}{{\varDelta }_{\rm{e}}-{\varDelta }_{\rm{rev}}}\right)+{F}_{\rm{rev}}$$ (11) 卸载
$${F_{{\rm{unload}}}}(\varDelta,{\varDelta _{{\rm{rev}}}}) = {F_{{\rm{rev}}}}\left[ {1 - {{\left(1 - \frac{\varDelta }{{{\varDelta _{{\rm{rev}}}}}}\right)}^\beta }} \right]$$ (12) 加载时, 若加载起点为原点(0, 0), 加载方程可简化为
$${F_{{\rm{load}}}}(\varDelta,{\varDelta _{{\rm{rev}}}}) = {F_{\rm{e}}}{\left(\frac{\varDelta }{{{\varDelta _{\rm{e}}}}}\right)^\beta }$$ (13) 式中Δ表示滞弹簧的位移变形, Δe是滞弹簧弹性变形量的极限值, F表示滞弹簧的恢复力值, Fe是弹簧弹性力的极限值.
滞弹簧耗散的能量可以表示为
$$ \begin{split} {E_{{\rm{sum}}}} = &\sum\limits_{k = 1}^n {\left[ {\int_{{\varDelta _{{{{\rm{rev}}(k - 1)}}}}}^{{\varDelta _{{{{\rm{rev}}(k - 1)}}}} + {\varDelta _{{k}}}} {{F_{{\rm{load}}}}(\delta,{\varDelta _{{{{\rm{rev}}}}}}){\rm{d}}\delta } } \right.}+\\ & \left. {{\rm{ }}\int_{{\varDelta _{{{{\rm{rev}}(k)}}}}}^{{\varDelta _{{{{\rm{rev}}(k)}}}} + {\varDelta _{{{k + 1}}}}} {{F_{{\rm{unload}}}}(\delta,{\varDelta _{{{{\rm{rev}}}}}}){\rm{d}}\delta } } \right] + \int_{{\varDelta _{{{{\rm{rev}}(n)}}}}}^\varDelta {F{\rm{d}}\delta } \end{split} $$ (14) 其中
$$\varDelta = {\varDelta _{{{{\rm{rev}}(k)}}}} - {\varDelta _{{{{\rm{rev}}(k - 1)}}}}$$ (15) Esum表示整个滞弹簧运动过程中由于弹性滞后效应而耗散的能量, Δrev表示加载和卸载过程中滞弹簧转折点的位移值, Frev表示加载和卸载过程中滞弹簧转折点的力值.
当滚动颗粒在平衡位置往复摆动时, 滚动过程可分为正向加载、正向卸载、负向加载、负向卸载四个阶段. 根据滞弹簧变形与滚动角速度, 式(11)和式(12)计算卸载与卸载时滞弹簧元件的弹性恢复力, 负向时弹性恢复力取负.
将滞弹簧与黏壶、摩擦器及非承拉节点等元件进行组合, 提出一种新的滞弹性滚动阻力离散元模型(HDEM), 如图4所示. 滞弹簧可以表征颗粒堆积稳定过程中与运动速度无关的滚动阻力耗能.
3. 验证
3.1 滞弹簧有效性验证
当一个在刚性平面上的滚动圆柱试件速度减小到一定程度后, 会在一个平衡位置往复摆动. 由于动能耗散, 摆动幅度逐渐减小直至为零, 如图5所示. 常规DEM滚动阻力模型参数可以通过测量试件的摆动来识别[36].
基于常规DEM模型, 自由滚动的圆柱试件的运动平衡方程为
$$\sum {{F_{{x}}} = 0} \qquad\qquad\qquad$$ $$m\ddot x + {c_{\rm{s}}}\dot x{\rm{ + }}{K_{\rm{s}}}x + \frac{{{M_{\rm{r}}}}}{R} = 0$$ (16) $${M_{\rm{r}}} = {J_{\rm{z}}}\ddot \theta + {c_{\rm{r}}}\dot \theta + {K_{\rm{r}}}\theta \qquad $$ (17) 或
$$\sum {{M_{\rm{z}}} = 0}\qquad\qquad\qquad\qquad $$ $$({J_{\rm{c}}} + m{R^2})\ddot \theta + {c_{\rm{r}}}\dot \theta + {K_{\rm{r}}}\theta + {F_{{x}}}R = 0$$ (18) $${J_{\rm{z}}} = {J_{\rm{c}}} + m{R^2}\qquad\qquad\qquad\qquad$$ (19) $${F_{{x}}} = m\ddot x\qquad\qquad\qquad\qquad\qquad$$ (20) 其中, θ为滚动角位移, Ks和Kr分别表示切向和滚动弹簧刚度系数, cs和cr分别表示切向和滚动方向阻尼系数, Jz是对接触点的转动惯量, Jc是对颗粒形心处的转动惯量, Fx是水平惯性力, Mr是惯性力矩, R为圆柱试件半径.
方程(18)为有阻尼振动方程, 可得到滚动角位移表达式为
$$\theta = {{\rm{e}}^{ - \xi \omega t}}A\sin \left( {\omega t + \alpha } \right)$$ (21) 其中ω为振动圆频率, ξ为阻尼比, A为最大幅值, α为相位角. 而阻尼系数和滚动刚度可表达为
$${c_{\rm{r}}} = 2{J_{\rm{z}}}\omega \xi $$ (22) $${K_{\rm{r}}} = {\omega ^2}{J_{\rm{z}}}$$ (23) 因此, 只需测得颗粒的往复摆动曲线, 按照式(21) ~ 式(23)可识别出振动圆频率ω和阻尼比ξ, 进而识别出阻尼系数cr和滚动刚度Kr.
采用激光位移传感器微振动位移检测试验装置, 可试验测得圆柱体试件往复摆动曲线, 试验装置如图6.
通过检测试验, 我们可以得到圆柱滚动的时间历程曲线, 如图7所示.
![]()
图 7 圆柱滚动角位移时程曲线Figure 7. Angular displacement variation versus time for a rolling cylinder为验证滞弹簧元件有效性, 分别使用常规DEM滚动阻力模型与HDEM模型对颗粒的纯滚动过程进行了离散元模拟, 如图8所示. 图9是圆柱体在平衡位置摆动过程的数值模拟结果. 图中常规DEM模型中弹簧角位移与弹性力的线性变化关系, 而HDEM模型中滞弹簧的角位移与弹性力形成滞回环, 与所建立的位移−荷载公式一致.
图10中黑色曲线是采用激光位移传感试验测得的滚动弹簧刚度值, 并使用常规DEM模型进行离散元数值模拟得到的动能衰减包络线. 紫色曲线是将常规DEM模型中阻尼系数放大1.5倍后的动能衰减包络线, 蓝色曲线是HDEM模型模拟出的动能衰减包络线. 在摆动起始时刻将黏壶的阻尼系数设置为0, 可以看出将阻尼系数放大后, 动能衰减曲线近似向下平移. 而HDEM模拟的动能衰减曲线在起始时刻7 s时与蓝色曲线交汇, 交汇前位于蓝色曲线上方, 交汇后位于下方, 说明在临近静止过程中滞弹簧的耗能大于黏壶.
参数β能反映滞弹簧的耗能能力, 其值越小耗能能力越强. 通过与试验数据进行对比, 经过试算拟合可得β值. 对10个不同材质圆柱形试件测量结果拟合得到的β值如图11所示. 聚氨酯圆柱β平均值为0.844, 铝圆柱β平均值为0.874. 聚氨酯的平均值要小于铝的平均值, 分析原因为聚氨酯材料质地较软, 在弹性滞后过程中会产生更多的能耗.
虽然材料的弹性滞后耗能同样存在法向与切向. 但由于离散元接触模型中法向与切向刚度比转动方向刚度大得多, 法向与切向振动频率高, 速度相关型的黏壶元件耗能要比滞弹簧耗能大得多, 因此HDEM模型中虑法向、切向滞弹簧与转动方向滞弹簧相比作用不显著, 在模型中可不考虑使用.
3.2 耗能分析
图12为橡胶圆柱自由滚动激光位移传感器试验、常规DEM模型和HDEM模型离散元数值模拟得到的时间−相对位移图. 从图中可以看出, 常规DEM和HDEM在振荡早期的动能衰减差异不大. 但从整体上, 常规DEM模型计算得到的滚动到静止时间较试验结果要长, 显示常规DEM模型不易达到静止稳定; HDEM模型与试验结果更为接近. 因此可以得出HDEM模型在接近静止时具有更强的能量耗散能力, 与试验结果吻合更好. 试验曲线不光滑是由于仪器采集误差造成的.
HDEM模型中有滞弹簧和黏壶两个耗能元件. 为分析速度无关滚动阻力与速度相关滚动阻力关系, 提取图12中数值模拟结果,得到体系总动能变化如图13(a)与黏壶作用耗能变化如图13(b). 可以看出, 体系动能不断衰减, 黏壶能量耗散逐渐减小.
图14(a)是HDEM模型中滞弹簧耗能与黏壶耗能变化. 图14(b)是他们的比值随时间变化趋势. 可以看出随着滚动速度降低, HDEM模型中滞弹簧元件所代表的与速度无关的能量耗散比例越来越大.
3.3 频率拟合结果
选取4组橡胶圆柱和铝圆柱试件的试验检测数据和常规DEM模型、HDEM模型计算结果数据进行对比(如图15和图16). 由于常规DEM模型表达的振动具有频率不变特性, 常规DEM模型与观测试验结果对比可以看出, 摆动速度接近0时, 试验中圆柱体摆动频率有增大现象; HDEM模型的数值模拟结果与试验结果一致, 在摆动速度接近0时同样表现出频率增大的现象. HDEM模型能够模拟滚动试件临近静止时刻摆动频率变高的现象.
![]()
图 15 橡胶圆柱数值模拟与试验对比结果Figure 15. Comparison of the numerical simulation and experimental results for the rubber cylinder![]()
图 16 铝圆柱数值模拟与试验对比结果Figure 16. Comparison of the numerical simulation and experimental results for the aluminum cylinder从模型计算结果与试验结果对比可以看出, 本文建立的滞弹簧滚动阻力模型能够很好地模拟颗粒滚动速度接近于零状态时的能量耗散过程, 能够对颗粒滚动阻力现象进行合理的滚动阻力机理解释. 建立的滞弹性滚动阻力可为颗粒材料堆积稳定问题研究提供方法.
4. 结 论
本文研究建立了滚动阻力滞弹性表达的HDEM模型, 与常规DEM模型的数值模拟结果进行了比较. 通过圆柱试件在平台上的自由滚动试验验证了该模型的有效性, 并得出以下结论:
(1) HDEM滚动阻力模型模拟结果与试验现象吻合, 能较好地解释颗粒在临近静止阶段的能量耗散特性;
(2) 弹性滞后引起的滚动阻力包括速度相关型的和速度无关型的两部分, 滚动试件临近静止时刻, 与速度无关的阻力成分占比越来越大;
(3) HDEM滚动阻力模型能较好地拟合橡胶材料与铝材料圆柱试件的摆动频率, 且能很好地模拟试验中临近静止时刻频率变高的现象.
-
![]()
图 7 圆柱滚动角位移时程曲线
Figure 7. Angular displacement variation versus time for a rolling cylinder
![]()
图 15 橡胶圆柱数值模拟与试验对比结果
Figure 15. Comparison of the numerical simulation and experimental results for the rubber cylinder
-
[1] Zhang J, Zou GY, Zhang N, et al. Dynamic analysis of a vehicle with leaf spring based on the hysteresis model. International Journal of Vehicle Performance, 2018, 4(3): 282-304 doi: 10.1504/IJVP.2018.095309
[2] Herrmann HJ, Luding S. Modeling granular media on the computer. Continuum Mechanics and Thermodynamics, 1998, 10(4): 189-231 doi: 10.1007/s001610050089
[3] 王云霞, 梁志杰, 张东兴等. 基于离散元的玉米种子颗粒模型种间接触参数标定. 农业工程学报, 2016, 32(22): 36-42 (Wang Yunxia, Liang Zhijie, Zhang Dongxing, et al. Calibration method of contact characteristic parameters for corn seeds based on EDEM. Transactions of the Chinese Society of Agricultural Engineering, 2016, 32(22): 36-42 (in Chinese) doi: 10.11975/j.issn.1002-6819.2016.22.005 [4] Höhner D, Wirtz S, Scherer V. Experimental and numerical investigation on the influence of particle shape and shape approximation on hopper discharge using the discrete element method. Power Technology, 2013, 235: 614-627 doi: 10.1016/j.powtec.2012.11.004
[5] 孙其诚, 王光谦. 流体动力学//中国科学院编. 自然界和工程中的颗粒物质与颗粒流. 北京: 科学出版社, 2014: 54-71 (Sun Qicheng, Wang Guangqian. Granular matter and granular flow in nature and engineering//Chinese Academy of Sciences ed. Fluid Dynamics. Beijing: Science Press, 2014: 54-71 (in Chinese))
[6] Bardet JP, Huang Q. Numerical modeling of micro-polar effects in idealized granular materials. American Society of Mechanical Engineers, 1992, 37: 85-92
[7] Morgan JK. Capturing physical phenomena in particle dynamics simulations of granular fault gouge//ACES Workshop Proceedings, 2003: 23-30
[8] Sakaguchi H, Ozaki E, Igarashi T. Plugging of the flow of granular materials during the discharge from a silo. International Journal of Modern Physics B, 1993, 7: 1949-1963 doi: 10.1142/S0217979293002705
[9] Jiang MJ, Yu HS, Harris D. A novel discrete model for granular material incorporating rolling resistance. Computers and Geotechnics, 2005, 32(5): 340-357 doi: 10.1016/j.compgeo.2005.05.001
[10] 魏新容, 段绍臻, 孙金龙等. 基于碰撞模型的斜坡滚石颗粒速度预测. 力学学报, 2020, 52(3): 707-715 (Wei Xinrong, Duan Shaozhen, Sun Jinlong, et al. Velocity prediction of slope rolling stone particle based on collision model. Chinese Journal of Theoretical and Applied Mechanics, 2020, 52(3): 707-715 (in Chinese) doi: 10.6052/0459-1879-20-039 [11] Zhang Y, Mollon G, Descartes S. Significance of third body rheology in friction at a dry sliding interface observed by a multibody meshfree model: Influence of cohesion between particles. Tribology International, 2020, 145: 106188 doi: 10.1016/j.triboint.2020.106188
[12] Zhou Y, Shi Z, Zhang Q, et al. Damming process and characteristics of landslide-debris avalanches. Soil Dynamics and Earthquake Engineering, 2019, 121: 252-261 doi: 10.1016/j.soildyn.2019.03.014
[13] Lopez RDF, Larsson S, Silfwerbrand J. A discrete element material model including particle degradation suitable for rockfill embankments. Computers and Geotechnics, 2019, 115: 103166 doi: 10.1016/j.compgeo.2019.103166
[14] 李刚, 李恩兴. 某堆积体斜坡变形特征及稳定性分析. 地质灾害与环境保护, 2020, 31(3): 66-71+76 (LI Gang, LI Enxing. A study of deformation features and stability of a talus landslide. Journal of Geologica Hazards and Environment Preservation, 2020, 31(3): 66-71+76 (in Chinese) doi: 10.3969/j.issn.1006-4362.2020.03.010 [15] 邹宇雄, 马刚, 李易奥等. 抗转动对颗粒材料组构特性的影响研究. 岩土力学, 2020, 41(8): 2829-2838 (Zou Yuxiong, MA Gang, LI Yiao, et al. Impact of rotation resistance on fabric of granular materials. Rock and Soil Mechanics, 2020, 41(8): 2829-2838 (in Chinese) [16] 祁原, 黄俊杰, 陈明祥. 可破碎颗粒体在动力载荷下的耗能特性. 力学学报, 2015, 47(2): 252-259 (Qi Yuan, Huang Junjie, Chen Mingxiang. Energy dissipation characteristics of crushable granules under dynamic excitations. Journal of Theoretical and Applied Mechanics, 2015, 47(2): 252-259 (in Chinese) doi: 10.6052/0459-1879-14-145 [17] 魏志刚, 陈海波. 一种新的橡胶材料弹性本构模型. 力学学报, 2019, 51(2): 473-483 (Wei Zhigang, Chen Haibo. A new elastic model for rubber-like materials. Chinese Jounal of Theoretical and Applied Mechanics, 2019, 51(2): 473-483 (in Chinese) doi: 10.6052/0459-1879-18-303 [18] 谈炳东, 许进升, 孙朝翔等. 短纤维增强三元乙丙橡胶横观各向同性黏—超弹性本构模型. 力学学报, 2017, 49(3): 677-684 (Tan Bingdong, Xu Jinsheng, Sun Chaoxiang, et al. A transversely isotropic visco-hyperelastic constitutive model for short fiber reinforced EPDM. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(3): 677-684 (in Chinese) doi: 10.6052/0459-1879-16-380 [19] Ai J, Chen JF, Rotter JM, et al. Assessment of rolling resistance models in discrete element simulations. Powder Technology, 2011, 206(3): 269-282 doi: 10.1016/j.powtec.2010.09.030
[20] Reynolds O. On rolling-friction. Philosophical Transactions of the Royal Society of London, 1876, 166(1): 155-174
[21] Tabor D. The mechanism of rolling friction. Philosophical Magazine Series 7, 1952, 43(345): 1055-1059 doi: 10.1080/14786441008520246
[22] Flom DG, Bueche AM. Theory of rolling friction for spheres. Journal of Applied Physics, 1960, 3(3): 1725-1730
[23] Greenwood JA, Minshall H, Tabor D. Hysteresis losses in rolling and sliding friction. Proceedings of the Royal Society of London, 1961, 259(1299): 480-507
[24] Brilliantov NV, Poschel T. Rolling friction of a viscous sphere on a hard plane. Europhysics Letters, 1998, 42(5): 511-516 doi: 10.1209/epl/i1998-00281-7
[25] Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Géotechnique, 1979, 29(1): 47-65
[26] 刘嘉英, 周伟, 马刚等. 颗粒材料三维应力路径下的接触组构特性. 力学学报, 2019, 51(1): 26-35 (Liu Jiaying, Zhou Wei, Ma Gang, et al. Contact fabric characteristics of granular materials under three dimensionalstress paths. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(1): 26-35 (in Chinese) doi: 10.6052/0459-1879-18-338 [27] 刘巨保, 王明, 王雪飞等. 颗粒群碰撞搜索及 CFD-DEM 耦合分域求解的推进算法研究. 力学学报, 2021: 1-16 (Liu Jubao, Wang Ming, Wang Xuefei, et al. Research on particle swarm collision search and advancement algorithm for CFD-DEM coupling domain solving. Chinese Journal of Theoretical and Applied Mechanics, 2021: 1-16 (in Chinese) [28] Iwashita K, Oda M. Rolling resistance at contacts in simulation of shear band development by DEM. Journal of Engineering Mechanics, 1998, 124(3): 285-192 doi: 10.1061/(ASCE)0733-9399(1998)124:3(285)
[29] Rorato R, Arroyo M, Gens A, et al. Image-based calibration of rolling resistance in discrete element models of sand. Computers and Geotechnics, 2021, 131: 103929 doi: 10.1016/j.compgeo.2020.103929
[30] Zamir S, Mehari T, David W. A coupled sliding and rolling friction model for DEM calibration. Journal of Terramechanics, 2017, 72: 9-20 doi: 10.1016/j.jterra.2017.03.003
[31] Zhu HP, Yu AB. A theoretical analysis of the force models in discrete element method. Powder Technology, 2005, 161(2): 122-129
[32] Zhou YC, Wright BD, Yang RY, et al. Rolling friction in the dynamic simulation of sand-pile formation. Physica A Statistical Mechanics & Its Applications, 1999, 269(2-4): 536-553
[33] Brilliantov NV, Poeschel T. Rolling friction of a viscous sphere on a hard plane. Europhys. Lett., 2007, 42(5): 511-516
[34] 孙珊珊, 苏勇, 季顺迎. 颗粒滚动-滑动转换机制及摩擦系数的试验研究. 岩土力学, 2009, 30(S1): 110-115 (Sun Shanshan, Su Yong, Ji Shunying. Experimental study of rolling-sliding transiton and friction coefficients of particles. Rock and Soil Mechanics, 2009, 30(S1): 110-115 (in Chinese) [35] 高政国, 王佃瑞, 张雅俊. 一种颗粒滚动阻力模型参数的动力试验测定装置. 中国专利, 202010913790.8. 2020-09-03 (Gao Zhengguo, Wang Dianrui, Zhang Yajun. A dynamic test method for parameter identification of particle rolling resistance model. Chinese Patend, 202010913790.8. 2020-09-03 (in Chinese))
[36] 王佃瑞. 颗粒滚动阻力与堆积稳定特性研究. [硕士论文]. 北京: 北京航空航天大学, 2020 (Wang Dianrui. Study on rolling resistance and stacking stability characteristics of particles. [Master Thesis]. Beijing: Beihang University, 2020 (in Chinese))
-
期刊类型引用(1)
1. 董朋昆，高政国. 颗粒滚动阻力本构参数标定试验研究. 上海理工大学学报. 2023(06): 574-583 . 
百度学术
其他类型引用(2)
下载:
计量
- 文章访问数:
- HTML全文浏览量:
- PDF下载量:
- 被引次数: 3
