一种数控机床双驱进给系统模态建模方法与流程

本发明属于机械振动分析
技术领域：
，特别是一种数控机床双驱进给系统模态建模方法。
背景技术：
：作为机床重要部件之一的进给系统，其动态性能是决定机床性能的重要因素，在设计阶段就需要考虑机床的动态性能，机床模态要远离共振区域，减小共振对加工的影响。目前主要有建模和有限元软件模拟两种方法。许向荣(许向荣.滚珠丝杠副直线导轨进给单元动态性能研究[d].济南：山东大学，2011)考虑滚动结合面的影响，对滚珠丝杠副直线导轨单元进行动力学建模，根据拉格朗日方程推导了进给单元的微分振动方程，获得了进给单元的质量矩阵、刚度矩阵及阻尼矩阵。张会端(张会端.机床进给系统的动力学分析[d].吉林：吉林大学，2009)根据机床进给系统的特点，建立其动力学模型，列出系统的动能、势能、耗能函数和外载荷所作的虚功，利用拉格朗日方程推导系统的运动学方程。上述研究由于均把系统当做连续系统研究，导致无法求解出具体的振型函数，所以无法求解出模态参数。王永强(王永强.滚珠丝杠进给系统自适应建模理论与方法研究[d].济南：山东大学，2013)运用拉个朗日方程建立了滚珠丝杠进给系统的二维动力学模型，求解出系统的一阶固有频率，由于建立的模型是在平面上的单自由度模型，所以不能求解二阶及以上的固有频率。在模态分析中，利用有限元软件进行分析也是常见的手段之一，在边界条件设置准确的情况下，有限元软件分析可以做到比较精准的分析。但是也存在明显的局限性：软件的通用性不强，每次分析都需要准确设置边界情况，操作比较繁琐；软件后处理只输出结果，对于每个不同的模型会有不同的输出结果，但是不能显示出对模态影响大的因素。技术实现要素：本发明所解决的技术问题在于提供一种数控机床双驱进给系统模态建模方法，以解决现有技术中连续系统振型方程无法求解和有限元软件模拟操作繁琐等问题；本发明提供一种能够具有可操作性的且求解效率更高的模态分析方法。实现本发明目的的技术解决方案为：一种数控机床双驱进给系统模态建模方法，包括以下步骤：步骤1、确定系统的广义坐标：判断系统的自由度数，并选取广义坐标qj来描述系统的运动状态，其坐标数目和自由度数相同；2、步骤2、建立机床进给系统的动能表达式：即分别建立丝杠的动能和工作台的动能表达式，进一步建立机床进给系统总的动能表达式；步骤3、建立进给系统的势能表达式：即分别建立导轨和滑块接触面的势能、轴承内外圈接触面的势能、丝杠和螺母接触面的势能、丝杠的轴向势能和扭转势能；进一步建立进给系统总势能表达式；步骤4、建立进给系统耗能表达式：即分别建立导轨和滑块接触面的耗能、轴承内外圈接触面的耗能、丝杠和螺母接触面的耗能、丝杠的耗能；进一步建立进给系统总耗能表达式；步骤5、建立滑块坐标与工作台坐标的的位置转化关系式：通过坐标转换得到滑块坐标与工作台坐标的数学关系式；步骤6、将能量表达式带入拉格朗日方程，得到进给系统的运动微分方程；将运动微分方程写成矩阵形式：式中m为质量矩阵，c为阻尼矩阵，k为刚度矩阵，f为激振力矩阵；步骤7、求解模态固有频率；简化方程(7)为判断矩阵行列式的值是否恒为零，分别采用矩阵迭代法和根据来求解模态固有频率。本发明与现有技术相比，其显著优点：(1)本发明的建模方法更加直观：变量可以控制，根据变量的改变可以观察模态的改变，变量与模态之间的关系更加直观。(2)本发明的建模方法改进更具有针对性：从质量矩阵m和刚度矩阵k的数值可以分析推导出对模态影响较大的参数，在改进的时候可以通过分析薄弱环节，有针对性的改进。(3)本发明的建模方法具有更好的可操作性：由能量法求解进给系统固有频率的方法，具有良好的可操作性，求解的效率比较高，借助于正交试验法可以实现机床进给系统的固有频率反向设计，对比有限元仿真具有一定的优势。下面结合附图对本发明作进一步详细描述。附图说明图1为本发明方法流程示意图。图2为进给系统位置关系和系统坐标图。图3(a-c)分别为滑块及工作台坐标的空间示意图、主视图、俯视图。具体实施方式为了说明本发明的技术方案及技术目的，下面结合附图及具体实施例对本发明做进一步的介绍。结合图1的流程示意图，本发明的一种数控机床双驱进给系统模态建模方法，包括以下步骤：步骤1、确定系统的广义坐标：判断系统的自由度数，并选取广义坐标来描述系统的运动状态，其坐标数目和自由度数相同。结合图2，以处于静止状态时的工作台重心位置为坐标原点，建立o-xyz主坐标系，x轴的正方向垂直于导轨，y轴的正方向垂直于工作台上表面向上，z轴的正方向为工作台远离电机的方向，z轴方向为工作台进给方向。建立广义坐标：qj＝{xw,yw,zw,α,β,γ,z1,z2,z3,θ1,θ2,θ3,x,y,…}j＝(1,2,3,....n)；其中j为广义坐标系qj的坐标编号，如q1＝{xw}；“…”为省略的其他自由度，为n为系统的自由度数目。其中，xw、yw、zw分别为工作台在x、y、z方向上的横向振动位移；α、β、γ分别为工作台绕x、y、z方向的扭转角度；z1为丝杠左端(远离电机一端)的轴向位移，θ1为丝杠绕z方向的转动角度；z2为螺母中心的轴向位移，θ2为螺母绕z方向的转动角度；z3为丝杠右端(靠近电机一端)的轴向位移，θ3为丝杠绕z方向的转动角度；x、y分别为丝杠的径向平行与x轴和y轴的位移；下面均以自由度数n＝14为例，即qj＝{xw,yw,zw,α,β,γ,z1,z2,z3,θ1,θ2,θ3,x,y}j＝(1,2,3,....n)。步骤2、建立机床进给系统的动能模型：即分别建立丝杠的动能和工作台的动能表达式，进一步建立机床进给系统总的动能表达式；步骤2.1、建立工作台的动能方程：工作台的动能包括在x、y、z方向上的横向振动动能和绕x、y、z方向的扭转动能，工作台的动能tw表达式为：式中mw为工作台的质量，jx为工作台绕x轴的转动惯量，jy工作台绕y轴的转动惯量，jz工作台绕z轴的转动惯量。步骤2.2建立丝杠的动能方程：丝杠的动能包括丝杠的扭转动能、轴向动能以及横向动能，(1)建立丝杠轴向动能：对于螺母左端丝杠，设距离丝杠左端为z的位置有长度为dz的极小一段，则该位置在z向的位移方程为：其中b为螺母左端丝杠长度；对上述公式进行积分得到螺母左端丝杠的轴向动能tl表达式：同理，对螺母右端丝杠进行积分运算，右端丝杠动能tr为mbs为丝杠质量，a为丝杆的有效长度(两个轴承之间丝杠的长度)；(2)建立丝杠横向动能表达式为：(3)建立丝杠扭转动能表达式为：jbsl为螺母左端丝杠的转动惯量，jbsr为螺母右端丝杠的转动惯量；(4)建立丝杠总动能tbs表达式为：2.3、机床进给系统总动能为：步骤3、建立进给系统的势能表达式：即分别建立导轨和滑块接触面的势能、轴承内外圈接触面的势能、丝杠和螺母接触面的势能、丝杠的轴向势能和扭转势能；进一步建立进给系统总势能表达式；步骤3.1、建立导轨和滑块接触面的势能方程：导轨和滑块接触面的势能us表达式为：其中，ksx为滑块与导轨侧向接触刚度，ksy为滑块与导轨法向接触刚度。步骤3.2、建立轴承内外圈结合面势能方程：轴承内外圈结合面势能ub表达式为：其中，kbx、kby为轴承径向刚度，kbz为轴承轴向刚度。步骤3.3、建立丝杠螺母结合面势能势能方程：杠螺母结合面势能势能un表达式为：其中，knx、kny为丝杠螺母副径向接触刚度，knz为丝杠螺母副轴向接触刚度。步骤3.4、建立丝杠的势能方程：丝杠的势能包括丝杠的轴向势能和扭转势能：ubs＝uz+uθ(1)丝杠轴向势能uz表达式：uz＝kbszl(z2-z1)2+kbszr(z3-z2)2其中，kbszl为螺母左侧丝杠轴向刚度，kbszr为螺母右侧丝杠轴向刚度。(2)丝杠扭转势能uθ表达式：uθ＝kθl(θ2-θ1)2+kθr(θ3-θ2)2其中，kθl为螺母左侧丝杠扭转刚度，kθr为螺母右侧丝杠扭转刚度。步骤3.5、建立进给系统总势能的表达式：步骤4、建立进给系统耗能表达式：即分别建立导轨和滑块接触面的耗能、轴承内外圈接触面的耗能、丝杠和螺母接触面的耗能、丝杠的耗能；进一步建立进给系统总耗能表达式；步骤4.1、建立导轨和滑块接触面耗能方程：轨和滑块接触面耗能ds表达式：其中，csx为滑块与导轨侧向接触阻尼，csy为滑块与导轨法向接触阻尼；xsi为第i个滑块在x方向的位移，ysi为第i个滑块在y方向的位移；步骤4.2、建立轴承内外圈结合面耗能方程：轴承内外圈结合面耗能db表达式：其中，cbx、cby为轴承径向阻尼，cbz为轴承轴向阻尼。步骤4.3、建立丝杠螺母结合面耗能方程：丝杠螺母结合面耗能dn表达式：其中，cnx、cny为丝杠螺母副径向接触阻尼，cnz为丝杠螺母副轴向接触阻尼。步骤4.4、建立丝杠的耗能方程：丝杠的耗能dbs表达式：其中，cbsl为螺母左侧丝杠轴向阻尼，cbsr为螺母右侧丝杠轴向阻尼。步骤4.5、建立进给系统总耗能的表达式：步骤5、建立滑块坐标与工作台坐标的的位置转化关系式：建立动态坐标系，通过坐标转换得到滑块坐标与工作台坐标在动态坐标系；结合图3滑块与工作台坐标示意图，以工作台任意时刻的重心为坐标原点，建立o1-x1y1z1动态坐标系，x1、y1、z1的方向分别与x、y、z方向一致。四个滑块中心点分别为s1,s2,s3,s4，第i个滑块中心与原点o在x方向上的距离为dsi，在y方向上的距离为hsi，在z方向上的距离为lsi。工作台在x、y、z三个方向上的横向振动位移分别为xw、yw、zw，绕x、y、z的扭转位移分别为α、β、γ。步骤5.1通过坐标转换得到滑块坐标与工作台坐标在动态坐标系下的关系为：5.2由于α、β、γ的值很小，所以tanα＝α，tanβ＝β，tanγ＝γ，所以上式可以化简为：步骤6：将能量表达式带入拉格朗日方程，得到进给系统的运动微分方程；将运动微分方程写成矩阵形式。其中qj为系统的坐标，为系统的速度，t为系统的动能，u为系统的势能，t为时间，q'j为对应的力，n为系统的自由度数。下面仍以自由度n＝14为例，分别计算广义坐标qj＝{xw,yw,zw,α,β,γ,z1,z2,z3,θ1,θ2,θ3,x,y}中每个坐标系下拉格朗日方程中各项导数；(1)以q1＝{xw}为例，联立方程式(1)-(6)对xw计算拉格朗日方程中各项导数：将上述各式整理得运动微分方程为：(2)对广义坐标qj＝{xw,yw,zw,α,β,γ,z1,z2,z3,θ1,θ2,θ3,x,y}每个坐标分别求导，微分方程为：式中f1-f14分别为每组方程对应的激振力，由于求解过程中令方程右端为零，所以f1-f14值未知对求解过程没有影响。(3)将运动微分方程写成矩阵形式式中m为质量矩阵，c为阻尼矩阵，k为刚度矩阵，f为激振力矩阵。其中：步骤7、求解模态固有频率：步骤7.1、矩阵行列式的值不恒为零时解法：由于有阻尼的固有频率ωd和无阻尼的固有频率ωn存在的关系，当进给系统阻尼比ξ＝10％时，ωd＝0.99499ωn，且一般情况下进给系统阻尼比ξ≤10％，所以求解进给系统固有频率时不考虑阻尼矩阵c，同时也不考虑作用力f。所以方程(7)化简为式(8)有的解为x为振幅，ωn为固有频率，为相位角。将式(9)解代入式(8)，得到如下振型方程上式存在非零解的条件是：当且仅当系数行列式为零，即：式中δ(ω2)为特征值行列式，将其展开后可以得到ωn的n次代数方程式：ωn2n+a1ωn2(n-1)+a2ωn2(n-2)+…+an-1ωn2+an＝0式中a1,a2…an为常系数。这一方程有n个根ωn2(r＝1,2,…n)，这些根称为特征值，特征值的平方根ωn称为进给系统的固有频率。步骤7.2、矩阵行列式的值恒为零时的解法：特别地：在机床进给系统中，当工作台位于丝杠中点时矩阵为半正定矩阵，所以矩阵行列式的值恒为零，此时用上述方法无法求解方程，需采用矩阵迭代法求解。同时由于刚度矩阵是半正定矩阵，所以用matlab中pinv语言求解刚度矩阵的伪逆矩阵，通过迭代求取工作台位于中点时进给系统的固有频率。本方法具有一定的通用性，适用于双导轨双丝杠四滑块进给系统，符合上述要求的进给系统都可以用此建模方法求解。例1.以mch63数控机床进给系统为例，整理得到矩阵中所需参数(具体参见表3)，通过matlab编程求解得到系统固有频率参数，整理之后得到进给系统前六阶固有频率计算值，如表1所示。表1.mch63数控机床进给系统参数ksx(n/m)ksy(n/m)kbx(n/m)kby(n/m)kbz(n/m)kbszr(n/m)kbs(n/m)7.7*1083.45*1092.4*1082.4*1083.2*1083.2*1081.6*108knz(n/m)kbszl(n/m)kθl(n/m)kθr(n/m)kθ(n/m)knx(n/m)kny(n/m)1*1093.2*1081.734*1041.734*1048.67*1038*1088*108ds(m)hs(m)jx(kgm2)jy(kgm2)jz(kgm2)a(m)b(m)0.327m0.17045.73285.22958.8751.20～1.2mw(kg)mbs(kg)jbs(kgm2)jbsl(kgm2)jbs(kgm2)eg727.6413.2672.540*10-3jbsb/a2.540*10-32.06*10117.938*1010ls1＝ls2(m)ls3＝ls4(m)ρbs(kg/m3)ρw(kg/m3)jbsr(kgm2)0.3110.27178507300jbs(a-b)/a表中ds为滑块与原点o在x方向上的距离，hs为滑块与原点o在y方向上的距离，e为丝杠材料弹性模量，g为丝杠材料剪切弹性模量(可由e、g计算出刚度值k，上表已算出),ls1为第一个滑块与原点在z方向的距离，ls2为第二个滑块与原点在z方向的距离，ls3为第三个滑块与原点在z方向的距离，ls4为第四个滑块与原点在z方向的距离，ρbs为丝杠材料的密度，ρw为工作台材料的密度。表2固有频率计算值阶数建模理论计算频率/hz仿真计算频率/hz误差145.248.67.5％291.6694.12.6％3105.5106.61.0％4227.4227.70.1％5269.5306.912.1％6615.1615.220.01％7615.1616.040.15％表2中建模理论计算结果与ansysy软件仿真结果误差最大为第五阶频率：12％，其他阶数的误差值都比较小，说明了建模方法具有一定的准确性。通过矩阵可以解耦出结构和固有频率的关系，从矩阵中的参数可以比较出对固有频率影响大的因素，在进行优化时，可以首先优化这些参数，使优化的目标更加明确。虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明。本发明所属
技术领域：
中具有通常知识者，在不脱离本发明的精神和范围内，当可作各种的更动与润饰。因此，本发明的保护范围当视权利要求书所界定者为准。当前第1页12