Mechanical dynamics modeling and analysis of the special machine for vertical inner cavity
-
摘要: 针对立式四轴舱体内腔加工专机，基于动态子结构法进行了机械动力学建模。求解了专机动结合部的滚珠丝杠副和直线导轨副的动态刚度特性。由建立的动力学模型分析整机的固有模态，应用解析法求解整机的固有频率。应用有限元软件建立专机有限元模型，求解结合部刚性接触与设定了动态刚度性的整机固有频率与振型及谐响应分析结果。结果显示建立的含动态刚度的动力学模型能够有效镜像专机动力学特性，可应用于专机的动力分析与数字孪生虚拟样机研究。Abstract: For the vertical four-axis inner cavity machining special machine,the mechanical dynamics modeling is carried out based on the dynamic sub-structure method. The dynamic stiffness characteristics of the ball screw pair and linear guide pair of the special motorized joint are solved. The natural mode of the whole machine is analyzed by the established dynamic mode, and the natural frequency is solved by the analytical method. The finite element mode of the machine is established by the finite element software, and the analysis results of the natural frequency, mode shape and harmonic response of the whole machine with the rigid contact of the joint and the dynamic stiffness are set. The results show that the established dynamic model with dynamic stiffness can effectively mirror the dynamic characteristics of the machine, and can be applied to the dynamic analysis of the machine and the research of the digital twin virtual prototype.
-
随着我国航空航天技术的发展，越来越多的新结构、新工艺正应用在各类回转体类航天产品中。例如航天薄壁类圆形舱体的内壁需要镗铣加工，舱体内腔存在许多凸台、零件安装面和通孔等特征需要加工[1]。内腔加工专机是能够满足此类机加工需求的专用设备，已经得到了广泛应用。
航天典型产品的内腔加工，刀具路径为复杂空间轨迹，加工中专机的动态误差远远超过几何误差等静态误差，是影响加工误差的主要因素[2]。对内腔加工专机开展机械动力学特性的研究与分析是加工动态误差辨识与补偿，提高内腔加工精度与产品质量的前提。目前工程中对内腔加工专机进行动力学建模的普遍方法是，设置专机CAD模型的材料数据与刚性接触，采用有限元分析软件求解固有频率与振型[3-4]，分析动力学特性。但是，这种简化假设没有考虑专机结合部的动态刚度特性，根据研究统计，机床整机60%~80%的刚度来自各种结合部[5]。因此无论是从专机设计、装配和动特性分析，还是从结构修改、建立虚拟样机来说，必须考虑专机结合部的动态刚度特性[6]。而且采用有限元分析软件进行的机械动力学建模分析为离线模型，不适用于数字孪生技术中虚拟实体镜像物理实体，数据实时交互的研究[7]。
本文以立式四轴舱体内腔加工专机为研究对象，基于动态子结构法建立机械动力学模型，求解了专机结合部的滚珠丝杠副、直线导轨副动态刚度。应用有限元软件建立专机有限元模型，求解结合部刚性接触与设定动态刚度性的整机固有频率与振型，最后对专机进行谐响应分析求解了加工振动激励下的专机响应。
1. 基于子结构法的机械动力学建模
立式四轴内腔加工用于加工大直径尺寸舱体，如图1所示主要由床身、立柱、滑枕、工作台和主轴部件组成，其中床身为T型底座与横梁的固连的整体。各部件由伺服电机通过丝杠螺母机构驱动，在直线导轨上往复运动，并通过光栅尺实现闭环控制，实现四轴联动加工。滑枕上安装有第二回转轴用于适当调节主轴角度，适应特定舱体内腔加工面。根据专机驱动形式，可划分为床身、立柱、滑枕、工作台和主轴5个子结构。
建立如图2所示坐标系，其中A轴不作为四轴联动轴。可以分析得知专机的刚性结合部有床身与机床垫铁的地脚螺栓结合部，滑枕与主轴的螺栓连接刚性结合部。专机的动结合部有X轴的工作台与床身的丝杠螺母副、直线导轨副结合部；Z轴的立柱与床身的丝杠螺母副、直线导轨副结合部；Y轴的立柱与滑枕的丝杠螺母副、直线导轨副结合部。
以机床床身底部的端点作为整个坐标系的原点，根据振动特点，将机床简化为具有20个自由度的集中参数动力学模型，各自由度的分配如下：
$ {q_1} $ 为床身沿$ y $ 轴的平动坐标$ {y_1} $ ；$ {q_2} $ 为床身绕${\textit{z}}$ 轴的转动坐标$ {\theta _{{\textit{z}}1}} $ ；$ {q_3} $ 为床身绕$ x $ 轴的转动坐标$ {\theta _{x1}} $ ；$ {q_4} $ 为立柱沿$ x $ 轴的平动坐标$ {x_1} $ ；$ {q_5} $ 为立柱沿$ y $ 轴的平动坐标$ {y_2} $ ；$ {q_6} $ 为立柱绕$ x $ 轴的转动坐标$ {\theta _{x2}} $ ；$ {q_7} $ 为立柱绕${\textit{z}} $ 轴的转动坐标$ {\theta _{{\textit{z}}2}} $ ；$ {q_8} $ 为滑枕沿$ x $ 轴的平动坐标$ {x_2} $ ；$ {q_9} $ 为滑枕沿$ {\textit{z}} $ 轴的平动坐标$ {{\textit{z}}_1} $ ；$ {q_{10}} $ 为滑枕绕$ y $ 轴的转动坐标$ {\theta _{y1}} $ ；$ {q_{11}} $ 为滑枕绕$ {\textit{z}} $ 轴的转动坐标$ {\theta _{{\textit{z}}3}} $ ；$ {q_{12}} $ 为主轴沿$ x $ 轴的平动坐标$ {x_3} $ ；$ {q_{13}} $ 为主轴沿$ y $ 轴的平动坐标$ {y_3} $ ；$ {q_{14}} $ 为主轴沿$ {\textit{z}} $ 轴的平动坐标$ {{\textit{z}}_2} $ ；$ {q_{15}} $ 为主轴绕$ x $ 轴的转动坐标$ {\theta _{x3}} $ ；$ {q_{16}} $ 为主轴绕$ {\textit{z}} $ 轴的转动坐标$ {\theta _{{\textit{z}}4}} $ ；$ {q_{17}} $ 为工作台沿$ y $ 轴的平动坐标$ {y_4} $ ；$ {q_{18}} $ 为工作台沿$ {\textit{z}} $ 轴的平动坐标$ {{\textit{z}}_3} $ ；$ {q_{19}} $ 为工作台绕$ x $ 轴的转动坐标$ {\theta _{x4}} $ ；$ {q_{20}} $ 为工作台绕$ {\textit{z}} $ 轴的转动坐标$ {\theta _{{\textit{z}}5}} $ 。依照划分的子结构、建立的坐标系，得到内腔加工专机机械动力学建模系统简图如图3所示。
图3中
$ {c_i} $ 为各子结构质心，$ {a_i} $ 、$ {b_i} $ 与$ {d_i} $ 为各子结构质心相对位置关系。用拉格朗日方程建立系统的运动方程，拉格朗日方程的一般表达式[8]为
$$ \frac{{\text{d}}}{{{\text{d}}t}}\left[ {\frac{{\partial T}}{{\partial {{\dot q}_j}}}} \right] - \frac{{\partial T}}{{\partial {q_j}}} = {Q_j} \quad j = 1,2, \cdots ,n $$ (1) 式中：
$ T $ 为系统的总动能；$ {q_j} $ 为系统的广义坐标；$ {Q_j} $ 为广义力。对于所研究的机械系统，$ {Q_j} $ 由广义势力$ - \dfrac{{\partial U}}{{\partial {q_j}}} $ 、广义线性阻尼力$ - \dfrac{{\partial D}}{{\partial {q_j}}} $ 和广义激振力$ {Q_j}' $ 组成，即：$$ {Q_j} = - \frac{{\partial U}}{{\partial {q_j}}} - \frac{{\partial D}}{{\partial {{\dot q}_j}}} + {Q_j}' $$ (2) 将式（2）代入式（1）中，得
$$ \frac{{\rm{d}}}{{{\rm{d}}t}}\left[ {\frac{{\partial T}}{{\partial {{\dot q}_j}}}} \right] - \frac{{\partial T}}{{\partial {q_j}}} + \frac{{\partial U}}{{\partial {q_j}}} + \frac{{\partial D}}{{\partial {{\dot q}_j}}} = {Q_j}' \quad j = 1,2, \cdots ,n $$ (3) 式中：
$ U $ 、$ D $ 分别为系统的总势能、瑞利耗能函数。根据动力学模型，算得
$ T $ 、$ U $ 、$ D $ ，系统的总动能为$$ \begin{split} T = &\frac{1}{2}{m_1}\dot q_1^2 + \frac{1}{2}{I_{{\textit{z}}1}}\dot q_2^2 + \frac{1}{2}{I_{x1}}\dot q_3^2 + \frac{1}{2}{m_2}({{\dot q}_4} - \\& {a_2}{{\dot q}_2}{)^2} + \frac{1}{2}{m_2}({{\dot q}_5} + {{\dot q}_1} - {d_3}{{\dot q}_2} - {d_2}{{\dot q}_3}) + \\& \frac{1}{2}{I_{x2}}\dot q_6^2 + \frac{1}{2}{I_{{\textit{z}}2}}\dot q_7^2 + \frac{1}{2}{m_3}[{{\dot q}_8} - \\& ({a_2} + {a_3}){{\dot q}_2} + {{\dot q}_4} - {a_3}{{\dot q}_7}{]^2} + \\& \frac{1}{2}{m_3}{[{{\dot q}_9} - ({a_2} + {a_3}){{\dot q}_3} + {{\dot q}_4} - {a_3}{{\dot q}_6}]^2} + \\& \frac{1}{2}{I_{y3}}\dot q_{10}^2 + \frac{1}{2}{I_{{\textit{z}}3}}\dot q_{11}^2 + \frac{1}{2}{m_4}[{{\dot q}_{12}} - \\& ({a_2} + {a_3}){{\dot q}_2} + {{\dot q}_4} - {a_3}{{\dot q}_7} + {{\dot q}_8} + ({a_8}){{\dot q}_{11}}{]^2} + \\& \frac{1}{2}{m_4}[{{\dot q}_{13}} + {{\dot q}_1} + ({d_6} + {d_{10}} - {d_3}){{\dot q}_2} - {d_2}{{\dot q}_3} + \\& {{\dot q}_5} + ({d_6} + {d_{10}}){{\dot q}_7} + {d_{10}}{{\dot q}_{11}}{]^2} + \frac{1}{2}{m_4}[{{\dot q}_{14}} + \\& ({a_2} + {a_3} + {a_5}){{\dot q}_3} + ({a_3} + {a_5}){{\dot q}_6} + {{\dot q}_9} - \\& {d_{10}}{{\dot q}_{10}}{]^2} + \frac{1}{2}{I_{x4}}\dot q_{15}^2 + \frac{1}{2}{I_{{\textit{z}}4}}\dot q_{16}^2 + \frac{1}{2}{m_5}[{{\dot q}_{17}} + \\& {{\dot q}_1} + ({d_3} + {d_7}){{\dot q}_2}{]^2} + \frac{1}{2}{m_5}{({{\dot q}_{18}} + {a_4}{{\dot q}_2})^2} + \\& \frac{1}{2}{I_{x5}}\dot q_{19}^2 + \frac{1}{2}{I_{{\textit{z}}5}}\dot q_{20}^2 \end{split} $$ (4) 系统的总势能为
$$ \begin{split} U =& \frac{1}{2}\{ 2K_y^{\text{①}}{({q_1} - {d_1}{q_2})^2} + 2K_y^{\text{①}}({q_1}+ \\& {d_2}{q_2}{)^2} + 2K_y^{\text{①}}{({q_1} - {b_1}{q_3})^2} + 2K_y^{\text{①}}({q_1}+ \\& {b_2}{q_3}{)^2} + 2K_{\theta {\textit{z}}}^{\text{①}}{q_2}^2 + 2K_{\theta x}^{\text{①}}{q_3}^2 + 2K_x^{\text{②}}{q_4}^2+ \\& 2K_y^{\text{②}}{({q_5} - {b_2}{q_6})^2} + 2K_y^{\text{②}}{({q_5} + {b_3}{q_6})^2}+ \\& 2K_y^{\text{②}}{({q_5} - {d_4}{q_7})^2} + 2K_y^{\text{②}}{({q_5} + {d_5}{q_7})^2}+ \\& 2K_{\theta x}^{\text{②}}{({q_6} - {q_3})^2} + 2K_{\theta {\textit{z}}}^{\text{②}}{({q_7} - {q_2})^2}+ \\& 2K_x^{\text{③}}{({q_8} - {a_6}{q_{11}})^2} + 2K_x^{\text{③}}{({q_8} + {a_7}{q_{11}})^2}+ \\& 2K_x^{\text{③}}{({q_8} - {b_6}{q_{10}})^2} + 2K_x^{\text{③}}{({q_8} + {b_7}{q_{10}})^2}+ \\& 2K_{\textit{z}}^{\text{③}}{q_9}^2 + 2K_{\theta y}^{\text{③}}{q_{10}}^2 + 2K_{\theta {\textit{z}}}^{\text{③}}{({q_{11}} - {q_7})^2}+ \\& 2K_x^{\text{④}}{({q_{12}} - {a_8}{q_{16}})^2} + 2K_x^{\text{④}}{({q_{12}} + {a_8}{q_{16}})^2}+ \\& 2K_y^{\text{④}}{q_{13}}^2 + 2K_{\textit{z}}^{\text{④}}{q_{14}}^2 + 2K_{\theta x}^{\text{④}}{q_{15}}^2+ \\& 2K_{\theta {\textit{z}}}^{\text{④}}{({q_{16}} - {q_{11}})^2} + 2K_y^{\text{⑤}}{({q_{17}} - {b_4}{q_{19}})^2}+ \\& 2K_y^{\text{⑤}}{({q_{17}} + {b_5}{q_{19}})^2} + 2K_y^{\text{⑤}}{({q_{17}} - {d_8}{q_{20}})^2}+ \\& 2K_y^{\text{⑤}}{({q_{17}} + {d_9}{q_{20}})^2} + 2K_{\textit{z}}^{\text{⑤}}{q_{17}}^2+ \\& 2K_{\theta x}^{\text{⑤}}{({q_{19}} - {q_3})^2} + 2K_{\theta {\textit{z}}}^{\text{⑤}}{({q_{20}} - {q_2})^2}\} \end{split} $$ (5) 系统的瑞利耗能函数为
$$ \begin{split}D=&\frac{1}{2}\{2{C}_{y}^{\text{①}}{(}{{\dot{q}}_{1}- {{d}_{1}}{{\dot {q}}_{2}}{)}^2}+2{C}_{y}^{\text{①}}({\dot{q}}_{1}+\\ & {d}_{2}{\dot{q}}_{2}{)}^{2}+2{C}_{y}^{\text{①}}{(}{{\dot{q}}_{1}- {{b}_{1}}{{\dot {q}}_{3}}{)}^2}+2{C}_{y}^{\text{①}}({\dot{q}}_{1}+\\ & {b}_{2}{\dot{q}}_{3}{)}^{2}+2{C}_{\theta {\textit{z}}}^{\text{①}}{\dot{q}}_{2}{}^{2}+2{C}_{\theta x}^{\text{①}}{\dot{q}}_{3}^{2}+2{C}_{x}^{\text{②}}{\dot{q}}_{4}^{2}+ \\ &2{C}_{y}^{{\text{②}}}{(}{{\dot{q}}_{5}}- {{b_2}{{\dot {q}}_6}{)}^2}+2{C}_{y}^{{\text{②}}}{(}{{\dot{q}}_{5}+ {b_3}{{\dot {q}}_6}{)}^2}+\\ & 2{C}_{y}^{\text{②}}{(}{{\dot{q}}_{5}}- {{d_4}{{\dot {q}}_7}{)}^2}+2{C}_{y}^{{\text{②}}}{(}{{\dot{q}}_{5}+ {d_5}{{\dot {q}}_7}{)}^2}+\\ & 2{C}_{\theta x}^{\text{②}}{(}{{\dot{q}}_{6}}- {{{\dot {q}}_3}{)}^2}+2{C}_{\theta {\textit{z}}}^{\text{②}}{(}{{\dot{q}}_{7}- {{\dot {q}}_2}{)}^2}+\\& 2{C}_{x}^{\text{③}}{(}{{\dot{q}}_{8} - {a_6}{{\dot {q}}_{11}}{)}^2}+2{C}_{x}^{\text{③}}{(}{{\dot{q}}_{8}+ {a_7}{{\dot {q}}_{11}}{)}^2}+\\& 2{C}_{x}^{\text{③}}{(}{{\dot{q}}_{8}- {b_6}{{\dot {q}}_{10}}{)}^2}+2{C}_{x}^{\text{③}}{(}{{\dot{q}}_{8}+ {b_7}{{\dot {q}}_{10}}{)}^2}+\\& 2{C}_{{\textit{z}}}^{\text{③}}{\dot{q}}_{9}^{2}+2{C}_{\theta y}^{\text{③}}{\dot{q}}_{10}^{2}+2{C}_{\theta {\textit{z}}}^{\text{③}}{(}{{\dot{q}}_{11}- {{\dot {q}}_7}{)}^2}+ \\&2{C}_{x}^{\text{④}}{(}{{\dot{q}}_{12}- {a_8}{{\dot {q}}_{16}}{)}^2}+2{C}_{x}^{\text{④}}{(}{{\dot{q}}_{12} + {a_8}{{\dot {q}}_{16}}{)}^2}+ \\&2{C}_{y}^{\text{④}}{\dot{q}}_{13}^{2}+2{C}_{{\textit{z}}}^{\text{④}}{\dot{q}}_{14}^{2}+2{C}_{\theta x}^{\text{④}}{\dot{q}}_{15}^{2}+ \\&2{C}_{\theta {\textit{z}}}^{\text{④}}{(}{{\dot{q}}_{16}- {{\dot {q}}_{11}}{)}^2}+2{C}_{y}^{\text{⑤}}{(}{{\dot{q}}_{17}- {b_4}{{\dot {q}}_{19}}{)}^2}+ \\&2{C}_{y}^{\text{⑤}}{(}{{\dot{q}}_{17}+ {b_5}{{\dot {q}}_{19}}{)}^2}+2{C}_{y}^{\text{⑤}}{(}{{\dot{q}}_{17}+ {b_5}{{\dot {q}}_{19}}{)}^2}+ \\&2{C}_{y}^{\text{⑤}}{(}{{\dot{q}}_{17}+ {d_9}{{\dot {q}}_{20}}{)}^2}+2{C}_{{\textit{z}}}^{\text{⑤}}{\dot{q}}_{17}^{2}+ \\&2{C}_{\theta x}^{\text{⑤}}{(}{{\dot{q}}_{19}- {q_3}{)}^2}+2{C}_{\theta {\textit{z}}}^{\text{⑤}}{(}{{\dot{q}}_{20}- {{\dot {q}}_2}{)}^2}\}\end{split} $$ (6) 式（4）~（6）中，
$ {m_i} $ （$ i = 1,2, \cdots ,5 $ ）依次为5个子结构：床身、立柱、滑枕、工作台和主轴的质量；$ {I_{\varphi i}} $ （$ \varphi = x,y,{\textit{z}} $ ，$ i = 1,2, \cdots ,5 $ ）为5个子结构绕各轴的转动惯量；$ K_\omega ^\alpha $ （$ \omega = x,y,{\textit{z}} $ ，$ \alpha ={\text{①}},{\text{②}},\cdots,{\text{⑤}} $ ）为5个动结合部的各向线刚度；$ K_{\theta \omega }^\alpha $ （$ \omega = x,y,{\textit{z}} $ ，$ \alpha ={\text{①}},{\text{②}},\cdots ,{\text{⑤}} $ ）为5个动结合部的各向角刚度；同理，$ C $ 为各结合部阻尼。将式（4）~（6）代入式（3）中，整理后得专机动力学方程
$$ {M_{20 \times 20}}{\ddot q_{20 \times 1}} + {C_{20 \times 20}}{\dot q_{20 \times 1}} + {K_{20 \times 20}}{q_{20 \times 1}} = {Q_{20 \times 1}} $$ (7) 2. 专机结合部的刚度求解
式（5）中专机的动力学方程中各项刚度参数为各结合部参数，专机的各种结合部可以分为子结构螺栓连接的刚性接触结合部、各移动轴由滚珠丝杠副和直线导轨副组成的动结合部，针对各结合部特点，开展求解。
2.1 刚性结合部的刚度
专机的刚性接触部即螺栓结合部刚度由吉村允孝法[9]求解，根据螺栓参数、等效元素个数由吉村允孝法得表1。
表 1 专机刚性结合部刚度刚性结合部 线刚度/(N/m)、角刚度/(N·m/rad) 床身与地面 $ {K}_{y}^{\text{①}}=7.642\times {10}^{8} $、$ {K}_{\theta {\textit{z}}}^{\text{①}}=3.629\times {10}^{7} $、$ {K}_{\theta x}^{\text{①}}=5.473\times {10}^{7} $ 滑枕与主轴 $ {K}_{x}^{\text{④}}=7.218\times {10}^{7} $、$ {K}_{y}^{\text{④}}=2.512\times {10}^{8} $、$ {K}_{{\textit{z}}}^{\text{④}}=2.512\times {10}^{8} $
$ {K}_{\theta x}^{\text{④}}=4.538\times {10}^{7} $、$ {K}_{\theta {\textit{z}}}^{\text{④}}=4.341\times {10}^{7} $2.2 滚珠丝杠副的刚度
对于动结合部进行动力学建模分析，滚珠丝杠副轴向刚度
$ {k_x} $ 为与滚珠丝杠副相关联的零部件刚度的串联总和[10]，丝杠副主要由轴承、丝杠和螺母组件组成，其动力学模型如图4所示。由此得到丝杠传动系统的轴向刚度
$ {k_x} $ 可表示为$$ \frac{1}{{{k_x}}} = \frac{1}{{{k_S}}} + \frac{1}{{{k_N}}} + \frac{1}{{{k_B}}} + \frac{1}{{{k_H}}} $$ (8) 式中：
$ {k_x} $ 为滚珠丝杠副的轴向刚度；$ {k_S} $ 为丝杠轴的轴向刚度；$ {k_N} $ 为螺母组件的轴向刚度；$ {k_B} $ 为支撑轴承的轴向刚度；$ {k_H} $ 为螺母座的刚性，N/μm。丝杠的安装方式为两端固定，丝杠轴的轴向刚度由下式求得
$$ {k_S} = \frac{{A \cdot E \cdot L}}{{1000 \cdot a \cdot b}} $$ (9) 式中：
$ A $ 为丝杠轴的断面面积，mm²；$ E $ 为丝杠轴的杨氏模数；$ L $ 为安装间距，mm；$ a $ 与$ b $ 为螺母距两端面的距离，mm。当在$ a = b = \dfrac{L}{2} $ 的位置时，$ {k_S} $ 的值为最小，轴向弹性位移量最大为$$ {k_S} = \frac{{4A \cdot E}}{{1000L}} $$ (10) 螺母的轴向刚性
$ {k_N} $ ，支撑轴承的轴向刚性$ {k_B} $ 与螺母座的轴向刚度$ {k_H} $ 可通过查询THK手册曲线结合机床螺母预压情况、使用情况得出[11-12]。综合以上丝杠螺母的各组成部分刚性，代入式中即可求解丝杠螺母动结合部的轴向刚性。丝杠螺母动结合部的法向刚性计算方法由基于赫兹基础理论法[13]求解，这里不再赘述。
得到丝杆螺母动结合部的刚度，其中：轴向刚度
$ {k_x} = 3.76 \times {10^5} $ N/mm，法向刚度为$ {k_y} = {k_{\textit{z}}} = 2.13 \times {10^6} $ N/mm。2.3 直线导轨副的刚度
建立直线导轨的动力学模型如图5所示。
滚动导轨承受负荷时，钢球、滑块在容许负荷范围产生弹性变形[14]。这时的形变量与外加负荷之比率就是刚性。由于滚动导轨的径向间隙影响刚性，因此机床选THK品牌标准的C0间隙即中预压，为适用于加工中心的预压种类。根据选定的间隙标准、滑块形式，查询手册可以得到直线导轨的刚性
$ {k_x}' $ ，直线导轨的法向刚度求解方法与滚珠丝杠相同。得到直线导轨动结合部的刚度为：轴向刚度
$ {k_x} = 2.53 \times {10^5} $ N/mm，法向刚度为$ {k_y} = {k_{\textit{z}}} = 4.78 \times {10^6} $ N/mm。将丝杠螺母刚度与直线导轨刚度代入式（5）即可得到专机考虑了动结合部刚度的动力学方程。
3. 内腔加工专机动力学特性分析
3.1 应用动态子结构法求解固有频率
多自由度系统的固有频率是通过求解系统的无阻尼自由振动方程得到的[15]。多自由度系统无阻尼自由振动的运动方程可由式（7）中令
$ {C_{20 \times 20}} = \left[ 0 \right] $ 、$ {Q_{20 \times 1}} = \left[ 0 \right] $ ，得$$ \left[ {\boldsymbol{M}} \right]\left\{ {\ddot q} \right\} + \left[ {\boldsymbol{K}} \right]\left\{ q \right\} = 0 $$ (11) 式中：
$ \left[ {\boldsymbol{M}} \right] $ 为专机的质量矩阵；$ \left[ {\boldsymbol{K}} \right] $ 为专机刚度矩阵，$ q $ 为系统的广义坐标。设方程的解为
$$ \left\{ q \right\} = \left\{ A \right\}{e^{i{\omega _n}t}} $$ (12) 式中：
$ \left\{ A \right\} $ 为系统自由振动时的振幅向量；$ {\omega _n} $ 为振动频率，$ t $ 为时间。将由式（12）及其对时间的二阶导数代入式（11）中，消去因子
$ {e^{i{\omega _n}t}} $ 后得$$ ({\boldsymbol{K}} - \omega _n^2{\boldsymbol{M}})\{ A\} = 0 $$ (13) 求解式（13）常称为特征值问题。要得到振动解（非零解），必须
$ \left\{ A \right\} $ 的系数行列式等于零，即$$ {\rm{det}}({\boldsymbol{K}} - \omega _n^2{\boldsymbol{M}}) = 0 $$ (14) 式（14）称为系统的特征方程，求解即可得到
$ n $ 个根：$ \omega _1^2 $ ，$ \omega _2^2 $ ，$ \omega _3^2 $ ，…，$ \omega _n^2 $ 称为特征值，将特征值分别开方后求得$ n $ 个$ {\omega _r} $ （$ r = 1,2, \cdots ,n $ ），称为系统的$ n $ 个固有频率，按大小顺序排列：$ {\omega _1} \leqslant {\omega _2} \leqslant \cdots \leqslant {\omega _n} $ ，分别为一阶固有频率、二阶固有频率、…、$ n $ 阶固有频率。式（14）中刚度矩阵
${\boldsymbol{K}}$ 各元素均已求得，质量矩阵${\boldsymbol{M}}$ 中子结构的质量及转动惯量由三维建模后设置密度，经CAD软件求解得到，如表2所示。表 2 机床子结构的质量及转动惯量子结构 质量/kg 转动惯量/(kg·m²) 床身 $ {m_1} = 6.316 \times {10^3} $ $ {I_{{\textit{z}}1}} = 6.37 \times {10^3} $、$ {I_{x1}} = 4.43 \times {10^3} $ 立柱 $ {m_2} = 834 $ $ {I_{x2}} = 77 $、$ {I_{{\textit{z}}2}} = 180 $ 滑枕 $ {m_3} = 275 $ $ {I_{y3}} = 46 $、$ {I_{{\textit{z}}3}} = 13 $ 主轴 $ {m_4} = 213 $ $ {I_{x4}} = 109 $、$ {I_{{\textit{z}}4}} = 99 $ 工作台 $ {m_5} = 2.052 \times {10^3} $ $ {I_{x5}} = 191 $、$ {I_{{\textit{z}}5}} = 238 $ 计算式（14）即可得到基于动态子结构法的专机前四阶固有频率为：35.1 Hz、36.9 Hz、76.4 Hz和73.7 Hz。
3.2 应用刚性简化的有限元分析求解固有频率
在实际机械设备设计中，通常不会考虑动结合部的刚性问题，将所有结合部都按照刚性固定方式处理[16]。采用刚性简化的方式，基于有限元软件求解专机的固有频率。
首先对专机进行合理简化，去除模型中的螺纹凸缘等特征，其中电机、转台、加长杆角度头、刀具和卡盘等部件采用实体建模、抽壳特征处理。设置材料属性如表3所示。
表 3 有限元仿真部件材料属性设置部件 材料 密度/(kg/m3) 弹性模量/GPa 泊松比 床身、立柱 HT200 7 340 120 0.25 滑枕、主轴、
工作台Q235 7 850 200 0.3 将三维模型导入有限元软件中建立有限元仿真模型，设置所有结合部刚性接触，划分网格如图6所示。
模态分析求解的边界条件设置为床身底面与垫铁的接触面固定。求得基于动态子结构法的专机前四阶固有频率为：40.2 Hz、42.9 Hz、90.6 Hz和91.7 Hz。
3.3 应用设置动结合部刚度的有限元法求解固有频率
为验证建立考虑动结合部刚性的子结构动力学模型正确性，利用有限元软件设置动结合部刚性，为简化计算处理，将各子结构接触面处刚性固连，在滚珠丝杠副及直线导轨处设置动结合部，求解专机固有频率。模态分析的材料设置、网格划分及边界条件与刚性简化计算一致。
设置专机的动结合部（X轴的工作台与床身的丝杠螺母副、直线导轨副结合部；Z轴的立柱与床身的丝杠螺母副、直线导轨副结合部；Y轴的立柱与滑枕的丝杠螺母副、直线导轨副结合部）的刚度进行模态分析。代入动结合部法向与切向的刚度，如图7a所示，在各丝杠螺母副的结合处设置3个弹簧单元模拟结合部刚度。
如图7b所示，以床身与工作台的动结合部为例，以床身丝杠螺母与工作台、床身导轨与工作台滑块的接触面法向建立X方向弹簧单元，以右手系建立X方向、Z方向弹簧单元。设置所有弹簧单元的刚度，完成动结合部刚度定义。
由此可以求解专机含动结合部刚度的固有频率为34.2 Hz，38.2 Hz、78.6 Hz和79.4 Hz。
专机的前四阶振型如图8所示。
从图8可知专机的机械动力特性，一阶振型为立柱与主轴头绕Y轴的偏摆，二阶振型为立柱与主轴头绕Z轴的转动，三阶振型为立柱与床身绕Z轴的扭转，四阶振型变化为主轴头滑枕绕X轴的扭转。根据固有频率与振型分析研究结果表明：（1）主轴部件的动力响应对刀尖点的动力响应贡献最大，应选用高刚性的镗铣动力头、加长杆角度头，并加强主轴部件与A轴回转台的刚性连接。（2）立柱组件与床身组件的动结合处中导轨滑块的跨距较小，致使二阶振型立柱组件绕Z轴转动幅度较大，应在该处选用高刚性大尺寸导轨滑块，并适当增大跨距。针对上述结构修改后，较大幅度提高了整机的动态性能。
3种求解方法得到的专机前四阶固有频率为如表4所示。
表 4 含动结合部刚度求解的专机固有频率固有频率/Hz 一阶$ {f_1} $ 二阶$ {f_2} $ 三阶$ {f_3} $ 四阶$ {f_4} $ 理论方法（基准） 35.1 36.9 76.4 73.7 刚性接触 40.2 42.9 90.6 91.7 设置动结合部刚性的
有限元法34.2 38.2 78.6 79.4 刚性接触与基准
相比误差14.5% 16.3% 18.6% 24.4% 设置动结合部刚性的有限元法与基准相比误差 2.6% 3.5% 2.9% 7.7% 可以看出，结合部动态参数的描述是影响机床整机动力学性能的重要因素，特别是对于高阶固有频率特性计算，影响更为明显。因此研究动结合部的动态刚度，对于求解机床动力学特性有着重要意义。由表4的数据可知，本文建立的考虑结合部动态刚度的子结构法可以正确求解专机动力学模型，求解固有频率与振型，建立的动力学方程可以应用到数字孪生动力学仿真应用中。
3.4 内腔加工专机的谐响应分析
模态分析只能求得设备本身的固有频率与振型，而谐响应分析可求得机床在不同频率简谐载荷作用下的位移响应，对于分析加工过程中刀具振动对专机加工精度的影响具有重要意义。
在刀具刀尖部添加大小为
$ F = {F_0}\sin \omega t $ ，$ {F_0} = 200 $ N，方向垂直于动平台的应力幅值；简谐力的频率变化范围为0~200 Hz（专机工作中刀具与工件振动激励的频率范围）；设置载荷子步数为50；其余边界条件按计算的动结合部刚性设置（理论方法）。得到刀尖点处X、Y、Z的3个方向上的位移响应结果（振幅与频率的关系）如图9所示。在激励频率为38 Hz、79 Hz时响应较大，并且刀尖点X、Y、Z三个方向均产生了较大响应，在其他频率处没有发生共振，并且38 Hz与机床前两阶固有频率吻合，79 Hz与第三阶固有频率相吻合，为该机床的敏感频率，因此在专机在机加工时应避开上述振动频率，从而保证加工精度和效率。
4. 结语
本文基于动态子结构法对立式四轴舱体内腔加工专机进行机械动力学建模，求解了专机动结合部滚珠丝杠副和直线导轨副的动态刚度特性。结果显示建立的含动态刚度的动力学模型能够有效镜像专机动力学特性，可以准确表达专机的真实动态性能，可应用于专机的动力分析与数字孪生虚拟样机研究。
-
表 1 专机刚性结合部刚度
刚性结合部 线刚度/(N/m)、角刚度/(N·m/rad) 床身与地面 $ {K}_{y}^{\text{①}}=7.642\times {10}^{8} $、$ {K}_{\theta {\textit{z}}}^{\text{①}}=3.629\times {10}^{7} $、$ {K}_{\theta x}^{\text{①}}=5.473\times {10}^{7} $ 滑枕与主轴 $ {K}_{x}^{\text{④}}=7.218\times {10}^{7} $、$ {K}_{y}^{\text{④}}=2.512\times {10}^{8} $、$ {K}_{{\textit{z}}}^{\text{④}}=2.512\times {10}^{8} $
$ {K}_{\theta x}^{\text{④}}=4.538\times {10}^{7} $、$ {K}_{\theta {\textit{z}}}^{\text{④}}=4.341\times {10}^{7} $
下载: 导出CSV
表 2 机床子结构的质量及转动惯量
子结构 质量/kg 转动惯量/(kg·m²) 床身 $ {m_1} = 6.316 \times {10^3} $ $ {I_{{\textit{z}}1}} = 6.37 \times {10^3} $、$ {I_{x1}} = 4.43 \times {10^3} $ 立柱 $ {m_2} = 834 $ $ {I_{x2}} = 77 $、$ {I_{{\textit{z}}2}} = 180 $ 滑枕 $ {m_3} = 275 $ $ {I_{y3}} = 46 $、$ {I_{{\textit{z}}3}} = 13 $ 主轴 $ {m_4} = 213 $ $ {I_{x4}} = 109 $、$ {I_{{\textit{z}}4}} = 99 $ 工作台 $ {m_5} = 2.052 \times {10^3} $ $ {I_{x5}} = 191 $、$ {I_{{\textit{z}}5}} = 238 $
下载: 导出CSV
表 3 有限元仿真部件材料属性设置
部件 材料 密度/(kg/m3) 弹性模量/GPa 泊松比 床身、立柱 HT200 7 340 120 0.25 滑枕、主轴、
工作台Q235 7 850 200 0.3
下载: 导出CSV
表 4 含动结合部刚度求解的专机固有频率
固有频率/Hz 一阶$ {f_1} $ 二阶$ {f_2} $ 三阶$ {f_3} $ 四阶$ {f_4} $ 理论方法（基准） 35.1 36.9 76.4 73.7 刚性接触 40.2 42.9 90.6 91.7 设置动结合部刚性的
有限元法34.2 38.2 78.6 79.4 刚性接触与基准
相比误差14.5% 16.3% 18.6% 24.4% 设置动结合部刚性的有限元法与基准相比误差 2.6% 3.5% 2.9% 7.7%
下载: 导出CSV
-
[1] 刘宏. 锥体内腔多轴数控铣削加工工艺技术[J]. 制造业自动化, 2019, 41(9): 144-148. DOI: 10.3969/j.issn.1009-0134.2019.09.033 [2] 刘彪. 薄壁件铣削稳定性与动态加工误差研究[D]. 北京: 北京理工大学, 2015. [3] Zhang G, Huang Y, Shi W. Predicting dynamic behaviours of a whole machine tool structure based on computer-aided engineering[J]. International Journal of Machine Tools and Manufacture, 2003(40): 699-706.
[4] 张广鹏, 史文浩, 黄玉美. 机床导轨结合部的动态特性解析方法及其应用[J]. 机械工程学报, 2013, 38(10): 114-117. [5] Liu H, Sun Y, Liu Z. Simulation and experiment of dynamic properties of joint surfaces based on fractal theory[J]. Shock and Vibration, 2015(9): 21-29.
[6] 王书亭, 李杰, 刘涛, 等. 机械固定结合面刚度特性建模[J]. 华中科技大学学报:自然科学版, 2011, 39(8): 26-31. [7] 周石恩. 基于数字孪生的复杂产品装配建模与精度分析方法[D]. 杭州: 浙江大学, 2019. [8] 廖伯瑜, 周新民, 尹志宏. 现代机械动力学及工程应用[M]. 北京: 机械工业出版社, 2004. [9] 毛宽民, 李斌, 雷声. 机床结合部动力学建模及应用[M]. 武汉: 武汉理工大学出版社, 2018. [10] 李景奎, 石宏, 郭建烨. 基于有限元机床导轨结合面参数的特性分析[J]. 航空制造技术, 2013(5): 92-94. [11] Huang Y M, Fu W P, Dong L X. Research on thenormal dynamic characteristic parameters of jointsurface[J]. Journal of Mechanical Engineering, 1993, 22(3): 54-55.
[12] Wu J S S, Chang J C, Hung J P. The effect of contactinterface on dynamic characteristics of compositestructures[J]. Mathematics and Computers in Simulation, 2007, 41(6): 13-19.
[13] 蒋书运, 祝书龙. 带滚珠丝杠副的直线导轨结合部动态刚度特性[J]. 机械工程学报, 2010, 46(1): 92-93. [14] 屈重年, 伍良生, 肖毅川, 等. 机床导轨技术研究综述[J]. 制造技术与机床, 2012(1): 30-36. DOI: 10.3969/j.issn.1005-2402.2012.01.010 [15] 张景. 动力学系统建模[M]. 北京: 国防工业出版社, 2000. [16] 王玉宝. 五轴联动数加工中心结构动特性分析与优化[D]. 天津: 天津大学, 2010. -
期刊类型引用(1)
1. 叶磊，成群林，雷怡凡，平昊，李凌霄，吕盾，周井文. 数控卧式内腔加工专机联合仿真建模与分析. 工具技术. 2023(12): 75-81 . 
百度学术
其他类型引用(0)
计量
- 文章访问数: 113
- HTML全文浏览量: 62
- PDF下载量: 194
- 被引次数: 1
