Dynamic error of multiaxis machine tools considering position dependent structural dynamics and axis coupling inertial forces

Abstract

During the working process of high-speed multiaxis machine tools, inertial forces can cause vibration and deformation of mechanical structure, which lead to the dynamic error of tool center point (TCP) relative to worktable and can adversely affect the machining performance. Considering the varying feed positions and accelerations during machining, a parameter-varying multi-rigid-body dynamic model of a 3-axis gantry machine tool is proposed. This model represents the position dependent structural dynamics and inertial forces, which can simulate the dynamic error of TCP relative to worktable within the entire workspace. The results show that the dynamic error in one direction is affected by the feed motions of multiple feed axes. The magnitudes of the dynamic error significantly vary with the position of Z-axis. And the dynamic errors in Y- and Z-direction show different varying trends. Then the theoretical model is used to discuss the dynamic error and position dependency. The expressions of TCP dynamic response and inertial forces reveal the reason why the dynamic errors in Y- and Z-direction show different varying trends.

Introduction

The increasing feed accelerations of feed axes enhance the inertial forces effect during working process of high-speed machine tools. The inertial forces can cause vibration and deformation of mechanical structure, which lead to dynamic error of TCP relative to worktable (DETW).^1,2 For multi-axis machine tools with kinematic coupling, axis coupling effect brings the coupling forces among multiple axes and varying structural dynamics.^3–6 These effects result in the complex characteristics of DETW during feed motions, which make the analysis and control of dynamic error of high-speed multiaxis machine tools more challenging.

The dynamic error induced by inertial forces during accelerating and decelerating stage has received great interest in the process of design and analysis of high-speed and high precision machine tools. In recent years, the related content has been added in international standards. Slocum, ⁷ Heisel and Gringel, ⁸ and Reinhart and Weissenberger ⁹ indicated that inertial forces can cause vibration and deformation of mechanical structure. Brecher et al. ¹⁰ defined the distance of the gravity center to the incidence point of drive force as Steiner Distance which is one of the crucial factors that influence acceleration forces induced error. This kind of error was considered as an important error source when designing high precision machine tools. Bringmann and Maglie, ¹¹ Nguyen et al., ¹² and Thoma et al. ¹³ proposed in-talk and cross-talk concepts to describe the dynamic deviation at TCP of machine tools with inertial effect. And Inertial Cross-Talk was presented in ISO 230-8:2009 to evaluate the vibrations of machine tools during accelerations and decelerations.

The dynamic error cannot be directly measured by rotary encoder and linear encoder,^14,15 which result in that the traditional feedback controller is incapable of controlling DETW. Therefore, it is necessary to develop a theoretical model which can estimate DETW. Parenti et al. ¹⁶ analyzed the effect of inertial deformation on machine tool accuracy by using finite element model. A simplified single degree of freedom (DOF) model was used to design a feed-forward compensation filter. Matsubara et al. ¹⁷ and Nagaoka and Matsubara ¹⁸ developed a transfer function model of machine tool structure using modal superposition method. A model-reference feed-forward controller was designed to control the load position of machine tool axes. Altintas and Khoshdarregi ¹⁹ characterized the structural dynamic mode of feed drive by using a second order transfer function model. The integrated vibration avoidance and contouring error compensation is developed using input shaping filter and pre-compensation. Uhlmann et al. ²⁰ presented a reduced finite element model of machine tool structure. Denkena et al. ²¹ built a multi-rigid-body model of 2-axis feed system. They both proposed a Kalman filter containing the model of mechanical structure to estimate the dynamic dislocation and velocity of TCP. The Kalman filter was applied in the cascaded control of machine tool feed axes. Ansoategui et al.^22,23 characterized the dynamics of feed drive and structure using 3-DOF lumped model. This model was used to estimate the maximum jerk of feed axes. Huang et al. ²⁴ constructed a 4-DOF dynamic model of a single axis feed drive system to predict the tracking error and the elastic deformation. This model included the flexibility of the transmission and the machine bed. However, all the theoretical models in the above researches are linear time invariant (LTI) model. This kind of model cannot correctly represent the varying structural dynamics and inertial forces under axis coupling effect.

In recent years, the researchers have taken the axis coupling effect seriously for the demands of high-speed and high precision machining. The coupling forces were considered in the dynamic analysis of multiaxis machine tools, including inertial force,^4,5 gravity force, ⁶ and friction force.^3,25 In high-speed working condition, the inertial force effect is dominant because of the high acceleration. Many researchers have focused on the position dependent dynamics of multi-axis machine tool structure,^26–29 including linear axis ³⁰ and rotary axis. ³¹ These researches showed that the multiaxis configuration varies with the feed positions. The changing configuration leads to the redistribution of the structural inertia, which causes the variations of structural dynamics and inertial forces. In order to deal with the effect of the varying dynamics and inertial forces, Hiramoto et al. ³² proposed drive at the center of gravity principle to reduce inertial vibration. The changing center of gravity induced by feed motion was considered. Steinlin et al. ³³ built a proportional model of cross-talk error and developed off-line compensation. The model was derived by linear fit of measured data. The proportion factor of the model was acceleration dependent and position dependent. In summary, limited open literature to date has reported modeling and discussion of DETW and its variation with feed positions. The existing models are too simple to characterize the structural dynamics and axis coupling inertial forces of a multiaxis machine tool. Therefore, it is necessary to develop a new theoretical model to study DETW and its position dependency. Detailed discussion of the position dependency of DETW is beneficial for the design and control of the dynamic error of high-speed multiaxis machine tools.

In this paper, a vertical gantry milling machine which is a typical 3-linear-axis machine tool with kinematic coupling is studied. Considering the varying feed positions and feed accelerations, the theoretical model is parameter-varying multi-rigid-body dynamic model. Dynamic response at TCP in the orthogonal direction of the feed direction is used to evaluate the DETW during accelerating stage. The variation of DETW with the feed positions is analyzed. By using the theoretical model, the effect of structural dynamics on the dynamic error and the mechanism of the position dependency of the dynamic error are discussed.

Dynamic modeling and solution process

Parameter-varying dynamic model considering the variation of feed positions

A gantry-type vertical milling machine which is a typical 3-linear-axis structure with kinematic coupling is studied. The equivalent dynamic model of the machine tool structure is multi-rigid-body model ³⁴ as shown in Figure 1. The structural components are considered as rigid bodies. The machine tool structure is composed of four rigid bodies. The machine bed is written as ${\mathrm{B}}_0$ . The crossbeam of Y-axis is ${\mathrm{B}}_1$ . The slider of X-axis is ${\mathrm{B}}_2$ . And ${\mathrm{B}}_3$ is the headstock of Z-axis and spindle. The flexibilities of kinematic joints are characterized by linear spring and damping elements (LSDE). The axial flexibility of lead screw transmission is modeled as 1-DOF LSDE. 2-DOF LSDE is used to represent the vertical and lateral flexibilities of linear guideway. And the machine support is modeled as 3-DOF LSDE.

Figure 1.

Equivalent dynamic model of machine tool structure.

The geometric parameters and feed motions of the machine tool are shown in Figure 2. Let the origin of the global reference frame RXYZ be located in the TCP position when all feed axes are at reference position. $L_{\mathrm{ij}\mathrm{d}}$ is the geometric parameters of the structure, which refers to the relative position from the center of gravity (COG) of the rigid body ${\mathrm{B}}_i$ to the COG of the rigid body ${\mathrm{B}}_j$ in d-direction. The feed motions of the linear feed axes are expressed by $x_{\mathrm{g}}(t)$ , $y_{\mathrm{g}}(t)$ , and $z_{\mathrm{g}}(t)$ . The varying offsets between the structural components can be derived by:

$$\begin{array}{c}{x}_{\mathrm{o}}\left(t\right)=x_{\mathrm{g}}\left(t\right)+L_{21\mathrm{x}}\hfill \\ \multicolumn{1}{c}{y_{\mathrm{o}}\left(t\right)=y_{\mathrm{g}}\left(t\right)+L_{10\mathrm{y}}}\\ \multicolumn{1}{c}{z_{\mathrm{o}}\left(t\right)=z_{\mathrm{g}}\left(t\right)+L_{32\mathrm{z}}}\end{array}$$

(1)

Figure 2.

Geometric parameters and feed motions of the machine tool.

The structural components is modeled as rigid body with spatial 6-DOF. The DOFs of a structural component are given by:

$${\mathbf{q}}_i={\left\{\begin{array}{cccccc}{u}_i\hfill & v_i\hfill & w_i\hfill & {\alpha }_i\hfill & {\beta }_i\hfill & {\gamma }_i\hfill \end{array}\right\}}^{\mathrm{T}}$$

(2)

where $u_i$ , $v_i$ , and $w_i$ represent the translational motions of the i th body along X-, Y-, and Z-axis, respectively. And ${\alpha }_i$ , ${\beta }_i$ , and ${\gamma }_i$ are the rotational motions of the i th body around X-, Y-, and Z-axis, respectively.

The structure of the machine tool is typically treelike multi-body structure. Using the relative coordinate system, the DOFs ${\mathbf{q}}_i$ are defined as the movements of rigid body ${\mathrm{B}}_i$ relative to its lower body ${\mathrm{B}}_{i-1}$ . The DOF of the four rigid bodies are chosen as the generalized coordinates of the entire dynamic system. So the generalized coordinates vector can be written as $\mathbf{q}={\left\{\begin{array}{cccc}{\mathbf{q}}_0\hfill & {\mathbf{q}}_1\hfill & {\mathbf{q}}_2\hfill & {\mathbf{q}}_3\hfill \end{array}\right\}}^{\mathrm{T}}$ .

The $4\times 4$ homogeneous transformation matrix ${{}^{i-1}\mathbf{T}}_i$ is used to describe the orientation and position of rigid body ${\mathrm{B}}_i$ relative to its lower body ${\mathrm{B}}_{i-1}$ . ³⁵ The transformation matrix ${{}^{i-1}\mathbf{T}}_i$ can be given by:

$${{}^{i-1}\mathbf{T}}_i={\mathbf{T}}_{\mathrm{l}i}{\mathbf{T}}_{\mathrm{q}i}$$

(3)

${\mathbf{T}}_{\mathrm{q}i}$ represents the homogeneous transformation matrix corresponding to the DOF ${\mathbf{q}}_i$ of rigid body ${\mathrm{B}}_i$ . Because ${\mathbf{q}}_i$ represent the elasic deformation, ${\mathbf{T}}_{\mathrm{q}i}$ can be simplified as:

$${\mathbf{T}}_{\mathrm{q}i}=\left[\begin{array}{cccc}1\hfill & -{\gamma }_i\hfill & {\beta }_i\hfill & u_i\hfill \\ \multicolumn{1}{c}{{\gamma }_i}& 1\hfill & -{\alpha }_i\hfill & v_i\hfill \\ \multicolumn{1}{c}{-{\beta }_i}& {\alpha }_i\hfill & 1\hfill & w_i\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]$$

(4)

${\mathbf{T}}_{\mathrm{l}i}$ is the homogeneous transformation matrix corresponding to the geometric parameters shown in Figure 2. If the relative position from the COG of ${\mathrm{B}}_i$ to the COG of ${\mathrm{B}}_{i-1}$ is $(L_{i\mathrm{x}},L_{i\mathrm{y}},L_{i\mathrm{z}})$ , ${\mathbf{T}}_{\mathrm{l}i}$ is:

$${\mathbf{T}}_{\mathrm{p}i}=\left[\begin{array}{cccc}1\hfill & 0\hfill & 0\hfill & L_{i\mathrm{x}}\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & L_{i\mathrm{y}}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & L_{i\mathrm{z}}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]$$

(5)

The transformation matrix of rigid body ${\mathrm{B}}_i$ with respect to the global reference frame RXYZ can be derived by:

$${{}^{\mathrm{R}}\mathbf{T}}_i={{}^{\mathrm{R}}\mathbf{T}}_0{{}^{\mathrm{0}}\mathbf{T}}_1\cdots {{}^{i-1}\mathbf{T}}_i$$

(6)

${{}^{\mathrm{R}}\mathbf{T}}_i$ can be divided into four parts as below:

$${{}^{\mathrm{R}}\mathbf{T}}_i=\left[\begin{array}{cc}\mathbf{Ro}{\mathbf{t}}_{3\times 3}\hfill & \mathbf{Po}{\mathbf{s}}_{3\times 1}\hfill \\ \multicolumn{1}{c}{\mathbf{Per}{\mathbf{s}}_{1\times 3}}& \mathrm{Scale}\hfill \end{array}\right]$$

(7)

where $\mathbf{Ro}{\mathbf{t}}_{3\times 3}$ is the rotational matrix, $\mathbf{Po}{\mathbf{s}}_{3\times 1}$ is the translational vector, $\mathbf{Per}{\mathbf{s}}_{1\times 3}$ represents projective transformation, and Scale is the scaling factor. The translational movements of the i th rigid body ${\mathrm{B}}_i$ with respect to the global reference frame RXYZ can be obtained from $\mathbf{Po}{\mathbf{s}}_{3\times 1}$ . The rotational movements of ${\mathrm{B}}_i$ with respect to RXYZ are the sum of the rotational movements of ${\mathrm{B}}_i$ and its lower bodies.

The deformations of LSDE can be derived by the position vector of the LSDE relative to the COG of structural components and the DOFs of structural components using homogeneous transformation. The position vector of LSDE relative to the COG of structural components is ${\mathbf{p}}_{\mathrm{sd}}={\left\{\begin{array}{ccc}{p}_{\mathrm{x}}\hfill & p_{\mathrm{y}}\hfill & p_{\mathrm{z}}\hfill \end{array}\right\}}^{\mathrm{T}}$ . The actual position ${\mathbf{q}}_{\mathrm{sde}}$ of the LSDE due to the motion of structural components can be calculated as:

$${\mathbf{q}}_{\mathrm{sde}}=\left[\begin{array}{cccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & p_z\hfill & -p_y\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & -p_z\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & p_y\hfill & -p_x\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]{\mathbf{q}}_i$$

(8)

Therefore, the deformation of the LSDE can be given by:

$$\begin{array}{c}{d}_{\mathrm{x}}=u_i+p_z{\beta }_i-p_y{\gamma }_i\hfill \\ d_{\mathrm{y}}=v_i-p_z{\alpha }_i+p_x{\gamma }_i\hfill \\ d_{\mathrm{z}}=w_i+p_y{\alpha }_i-p_x{\beta }_i\hfill \end{array}$$

(9)

The movements of all rigid bodies and the deformations of all LSDE have been derived. The kinetic energy $E_{\mathrm{k}}$ , the potential energy $E_{\mathrm{p}}$ , and the dissipative energy $E_{\mathrm{d}}$ of the entire dynamic system can be expressed. The dynamic equations can be derived by using Lagrange’s equations with dissipation forces ³⁶ which is written as:

$$\frac{d}{\mathrm{dt}}\left(\frac{\partial E_{\mathrm{k}}}{\partial {\stackrel{·}{q}}_j}\right)-\frac{\partial E_{\mathrm{k}}}{\partial q_j}+\frac{\partial E_{\mathrm{p}}}{\partial q_j}+\frac{\partial E_{\mathrm{d}}}{\partial {\stackrel{·}{q}}_j}=Q_j$$

(10)

here j is the number of the generalized coordinates. $Q_j$ is the external load.

The dynamic stiffening terms in the dynamic equations are neglected. The terms corresponding to the feed accelerations ${\stackrel{··}{x}}_{\mathrm{g}}$ , ${\stackrel{··}{y}}_{\mathrm{g}}$ , and ${\stackrel{··}{z}}_{\mathrm{g}}$ are collected. These collected terms form the inertial forces vector ${\mathbf{F}}_{\mathrm{i}\mathrm{nertia}}$ . All the terms belong to ${\mathbf{F}}_{\mathrm{i}\mathrm{nertia}}$ are moved to the right side of the dynamic equations and considered as the input forces. The dynamic equations can be expressed as:

$$\begin{array}{c}\mathbf{M}(x_{\mathrm{g}},y_{\mathrm{g}},z_{\mathrm{g}})\stackrel{··}{\mathbf{q}}\left(t\right)+\mathbf{C}\left(z_{\mathrm{g}}\right)\stackrel{·}{\mathbf{q}}\left(t\right)+\mathbf{K}\left(z_{\mathrm{g}}\right)\mathbf{q}\left(t\right)\hfill \\ \multicolumn{1}{c}{={\mathbf{F}}_{\mathrm{inertia}}(x_{\mathrm{g}},y_{\mathrm{g}},z_{\mathrm{g}},{\stackrel{··}{x}}_{\mathrm{g}},{\stackrel{··}{y}}_{\mathrm{g}},{\stackrel{··}{z}}_{\mathrm{g}})}\end{array}$$

(11)

where $\mathbf{M}$ , $\mathbf{C}$ , and $\mathbf{K}$ are mass matrix, damping matrix, and stiffness matrix, respectively.

It can be seen that the derived dynamic equations are parameter-varying equations. The dynamic matrices vary with the feed positions, which reflect the position dependent structural dynamics. The input forces are a function of the feed positions and accelerations, which characterize the axis coupling inertial forces determined by the feed positions and accelerations.

Solution process of dynamic error of TCP relative to worktable

The DETW during the accelerating stage of one feed axis is simulated in this study. Let the feedrate is 5000 mm/min and the acceleration is 2  $\mathrm{m}/{\mathrm{s}}^2$ . The ideal displacement during the accelerating stage is 20 mm, which is compared to the motion range of the gantry machine tool. If the simulation is carried out in a small motion range, a simplified linear time-invariant dynamic system can be obtained by neglecting the variation of feed position during the accelerating stage.

In order to calculate the feed acceleration, a simplified servo model is proposed, as shown in Figure 3. The mechanical system of feed drive is modeled using 2-DOF lumped model. Here $j_{\mathrm{m}}$ is motor inertia, $j_{\mathrm{l}}$ is load inertia, c is damping, k is stiffness, $\zeta$ is damping ratio, ${\omega }_n$ is undamped natural frequency. The position controller is P controller. The position loop gain is $k_{\mathrm{p}}$ . PI controller is used for speed control. Here $k_{\mathrm{v}}$ is speed loop gain, $T_{\mathrm{v}}$ is speed loop time constant.

Figure 3.

Simplified servo model.

A multi-input-multi-output (MIMO) state space model of the system is proposed to calculate DETW. The input signals of the state space model are the feed accelerations of the three feed axes. The output signals of the state space model are DETW in X-, Y- and Z-direction. The state-dynamic matrix $\mathbf{A}$ is given by:

$$\mathbf{A}=\left[\begin{array}{cc}0\hfill & \mathbf{E}\hfill \\ \multicolumn{1}{c}{-{\mathbf{M}}^{-1}\mathbf{K}}& -{\mathbf{M}}^{-1}\mathbf{C}\hfill \end{array}\right]$$

(12)

The input matrix $\mathbf{B}$ is:

$$\mathbf{B}=\left(\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{{\mathbf{M}}^{-1}}\end{array}\right)$$

(13)

The input vector $\mathbf{u}$ can be derived from the inertial forces ${\mathbf{F}}_{\mathrm{i}\mathrm{nertia}}$ in equation (11). The input vector $\mathbf{u}$ is expressed as:

$$\mathbf{u}=\left(\begin{array}{ccc}{\mathbf{F}}_{\mathrm{x}}\hfill & {\mathbf{F}}_{\mathrm{y}}\hfill & {\mathbf{F}}_{\mathrm{z}}\hfill \end{array}\right)$$

(14)

${\mathbf{F}}_{\mathrm{x}}$ , ${\mathbf{F}}_{\mathrm{y}}$ , and ${\mathbf{F}}_{\mathrm{z}}$ can be calculated by equation (15). Detail expressions of the input vector can be found in Appendix.

$$\begin{array}{c}{\mathbf{F}}_{\mathrm{x}}={\mathbf{F}}_{\mathrm{inertia}}{|}_{{\stackrel{··}{x}}_{\mathrm{g}}=1,{\stackrel{··}{y}}_{\mathrm{g}}=0,{\stackrel{··}{z}}_{\mathrm{g}}=0}\hfill \\ \multicolumn{1}{c}{{\mathbf{F}}_{\mathrm{y}}={\mathbf{F}}_{\mathrm{inertia}}{|}_{{\stackrel{··}{x}}_{\mathrm{g}}=0,{\stackrel{··}{y}}_{\mathrm{g}}=1,{\stackrel{··}{z}}_{\mathrm{g}}=0}}\\ \multicolumn{1}{c}{{\mathbf{F}}_{\mathrm{z}}={\mathbf{F}}_{\mathrm{inertia}}{|}_{{\stackrel{··}{x}}_{\mathrm{g}}=0,{\stackrel{··}{y}}_{\mathrm{g}}=0,{\stackrel{··}{z}}_{\mathrm{g}}=1}}\end{array}$$

(15)

The transformation matrix of TCP with respect to the global reference frame RXYZ is ${{}^{\mathrm{R}}\mathbf{T}}_{\mathrm{tcp}}$ . The transformation matrix of worktable with respect to RXYZ is ${{}^{\mathrm{R}}\mathbf{T}}_{\mathrm{wtp}}$ . According to equation (6), ${{}^{\mathrm{R}}\mathbf{T}}_{\mathrm{tcp}}$ and ${{}^{\mathrm{R}}\mathbf{T}}_{\mathrm{wtp}}$ can be given by:

$$\begin{array}{c}{{}^{\mathrm{R}}\mathbf{T}}_{\mathrm{tcp}}={{}^{\mathrm{R}}\mathbf{T}}_0{{}^0\mathbf{T}}_1{{}^1\mathbf{T}}_2{{}^2\mathbf{T}}_3{{}^3\mathbf{T}}_{\mathrm{tcp}}\hfill \\ {{}^{\mathrm{R}}\mathbf{T}}_{\mathrm{wtp}}={{}^{\mathrm{R}}\mathbf{T}}_0{{}^0\mathbf{T}}_{\mathrm{wtp}}\hfill \end{array}$$

(16)

The translational vector ${\mathbf{P}}_{\mathrm{tcp}}$ of ${{}^{\mathrm{R}}\mathbf{T}}_{\mathrm{tcp}}$ represents the position of the TCP in global reference frame. And the translational vector ${\mathbf{P}}_{\mathrm{wtp}}$ of ${{}^{\mathrm{R}}\mathbf{T}}_{\mathrm{wtp}}$ represents the position of the worktable in global reference frame. Here ${\mathbf{P}}_{\mathrm{tcp}}$ and ${\mathbf{P}}_{\mathrm{wtp}}$ are the linear combination of the generalized coordinates $\mathbf{q}$ . The matrix form of ${\mathbf{P}}_{\mathrm{tcp}}$ and ${\mathbf{P}}_{\mathrm{wtp}}$ can be written as:

$$\begin{array}{c}{\mathbf{p}}_{\mathrm{tcp}}={\mathbf{C}}_{\mathrm{tcp}}\mathbf{q}\hfill \\ {\mathbf{p}}_{\mathrm{wtp}}={\mathbf{C}}_{\mathrm{wtp}}\mathbf{q}\hfill \end{array}$$

(17)

The DETW in X-, Y-, and Z-direction can be derived by equation (18). The expressions of $x_{\mathrm{e}}$ , $y_{\mathrm{e}}$ , and $z_{\mathrm{e}}$ are shown in Appendix.

$$\left\{\begin{array}{c}{x}_{\mathrm{e}}\hfill \\ \multicolumn{1}{c}{y_{\mathrm{e}}}\\ \multicolumn{1}{c}{z_{\mathrm{e}}}\end{array}\right\}=({\mathbf{C}}_{\mathrm{tcp}}-{\mathbf{C}}_{\mathrm{wtp}})\mathbf{q}$$

(18)

The output matrix $\mathbf{C}$ is a $3\times 48$ matrix:

$$\mathbf{C}=\left(\begin{array}{cc}{\mathbf{C}}_{\mathrm{tcp}}-{\mathbf{C}}_{\mathrm{wtp}}\hfill & 0_{3\times 24}\hfill \end{array}\right)$$

(19)

The above equations can be stiff ordinary differential equations (ODEs). Therefore, the solution algorithm for solving stiff ODEs was used in this study, such as ode15s and Gear method. The block diagram of the solution process is shown in Figure 4. $x_{\mathrm{c}}$ , $y_{\mathrm{c}}$ , and $z_{\mathrm{c}}$ are the setpoints generated by the interpolator of CNC system. $G_{\mathrm{x}}(s)$ , $G_{\mathrm{y}}(s)$ , and $G_{\mathrm{z}}(s)$ are the transfer functions transformed from the servo model in Figure 3. The output of the transfer functions are feed positions $x_{\mathrm{g}}$ , $y_{\mathrm{g}}$ , and $z_{\mathrm{g}}$ . Matrices $\mathbf{A}$ , $\mathbf{B}$ , and $\mathbf{C}$ are calculated according to the feed positions. The feed accelerations are obtained by differentiating the feed positions. The input vector $\mathbf{u}$ is calculated according to the feed positions and feed accelerations. Then the DETW can be calculated using the MIMO state space model. Here $x_{\mathrm{e}}$ , $y_{\mathrm{e}}$ , and $z_{\mathrm{e}}$ refer to the DETW in X-, Y-, and Z-direction, respectively. The actual TCP positions $x_{\mathrm{t}\mathrm{cp}}$ , $y_{\mathrm{t}\mathrm{cp}}$ , and $z_{\mathrm{t}\mathrm{cp}}$ are calculated as the sum of feed positions and DETW.

Figure 4.

Block diagram of the solution process of dynamic error of TCP relative to worktable.

Experimental setup

Etalon LaserTracer MultiTrace system was used for measuring DETW. The experimental setup is shown in Figure 5. The sampling frequency was set to 1000 Hz. The maximum trace velocity of the measuring system was 10  $\mathrm{m}/\min$ . The maximum trace acceleration was 0.5 g. A temperature measuring system was used to eliminate the influence of temperature varying on the measured results. The reflector was 120° Cateye and was mounted on the end of the spindle. Four laser tracers were located on the worktable. Therefore, the measurement signals were dynamic displacement response of tool relative to worktable. The machine tool was working under G01 mode. The testing trajectory was simple straight line. The dynamic response in the orthogonal direction of the feed direction was used to evaluate DETW. At last, the measurement signals were filtered using a low-pass filter to reduce the measurement noise.

Figure 5.

Experimental setup for measurement of dynamic error of TCP relative to worktable.

The feed rate was set to 5000  $\mathrm{m}\mathrm{m}/\min$ . The maximum acceleration was 2  $\mathrm{m}/{\mathrm{s}}^2$ . The feed axis started to accelerate from 0 to 5000  $\mathrm{m}\mathrm{m}/\min$ at 0.1 s. Figure 6 shows the experimental and simulated feed accelerations of X-, Y-, and Z-axis. The results show that the servo model can calculate the feed accelerations with high precision.

Figure 6.

Experimental and simulated feed accelerations of X-, Y-, and Z-axis.

Results and discussion

Dynamic error of TCP relative to worktable

Figure 7 shows theoretical and experimental DETW in the orthogonal directions of feed direction during the accelerating stage of each feed axis. The feed positions were $x_{\mathrm{g}}=-400\mathrm{mm}$ , $y_{\mathrm{g}}=-500\mathrm{mm}$ , and $z_{\mathrm{g}}=-100\mathrm{mm}$ . Let the DETW in Y-direction during the accelerating stage of X-axis be denoted by $e_{\mathrm{y}}^{\mathrm{x}}$ . The subscript refers to the direction of DETW. The superscript refers to the feed axis. The same notation is used for DETW and frequency response hereinafter.

Figure 7.

Dynamic error of TCP relative to worktable: (a) dynamic error during the accelerating stage of X-axis, (b) dynamic error during the accelerating stage of Y-axis, and (c) dynamic error during the accelerating stage of Z-axis.

The DETW in one direction is not only influenced by the feed axis in this direction but also influenced by other feed axes. Figure 7(a) shows that $e_{\mathrm{y}}^{\mathrm{x}}$ and $e_{\mathrm{z}}^{\mathrm{x}}$ are less than 2  $\mu \mathrm{m}$ . Therefore, the feed motion of X-axis has no effect on the DETW in Y- and Z-direction. In Figure 7(b), $e_{\mathrm{x}}^{\mathrm{y}}$ can be neglected. The maximum magnitude of $e_{\mathrm{z}}^{\mathrm{y}}$ is 75  $\mu \mathrm{m}$ . The maximum magnitude of $e_{\mathrm{y}}^{\mathrm{y}}$ is −41  $\mu \mathrm{m}$ . In Figure 7(c), $e_{\mathrm{x}}^{\mathrm{z}}$ can be neglected. The maximum magnitude of $e_{\mathrm{y}}^{\mathrm{z}}$ is 14  $\mu \mathrm{m}$ . The maximum magnitude of $e_{\mathrm{z}}^{\mathrm{z}}$ is 25  $\mu \mathrm{m}$ . In summary, the DETW in Y- and Z-direction are affected by both of the feed motions of Y- and Z-axis. The comparisons of theoretical and experimental results verify the effectiveness of the proposed model.

In order to discuss to influence of structural dynamics on DETW, frequency domain analysis and numerical modal analysis are undertaken. Figure 8 shows the frequency responses of DETW. The excitation is the acceleration of one feed axis. Figure 9 shows the first 9 mode shapes. It is obvious that the second and third modes are the relevant modes to the DETW during the accelerating stage of Y- and Z-axis. The second mode is the most dominant mode for $G_{\mathrm{y}}^{\mathrm{y}}$ and $G_{\mathrm{z}}^{\mathrm{y}}$ . While the third mode is more dominant than the second mode for $G_{\mathrm{y}}^{\mathrm{z}}$ and $G_{\mathrm{z}}^{\mathrm{z}}$ . The mode shapes of the second and third modes correspond to the nodding vibration of the moving structural component according to Figure 9.

Figure 8.

Frequency response of dynamic error: (a) frequency response of dynamic error during the accelerating stage of Y-axis and (b) frequency response of dynamic error during the accelerating stage of Z-axis.

Figure 9.

First 9 mode shapes of the machine tool structure.

The expressions of DETW is written as equations (24)–(26) in Appendix. The DETW in Y-direction is composed of the terms corresponding to $v_1$ , ${\alpha }_1$ , ${\gamma }_1$ , $v_2$ , ${\alpha }_2$ , ${\gamma }_2$ , $v_3$ , ${\alpha }_3$ , and ${\gamma }_3$ . The DETW in Z-direction is composed of the terms corresponding to $w_1$ , ${\alpha }_1$ , ${\beta }_1$ , $w_2$ , ${\alpha }_2$ , ${\beta }_2$ , $w_3$ , ${\alpha }_3$ , and ${\beta }_3$ . Figure 10 shows the compositions of $e_{\mathrm{z}}^{\mathrm{y}}$ and $e_{\mathrm{y}}^{\mathrm{z}}$ . The errors corresponding to each generalized coordinate are divided. It can be found that the errors corresponding to ${\alpha }_1$ and ${\alpha }_2$ are the main error sources. In Figure 10(a), more than 80% of $e_{\mathrm{z}}^{\mathrm{y}}$ consists of the error corresponding to ${\alpha }_1$ . Figure 10(b) shows that the errors corresponding to ${\alpha }_1$ and ${\alpha }_2$ have similar magnitude. The errors corresponding to the other generalized coordinates are much smaller. The compositions of DETW can be used to analyze the weakness of machine tool structure. For the milling machine studied in this paper, the linear guideways of crossbeam and slider are the weak components.

Figure 10.

Composition of the dynamic error: (a) dynamic error in Z-direction during the accelerating stage of Y-axis ( $e_{\mathrm{z}}^{\mathrm{y}}$ ) and (b) dynamic error in Y-direction during the accelerating stage of Z-axis ( $e_{\mathrm{y}}^{\mathrm{z}}$ ).

Position dependency of dynamic error

$e_{\mathrm{z}}^{\mathrm{y}}$ and $e_{\mathrm{y}}^{\mathrm{z}}$ vary with the feed position of Z-axis and show different trends. Let the feed position of Z-axis vary from 0 mm to $-300$  mm. Figure 11 shows the variation of $e_{\mathrm{z}}^{\mathrm{y}}$ with Z-position. And the variation of $e_{\mathrm{y}}^{\mathrm{z}}$ with Z-position is shown in Figure 12. $e_{\mathrm{z}}^{\mathrm{y}}$ shows a downward trend from $z_{\mathrm{g}}=0\mathrm{mm}$ to $z_{\mathrm{g}}=-300\mathrm{mm}$ . The maximum magnitude of $e_{\mathrm{z}}^{\mathrm{y}}$ decreases from $89\mu \mathrm{m}$ to $58\mu \mathrm{m}$ . The relative change of $e_{\mathrm{z}}^{\mathrm{y}}$ is 35%. While $e_{\mathrm{y}}^{\mathrm{z}}$ shows a reverse trend compared with $e_{\mathrm{z}}^{\mathrm{y}}$ . The maximum magnitude of $e_{\mathrm{y}}^{\mathrm{z}}$ increases from $-12\mu \mathrm{m}$ to $-20\mu \mathrm{m}$ . The relative change of $e_{\mathrm{y}}^{\mathrm{z}}$ is up to 67%. The simulated results of $e_{\mathrm{x}}^{\mathrm{x}}$ , $e_{\mathrm{y}}^{\mathrm{y}}$ , and $e_{\mathrm{z}}^{\mathrm{z}}$ at different Z-position are shown in Figure 13. It can be found that the relative change of the DETW in feed direction is smaller than that of the DETW in orthogonal direction. Comparison of theoretical and experimental results of the DETW in orthogonal direction verify that the proposed model can characterize the position dependency of DETW.

Figure 11.

Position dependency of the dynamic error in Z-direction during the accelerating stage of Y-axis.

Figure 12.

Position dependency of the dynamic error in Y-direction during the accelerating stage of Z-axis.

Figure 13.

Position dependency of the dynamic error in feed direction: (a) dynamic error in X-direction during the accelerating stage of X-axis, (b) dynamic error in Y-direction during the accelerating stage of Y-axis, and (c) dynamic error in Z-direction during the accelerating stage of Z-axis.

The variations of inertial forces and structural dynamics with feed positions determine the position dependency of DETW. The inertial forces generated by Y-acceleration vary with Z-position. While the inertial forces generated by Z-acceleration do not vary with Z-position. The inertial forces can be expressed by equation (11) and can be calculated by multiplying ${\mathbf{F}}_{\mathrm{x}}$ , ${\mathbf{F}}_{\mathrm{y}}$ , and ${\mathbf{F}}_{\mathrm{z}}$ by the feed accelerations. The expressions of ${\mathbf{F}}_{\mathrm{x}}$ , ${\mathbf{F}}_{\mathrm{y}}$ , and ${\mathbf{F}}_{\mathrm{z}}$ are given by equations (27), (29), and (31) in Appendix. It can be found that ${\mathbf{F}}_{\mathrm{y}}$ is a function of $z_{\mathrm{g}}$ . While ${\mathbf{F}}_{\mathrm{y}}$ is independent of $z_{\mathrm{g}}$ . Because ${\alpha }_1$ and ${\alpha }_2$ are the most important generalized coordinates for DETW, the coefficients corresponding to ${\alpha }_1$ and ${\alpha }_2$ in ${\mathbf{F}}_{\mathrm{y}}$ and ${\mathbf{F}}_{\mathrm{z}}$ can be written as equation (20). From these equations, it is clear that there is a positive correlation between the inertial forces generated by Y-acceleration and Z-position.

$$\begin{array}{c}{c}_{\alpha 1}^{\mathrm{y}}=196.38+800z_{\mathrm{g}}\hfill \\ c_{\alpha 2}^{\mathrm{y}}=178.40+800z_{\mathrm{g}}\hfill \\ c_{\alpha 1}^{\mathrm{z}}=331.20\hfill \\ c_{\alpha 2}^{\mathrm{z}}=189.60\hfill \end{array}$$

(20)

The dynamic response at TCP of the machine tool structure vary with the feed positions. And the dynamic response at TCP in different direction shows different trend. The expressions of DETW in X-, Y-, and Z-direction are given by equations (24)–(26) in Appendix. Considering the geometric parameters of the machine tool, these equations can be written as equations (21)–(23). These equations show that there is a negative correlation between $y_{\mathrm{e}}$ and Z-position. While $z_{\mathrm{e}}$ is independent of Z-position.

$$\begin{array}{c}{x}_{\mathrm{e}}=u_1+(-0.777+z_{\mathrm{g}}){\beta }_1+0.659{\gamma }_1+\hfill \\ u_2+(-1.00+z_{\mathrm{g}}){\beta }_2+0.482{\gamma }_2+\hfill \\ u_3-0.95{\beta }_3+0.245{\gamma }_3\hfill \end{array}$$

(21)

$$\begin{array}{c}{y}_{\mathrm{e}}=v_1-(0.777+z_{\mathrm{g}}){\alpha }_1+(0.377+x_{\mathrm{g}}){\gamma }_1+\hfill \\ v_2-(1+z_{\mathrm{g}}){\alpha }_2+v_3+0.95{\alpha }_3\hfill \end{array}$$

(22)

$$\begin{array}{c}{z}_{\mathrm{e}}=w_1-0.659{\alpha }_1-(0.377+x_{\mathrm{g}}){\beta }_1+\hfill \\ w_2-0.482{\alpha }_2+w_3-0.245{\alpha }_3\hfill \end{array}$$

(23)

Figure 14 shows the variation of $G_{\mathrm{z}}^{\mathrm{y}}$ and $G_{\mathrm{y}}^{\mathrm{z}}$ at different Z-position. It can be found that the magnitude of the second mode decreases if Z-position vary from 0 to $-500$  mm. On the contrary, the magnitude of the third mode increases if Z-position vary from 0 to $-500$  mm. In summary, the position dependency of $e_{\mathrm{z}}^{\mathrm{y}}$ is only affected by the variation of dynamic response at TCP. While the position dependency of $e_{\mathrm{y}}^{\mathrm{z}}$ is only affected by the variation of inertial forces. Therefore, the varying trend of $e_{\mathrm{z}}^{\mathrm{y}}$ with Z-position is different from that of $e_{\mathrm{z}}^{\mathrm{y}}$ .

Figure 14.

Frequency response for different Z-positions: (a) dynamic error in Z-direction during the accelerating stage of Y-axis and (b) dynamic error in Y-direction during the accelerating stage of Z-axis

Conclusions

This research was undertaken to model and analyze the dynamic error of TCP relative to worktable of multiaxis machine tool. The variations of structural dynamics and axis coupling inertial forces with feed positions are considered. The position dependency of dynamic error is discussed. The main conclusions are as follows:

1. The parameter-varying dynamic model of machine tool structure is proposed considering the variation of feed positions. The proposed model can describe the position dependent structural dynamics and inertial forces, and can calculate the dynamic error when the feed axes are in different positions.

2. The dynamic error in one direction is affected by the feed motions of multiple feed axes. The magnitudes of dynamic errors in Y- and Z-direction considerably vary with the feed position of Z-axis and have different varying trend. The relative decrease of the dynamic error in Y-direction is 35% when the feed position of Z-axis changes from 0 to $-300$  mm. The relative increase of the dynamic error in Z-direction is up to 67%.

3. The frequency domain analysis, numerical modal analysis, and the compositions of dynamic error reflect the effect of the structural dynamics on the dynamic error. The results indicate that the linear guideways of crossbeam and slider are the weak components of the machine tool structure.

4. The expressions of TCP dynamic response and inertial forces reveal the reason why the dynamic errors in Y- and Z-direction have different varying trend. The position dependency of the dynamic error in Y-direction is only determined by the variation of inertial forces. While the position dependency of the dynamic error in Z-direction is only determined by the variation of TCP dynamic response.