应用有限单元法的机床进给系统轴扭特性分析

陈 晔，熊晓航，赵春雨，张义民

（1.辽宁工业大学机械工程与自动化学院，辽宁 锦州 121001；2.东北大学机械工程与自动化学院，辽宁 沈阳 110819；3.杭州电子科技大学机械工程学院，浙江 杭州 310018）

1 引言

在机床的滚珠丝杠进给系统中，工作台和螺母与丝杠相连接，在受到工作载荷的作用时，丝杠和螺母会产生轴向和扭转的接触变形，从而导致丝杠产生轴向振动和扭转振动，进而使工作台产生轴向位移，影响进给系统的定位精度[1-2]。因此，研究机床传动系统的动力学是提高工作台定位精度的基础，国内外许多学者都对此进行了大量的研究工作。文献[3-5]对丝杠螺母结合面分别给出了有限元模型，并对动力学特性进行了分析；文献[6-9]对进给系统应用集中质量法建模，并对系统的振动响应进行了分析；文献[10]提出一种基于时间序列的进给系统建模方法，并用于加速度响应特性的预测；文献[11-12]分别对传动丝杠建立了振型方程，并对相应振型进行了分析，文献[13]考虑了摩擦和间隙因素对滚珠丝杠进给系统建模与分析。工作台与丝杠通过螺母接触相连，国内外对其接触特性与刚度进行了很多研究[14-20]，如应用Frenet-Serret坐标系研究滚珠丝杠副的动力学特性、和接触刚度求解；对滚珠丝杠副的滚动体承载、蠕滑以及摩擦特性进行研究。基于上述学者的研究，可较为准确的求取丝杠螺母副的接触刚度。

机床工作台位置在工作时是沿丝杠轴向不断变化的，这会导致机床进给系统的动态性能、固有频率以及振动响应等均随之发生变化，而加工过程中工作台所受的多为变载荷，所带来的振动问题将更为突出，将会影响加工精度。而上述文献中，多将丝杠作为集中质量处理建立数学模型或仅仅应用有限元分析软件建立进给系统的物理模型，影响计算；精度；同时工作台位置对系统，特别是工作台自身的动态特性影响的分析并未提及。因此，对于机床进给系统精细的动力学建模以及对工作台的动态特性分析还有必要进行深入研究。

基于有限元法，以某数控车床为例，将丝杠进行单元划分，建立考虑轴向和扭转耦合的进给系统的动力学方程，讨论了单元体数目对提高运算精度和计算时间的影响，对于机床工作台在不同位置下影响其轴向振动模态所对应的固有频率、工作台动态响应进行研究，并进行了试验验证，证明在合理选取单元体数目的情况下该模型具有较高准确性以及运算效率。

2 基于有限元法的系统动力学方程

2.1 进给系统模型与单元体划分

以丝杠螺母副以及工作台刀具为研究对象，考虑轴向、扭转变形，建立的进给单元的动力学模型，如图1所示。

图1 滚珠丝杠进给系统模型Fig.1 Simple Model of Ball Screw Drive System

工作台简化为质量为ms的质量块，以丝杠左端为坐标原点建立坐标系，设工作台沿z方向的位移为u s(t)，扭转位角位移为φs(t)，丝杠螺母副沿z方向的轴向接触刚度和阻尼分别记为k u、cu；绕z轴的扭转刚度和阻尼记为kφ、cφ。前、后支承轴承沿z方向移动的轴向刚度和阻尼分别记为；两轴承绕z轴转动的扭转刚度和阻尼分别记为

由于丝杠为连续体，其振型函数和运动方程应是关于时间t和空间z的函数，应用有限元法，将丝杠划分为n个单元体，图2所示为丝杠第i个单元体（i=1，2…n），其中，节点i和节点i+1上沿z方向的位移分别为u i(t)和u i+1(t)，绕z轴扭转的角位移分别为φi(t)和φi+1(t)。设单元广义坐标自由度顺序为：

图2 丝杠单元体Fig.2 Element of Screw

对节点单元应用Lagrange方程，得到单元i的质量矩阵m i和刚度矩阵k i，i=2，3，…，n-1。

式中：ρ—丝杠密度；l为丝杠单元长度；A—丝杠截面积；E—丝杠的弹性模量；I—丝杠截面极惯性矩；J—丝杠截面中性轴极惯性矩；G—剪切弹性模量。

2.2 特殊单元

同时，考虑特殊单元体：和左右轴承接触的两个边界单元，以及与工作台螺母接触的单元体，并分别求出其各自的质量矩阵和刚度矩阵以及相关耦合项。

（1）与左右轴承接触的边界单元

边界单元的质量矩阵与中间单元的相同，即：

以左侧边界单元为例，如图3所示。对节点单元应用La⁃grange方程，可得其单元刚度矩阵：

图3 丝杠边界单元体Fig.3 Boundary Element of Screw

同理，可得到右侧单元的刚度矩阵：

（2）与工作台连接的第j个丝杠单元

设丝杠第j个节点与螺母工作台接触，如图4所示，节点j和节点j+1上沿z方向的位移分别为u j(t)和u j+1(t)，绕z轴扭转的角位移分别为。计算工作台的动能：

图4 工作台与丝杠单元体Fig.4 The Worktable and Screw Element

由式（7）可知，丝杠单元体j的质量矩阵中不存在与工作台质量的耦合。

计算工作台的势能：

将V s分别对u j、φj、u s和φs求偏导，得：

从上式可以看出，第j个丝杠单元刚度矩阵中存在关于u s和φs的耦合项。

2.3 单元组集

丝杠单元体的单元质量矩阵和刚度矩阵按丝杠位移向量通过矩阵变换（式（10）），可得到丝杠的质量和刚度矩阵。

式中：A i—变换矩阵。

考虑工作台轴向位移的自由度，整个系统的自由度为：

2.4 进给系统动力学方程

结合丝杠的质量矩阵和刚度矩阵以及相关耦合项，可得系统的质量和刚度矩阵，进而求得进给系统运动微分方程为：

式中：M—系统质量矩阵；K—系统刚度矩阵；F—激励向量；q—系统广义坐标。

这里研究的对象为典型的数控车床z轴进给系统，其相关的参数，如表1所示。

表1 进给系统的参数Tab.1 Technical Data of Feed Drive System

3 机床工作台动态特性分析

令F=0，方程（3）的特征方程为：

由方程（13）可求解出系统的所有固有频率和固有振型，振动系统的固有频率数目与系统自由度数目相同，由于这里系统的动力学方程是基于有限单元法建立，故系统的固有频率数目众多，由于篇幅关系，这里只研究对系统影响较大的前3阶固有频率进行分析。由于丝杠单元体n的数目影响方程的运算精度和运算效率，数目过少运算速度快但精度低，数目过高则反之，故取不同单元体试算z∕L=0.5处的前3阶固有频率，发现当n取120以上时，前三阶固有频率值均趋于稳定收敛，而当n超过90之后，由于迭代次数增多，程序运行的时间损耗也大幅度提高。故综合考虑运算精度和速度这里仿真运算中均取丝杠单元体数目n=120。

取工作台位置为z∕L=0.5处，应用MATLAB求解方程（13），得到系统的前3阶固有频率及其对应的模态，如表2所示。由表2可以看出，工作台轴向振动模态所对应的是系统的第1、3阶固有频率。

表2 z/L=0.5时系统前3阶固有频率Tab.2 The First 3 Natural Frequencies as z/L=0.5

但机床在工作时，工作台位置并非是固定不变的，这也就意味着在动力学方程（3）中，与工作台发生耦合的丝杠单元是不断变化的，故在工作过程中系统为时变系统，其质量矩阵和刚度矩阵均会因为工作台位置变化而发生变化，进而系统的固有特性也会发变化。工作台位置在z∕L=0.1～0.9范围内变化时，影响工作台轴向振动模态的系统第1、3阶固有频率的变化规律，如图5所示。从图5（a）中可以看出，系统的第1阶固有频率呈现先减后增的变化趋势，其中固有频率最大的工作台位置出现在丝杠z∕L=0.1和0.9处，即丝杠端部；而固有频率最大的工作台位置出现在丝杠z∕L=0.5处，即丝杠中部处，固有频率的变化率达到了50.58%。

从图5（b）中可以看出，系统的第2阶固有频率只在工作台处于两端处突然呈现较大的变化，变化率为8.48%，而沿丝杠其他位置处基本保持不变。由于工作时，工作台位置不会达到两端部，故可以认为工作台在z轴位置变化对系统的第2阶固有频率影响不大。

从图5（c）中可见，系统的第3阶固有频率的变化趋势与第1阶恰好相反，呈现先增大后减小的变化规律，其中固有频率最大的工作台位置均出现在丝杠z∕L=0.5处，而固有频率最小的工作台位置亦都出现在丝杠z∕L=0.1和0.9处。固有频率的变化率为10.3 %。

图5 系统前3阶固有频率随工作台位置的变化规律Fig.5 Variation of the First Three Natural Frequencies with Positions of Worktable

4 试验验证

为验证这一规律，这里采用一台典型的数控车床作为试验对象，采用锤击法对机床工作台特殊点进行了动态特性试验，如图6所示。

图6 试验系统Fig.6 Experiment System

为使测试更加准确，应用086M92ES型电子冲击力锤，可提供可重复的稳定冲击力，并在工作台的关键点处布置压电式加速度传感器的测振点，加速度传感器主要拾取在敲击点处施加沿z轴的脉冲激励后工作台的速度振动响应信号，经过滤波、快速傅里叶变换等处理，便可以得到工作台在频域上的速度振动响应。试验选取工作台处在z∕L=0.2、z∕L=0.35、z∕L=0.5、z∕L=0.65以及z∕L=0.8五个位置作为试验点，拾取速度响应信号。通过测试得到的频率成分，发现在有两个明显的频率分别与理论解中的第1和3阶固有频率相对应，因此将测试值与理论值一起绘出，如图5a、c所示，第1阶固有频率的仿真结果与试验的误差在4.5%以内，第3阶固有频率的仿真结果与试验结果的误差在2%以内，且变化趋势与仿真计算得出的结论相吻合。由此可见，这里的进给系统动力学方程在求解系统固有频率，特别是求系统中、低段固有频率时均具有较高的精度。

5 结论

（1）应用有限元法，将丝杠作为连续体划分单元体，建立机床进给系统轴向和扭转振动的动力学方程，可求解系统、特别是丝杠连续体各个位置处的动态特性。通过讨论单元体数目对提高运算精度和计算时间的影响，得到当丝杠单元体数目取120时，系统前三阶固有频率的仿真结果均趋于稳定收敛，此时仿真结果与试验结果的误差低于4.5%。（2）轴向振动模态频率随工作台位置呈现非线性的变化规律：第1阶固有频率呈现先减后增的变化趋势，其中固有频率最大的工作台位置出现在丝杠端部；第3阶固有频率则为先增后减的变化规律，其中固有频率最大的工作台位置均出现在丝杠中部。第2阶固有频率受工作台位置的影响并不显著。