基于LuGre 摩擦模型的接触约束法旋转柔性梁斜碰撞研究1)

范纪华 谷通顺 王明强 谌 宏 陈立威

∗(江苏科技大学机械工程学院,江苏镇江 212003)

†(江苏科技大学苏州理工学院,江苏张家港 215600)

∗∗(无锡交通高等职业技术学校机械工程学院,江苏无锡 214000)

引言

接触碰撞是工程界普遍存在的现象,例如航天器的对接和抓取[1-3]、太阳能帆板的展开[4-5]、旋转关节的间隙碰撞[6]、车辆碰撞[7]、火箭的回收着陆[8]、工程机械作业[9-10].由于碰撞过程存在法向碰撞和切向摩擦作用,具有非光滑、不连续的特性[10-11].因此,有必要对高效精确的接触碰撞建模方法进行研究.

在接触碰撞过程中存在法向接触速度跳跃、能量突变和切向摩擦的滑动−黏滞互相切换等问题.对于解决柔性体接触碰撞过程中含切向黏滞/滑动状态摩擦的碰撞非光滑问题,国内外学者提出了光滑化方法[11-13],该方法采用光滑曲线来描述Coulomb 摩擦模型滑动−黏滞的切换过程,有效解决了切换引起的稳定性问题,然而该方法存在失去黏滞时动力学特征的缺陷,针对光滑化方法存在的缺陷,国内外学者提出建立互补关系来求解接触碰撞非光滑问题[14-17].互补关系算法通过Hertz 接触理论描述法向接触力,在切向上建立正负加速和正负摩擦余量的互补关系从而通过求解出摩擦余量得出静摩擦力[14].然而互补算法需要计算互补方程同时还需要进行状态切换条件的判断,相对增加了计算求解时间[18],有学者[18-21]提出直接根据接触碰撞过程黏滞/滑移条件进行程序判断计算的方法.文献[18]对柔性滑块接触问题研究中,采用弹簧阻尼模型计算碰撞力,以微小速度差作为滑动和黏滞的切换判断条件,黏滞时直接以滑块受到的外力作为静摩擦力,并且对于多节点接触问题求解不需要建立和计算互补方程,通过两个仿真案例说明该方法有效可行.文献[21]基于时间步进法直接引入预碰撞速度二次函数判别碰撞前黏滞/滑动状态,并对Kane[22]案例进行研究证明其方法的正确性.以上文献对于接触碰撞过程通常采用Coulomb 摩擦模型进行描述,当对柔性梁接触碰撞的非光滑问题进行研究时,为了对柔性梁碰撞中黏滞时摩擦力的准确描述需要引入约束方程进行计算.而LuGre 摩擦模型可通过平均鬃毛变形来反应实时摩擦系数变化,其包含了运动过程中最大静摩擦、预滑动、黏滞等现象,直接避免黏滞/滑动切换的判断,减少建模过程中动力学方程的推导,且能完整描述预滑动阶段的摩擦特性,有学者[23-27]将LuGre 模型用于多体系统接触碰撞中,证明该模型能够解决接触碰撞过程中的非光滑问题.以上文献对于接触碰撞过程中法向接触力的计算通常采用连续接触法,该方法需要进行嵌入假设,而接触碰撞过程的实际情况应该是非嵌入的,此时更符合工程实际情况的附加约束法被国内外学者提出[28-32].段玥晨等[30]证明该方法能够方便准确求解出法向碰撞力并与连续接触法的作对比.Yuan 等[32]针对柔性梁和刚性斜面斜碰撞问题进行研究,在法向碰撞时引入约束方程求解碰撞力,在切向上根据滑动/黏滞判断条件分别采用切向接触约束法直接求解静摩擦力和Coulomb 模型计算滑动摩擦力的混合建模方法,并与采用修正库伦模型和无摩擦模型作对比证明该方法能够相对准确描述碰撞时的摩擦力.

如上述文献研究可知,接触约束法和LuGre 模型在碰撞问题研究中有各自的优势,然而目前柔性梁采用接触约束法结合LuGre 摩擦模型对斜碰撞问题研究的文献相对尚少.本文首先以刚柔耦合多体系统动力学理论为基础建立无碰撞时重力作用下旋转柔性梁的动力学方程;其次在碰撞时引入法向加速度约束计算出法向碰撞力建立无摩擦的碰撞动力学方程,在切向摩擦力计算上对于LuGre 模型提出两种不同处理方法,第一种是对LuGre 摩擦模型在黏滞时建立接触约束方程引入拉格朗日乘子计算出黏滞时切向摩擦力,滑动时根据摩擦系数和法向碰撞力计算出相应摩擦力(与Coulomb 摩擦模型处理切向摩擦力的方法一致).第二种是无论滑动还是黏滞状态直接使用LuGre 摩擦模型的摩擦力计算方法,碰撞全部过程采用相同的摩擦力计算方程.最后基于系统动力学方程编写含摩擦的碰撞动力学仿真程序.通过动力学仿真发现LuGre 模型在黏滞状态下两种不同建模方式的动力学特性无差别,Coulomb 摩擦模型建模方式的动力学特性随着碰撞次数增加与LuGre 摩擦模型的逐渐产生区别.

1 力学模型描述

如图1 所述,柔性梁绕铰O在重力场作用下作大范围定轴旋转运动,柔性梁上任意一点P与刚性坡面上点P′发生碰撞.以铰O为惯性坐标系原点,水平方向为惯性坐标系X轴,建立惯性坐标系OXY.在未变形时的柔性梁的中轴线为浮动坐标系的x轴,铰O为原点建立单位矢量浮动坐标系Oxy.惯性坐标系OX的轴和浮动坐标系Ox的夹角为柔性梁旋转运动时的转角θ.取刚性斜坡与水平面倾角斜率为α,设O1X1轴与斜坡重合建立碰撞面坐标系O1X1Y1.假设柔性梁为欧拉−伯努利梁,忽略梁的轴向变形,考虑横向弯曲变形引起的纵向缩短量(即非线性耦合变形量).设柔性梁的物理参数分别为:长度为L,横截面为S,密度为ρ,弹性模量为E.

图1 柔性梁碰撞模型图Fig.1 Schematic of impact system of flexible beam

柔性梁在浮动坐标上任意一点Q变形后位移矢量表示为

式(1)中是uy表示任意点Q横向变形位移、ux表示为横向变形引起的纵向缩短量,采用假设模态法将柔性变形位移离散表示为

式中,Φy(x)表示为柔性梁的横向振动模态函数的行阵,B(t)表示横向振动模态坐标列阵,H(x)表示耦合形函数矩阵.

2 碰撞动力学建模

2.1 无碰撞时刚柔耦合动力学方程

柔性梁上任意一点Q变形后位置矢量在惯性坐标下表示为

式中,W表示惯性坐标到浮动坐标系的方向余弦矩阵,ρ0表示任意点Q未变形前在浮动坐标系的位置.柔性动能系统的动能为

不考虑剪切和扭转,柔性梁系统的势能为

式中,VE和VG分别为柔性梁的弹性势能和重力势能,g为惯性系下的重力加速度矢量阵,h(x) 为任意点Q变形后在浮动坐标Oxy上的位置矢量阵.将式(8)和式(9)代入第二类拉格朗日方程

式中,F为非保守力对应的广义力,其中F=[τ,0]T,τ为转动铰中的驱动力矩,通过推导可得到系统的动力学方程

2.2 无摩擦系统碰撞动力学方程

梁上一点P与刚性斜坡上点P′会发生碰撞,碰撞面的单位法向量和切向量在惯性坐标系分别为n=[nx,ny]T和τ=[τx,τy]T,WP为惯性坐标系到刚性斜坡坐标系的转换矩阵.在发生碰撞时,梁上碰撞点P相对于斜面接触点P′的在斜坡法向上位移相同,因此法向位移约束方程为

2.3 碰撞的初始条件求解

在接触约束法中,初始接触瞬间速度是不连续的,在碰撞接触时直接添加约束会导致该时刻运动不协调从而导致违约计算困难的情况.为了实现运动的协调,在初始碰撞时刻需要实现无碰撞到碰撞的速度跳跃,得出碰撞时刻添加约束的速度.对于得出速度跳跃值,许多学者进行了研究.Dong 等[29]根据连续介质力学的间断面理论得到碰撞初始速度,段玥晨等[30]利用冲量−动量法求出速度的跳跃值.本文采用冲量动量法求出碰撞的初始速度

2.4 黏滞状态时碰撞动力学方程

柔性梁上碰撞点P在碰撞过程中还包含切向的摩擦力作用.切向运动可以分为黏滞状态和滑移状态,在碰撞接触过程中满足一定条件可以相互切换.在黏滞的状态下,因为切向没有发生运动,需要建立切向的约束方程求解摩擦力,在滑移状态下,需要根据摩擦系数来确定,通过引入摩擦力势能求解.

当系统处于黏滞状态时,柔性梁上碰撞点P和斜坡上接触点P′在斜坡上切向位移相同,由此切向的约束方程为

研究发现斜向碰撞过程和法向碰撞力和切向摩擦力变化同步[33],目前大多数学者认为切向碰撞力和法向摩擦力相互独立[10,34].可以将无摩擦动力学和切向接触约束方程结合,黏滞状态的动力学方程为

为了同时满足位移、速度和加速度的约束方程,分别对法向和切向接触约束方程和切向约束方程进行违约修正

修正后,黏滞状态下碰撞动力学方程为

由黏滞状态下碰撞动力学方程求出,即法向和切向的接触约束反力为

2.5 滑动状态时碰撞动力学方程

在碰撞中发生滑移时,碰撞点法向间距、速度、加速度始终为零,而切向产生碰撞点的位移、速度、加速度不为零,本文分别采用Coulomb 摩擦模型和LuGre 摩擦模型用来描述滑动过程中的实际摩擦系数.Coulomb 摩擦模型取动摩擦系数µd,LuGre 模型摩擦系数为

式中,z和表示梁上碰撞点P与斜坡接触点P′的平均鬃毛变形和变形速率;σ0,σ1,σ2分别表示鬃毛刚度系数、鬃毛阻尼系数、黏性摩擦系数;vs为stribeck速率;µd和µs分别表示动摩擦系数、静摩擦系数.

柔性梁上碰撞点摩擦力如下

2.6 碰撞时状态的切换准则

接触/分离切换的准则是通过与斜坡的法向间距矢量和法向相对速度,判断是否是碰撞接触状态需要同时满足:

(1)柔性梁与斜面接触;

(2)梁上碰撞点相对于斜坡的法向速度为负.

黏滞/滑移状态切换准则:在切向接触状态中存在黏滞、正滑移、逆滑移现象,在切向接触过程其滑动和黏滞的判断条件为

其中,Coulomb 模型中µd和µ相同,因为零时刻速度程序难以获得,当切向速度小于微小速度Dv时,本文认为发生黏滞.

基于接触约束法的3 种碰撞动力学建模方式如下:

(1)Coulomb 模型黏滞时考虑切向约束混合建模方式(即hybrid Coulomb model)和LuGre 模型黏滞时考虑切向约束混合建模方式(即hybrid LuGre model).

黏滞状态时动力学方程

滑动状态时动力学方程

(2)LuGre 模型黏滞时不考虑切向约束采用统一摩擦力公式建模方式(即LuGre model).

碰撞时动力学方程

3 动力学仿真

3.1 与现有文献对比

取柔性梁的参数为:L=1.0 m,S=3.14×10−4m2,ρ=27 667 kg/m3,E=68.952 GPa,I=7.85×10−9.假设柔性梁在初始状态处于静止和未变形状态;柔性梁在初始角度θ=π/4 处在重力场下作自由下落.旋转过程中柔性梁和原点位于(0.3,0)处倾角α 为π/10的刚性斜坡发生;取动摩擦系数µd为0.3,如图2 和图3 所示柔性梁的角位移和末端Y方向位移与文献[32]仿真结果对比相同.

图2 角位移Fig.2 Angular displacement

图3 梁末端末端Y 方向位移Fig.3 Y-direction displacement of beam tip

3.2 仿真结果分析

取柔性梁尺寸与材料参数如3.1 节所示;取违约数χ 和β 为10 和1000.在惯性坐标系(0.5,0) 处取刚性斜坡坐标原点,斜坡坡倾角为π/10;柔性梁在初始角度θ=π/4 处在重力场下作自由下落与斜坡发生碰撞.取动摩擦系数µd为0.4、最大静摩擦系数µs为0.5,通过对比取σ0为1.0×105N/m、σ1为1.0×102N·s/m;假设初始时刻柔性梁未发生变形,未碰撞时和碰撞结束后鬃毛平均变形z和平均变形速率为0.

图4 和图5 分别表示旋转柔性梁在碰撞过程中角位移和角速度变化.如图 4 和图 5 所示LuGre 摩擦模型在仿真时考虑约束(黏滞时考虑切向约束方程混合建模方式——hybrid LuGre model) 和不考虑约束(黏滞时不考虑切向约束方程建模方式——LuGre model) 的角位移和角速度相同;如图4 所示Coulomb 摩擦模型(即hybrid Coulomb model) 的角位移与LuGre 摩擦模型考虑约束和不考虑约束的变化趋势相同,但变化并不完全相同,在第一次碰撞时和碰撞后3 种摩擦模型基本相同,第二次碰撞后Coulomb 摩擦模型和LuGre摩擦模型(hybrid LuGre model 和LuGre model) 的角位移存在不同,其反弹最大角位移相对偏大.如图5 所示3 种摩擦模型在碰撞后角速度曲线变化趋势相同,在第一次碰撞后LuGre 模型考虑约束和不考虑约束两者没有区别,Coulomb 模型与LuGre 模型考虑约束和不考虑约束的角速度有细微差别;但是第二次以后碰撞Coulomb 模型的速度变化范围相对于LuGre 模型的区别较大.碰撞过程黏滞状态时,LuGre模型考虑切向约束和不考虑的两种建模方式对旋转柔性梁的角位移和角速度无影响.

图4 柔性梁角位移Fig.4 Angular displacement of flexible beam

图5 柔性梁角速度Fig.5 Angular velocity of flexible beam

图6 所示为柔性梁末端横向变形,由图6(a)可知LuGre 模型两种方式(即hybrid LuGre model 和LuGre model)与Coulomb 模型碰撞后的横向变形的曲线变化趋势相同.由图6(a) 和图6(b) 看出在第一次碰撞后3 种模型的末端变形曲线基本重合;在第二次碰撞后,Coulomb 模型相对LuGre 模型两种方式柔性梁末端变形的变化基本相同但是变形量相对LuGre 模型两种方式稍微偏大;第三次碰撞后Coulomb 模型变形量相对于LuGre 模型有明显不同,碰撞后末端横向变形量明显偏大.LuGre 模型两种方式在多次碰撞中末端变形仍然相同.图7 和图8 分别表示柔性梁末端Y方向的位移和速度,从图7 和图8 可以看出,3 种模型方式(即hybrid Coulomb model,hybrid LuGre model 和LuGre model) 柔性梁末端在Y方向位移和速度变化基本相同.在仿真时间内,LuGre 模型两种方式Y方向位移和速度变化曲线重合,Coulomb 模型随着碰撞次数增多和LuGre 模型两种方式Y方向位移和速度曲线差别先变大后变小.LuGre 模型考虑约束和不考虑约束柔性梁末端变形、位移、速度变化相同.

图6 柔性梁末端横向变形图Fig.6 Transverse deformation of flexible beam tip

图6 柔性梁末端横向变形图(续)Fig.6 Transverse deformation of flexible beam tip(continued)

图7 柔性梁末端Y 方向位移Fig.7 Y-direction displacement of flexible beam tip

图8 柔性梁末端Y 方向速度Fig.8 Y-direction velocity of flexible beam tip

图9 表示3 种摩擦模型方式柔性梁碰撞过程中法向碰撞力的变化图,如图9 所示,在第一次碰撞中法向的碰撞力完全相同,随着碰撞次数增加3 种模型碰撞时间和法向碰撞力大小出现明显区别.在仿真时间内,LuGre 模型黏滞时考虑约束和不考虑约束的碰撞力大小无明显区别,而Coulomb 模型随着碰撞次数增加碰撞力相对LuGre 模型两种方式差异逐渐明显.图10 表示3 种摩擦模型柔性梁碰撞过程中切向摩擦力的变化,由图10(b)放大图可知,切向摩擦力在第一次碰撞过程中变化曲线基本相同,在碰撞开始阶段LuGre 模型两种方式相对Coulomb 模型摩擦力相对略微偏大,随后3 种方式曲线重合.LuGre模型的两种方式在多次碰撞中摩擦力没有明显区别,而Coulomb 摩擦力在多次碰撞后和LuGre 模型出现明显差异.说明LuGre 模型黏滞状态考虑约束和不考虑约束方程两种碰撞建模方式不影响系统的碰撞力和摩擦力.

图9 碰撞力变化图Fig.9 Variation of impact force

图10 摩擦力变化图Fig.10 Variation of friction force

图11 表示柔性梁3 种摩擦模型方式的能量变化,从图11(a)可以看出3 种摩擦模型方式机械能整体变化的趋势相同,随着碰撞次数增加,每次碰撞过程损耗的机械能逐渐减少.第一次碰撞过程中LuGre模型的两种方式消耗的能量比Coulomb 模型略微多点但是不明显,由于第一次碰撞过程中3 种模型法向碰撞力相同,能量差异只能是LuGre 模型切向摩擦力和Coulomb 模型出现略微不同而引起的.随着碰撞次数增加,柔性梁系统能量的大小变化差异呈现先增大后减小趋势,但在仿真过程中LuGre 模型考虑黏滞约束和不考虑黏滞约束方程机械能无明显区别.从图11(b) 可以看出,在碰撞过程中会引起动能的突变,在第一碰撞后动能变化基本相同,但是后面几次碰撞后逐渐可以看出Coulomb 模型相对其他两种模型明显不同.由图11(c) 看出,3 种模型方式的柔性梁弹性势能曲线趋势相同,对比图6 可以看出在碰撞时柔性梁会出现大的变形随后引起高频振荡,其弹性势能体现为在碰撞时产生大的突变、随后数值发生大范围周期变化.随着能量的损耗,弹性势能数值变化幅度能随着碰撞次数也逐渐减小、梁末端变形幅度也随着减小.Coulomb 模型的弹性势能随着碰撞次数增加和其他两种方式逐渐产生区别后随着能量消耗区别减小,在仿真过程中LuGre 模型两种方式的弹性势能没有区别.由图12(d) 看出,重力势能变化趋势和角位移变化趋势基本相同,重力势能的变化主要受刚性的转动影响,在第二次碰撞后Coulomb 模型和LuGre 模型两种方式有明显区别,柔性梁碰撞后Coulomb 模型反弹高度相对较高,机械能在反弹重力势能最大时刻能量以重力势能为主,还有少部分机械能以弹性势能存在,重力势能碰撞最大反弹高度区别能够反应出相同条件下Coulomb能量消耗略低而LuGre 模型考虑约束和不考虑约束能量消耗相同.反弹高度在第二次碰撞反弹同样影响碰撞力区别,第三次碰撞由于前两次LuGre 模型两种方式消耗能量大,出现反弹高度小再次碰撞就出现碰撞力和摩擦力减小,LuGre 模型两种方式能量损耗减小,从而呈现3 种摩擦模型方式能量之间差异减小的现象.在仿真时间内LuGre 模型黏滞状态下考虑约束和不考虑约束方程建模的两种方法所有能量曲线重合,说明LuGre 模型黏滞状态下引入计算切向约束反力和LuGre 模型计算黏滞状态下摩擦力的能量损耗相同.

图11 柔性梁能量变化图Fig.11 Energy variation of flexible beam

图11 柔性梁能量变化图(续)Fig.11 Energy variation of flexible beam(continued)

图12 表示第一次碰撞梁上碰撞点相对于斜坡原点的切向位移,从图12 看出在第一次碰撞过程中3 种模型方式的位移基本相同,且在碰撞过程中都反映出正反滑移和黏滞现象.其中LuGre 无论是否考虑约束的滑动位移完全相同,说明LuGre 模型黏滞状态下不考虑约束同样能准确表现出该模型黏滞状态的位移.图13 表示梁上碰撞点碰撞时在斜坡上切向滑动速度,可以看出随着碰撞次数增加,Coulomb 摩擦模型的速度相对于LuGre 模型的两种方式逐渐产生区别.图14 表示LuGre 摩擦模型碰撞过程摩擦系数变化,从图14(a) 和图14(b) 摩擦系数变化可知LuGre 摩擦模型黏滞考虑约束和不考虑约束的建模的摩擦系数变化相同,由图14(b)第一次碰撞时摩擦系数变化可以看出,在碰撞时摩擦系数会出现高于和低于动摩擦系数的值,对比切向滑动位移和速度变化,得出在开始滑动的阶段摩擦系数出现超过滑动摩擦系数(0.4)的现象,这与实际情况相符,预滑动到宏观滑动时推力必须大于临界摩擦力.这个现象反映在摩擦力上为碰撞开始阶段,LuGre模型的两种方式的摩擦力相对于Coulomb 模型偏大.图14(b)和图12 结合同时看出,在切向位移发生黏滞时(即位移不变时),摩擦系数值出现小于或大于滑动摩擦系数(0.4)的现象,这说明黏滞时摩擦力是小于当时状态最大静摩擦力而不是小于滑动摩擦力,反应在摩擦力上在黏滞时出现LuGre 模型的两种方式会出现大于Coulomb 摩擦力.因为摩擦系数不同会导致在碰撞开始阶段和黏滞状态下不同模型摩擦力不同,因此不同模型能量损耗等动力学特性也不同.采用Coulomb 模型虽然第一次碰撞时滑动位移、速度等动力学特性与LuGre 摩擦模型基本相同,但是随着碰撞次数增加,期间能量损耗逐渐累积的差异逐渐体现出来,Coulomb 模型相对LuGre 模型的两种方式动力特性区别也逐渐明显,从而明显看出柔性梁末端变形差异、Coulomb 模型能量损耗小、再次碰撞时间间隔较长等现象,不过这种现象随着能量消耗趋于稳定值差异也逐渐减小.如图5 柔性梁角速度的变化上亦能看出Coulomb 模型和LuGre 模型的区别,第一次接触碰撞因切向摩擦力微小区别导致柔性梁角速度变化微小;在第二次碰撞时,Coulomb 模型的摩擦力相对于LuGre 模型两种方式的有明显区别,在碰撞反弹后柔性梁角速度产生差异;在第三次及接下来的碰撞中由于多次碰撞累积效应,Coulomb 模型碰撞力明显偏大,导致碰撞时其柔性梁末端变形变大,进而柔性梁的柔性作用反向影响柔性梁大范围旋转运动,体现在柔性梁角速度变化上.在切向摩擦力上考虑和未考虑黏滞的LuGre 摩擦建模时摩擦系数、摩擦力和切向位移、速度并未出现差异,说明LuGre摩擦模型能够反应出该模型的黏滞时摩擦力、位移和速度的变化.在柔性梁斜向碰撞动力学建模时,法向接触碰撞使用接触约束法,符合碰撞过程中非嵌入的实际情况;切向使用LuGre 摩擦模型引入摩擦势能,能够反应碰撞过程中切向运动的实际情况.这样斜向碰撞动力建模相对于采用Hybrid LuGre model建模减少了计算程序中黏滞滑动的切换,同样推导公式相对方便简洁.

图13 碰撞时切向速度变化Fig.13 Variation of tangential velocity during impact

图14 碰撞时摩擦系数变化Fig.14 Variation of friction coefficient during impact

4 结论

(1) 在碰撞接触过程中,3 种模型的动力学特性相似,但随着碰撞次数的增加Coulomb 摩擦模型和LuGre 摩擦模型两种方式动力学特性差异逐渐明显,后随碰撞时能量消耗柔性梁系统趋于稳定,差异又逐渐减小.

(2)LuGre 摩擦模型相对于Coulomb 摩擦模型对于碰撞过程中的斜向摩擦描述更加准确,可以描述预滑动到宏观滑动之间的临界摩擦力及在黏滞状态下介于正负最大静摩擦力之间的摩擦力变化.

(3)LuGre 摩擦模型在黏滞状态下考虑约束和不考虑约束在仿真时间内动力学特性没有区别,摩擦力、摩擦系数、碰撞时切向位移等均相同.这些现象说明在黏滞状态下LuGre 摩擦模型的动力学特性和切向约束描述的相同且减少了黏滞状态动力学方程的推导,避免了计算程序中黏滞/滑动状态判断.