在多体系统中, 由于制造误差、材料变形、摩擦磨损等因素的影响, 间隙普遍存在于机构的运动副中。而且, 间隙处产生的碰撞力会增大部件间的磨损率及系统的振动噪声水平^[1]。另外, 现代的机械产品都在朝着高速高精度的方向发展, 关节间隙、磨损、部件柔性和润滑等因素严重影响着精密机构的动力学性能, 因此, 为精确预测机械系统的动力学行为, 动力学建模中尽可能多地考虑像部件柔性、多个关节间隙、磨损及润滑等的影响因素是很有必要的^[2]。

目前, 国内外学者对考虑间隙、部件柔性、磨损和润滑机构的动力学建模及分析进行了大量的研究。基于连续碰撞力模型^[3], Flores等^[4-6]提出了一种考虑摩擦和润滑效应的平面及空间含间隙机构动力学建模方法, 理论和实验研究了不同曲柄转速、间隙大小对机构动力学性能的影响, 结果表明润滑铰情况下的动力学响应曲线和理想铰十分接近。在此基础上, Tian等^[7-8]对含间隙润滑铰的平面及空间柔性多体系统分别进行了动力学分析。Khemili等^[9]研究了考虑部件柔性的平面多体系统动力学行为, 数值仿真结果显示部件柔性能够削弱碰撞力峰值。为提高动力学分析的精度, Flores^[10]对考虑多个间隙关节的刚性曲柄滑块机构进行了研究, 分析了间隙关节数量对机构动力学行为的影响, 发现两个间隙关节会明显增加系统的响应峰值。近来, Zheng等^[11]利用动力学分析软件ADAMS对一种封闭式高速压力系统进行了分析, 分析过程中作者将系统简化为含多个间隙的刚柔耦合曲柄滑块机构。蒋万松等^[12]针对翼伞-载荷多体系统进行了动力学分析。然而, 需要指出的是, 上述研究基本都是以曲柄滑块机构为模型进行的分析, 而对四连杆机构的研究较少, 白争锋等^[1]将一种改进的连续接触力混合模型嵌入到ADAMS中, 对含间隙平面四连杆机构的动力学行为进行了仿真分析。Flores^[13]在他的博士论文中则对空间含间隙四连杆机构进行了动力学分析。Erkaya等^[14]通过将关节间隙等效为一根无质量连杆对含两个间隙的平面四连杆机构进行了运动学和动力学分析。

然而, 对比上述含间隙机构的研究发现, 考察部件柔性或润滑对动力学响应的影响时, 动力学建模中通常只考虑了一个间隙关节, 而当动力学建模考虑多个间隙关节时, 部件柔性或润滑的影响又没有被考虑, 或是在动力学分析软件ADAMS中进行的, 这是因为每增加一个影响因素都会使多体系统的建模更加复杂, 然而可以提高动力学分析的精度。

鉴于此, 本文在上述研究的基础上, 以四连杆机构为研究对象, 建立了考虑关节润滑和两个间隙的动力学数值仿真模型, 对单个间隙关节的位置、两个间隙关节及润滑条件对曲柄驱动力矩和随动杆加速度的影响进行了分析和对比, 得出了一些有意义的结论, 本文的研究为以四连杆机构为简化模型的机械系统分析提供了理论基础。

1 多体系统动力学建模

一个受完整约束的刚性多体系统运动方程本质上是一组微分代数方程组(DAEs), 基于拉格朗日方程, 该DAEs具有如下的表达式^[4]

$ \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{M}}&{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_q^{\rm{T}}}\\ {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_q}}&{\bf{0}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\ddot q}}}\\ \mathit{\boldsymbol{\lambda }} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{Q}}\\ \mathit{\boldsymbol{\gamma }} \end{array}} \right] $

(1)

式中:M为整体的质量矩阵, 由系统部件的质量和转动惯量组成; Φ_q为约束方程的Jacobi矩阵; $\mathit{\boldsymbol{\ddot{q}}}$为广义坐标的加速度矢量; λ为和约束有关的拉格朗日乘子向量; Q为广义力向量; γ为和广义坐标的速度有关的向量。其表达式如下

$ \mathit{\boldsymbol{\gamma }} = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_q}\mathit{\boldsymbol{\ddot q}} = - {\left( {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_q}\mathit{\boldsymbol{\dot q}}} \right)_q} - 2{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{qt}}\mathit{\boldsymbol{\dot q}} - {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{tt}} $

(2)

式中:Φ_qt为约束方程关于坐标向量q和时间t的偏导数, 而Φ_tt为约束方程关于时间t的二阶偏导数。

为避免式(1) 在数值求解过程中出现约束违约, 引入Baumgarte^[15]稳定技术来抑制求解过程中误差的增长, 此时式(1) 被改写为

$ \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{M}}&{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_q^{\rm{T}}}\\ {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_q}}&{\bf{0}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\ddot q}}}\\ \mathit{\boldsymbol{\lambda }} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{Q}}\\ {\mathit{\boldsymbol{\gamma }} - 2\alpha \mathit{\boldsymbol{ \boldsymbol{\dot \varPhi} }} - {\beta ^2}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}} \end{array}} \right] $

(3)

式中:α和β为稳定系数。一般情况下, α和β均为正数时, 式(3) 的稳定性才可以保证, 而且, α=β能更快地稳定系统响应。

2 间隙铰接副建模 2.1 旋转间隙关节数学模型

轴和轴承构成的旋转间隙关节如图 1所示, 轴和轴承分别固联于体i和体j上, 局部坐标系x_io_iy_i和x_jo_jy_j分别位于体i和体j的质心处, e_ij为连接轴和轴承中心的中心距矢量, 表达式如下^[16]

$ {\mathit{\boldsymbol{e}}_{ij}} = \mathit{\boldsymbol{r}}_i^p - \mathit{\boldsymbol{r}}_j^p = \left( {{\mathit{\boldsymbol{r}}_i} + {\mathit{\boldsymbol{A}}_i}{{\mathit{\boldsymbol{s'}}}_i}} \right) - \left( {{\mathit{\boldsymbol{r}}_j} + {\mathit{\boldsymbol{A}}_j}{{\mathit{\boldsymbol{s'}}}_j}} \right) $

(4)

$ \mathit{\boldsymbol{n}} = \frac{{{\mathit{\boldsymbol{e}}_{ij}}}}{{\left\| {{\mathit{\boldsymbol{e}}_{ij}}} \right\|}} $

(5)

式中:n为中心距矢量的单位向量, 决定间隙关节处法向碰撞力的方向; r_i^p和r_j^p分别为轴和轴承中心在全局坐标系下的位置矢量; r_i和r_j分别为全局坐标到体i和体j的局部坐标的位置矢量; s′_i和s′_j分别为轴和轴承中心在局部坐标系下的位置矢量; 而A_i, A_j分别表示体i和体j的局部坐标向全局坐标的转换矩阵。

图 1表示的轴和轴承的碰撞点Q_m(m=i, j)在全局坐标系下的位置矢量表示为

$ \mathit{\boldsymbol{r}}_m^Q = {\mathit{\boldsymbol{r}}_m} + {\mathit{\boldsymbol{A}}_m}{{\mathit{\boldsymbol{s'}}}_m} + {\mathit{\boldsymbol{R}}_m}\mathit{\boldsymbol{n}} $

(6)

式中R_m(m=i, j)表示轴和轴承的半径。当轴和轴承发生碰撞时, 由于局部变形, 轴嵌入轴承的深度为

$ \delta = {e_{ij}} - c $

(7)

式中:δ表示碰撞点Q_i和Q_j之间的距离, e_ij=‖e_ij‖为中心距矢量的幅值, c=R_j-R_i为径向间隙。

2.2 接触碰撞力模型

间隙旋转铰碰撞力模型的建立是含间隙多体系统动力学建模的重要内容。Lankarani-Nikravesh模型由于考虑了碰撞过程中的能量损失而被广泛使用, 其表达式如下^[16]

$ {F_N} = K{\delta ^n} + D\dot \delta $

(8)

式中:Kδⁿ项表示关节部件间的弹性碰撞力; D$\dot{\delta }$项表示碰撞过程中的能量耗散; K, D, δ和$\dot{\delta }$分别表示刚度系数, 阻尼系数, 相对渗透深度和相对渗透速度。指数n通常取1.5。刚度系数K由碰撞体的几何和材料属性决定, 表达式如下

$ K = \frac{4}{{3\left( {{h_i} + {h_j}} \right)}}{\left( {\frac{{{R_i}{R_j}}}{{{R_i} - {R_j}}}} \right)^{\frac{1}{2}}} $

(9)

式中:h_i和h_j为碰撞体的材料属性, 他们依赖于各自的弹性模量和泊松比, 表达式为

$ {h_m} = \frac{{1 - v_m^2}}{{{E_m}}}\;\;\;m = i,j $

(10)

而阻尼系数D为

$ D = \mu {\delta ^n} $

(11)

式中:μ表示黏性阻尼因子, 依赖于刚度系数K、初始渗透速度$\dot{\delta }$^(-)和恢复系数c_e, 其表达式为

$ \mu = \frac{{3K\left( {1 - c_e^2} \right)}}{{4{{\dot \delta }^{\left( - \right)}}}} $

(12)

将式(9~12) 代入式(8), 得到含黏性阻尼项的碰撞力模型为

$ {F_N} = K{\delta ^n}\left( {1 + \frac{{3\left( {1 - c_e^2} \right)}}{4}\frac{{\dot \delta }}{{{{\dot \delta }^{\left( - \right)}}}}} \right) $

(13)

此时的法向碰撞力为

$ {\mathit{\boldsymbol{f}}_N} = {F_N}\mathit{\boldsymbol{n}} $

(14)

2.3 摩擦力模型

当轴和轴承干碰撞(无润滑)存在相对切向运动时, 接触表面除了法向碰撞力外还有摩擦力的作用, Coulomb摩擦力模型常用于计算接触碰撞或冲击过程中的摩擦行为, 但对低速时出现的粘质和微滑现象, 该模型计算的摩擦力容易发生突变, 导致计算不收敛。而改进的Coulomb摩擦力模型则具有较好的数值计算特性, 其表达式如下^[13]

$ {\mathit{\boldsymbol{f}}_T} = - {c_f}{c_d}{F_N}\frac{{{\mathit{\boldsymbol{v}}_t}}}{{\left\| {{\mathit{\boldsymbol{v}}_t}} \right\|}} $

(15)

式中:${{\mathit{\boldsymbol{v}}}_{t}}={{({{\mathit{\boldsymbol{\dot{r}}}}_{j}}^{Q}-{{\mathit{\boldsymbol{\dot{r}}}}_{i}}^{Q})}^{\rm{T}}}\mathit{\boldsymbol{t}}$为轴与轴承之间的相对切向速度, 单位切向量t为单位法向量n逆时针旋转90°得到的, c_f和c_d分别为摩擦因数和动力学修正系数, c_d按如下表达式进行修正

$ {c_d} = \left\{ {\begin{array}{*{20}{c}} 0&{\left\| {{\mathit{\boldsymbol{v}}_t}} \right\| \le {v_0}}\\ {\frac{{\left\| {{\mathit{\boldsymbol{v}}_t}} \right\| - {v_0}}}{{{v_1} - {v_0}}}}&{{v_0} \le \left\| {{\mathit{\boldsymbol{v}}_t}} \right\| \le {v_1}}\\ 1&{\left\| {{\mathit{\boldsymbol{v}}_t}} \right\| \ge {v_1}} \end{array}} \right. $

(16)

式中v₀和v₁为给定的相对切向速度的容差。

此时, 将式(14) 计算的法向碰撞力和式(15) 计算的摩擦力等效到体i和体j的质心O_i和O_j, 如图 2所示, 则作用在体i质心上的力和力矩分别为

$ {\mathit{\boldsymbol{f}}_i} = {\mathit{\boldsymbol{f}}_N} + {\mathit{\boldsymbol{f}}_T} $

(17)

$ {m_i} = - \left( {y_i^Q - {y_i}} \right)f_i^x + \left( {x_i^Q - {x_i}} \right)f_i^y $

(18)

作用在体j质心上的力和力矩分别为

$ {\mathit{\boldsymbol{f}}_j} = - {\mathit{\boldsymbol{f}}_i} $

(19)

$ {m_j} = - \left( {y_j^Q - {y_j}} \right)f_j^x + \left( {x_j^Q - {x_j}} \right)f_j^y $

(20)

将式(17~20) 计算的力和力矩作为广义外力组装到动力学方程式(3) 中, 运用四阶Runge-Kutta法就可以求解含间隙多体系统的动力学响应。

3 润滑模型

在大多数机械系统中, 为降低摩擦磨损, 提高系统性能, 使用润滑剂是行之有效的方法之一。因此, 间隙旋转副中润滑模型的建立极其重要。通常, 润滑油膜效应分为纯挤压油膜效应和楔形油膜效应, 纯挤压油膜效应用于轴颈和轴承间的相对径向速度远大于相对旋转速度的情况, 而楔形油膜效应的应用相反^[5]。本文假设润滑间隙旋转副的油膜效应为纯挤压油膜效应。图 3显示了一个润滑旋转关节的横截面图。F_r和F_t分别为润滑力沿径向和切向的分量, 其中径向为轴和轴承的偏心距方向, 切向为径向沿逆时针转动90°得到。γ为轴与轴承的连心线和X轴方向的夹角。

在Sommerfeld润滑条件下, 无限长轴颈-轴承模型的润滑力可以表示为^[5]

$ {F_r} = - \frac{{12{\rm{ \mathsf{ π} }}\mu LR_j^3\dot \varepsilon }}{{{c^2}{{\left( {1 - {\varepsilon ^2}} \right)}^{3/2}}}} $

(21)

$ {F_t} = \frac{{12{\rm{ \mathsf{ π} }}\mu LR_j^3\varepsilon \left( {\omega - 2\dot \gamma } \right)}}{{{c^2}\left( {2 + \dot \varepsilon } \right){{\left( {1 - {\varepsilon ^2}} \right)}^{1/2}}}} $

(22)

式中:L为轴-轴承长度; ω为轴和轴承之间的相对角速度; μ为润滑油的黏度; ε=e/c为偏心率, 即轴和轴承的偏心距与间隙的比值; $\dot{\gamma }=({{e}_{x}}{{\dot{e}}_{y}}-{{\dot{e}}_{x}}{{e}_{y}})/{{e}^{2}}$, 而e_x, e_y为式(4) 中偏心距矢量的分量。

由式(21) 可以看出, 当偏心率ε接近1, 即轴和轴承开始接触时, 润滑力将趋于无限大, 此时的润滑力计算将不收敛。为此, Flores^[5]提出了转换模态的概念, 即存在一定厚度不可压缩的润滑油膜, 当润滑铰处的油膜厚度小于该指定厚度时, 轴和轴承的碰撞为干碰撞, 而当油膜厚度大于该指定厚度时, 润滑力模型起作用, 此时, 含润滑旋转铰的润滑力表示为

$ {F_r} = \left\{ {\begin{array}{*{20}{c}} {{F_{{\rm{squeeze}}}}}&{{e_{ij}} \le c}\\ {\frac{{\left( {c + {e_0}} \right) - {e_{ij}}}}{{{e_0}}}{F_{{\rm{squeeze}}}} + \frac{{{e_{ij}} - c}}{{{e_0}}}{F_{{\rm{dry}}}}}&{c \le {e_{ij}} \le c + {e_0}}\\ {{F_{{\rm{dry}}}}}&{{e_{ij}} \ge c + {e_0}} \end{array}} \right. $

(23)

式中:e₀为指定的偏心距容差, F_dry为轴-轴承干碰撞时计算的碰撞力, Fsqueeze为润滑力。该模型保证了润滑力仿真求解过程中的连续性。

当轴和轴承之间含润滑剂时, 将式(23) 计算的润滑力向X轴和Y轴投影, 然后等效到体i和体j的质心O_i和O_j处, 将计算的力和力矩作为广义外力整合到动力学方程式(3) 中, 运用四阶Runge-Kutta法即可求解含润滑旋转铰多体系统的动力学响应。

4 算例

本文以含间隙铰接副的四连杆机构为研究对象, 分析不同位置处的单个间隙关节、2个间隙关节及润滑关节对机构动力学响应的影响。含间隙四连杆机构简图如图 4所示, 由机架1, 连架杆2, 连杆3和连架杆4构成。间隙关节1位于连架杆2和连杆3之间, 间隙关节2位于连杆3和连架杆4之间。本文从以下几个方面对该四连杆机构进行了动力学研究:(1) 假设旋转铰不含间隙, 即理想情况, (2) 假设仅间隙关节1存在, (3) 假设仅间隙关节2存在, (4) 假设间隙关节1和间隙关节2都存在, (5) 在各间隙关节处存在润滑。机构的结构参数和仿真参数如表 1和表 2所示^[13]。

由表 1中杆的长度可以看出, 该机构为曲柄摇杆机构^[17], 图 3中的随动杆为摇杆。动力学仿真过程中假设曲柄以20π rad/s的恒定角速度转动, 稳定系数α和β的值为250, 恢复系数c_e为0.9, 干碰撞时的摩擦因数为0.001, 在40 ℃时润滑油的黏度为400 cP。上述数据均参考文献[12]。初始时, 杆2和杆3共线, 杆4与水平方向的夹角为锐角。图 5~10记录了系统达到稳定状态时2个运动周期内的动力学响应曲线。

对比图 5~7的摇杆角加速度响应曲线可以看出, 间隙存在时, 响应曲线有很多振荡峰值, 且间隙关节2存在时的振动峰值比间隙关节1存在时的振荡峰值要大, 两个间隙关节(间隙关节12) 同时存在时的振荡峰值最大, 而在润滑条件下, 摇杆角加速度的振荡峰值明显减弱, 与理想情况下的响应曲线基本一致, 只是在间隙关节2存在时有一些小的波动峰值(图 6, 7黄色实线所示)。产生上述这些现象的原因是, 位于连杆和摇杆之间的间隙关节2使得摇杆在运动过程中受到脉冲形式的碰撞力, 加大摇杆的振动与冲击, 从而影响其角加速度。而间隙关节1远离摇杆, 使得间隙关节1处产生的碰撞力在经连杆传递过程中逐渐被削弱, 从而对摇杆的动力学响应影响相对较弱。间隙关节数量增加时, 关节之间的碰撞频率和碰撞强度加大, 从而增大了动力学响应的振动峰值。

对比图 8~10的曲柄驱动力矩响应曲线可以看出, 间隙关节1存在时曲柄驱动力矩响应曲线的振动峰值比间隙关节2存在时要大, 两个间隙关节都存在时的振动峰值最大, 而考虑润滑时的振动峰值明显减弱, 与理想曲线基本一致。间隙关节1比间隙关节2对曲柄驱动力矩的影响要大是因为间隙关节1位于曲柄和连杆之间, 该间隙产生的冲击震荡对曲柄产生直接的影响, 而远离曲柄的间隙关节2产生的脉冲力在传递过程中逐渐被削弱, 从而对曲柄驱动力矩产生较小的影响。

5 结束语

为提高多体系统动力学分析的精度及运动可靠性, 本文以机械工业中常见的四连杆机构为研究对象, 建立了考虑多个间隙关节和润滑的多体系统动力学数值模型, 分析了间隙关节1、间隙关节2、间隙关节12及润滑关节对曲柄驱动力矩及摇杆角加速的影响, 发现间隙的存在使间隙关节处产生的碰撞力会使系统产生高频振荡, 引起各动力学响应曲线的波动峰值, 但间隙关节的位置和数量及润滑会使波动峰值不同, 其中, 摇杆角加速度的振荡峰值在间隙关节2存在比间隙关节1存在时要大, 曲柄驱动力矩的振荡峰值在间隙关节2存在比间隙关节1存在时要小, 然而, 无论哪种响应输出, 当间隙关节12存在时, 其波动峰值都有很大的增加, 而润滑能削弱每种情况下的响应输出, 起到很好的缓冲效果。