高速曲柄压力机作为一种高效率、高精度的锻压设备，广泛应用于工业生产中^[1-2]。曲柄滑块机构是高速曲柄压力机主要的传动结构，用于传递运动和动力^[3], 其中曲轴的受力最为复杂，其承受着巨大的周期性冲击载荷和惯性力载荷，且支承颈和曲柄颈的轴肩处截面尺寸变化剧烈，导致截面处出现应力集中，而应力集中往往是零件损坏的根源^[4-5]。针对高速曲柄压力机的曲轴，衣锐对其进行了受力分析，根据相应载荷谱，用Workbench软件进行了静力学分析和疲劳分析^[6]。牛瑞霞等建立了曲轴的数学模型，进行了理论计算，并对其进行静态和模态分析，得到了曲轴的固有频率^[7]。张德峰等分析了曲轴在公称力冲击下的应力分布^[8]。王磊等建立了压力机转子-轴承系统的振动模型并采用Newmark方法进行求解，通过模拟仿真分析了转子-轴承系统的振动响应^[9]。

在目前对高速曲柄压力机曲轴的研究中，往往忽略了曲柄滑块机构高速运动时不平衡惯性力对曲轴的影响。经研究发现，随着曲柄转速的提高，惯性力所引起的附加转矩在一定速度下可超过由公称力产生的转矩^[10]。因此，在对曲轴进行分析时，忽略不平衡惯性力对曲轴的影响，会导致分析结果不符合实际。

本文基于联立约束法和达朗贝尔原理建立了曲柄滑块机构的动力学模型。首先采用达朗贝尔原理建立机构的动力学方程，然后，根据机构的位移闭环矢量方程，建立各构件的运动学方程及加速度方程，最后通过联立约束法联立所建立的运动学方程与动力学方程。该方法可以准确地反映出速度、加速度和惯性载荷之间的关系^[11]。将动力学方程组装成稀疏矩阵，并将其作为Function模块嵌入Simulink仿真模型中，通过模型对矩阵进行迭代求解。通过任意改变曲轴的转速和冲压力的大小，快速、准确地求出曲轴各个轴颈处载荷随时间的变化规律。根据仿真求得的载荷大小，用Workbench软件对压力机曲轴施加载荷，进行疲劳分析。由于考虑了惯性力对曲轴的影响，得到的分析结果更加接近实际，其为高速曲柄压力机曲柄和曲轴支撑结构的设计提供可靠的理论依据。

1 曲柄滑块机构动力学模型的建立

高速曲柄压力机高速运转时，曲柄滑块机构自身会产生周期性变化的不平衡惯性力，不平衡惯性力超过一定范围时会导致滑块下死点的动态性能急剧恶化，引起机床的强烈振动和摆动。因此为了保证冲压性能和机身的稳定性，部分高速压力机采用反向配置平衡副滑块的方法以抵消冲压时所产生的不平衡惯性力，该方法是一种较理想的不完全动平衡方法^[12-14]。

高速压力机中的曲柄滑块机构如图 1(a)所示放置，x方向为实际竖直方向, y方向为实际水平方向。其中连杆3和滑块4为压力机冲压部分，连杆5和副滑块6为平衡部分。曲轴的结构和其承受的载荷分布如图 1(b)所示：平衡曲柄颈驱动平衡副滑块，受到连杆5和滑块6产生的不平衡惯性力；压力曲柄颈驱动冲压滑块，承受冲压力和连杆3与滑块4产生的不平衡惯性力。

平衡部分和冲压部分均为曲柄滑块机构，由同一曲轴驱动。平衡曲柄颈和压力曲柄颈之间存在支承颈，忽略运动副之间的间隙，根据刚化原理，分析曲柄颈时，将支承颈视为刚体，其平衡状态不变。因此只需对一侧曲柄滑块机构进行受力分析，然后在相同的仿真条件下，通过改变曲柄连杆的参数分别求得两侧曲柄滑块机构动力学仿真结果。

图 2为各连杆和滑块的受力简图。

根据图 2所示的受力图，基于达朗贝尔原理建立机构动力学方程：

$ {F_{ji,x}} + {F_{ki,x}} = {m_i}{a_{{\rm{c}}i,x}} $

(1)

$ {F_{ji,y}} + {F_{ki,y}} = {m_i}{a_{{\rm{c}}i,y}} $

(2)

$ \begin{array}{*{20}{c}} { - {F_{ki,x}}\left( {{r_i} - {r_{{\rm{c}}i}}} \right)\sin {\theta _i} + {F_{ki,y}}\left( {{r_i} - {r_{{\rm{c}}i}}} \right)\cos {\theta _i} + }\\ {{F_{ij,x}}{r_{{\rm{c}}i}}\sin {\theta _i} - {F_{ij,y}}{r_{{\rm{c}}i}}\cos {\theta _i} = {I_i}{\alpha _i}} \end{array} $

(3)

$ {F_{34,x}} + F = {m_4}{{\ddot r}_1} $

(4)

$ {F_{43,y}} + {F_{14,y}} = 0 $

(5)

式中：F_{ji, x}，F_{ki, x}为x方向上连杆j和k对连杆i的作用力；F_{ji, y}，F_{ki, y}为y方向上连杆j和k对连杆i的作用力；a_{ci, x}，a_{ci, y}为连杆i在x和y方向上的加速度；r_i为连杆i的长度; r_ci为第i个连杆质心到连杆旋转中心的距离；${\ddot r}$为第i个连杆的质心加速度；m_i为第i个连杆的质量；α_i为第i个连杆的角加速度；I_i为第i个连杆的转动惯量。

为了建立连杆的力和加速度关系，对闭环矢量方程进行二阶求导并建立质心加速度方程，即方程(6)至(11)。然后将运动学方程与动力学方程组装成稀疏矩阵，编入MATLAB的程序文件中，稀疏矩阵如(12)所示。

$ {{\ddot r}_1} + {r_3}{\alpha _3}\sin {\theta _3} + {r_2}{\alpha _2}\sin {\theta _2} = - {r_2}\omega _2^2\cos {\theta _2} - {r_3}\omega _3^2\cos {\theta _3} $

(6)

$ - {r_3}{\alpha _3}\cos {\theta _3} + {r_2}{\alpha _2}\cos {\theta _2} = - {r_2}\omega _2^2\sin {\theta _2} - {r_3}\omega _3^2\sin {\theta _3} $

(7)

$ {a_{{\rm{c}}2,x}} = - {r_{{\rm{c}}2}}{\alpha _2}\sin {\theta _2} - {r_{{\rm{c}}2}}\omega _2^2\cos {\theta _2} $

(8)

$ {a_{{\rm{c}}2,y}} = {r_{{\rm{c}}2}}{\alpha _2}{\rm{con}}{\theta _2} - {r_{{\rm{c}}2}}\omega _2^2\sin {\theta _2} $

(9)

$ \begin{array}{l} {a_{{\rm{c}}3,x}} = - {r_2}{a_2}\sin {\theta _2} - {r_2}\omega _2^2\cos {\theta _2} - \\ \;\;\;\;\;\;\;\;\;{r_{{\rm{c3}}}}{\alpha _3}\sin {\theta _3} - {r_{{\rm{c3}}}}\omega _3^2\cos {\theta _3} \end{array} $

(10)

$ \begin{array}{l} {a_{{\rm{c}}3,y}} = - {r_2}{a_2}{\rm{con}}{\theta _2} - {r_2}\omega _2^2\sin {\theta _2} + \\ \;\;\;\;\;\;\;\;\;{r_{{\rm{c3}}}}{\alpha _3}{\rm{con}}{\theta _3} - {r_{{\rm{c3}}}}\omega _3^2\sin {\theta _3} \end{array} $

(11)

$ \left[ {\begin{array}{*{20}{c}} 1&0&1&0&0&0&0&0&0&0&{ - {m_2}}&0&0&0\\ 0&1&0&1&0&0&0&0&0&0&0&{ - {m_2}}&0&0\\ 0&0&{ - {r_2}{S_2}}&{{r_2}{C_2}}&0&0&0&1&0&0&0&0&0&0\\ 0&0&{ - 1}&0&1&0&0&0&0&0&0&0&{ - {m_3}}&0\\ 0&0&0&{ - 1}&0&1&0&0&0&0&0&0&0&{ - {m_3}}\\ 0&0&{ - {r_{{\rm{c}}3}}{C_3}}&{{r_3}{S_3}}&{\left( {{r_{{\rm{c}}3}} - {r_3}} \right){S_3}}&{\left( {{r_3} - {r_{{\rm{c}}3}}} \right){C_3}}&0&0&0&{ - {I_3}}&0&0&0&0\\ 0&0&0&0&1&0&0&0&{ - {m_4}}&0&0&0&0&0\\ 0&0&0&0&0&1&1&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&1&{{r_3}{S_3}}&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&{ - {r_3}{C_3}}&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0&1&0&0&0\\ 0&0&0&0&0&0&0&0&0&0&0&1&0&0\\ 0&0&0&0&0&0&0&0&0&{{r_{{\rm{c3}}}}{S_3}}&0&0&1&0\\ 0&0&0&0&0&0&0&0&0&{ - {r_{{\rm{c3}}}}{C_3}}&0&0&0&1 \end{array}} \right] \\ \left[ {\begin{array}{*{20}{c}} {{F_{12,x}}}\\ {{F_{12,y}}}\\ {{F_{32,x}}}\\ {{F_{32,y}}}\\ {{F_{43,x}}}\\ {{F_{43,y}}}\\ {{F_{14,y}}}\\ {{M_{12}}}\\ {\ddot r}\\ {{\alpha _3}}\\ {{a_{{\rm{c2}},x}}}\\ {{a_{{\rm{c2}},y}}}\\ {{a_{{\rm{c3}},x}}}\\ {{a_{{\rm{c3}},y}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ F\\ 0\\ { - {r_2}{C_2}\omega _2^2 - {r_3}{C_3}\omega _3^2}\\ { - {r_2}{S_2}\omega _2^2 - {r_3}{S_3}\omega _3^2}\\ { - {r_{{\rm{c2}}}}{C_2}\omega _2^2}\\ { - {r_{{\rm{c2}}}}{S_2}\omega _2^2}\\ { - {r_2}{C_2}\omega _2^2 - {r_{{\rm{c3}}}}{C_3}\omega _3^2}\\ { - {r_2}{S_2}\omega _2^2 - {r_{{\rm{c3}}}}{S_3}\omega _3^2} \end{array}} \right] $

(12)

式中:ω₂，ω₃为连杆2和连杆3的角速度；C₂=cos θ₂，S₂=sin θ₂，C₃=cos θ₃，S₃=sin θ₃。

2 曲柄滑块机构仿真模型及仿真结果 2.1 仿真初始条件

某型号压力机曲柄滑块构件参数如表 1所示。

压力机的公称力F=800 kN，滑块行程为30 mm，行程次数为1 250次/min。仿真起始位置为压力机的下死点位置，即θ₂=θ₃=0。曲柄滑块共线时的速度方程为：

$ \left\{ \begin{array}{l} - {r_2}{\omega _2}\sin {\theta _2} - {r_3}{\omega _3}\sin {\theta _3} = {{\dot r}_1}\\ {r_2}{\omega _2}\cos {\theta _2} + {r_3}{\omega _3}\cos {\theta _3} = 0 \end{array} \right. $

(13)

根据式(13)，求出行程次数为1 250次/min时，ω₃=-1.03 rad/s。仿真的初始条件为：ω₂=130.83 rad/s，ω₃=-1.03 rad/s，θ₂=θ₃=0。

2.2 Simulink仿真模型

曲柄滑块机构Simulink仿真模型如图 3所示。选择ode45变步长求解器，相对误差设为10^-3。

2.3 仿真结果

根据建立的Simulink仿真模型进行动力学仿真，取仿真时间为0.2 s，冲压力为800 kN。得到压力曲柄颈处载荷随时间的变化规律，如图 4所示。

图 4(a)为压力曲柄颈处载荷大小的变化。1个行程周期为0.048 s。其中0-0.04 s为空载阶段，压力曲柄受到连杆和滑块的不平衡惯性力，最大为25 kN; 0.04-0.048 s为冲压阶段，受到连杆和滑块的不平衡惯性力和冲压力, 最大为825 kN。图 4(b)为以x轴正方向为基准的压力曲柄颈处载荷方向的变化。

通过改变曲柄连杆参数及冲压力，在同样的仿真条件下得到平衡曲柄颈处载荷随时间的变化规律, 如图(5)。

图 5(a)为平衡曲柄处载荷大小的变化。1个行程周期为0.048 s，在第0秒和第0.048秒时, 滑块处于右极限位置，第0.024秒时滑块处于左极限位置。由图可以看出滑块处于极限位置时，平衡曲柄颈处的载荷最大，滑块位于行程中间位置时，平衡曲柄颈处的载荷最小。图 5(b)为以x轴正方向为基准的平衡曲柄颈处载荷方向的变化。

由于压力机反置了平衡副滑块，抵消了机构绝大部分的不平衡惯性力，所以支承颈处的不平衡惯性力相比于冲压力可以忽略不计。根据刚化原理，视曲柄颈为刚体, 在同样的仿真条件下得到支承颈处载荷随时间的变化规律, 如图 6所示。

由图 6(a)可以看出，曲轴支承颈处的惯性力载荷经平衡副滑块抵消后，最大的不平衡惯性力为6 kN，远远小于冲压力。所以曲轴支承颈主要承受冲压力的作用，承受的最大载荷为806 kN。

4 曲轴疲劳分析

根据现有研究得知，利用有限元分析法分析曲轴时，按照均布载荷形式对其加载, 得到的分析结果才更可靠, 更接近曲轴变形的真实情况^[15]。曲轴滑块的受力情况如图 1(a)所示，根据Simulink仿真求得的曲轴支承颈处和曲柄颈处载荷谱，采用Workbench软件对其加载，并进行静态分析和疲劳寿命分析^[16-19]。

图 7为曲轴的疲劳寿命分析结果。图 7(a)为曲轴的应力云图，最大应力为34.8 MPa，位于压力曲柄颈和支承颈的轴肩处，这是因为轴肩处截面尺寸发生突变，导致应力集中严重。图 7(b)为曲轴的疲劳寿命云图，最小疲劳寿命达到1.36×10⁹次，根据工程设计中的规定：对于应力疲劳曲线有水平渐近线的材料，若经历10⁷次应力循环后不损坏，即认为其可承受无穷多次应力循环，因此认为该曲轴已达到无限寿命。图 7(c)为曲轴的疲劳损伤云图，可知曲轴的损伤最大值为0.73，损伤数值小于1，即达到无限寿命时没有出现明显损伤。图 7(d)为曲轴的疲劳安全系数云图，曲轴的最低安全系数为1.09，安全系数大于1，因此，曲轴的安全系数满足设计要求。

5 结论

1) 采用联立约束法建立了曲柄滑块机构动力学模型，并用MATLAB中的Simulink模块进行了仿真分析, 快速、准确地求出了曲轴各个轴颈处载荷随时间的变化规律。

2) 仿真得到的各曲柄颈和支承颈处载荷的变化规律是根据运动学参数求得的，可以真实地反映出速度、加速度和惯性载荷之间的关系，为压力机曲轴、曲轴支撑结构和动平衡结构设计及压力机下死点精度的研究提供重要的理论依据。

3) 根据Simulink仿真求得的曲轴支承颈处和曲柄颈处载荷谱，采用Workbench软件对其进行加载，并进行静态分析和疲劳寿命分析。疲劳分析结果验证了曲轴的疲劳寿命符合设计要求，分析结果更加接近实际，为压力机曲轴的实际应用提供参考。