APP下载

一种滞弹簧耗能的新型离散元滚动阻力模型研究1)

2021-11-10高政国董朋昆张雅俊孙卉竹

力学学报 2021年9期
关键词:圆柱弹簧阻力

高政国 董朋昆 张雅俊 孙卉竹 迪 亚

* (北京航空航天大学交通科学与工程学院,北京 100191)

† (中国矿业大学(北京)力学与建筑工程学院,北京 100083)

引言

滚动阻力问题涉及车辆工程[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)表示弹性极限点.因此,可以通过如下幂指数表达式定义应力和应变之间的关系

图1 弹性滞后示意图Fig.1 Schematic diagram of elastic hysteresis

加载

当应变ε∈ (0,εe),表示材料处在弹性范围.定义参数β∈ (0,1),表示应变滞后于应力的程度.β值越接近0,加载曲线和卸载曲线所包围的面积越大,弹性滞后所引起的能量耗散也就越大.

为定义在复杂应力下的弹性滞后的能量耗散过程,可根据加载和卸载构造相应的能量耗散过程.以单个完整加载或单个完整卸载定义为一个子过程,将整个荷载过程分为多个子过程组合,可表达为

Δε是不足一个完整子过程的多余应变.子过程的能量密度表达式可写为

加载e表示每个荷载过程的能量密度,表示不能构成一个子过程的多余应变的能量密度.

2 滞弹簧与HDEM 模型

根据接触方向,常规DEM 接触模型[28]可分为法向接触模型,切向接触模型和滚动阻力模型,如图2所示.滚动方向阻力模型由滚动弹簧、滚动黏壶、摩擦器和非承拉节点组成.滚动模型所提供的滚动阻力可表示为

图2 常规DEM 模型Fig.2 DEM model

式中Kr是弹簧刚度系数,cr是滚动阻尼系数,μr是滚动摩擦系数.

从上式可看出,颗粒发生持续滚动时,滚动力矩总是等于最大摩擦力矩μrFn,当颗粒滚动不能持续时,滚动阻力小于最大摩擦力矩μrFn,常规DEM 模型所提供的阻力值与滚动角速度有关.滚动角速度越大,滚动模型所提供的阻力值越大.滚动弹簧是线弹性的,不能耗散能量,动能的耗散仅依靠黏壶黏滞力做功实现.为表征滚动摩擦中速度无关的摩擦力,根据上述建立弹性滞后的应力−应变关系,提出图3所示的滞弹簧元件.

滞弹簧的角位移与弹性力不同于常规DEM 滚动模型中滚动弹簧的线弹性关系.参照弹性滞后应力应变表达,滞弹簧的加载和卸载遵循以下关系

加载

式中Δ表示滞弹簧的位移变形,Δe是滞弹簧弹性变形量的极限值,F表示滞弹簧的恢复力值,Fe是弹簧弹性力的极限值.

滞弹簧耗散的能量可以表示为

Esum表示整个滞弹簧运动过程中由于弹性滞后效应而耗散的能量,Δrev表示加载和卸载过程中滞弹簧转折点的位移值,Frev表示加载和卸载过程中滞弹簧转折点的力值.

当滚动颗粒在平衡位置往复摆动时,滚动过程可分为正向加载、正向卸载、负向加载、负向卸载四个阶段.根据滞弹簧变形与滚动角速度,式(11)和式(12)计算卸载与卸载时滞弹簧元件的弹性恢复力,负向时弹性恢复力取负.

将滞弹簧与黏壶、摩擦器及非承拉节点等元件进行组合,提出一种新的滞弹性滚动阻力离散元模型(HDEM),如图4 所示.滞弹簧可以表征颗粒堆积稳定过程中与运动速度无关的滚动阻力耗能.

图4 HDEM 模型Fig.4 HDEM model

3 验证

3.1 滞弹簧有效性验证

当一个在刚性平面上的滚动圆柱试件速度减小到一定程度后,会在一个平衡位置往复摆动.由于动能耗散,摆动幅度逐渐减小直至为零,如图5 所示.常规DEM 滚动阻力模型参数可以通过测量试件的摆动来识别[36].

图5 颗粒自由滚动示意图Fig.5 Particle free-rolling on a flat surface

基于常规DEM 模型,自由滚动的圆柱试件的运动平衡方程为

其中,θ为滚动角位移,Ks和Kr分别表示切向和滚动弹簧刚度系数,cs和cr分别表示切向和滚动方向阻尼系数,Jz是对接触点的转动惯量,Jc是对颗粒形心处的转动惯量,Fx是水平惯性力,Mr是惯性力矩,R为圆柱试件半径.

方程(18)为有阻尼振动方程,可得到滚动角位移表达式为

其中ω为振动圆频率,ξ为阻尼比,A为最大幅值,α为相位角.而阻尼系数和滚动刚度可表达为

因此,只需测得颗粒的往复摆动曲线,按照式(21)~式(23)可识别出振动圆频率ω和阻尼比ξ,进而识别出阻尼系数cr和滚动刚度Kr.

采用激光位移传感器微振动位移检测试验装置,可试验测得圆柱体试件往复摆动曲线,试验装置如图6.

图6 试验装置图Fig.6 Experimental device diagram

通过检测试验,我们可以得到圆柱滚动的时间历程曲线,如图7 所示.

图7 圆柱滚动角位移时程曲线Fig.7 Angular displacement variation versus time for a rolling cylinder

为验证滞弹簧元件有效性,分别使用常规DEM 滚动阻力模型与HDEM 模型对颗粒的纯滚动过程进行了离散元模拟,如图8 所示.图9 是圆柱体在平衡位置摆动过程的数值模拟结果.图中常规DEM 模型中弹簧角位移与弹性力的线性变化关系,而HDEM 模型中滞弹簧的角位移与弹性力形成滞回环,与所建立的位移−荷载公式一致.

图8 颗粒转动模型示意图Fig.8 Particle rotation model

图9 滚动过程Fig.9 Rolling process

图10 中黑色曲线是采用激光位移传感试验测得的滚动弹簧刚度值,并使用常规DEM 模型进行离散元数值模拟得到的动能衰减包络线.紫色曲线是将常规DEM 模型中阻尼系数放大1.5 倍后的动能衰减包络线,蓝色曲线是HDEM 模型模拟出的动能衰减包络线.在摆动起始时刻将黏壶的阻尼系数设置为0,可以看出将阻尼系数放大后,动能衰减曲线近似向下平移.而HDEM 模拟的动能衰减曲线在起始时刻7 s 时与蓝色曲线交汇,交汇前位于蓝色曲线上方,交汇后位于下方,说明在临近静止过程中滞弹簧的耗能大于黏壶.

图10 HDEM 与DEM 模拟动能变化Fig.10 Kinetic energy evolution with HDEM and DEM simulation

参数β能反映滞弹簧的耗能能力,其值越小耗能能力越强.通过与试验数据进行对比,经过试算拟合可得β值.对10 个不同材质圆柱形试件测量结果拟合得到的β值如图11 所示.聚氨酯圆柱β平均值为0.844,铝圆柱β平均值为0.874.聚氨酯的平均值要小于铝的平均值,分析原因为聚氨酯材料质地较软,在弹性滞后过程中会产生更多的能耗.

图11 橡胶圆柱与铝圆柱的β 值Fig.11 β values for the rubber and aluminum cylinder

虽然材料的弹性滞后耗能同样存在法向与切向.但由于离散元接触模型中法向与切向刚度比转动方向刚度大得多,法向与切向振动频率高,速度相关型的黏壶元件耗能要比滞弹簧耗能大得多,因此HDEM 模型中虑法向、切向滞弹簧与转动方向滞弹簧相比作用不显著,在模型中可不考虑使用.

3.2 耗能分析

图12 为橡胶圆柱自由滚动激光位移传感器试验、常规DEM 模型和HDEM 模型离散元数值模拟得到的时间−相对位移图.从图中可以看出,常规DEM 和HDEM 在振荡早期的动能衰减差异不大.但从整体上,常规DEM 模型计算得到的滚动到静止时间较试验结果要长,显示常规DEM 模型不易达到静止稳定;HDEM 模型与试验结果更为接近.因此可以得出HDEM 模型在接近静止时具有更强的能量耗散能力,与试验结果吻合更好.试验曲线不光滑是由于仪器采集误差造成的.

图12 DEM 滚动模型与HDEM 模型模拟结果Fig.12 Rolling angle versus time for DEM and HDEM model

HDEM 模型中有滞弹簧和黏壶两个耗能元件.为分析速度无关滚动阻力与速度相关滚动阻力关系,提取图12 中数值模拟结果,得到体系总动能变化如图13(a)与黏壶作用耗能变化如图13(b).可以看出,体系动能不断衰减,黏壶能量耗散逐渐减小.

图13 体系动能与黏壶耗能Fig.13 Kinetic energy and energy dissipated by the damper

图14(a)是HDEM 模型中滞弹簧耗能与黏壶耗能变化.图14(b)是他们的比值随时间变化趋势.可以看出随着滚动速度降低,HDEM 模型中滞弹簧元件所代表的与速度无关的能量耗散比例越来越大.

图14 滞弹簧与黏壶耗能Fig.14 Energy consumption ratio of hysteresis spring and damper

3.3 频率拟合结果

选取4 组橡胶圆柱和铝圆柱试件的试验检测数据和常规DEM 模型、HDEM 模型计算结果数据进行对比(如图15 和图16).由于常规DEM 模型表达的振动具有频率不变特性,常规DEM 模型与观测试验结果对比可以看出,摆动速度接近0 时,试验中圆柱体摆动频率有增大现象;HDEM 模型的数值模拟结果与试验结果一致,在摆动速度接近0 时同样表现出频率增大的现象.HDEM 模型能够模拟滚动试件临近静止时刻摆动频率变高的现象.

图15 橡胶圆柱数值模拟与试验对比结果Fig.15 Comparison of the numerical simulation and experimental results for the rubber cylinder

图16 铝圆柱数值模拟与试验对比结果Fig.16 Comparison of the numerical simulation and experimental results for the aluminum cylinder

从模型计算结果与试验结果对比可以看出,本文建立的滞弹簧滚动阻力模型能够很好地模拟颗粒滚动速度接近于零状态时的能量耗散过程,能够对颗粒滚动阻力现象进行合理的滚动阻力机理解释.建立的滞弹性滚动阻力可为颗粒材料堆积稳定问题研究提供方法.

4 结论

本文研究建立了滚动阻力滞弹性表达的HDEM模型,与常规DEM 模型的数值模拟结果进行了比较.通过圆柱试件在平台上的自由滚动试验验证了该模型的有效性,并得出以下结论:

(1) HDEM 滚动阻力模型模拟结果与试验现象吻合,能较好地解释颗粒在临近静止阶段的能量耗散特性;

(2) 弹性滞后引起的滚动阻力包括速度相关型的和速度无关型的两部分,滚动试件临近静止时刻,与速度无关的阻力成分占比越来越大;

(3) HDEM 滚动阻力模型能较好地拟合橡胶材料与铝材料圆柱试件的摆动频率,且能很好地模拟试验中临近静止时刻频率变高的现象.

猜你喜欢

圆柱弹簧阻力
圆柱的体积计算
鼻阻力测定在儿童OSA诊疗中的临床作用
“圆柱与圆锥”复习指导
析弹簧模型 悟三个性质
零阻力
别让摩擦成为学习的阻力
如何求串联弹簧和并联弹簧的劲度系数
时间弹簧
阻力不小 推进当循序渐进
圆柱壳的声辐射特性分析