Characterizing machine tool spatial stiffness and predicting accuracy

Abstract

As a typical precision machine, the machine tool has complex structure and changeable working conditions, which leads to huge calculation data and low analysis efficiency when using traditional finite element method to obtain stiffness information. General test methods often have discrete values and cannot obtain stiffness information comprehensively and systematically. In this paper, a method for characterizing the stiffness of machine tool processing space based on multi-body system topology is proposed. Based on the matrix condensation theory, a method for constructing the equivalent condensation model of the stiffness of machine tool components is proposed. The stiffness characterization model of the components is established by using the stiffness equivalent condensation method. Based on the multi-body system theory, the stiffness distribution model of machine tool machining space is established, which can quickly and intuitively characterize the machining stiffness information of machine tool. Finally, based on the stiffness distribution model of the machining space of the machine tool, the machining accuracy of the S-shaped test specimen is predicted. Combined with the experimental analysis, the measurement results are consistent with the prediction results. The accuracy, feasibility, and effectiveness of the stiffness distribution model of machine tool machining space are verified.

Introduction

Machine tool stiffness is an important factor affecting the machining accuracy of a machine tool. Machine tool stiffness is affected by its topological layout, structural configuration, and other internal design factors, which are costly to adjust in the design and manufacturing process of machine tools and difficult to fully characterize.^1,2 In recent years, domestic and foreign experts and scholars on the stiffness information acquisition methods and stiffness characterization methods can be divided into stiffness numerical analysis methods and stiffness test methods.^3–5 The stiffness numerical analysis method uses the finite element to establish the digital model of the research object, and obtains the stiffness information through the numerical simulation and analysis of the model. The scholar Ying Jun et al. used the motorized spindle bearing of the CNC machine tool as the research object. Aiming at the problem that it is difficult to express the dynamic load of the motorized spindle, the bearing deformation and the bearing performance, and the bearing life, the nonlinear dynamic model of the motorized spindle bearing is established. Considering the influence of bearing defect shape, bearing deformation, and other factors on the performance degradation of the motorized spindle bearing, the stiffness deformation, performance degradation state, and theoretical life of the motorized spindle bearing are predicted under variable load conditions.^6–8 In order to effectively identify the weak parts of the cantilever structure in the tail frame and bed of the machine tool, Guo et al. used the finite element method to numerically simulate the cantilever structure of the machine tool, and established the stiffness reduction model of the cantilever structure. At the same time, the sensor was used to collect the stiffness deformation data of the cantilever structure of the machine tool. The results show that the finite element simulation results are in good agreement with the experimental structure, and the weak areas of the cantilever structure in the tail frame and bed are obtained.^9,10 The above stiffness numerical analysis method can not adapt to the complex and changeable working conditions in practical engineering, and the application range is greatly limited. The test and experimental methods of stiffness characterization are mainly to obtain stiffness information by direct measurement of the research object. Many scholars have designed many measurement methods to obtain stiffness information of different measurement objects. For different measurement objects, many scholars have designed many measurement methods^11–16; some scholars have introduced the shrinkage theory to simplify the stiffness information, and established a normalized stiffness characterization model for machine tool components,^17,18 which effectively reduces the amount of calculation. At the same time, the machine tool as a typical multi-body system, when the constituent parts of the machine tool is loaded, the deformation generated by itself and the external load will be transmitted along the tandem kinematic chain of the machine tool, and ultimately superimposed on the end of the machine tool actuator.^19,20 In terms of machine tool machining accuracy inspection and prediction, related scholars^21,22 designed relevant experiments to collect data such as contour error and surface quality of machine tool machining, so as to evaluate the advantages and disadvantages of machine tool machining accuracy.

In order to study the influence of the structural design of the machine tool on the accuracy of the machine tool, Cheng et al. ²³ used the multi-axis machine tool as the research object, and used the multi-body system theory to establish the volume error model of the machine tool. Based on this model, combined with the random process theory, the equivalent error of the machine tool processing space and the fluctuation of the processing accuracy were predicted. Scholars Okafor and Ertekin ²⁴ used the vertical machining center as the research object, used the rigid body kinematics theory to perform homogeneous coordinate transformation, and modeled the attitude error of the machining center. The model can predict the linear positioning accuracy of the machining center under multiple attitudes including roll and pitch.

In order to study the influence of cutting load on the accuracy of the machine tool during the machining process, Fan et al. ²⁵ proposed a mathematical model. The model can calculate the reaction force transmitted to the guide rail-slider of the machine tool according to the cutting load. According to the static equilibrium equation of the slip body deformation-free body, the contact deformation caused by the wear of the guide rail and the geometric error of the slider are solved, and then the position error of the machine tool after long-term operation can be predicted. Archenti and Nicolescu ²¹ established the reaction-deformation deviation calculation model of each joint in the machine tool space by arranging sensors to measure the load at the tool holder and workbench during the milling process. Through this model, the machining error caused by the joint stiffness and joint deformation of the machine tool can be solved during the whole machining process.

In order to study the influence of process parameters on the accuracy of machine tools during machining, scholars Chiu and Lee ²² used CNC machine tools as the research object. Firstly, relevant experiments were designed to collect data such as contour error and surface quality of machine tools. Based on fuzzy rules, a prediction model for these two indicators was established. Through this model, the contour and surface quality of machining can be predicted based on CNC machining parameters, so as to evaluate the machining accuracy of machine tools. Scholars Xie et al.^26,27 focused on the influence of feed speed on the position accuracy of machine tools, established a position accuracy prediction model through a large number of previous experimental data, and proposed a position accuracy control method combining real-time dynamic feedforward and feed speed adjustment. According to the position accuracy results obtained from real-time prediction, the adjustment of feed speed can be dynamically controlled, thereby limiting the machining error caused by the input signal.

Based on the above research status of machine tool machining accuracy analysis methods, it can be found that the machining accuracy of machine tools is affected by multiple factors such as its own structural design, cutting load, and processing parameters. At this stage, the prediction method of machine tool machining accuracy often only considers one of the factors, which is not consistent with the actual working conditions. Therefore, this paper proposes a universal method that can accurately predict the machining accuracy error of the machine tool based on the internal cause of the spatial stiffness of the machine tool and the external cause of the process system.

Based on the current research status of spatial stiffness characterization of machine tools, there are still problems such as incomplete characterization of machine tool stiffness, inability to normalize characterization, and difficulty in predicting machining accuracy. In this paper, a method for characterizing the stiffness of machine tool processing space based on multi-body system topology is proposed. Based on the matrix condensation theory, it is deduced that the stiffness characteristics of the component system can be characterized by the key geometric feature set that only bears the external load when the stiffness is equivalently condensed. The stiffness characterization model of machine tool components such as multi-constraint components and multi-position components is established by using the stiffness equivalent condensation method. Based on the multi-body system theory, the machine tool system is described as two series motion chains and the deformation transfer relationship is analyzed. The stiffness distribution model of the machine tool machining space is established, which can quickly and intuitively characterize the machining stiffness information of the machine tool. Finally, based on the stiffness distribution model of the machining space of the machine tool, the machining accuracy of the S-shaped test specimen is predicted. Combined with the experimental analysis, the measurement results are consistent with the prediction results. The accuracy, feasibility, and effectiveness of the stiffness distribution model of machine tool machining space are verified.

The spatial stiffness distribution of machine tools based on the theory of multibody systems

Topology configuration of machine tool multibody systems

A certain model of high-precision CNC (CNC machining center) machining center in the working state, can be abstractly viewed as a multi-body system consisting of five components: bed, column, table, spindle box, slide saddle.

In the multibody system of the machine tool, the inertial reference frame R is selected and labeled as D₀, the bed structure is labeled as D₁, and then the other constituent structures are numbered based on the direction away from the bed. Subsequently, the low-order body array method is applied to describe the topological structure of the machine tool multibody system, as shown in Figure 1:

Figure 1.

Structure and topology of machine tool multibody system.

The elements that constitute the topological structure of multi-body systems are called typical bodies of multi-body systems. In Figure 1, 1 is the bed, 2 is the column, 3 is the sliding saddle, 4 is the main axle box, 5 is the workbench, and 6 is the workpiece to be processed. Any component of the machine tool is a typical body in the multi-body system topology of the machine tool, and the order number definition of the nth-order low-order body D_S of the typical body D_K can be expressed using the following formula:

$$D^n\left(K\right)=S$$

(1)

In the formula: D– The low-order body operator in the topological structure of machine tool multi-body system; K, S– The code of the typical body in the multi-body system topology of the machine tool, where S is the n-order low-order body of K, and K is the n-order high-order body of S.

In the multi-body system topology of machine tools, two adjacent typical bodies should meet the following conditions:

$$D^n\left(K\right)=D\left(D^{n-1}\left(K\right)\right)$$

(2)

The relationship between the pose change of the key geometric features of each component in the multi-body system of the machine tool in its own local coordinate system and the spatial deformation error in the inertial coordinate system can be expressed as follows:

$$\left\{\begin{array}{c}4\hfill \\ \multicolumn{1}{c}{3}\\ \multicolumn{1}{c}{2}\\ \multicolumn{1}{c}{1}\\ \multicolumn{1}{c}{5}\\ \multicolumn{1}{c}{\prime 1}\end{array}\right\}=\left\{\begin{array}{cccccc}\mathit{I}\hfill & {\mathit{T}}_{43}\hfill & {\mathit{T}}_{42}\hfill & {\mathit{T}}_{41}\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \mathit{I}\hfill & {\mathit{T}}_{32}\hfill & {\mathit{T}}_{31}\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \mathit{I}\hfill & {\mathit{T}}_{21}\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \mathit{I}\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \mathit{I}\hfill & {\mathit{T}}_{51}\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \hfill & \mathit{I}\hfill \end{array}\right\}\left\{\begin{array}{c}4\hfill \\ \multicolumn{1}{c}{3}\\ \multicolumn{1}{c}{2}\\ \multicolumn{1}{c}{1}\\ \multicolumn{1}{c}{5}\\ \multicolumn{1}{c}{\prime 1}\end{array}\right\}$$

(3)

In the formula: Δ _K – The pose change of the key geometric features of each component of the whole machine in the inertial coordinate system; ε _K – The pose change of the key geometric features of each component in the local coordinate system; T – The pose transfer matrix between components, T _KS is the pose transfer matrix from component S to K, which represents the deformation error transfer relationship between the pose change ε _S of the key geometric features of component S and the pose change ε _K of component K.

Deformation transfer relationship of machine tool multi-body system

By analyzing the topological structure of the multi-body system of the machine tool, the motion form of the machine tool can be decomposed into two series motion chains shown in Figure 2. One is the tool chain i chain composed of “bed → column → saddle → spindle box → tool.” The other is the workpiece chain j chain composed of “bed → workbench → workpiece.” Therefore, the deformation error of the machining space at the end of the machine tool is mainly affected by the two motion chains of the tool chain and the workpiece chain.

Figure 2.

Kinematic chain model of machine tool.

Assuming that there are n components on the tool chain and m components on the workpiece chain, the deformation error of the final machining space at the end of the machine tool can be superimposed and expressed as the error subtraction of the two kinematic chains:

$$\left[f(x,y,z)\right]{=}_n-{\mathit{M}}_{m+1m}$$

(4)

In the formula: $\left[f(x,y,z)\right]$ – Deformation error in the machining space at the end of the machine; _n , _m – The deformation error of component n and component m; ${\mathit{M}}_{m+1}$ – The deformation error transfer matrix of the workpiece chain.

For any component K in a kinematic chain of a machine tool multi-body system, the factors affecting its error Δ _K mainly include the deformation of the component itself under load, and the deformation transfer error caused by the change of the position and posture of the adjacent low-order body components at the front end. The final synthesis error of the tool chain with n components and the workpiece chain with m components can finally obtain the deformation error of the machining space at the end of the machine tool:

$$\begin{array}{c}\multicolumn{1}{c}{\left[f(x,y,z)\right]=\sum _{i=1}^n[\left(\underset{p=i+1}{\overset{n}{\Pi }}{\mathit{M}}_p\right)·\prime i]}\\ -{\mathit{M}}_{m+1}\left\{\sum _{j=1}^m[\left(\underset{q=j+1}{\overset{m}{\Pi }}{\mathit{M}}_q\right)·\prime j]\right\}\hfill \end{array}$$

(5)

In the formula: M _p , M _q – They are the deformation error transfer matrices of the pth and qth components on the tool chain and the workpiece chain, respectively. $\prime i$ , $\prime j$ – is the deformation of the i and j components on the tool chain and the workpiece chain after loading, which is mainly composed of two parts: displacement deformation and rotation deformation.

Construction of the machine tool processing space stiffness distribution model

Based on the multi-body system theory, the analysis of the deformation transfer relationship of the machine tool components is completed, and the spatial stiffness of the machine tool is further characterized, and the spatial stiffness distribution model of the machine tool end processing is established.

It is assumed that there is an input equivalent combined load F _i ⁱⁿ on the load surface of the key geometric feature of a component i in the multi-body system of the machine tool, and an output equivalent combined load F _i ^out on the constraint surface of the key geometric feature. Because the constraint end of the component i is fixed, there is $\prime i-1$  = 0. In the multi-body system of the machine tool, the constraint end of the i component is also the load input end of the i−1 component. By deriving the reaction load transfer matrix from component i to component i−1, combined with the equivalent stiffness matrix K _i of each constituent component i of the machine tool, the deformation error of the end machining space of the machine tool multi-body system can be converted to , and thus the stiffness characterization model K of the end machining space of the machine tool multi-body system can be obtained:

$$\begin{array}{c}=\sum _{i=1}^n\left[{\mathit{M}}_{n+1}^{-1}{\mathit{C}}_{n+1}^{-1}{\mathit{K}}_i^{-1}\left(\underset{p=i+1}{\overset{n+1}{\Pi }}{\mathit{M}}_p{\mathit{C}}_p\right)\right]\hfill \\ -{\mathit{M}}_{m+1}\sum _{j=1}^m\left({\mathit{M}}_{m+1}^{-1}{\mathit{C}}_{m+1}^{-1}{\mathit{K}}_j^{-1}\left(\underset{q=j+1}{\overset{m+1}{\Pi }}{\mathit{M}}_q{\mathit{C}}_q\right)\right)\hfill \end{array}$$

(6)

$$\mathit{K}=\frac{1}{\sum _{i=1}^n\left[{\mathit{M}}_{n+1}^{-1}{\mathit{C}}_{n+1}^{-1}{\mathit{K}}_i^{-1}\left(\underset{p=i+1}{\overset{n+1}{\Pi }}{\mathit{M}}_p{\mathit{C}}_p\right)\right]-{\mathit{M}}_{m+1}\sum _{j=1}^m\left({\mathit{M}}_{m+1}^{-1}{\mathit{C}}_{m+1}^{-1}{\mathit{K}}_j^{-1}\left(\underset{q=j+1}{\overset{m+1}{\Pi }}{\mathit{M}}_q{\mathit{C}}_q\right)\right)}$$

(7)

In the formula: K _i , K _j – They are the equivalent stiffness matrix of the components i and j in the multi-body system of the machine tool; M _p , M _q – They are the deformation error transfer matrix of the components p and q in the multi-body system of the machine tool, which can be solved by the formula; C _p , C _q – These are the interaction load transfer matrices for components p and q in the machine tool multi-body system.

Solving characterization of spatial stiffness distribution model of machine tool machining

Machine tool component stiffness aggregation representation method

Taking the workbench components of a high-precision Computerized Numerical Control machining center shown in Figure 3 as an example, the EFGH surface is the mounting surface of the workpiece and its supporting fixture in contact with the workbench, which belongs to the bearing end of the workbench components. The load generated by cutting the workpiece will be transmitted through this way and input to the workbench components. The ABCD surface is the installation surface where the slider is in contact with the workbench. The workbench is installed on the bed part through this surface, which belongs to the constraint end of the workbench part. The load is transmitted and output to the bed part through this way. Therefore, A–H is the key geometric feature set that needs to be considered in the workbench components.

Figure 3.

Schematic diagram of equivalent aggregation principle for machine tool components.

For any multi-constrained component i with fixed position in the stiffness system of the machine tool, there are n key geometric features ABCDEFGH subjected to external loads, so the stiffness model of the component i can be represented by the stiffness matrix of A–H:

$${\mathit{K}}_i={\left[\begin{array}{ccc}{\mathit{K}}_{\mathrm{AA}}\hfill & \cdots \hfill & {\mathit{K}}_{\mathrm{AH}}\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\mathit{K}}_{\mathrm{HA}}}& \cdots \hfill & {\mathit{K}}_{\mathrm{HH}}\hfill \end{array}\right]}_{6n\times 6n}$$

(8)

At this time, the size of the stiffness matrix of component i is 6n × 6n.In order to realize the synthesis operation of the spatial stiffness of the whole machine in the later stage, it is necessary to carry out equivalent condensation, and the stiffness characteristics are reflected by the normalized stiffness matrix of 6 × 6. It is assumed that there is an input equivalent resultant force F _i ⁱⁿ on the loading surface EFGH, and there is an output equivalent resultant force F _i ^out on the constraint surface ABCD. The deformation effect of the input equivalent resultant force F _i ⁱⁿ and the output equivalent resultant force F _i ^out on component i is exactly the same as that of A–H multi-point loading. Then there is:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{K}}_{\mathit{i}}=\frac{1}{{\mathit{R}}_{\mathit{i}}}==\left[\begin{array}{cccccc}\frac{{\delta }_{\mathrm{xx}}}{f_x}\hfill & \frac{{\delta }_{\mathrm{xy}}}{f_y}\hfill & \frac{{\delta }_{\mathrm{xz}}}{f_z}\hfill & \frac{{\delta }_{\mathrm{xa}}}{m_a}\hfill & \frac{{\delta }_{\mathrm{xb}}}{m_b}\hfill & \frac{{\delta }_{\mathrm{xc}}}{m_c}\hfill \\ \multicolumn{1}{c}{\frac{{\delta }_{\mathrm{yx}}}{f_x}}& \frac{{\delta }_{\mathrm{yy}}}{f_y}\hfill & \frac{{\delta }_{\mathrm{yz}}}{f_z}\hfill & \frac{{\delta }_{\mathrm{ya}}}{m_a}\hfill & \frac{{\delta }_{\mathrm{yb}}}{m_b}\hfill & \frac{{\delta }_{\mathrm{yc}}}{m_c}\hfill \\ \multicolumn{1}{c}{\frac{{\delta }_{\mathrm{zx}}}{f_x}}& \frac{{\delta }_{\mathrm{zy}}}{f_y}\hfill & \frac{{\delta }_{\mathrm{zz}}}{f_z}\hfill & \frac{{\delta }_{\mathrm{za}}}{m_a}\hfill & \frac{{\delta }_{\mathrm{zb}}}{m_b}\hfill & \frac{{\delta }_{\mathrm{zc}}}{m_c}\hfill \\ \multicolumn{1}{c}{\frac{{\theta }_{\mathrm{ax}}}{f_x}}& \frac{{\theta }_{\mathrm{ay}}}{f_y}\hfill & \frac{{\theta }_{\mathrm{az}}}{f_z}\hfill & \frac{{\theta }_{\mathrm{aa}}}{m_a}\hfill & \frac{{\theta }_{\mathrm{ab}}}{m_b}\hfill & \frac{{\theta }_{\mathrm{ac}}}{m_c}\hfill \\ \multicolumn{1}{c}{\frac{{\theta }_{\mathrm{bx}}}{f_x}}& \frac{{\theta }_{\mathrm{by}}}{f_y}\hfill & \frac{{\theta }_{\mathrm{bz}}}{f_z}\hfill & \frac{{\theta }_{\mathrm{ba}}}{m_a}\hfill & \frac{{\theta }_{\mathrm{bb}}}{m_b}\hfill & \frac{{\theta }_{\mathrm{bc}}}{m_c}\hfill \\ \multicolumn{1}{c}{\frac{{\theta }_{\mathrm{bx}}}{f_x}}& \frac{{\theta }_{\mathrm{cy}}}{f_y}\hfill & \frac{{\theta }_{\mathrm{cz}}}{f_z}\hfill & \frac{{\theta }_{\mathrm{ca}}}{m_a}\hfill & \frac{{\theta }_{\mathrm{cb}}}{m_b}\hfill & \frac{{\theta }_{\mathrm{cc}}}{m_c}\hfill \end{array}\right]}\end{array}$$

(9)

In the formula: K _i – The equivalent aggregated stiffness matrix for component i under loading, with a scale of 6 × 6; R _i – The flexibility matrix that characterizes the relationship between the deformation of machine tool component i and the load; δ _ij – The displacement deformation (mm) of the bearing end of the component in the i (i = x, y, or z) direction caused by the force f _j (j = x, y, or z) or the moment m _j (j = a, b, or c) acting in the j direction. θ _ij – The rotational deformation (rad) of the bearing end of the component in the i (i = a, b, or c) direction caused by the force f _j (j = x, y, or z) or the moment m _j (j = a, b, or c) acting in the j direction.

Stiffness data acquisition of machine tool parts

Based on the characteristics of the spatial structure of the machine tool and the need for stiffness characterization, a unified machine tool coordinate system should be established first. As shown in Figure 4, it is specified that the x direction of the machine tool coordinate system is along the column stroke direction, the y direction of the machine tool system is along the sliding saddle stroke direction, and the z direction of the machine tool system is along the bed stroke direction. The rotation directions around the x axis, the y axis, and the z axis are a, b, and c, respectively, and the machine tool coordinate system satisfies the right-handed spiral rule.

Figure 4.

Schematic diagram of coordinate system for machine tool machining space.

A certain model of high-precision CNC (CNC machining center) machining center in the working state, can be abstractly viewed as a multi-body system consisting of five components: bed, column, table, spindle box, slide saddle. For each component, it can be divided into two categories according to its travel characteristics and constraint conditions. The stiffness modeling methods of the two types of components are generally the same, but the details are slightly different. One kind of multi-constraint parts, headstock, worktable, etc., its main feature is that the load and constraint conditions are complex, but the position of the parts and the load surface are fixed; The other type of multi-position parts, columns, sliding saddles, etc., whose main feature is that the parts contain a structure with variable stroke position, and the external load position of the parts will also move instantaneously along the stroke direction. Compared with the multi-constraint parts, this type of parts also need to discretization continuous stroke when establishing the stiffness model. It should be noted that the bed is connected with the column and the table at the same time, and the column is installed on the bed and fixed. Therefore, the part connected by the bed and the column should be analyzed according to the stiffness characterization method of the multi-constrained parts. The table has a stroke movement on the slide rail of the bed, so the part of the bed and the table should be analyzed according to the stiffness characterization method of the multi-position parts.

Combined with the finite element calculation, the deformation simulation data of the key geometric features of the machine tool components after loading are obtained, and the data are processed and calculated. Finally, the stiffness characterization model of the multi-constrained components is established in the form of matrix.

Taking the machine tool table component as an example, the stiffness characterization model of the multi-constraint component was established. The key geometric features of the workbench components are identified according to the method of stiffness equivalent condensation. Through analysis, it can be seen that in the workbench components, the mounting contact surface of the workpiece and its supporting fixture on the workbench is the bearing end of the component, and the mounting contact surface of the slider on the workbench is the constraint end of the component. Therefore, these two planes are the key geometric features of the table components. After establishing the finite element model of the workbench components, fixed constraints should be applied to the installation surface of the slider at the constraint end of the components, and the unit force (F _x , F _y , F _z ) in the three directions of X, Y, and Z and the unit torque (M _a , M _b , M _c ) in the three directions of A, B, and C should be applied to the working table at the bearing end in turn. Among them, F _x , F _y , F _z are 10,000 N, M _a , M _b , M _c are 10,000 N mm (Figure 5).

Figure 5.

Selection of key feature set for workbench.

According to the finite element simulation results shown in Figure 6, the deformation data of the feature points after the workbench is loaded can be obtained, and the workbench flexibility matrix shown in Formula (9) can be obtained. Similarly, the flexibility matrix of other parts of the machine tool can also be obtained by using the same method.

$${\mathit{R}}_{\mathrm{gzt}}=\left[\begin{array}{cccccc}8.4\times 10^{-7}\hfill & 6.7\times 10^{-8}\hfill & -5.3\times 10^{-8}\hfill & -1.2\times 10^{-10}\hfill & -1.7\times 10^{-10}\hfill & -1.9\times 10^{-9}\hfill \\ \multicolumn{1}{c}{2.1\times 10^{-7}}& 3.3\times 10^{-7}\hfill & 1.5\times 10^{-7}\hfill & 4.7\times 10^{-10}\hfill & -8.1\times 10^{-12}\hfill & -5.9\times 10^{-10}\hfill \\ \multicolumn{1}{c}{-5.6\times 10^{-8}}& 4.6\times 10^{-8}\hfill & 6.1\times 10^{-7}\hfill & 1.4\times 10^{-9}\hfill & 1.4\times 10^{-10}\hfill & 1.2\times 10^{-10}\hfill \\ \multicolumn{1}{c}{4.8\times 10^{-11}}& 1.7\times 10^{-10}\hfill & 1.1\times 10^{-9}\hfill & 3.0\times 10^{-12}\hfill & 2.1\times 10^{-12}\hfill & 2.1\times 10^{-12}\hfill \\ \multicolumn{1}{c}{5.2\times 10^{-9}}& 7.4\times 10^{-10}\hfill & 3.6\times 10^{-9}\hfill & 8.2\times 10^{-12}\hfill & 0.0\times 10^{-0}\hfill & 1.2\times 10^{-11}\hfill \\ \multicolumn{1}{c}{1.5\times 10^{-9}}& 3.7\times 10^{-11}\hfill & 1.3\times 10^{-10}\hfill & 2.1\times 10^{-12}\hfill & 0.0\times 10^{-0}\hfill & 3.7\times 10^{-12}\hfill \end{array}\right]$$

(10)

Figure 6.

Finite element simulation deformation results of the workbench: (a) deformation in X direction, (b) deformation in Y direction, (c) deformation in Z direction, (d) deformation in A direction, (e) deformation in B direction, and (f) deformation in C direction.

Similarly, stiffness characterization models of headstock, bed (connected to the column) and bed (connected to the table) can also be established in accordance with the methods described in the establishment of workbench stiffness characterization models, as shown in Figure 7.

Figure 7.

Flexibility matrix diagram.

Analysis of stiffness characteristics of machine tool processing space

The component stiffness matrix and the deformation transfer matrix association, can be obtained by the end of the machine tool machining space stiffness distribution model K, in order to more intuitively reflect the machine tool machining space stiffness of the specific distribution, can be established by the machining space coordinate system, will be obtained from the spatial stiffness distribution model results in the system in the form of the stiffness distribution of the form of the image to be expressed.

It is known that the machining space of the machine tool is a 2000 × 1850 × 2150 mm rectangular body, for the establishment of the machine tool machining space coordinate system for this range: define the slide saddle in the x travel range of the leftmost, the location of the tip of the tool point for the machine tool machining space coordinate system of the X-axis of the zero point; definition of the main spindle box in the y travel range of the lowest, the tip of the tool point of the location of the Y-axis of the zero point; define the worktable in the z travel range, the location of the tool tip of the zero point. When the table is closest to the end of the column in the z travel range, the position of the tool tip is the zero point of the Z-axis.

As shown in Figure 8, the ability of the machine tool to resist x-direction displacement deformation under the action of F _x is more sensitive in the direction of z-direction travel, while the change in x- and y-direction travel is not obvious. That is, under the action of F _x , if the workpiece clamped to the table moves back and forth on the z-travel of the bed, when the workpiece moves with the table to the zero point of the Z-axis of the machining space coordinate system, the displacement deformation of the machine tool is the largest compared to other positions, when the workpiece moves with the table to the furthest away from the end of the column, the displacement deformation of the machine tool is the smallest compared to other positions.

Figure 8.

Sliced diagram of stiffness distribution in machine tool machining space under external loads: (a) x-direction, (b) y-direction, and (c) z-direction.

Similarly, the same method can be used to analyze the stiffness characteristics of the machine tool by plotting slices of the stiffness distribution under other loads.

Machine tool machining accuracy prediction

Machine tool machining accuracy prediction model

In this paper, the milling process of the S-shaped test specimen shown in Figure 9 is selected as the application object. Using the spatial stiffness distribution model of the machine tool established above, the machining accuracy of the component surface of the S-shaped test specimen is predicted by considering the spatial stiffness of the machine tool and the influencing factors of the process system.

Figure 9.

The three-dimensional model of the S-shaped test piece.

If the mathematical characterization model of the geometric characteristics of the surface to be processed of the S-shaped test specimen is assumed to be $s(x,y,z)$ , the equivalent cutting load under the dynamic machining trajectory of the machine tool is assumed to be $F(f_x,f_y,f_z,m_a,m_b,m_c)$ , and the established machining space stiffness distribution model $\mathit{K}(x,y,z,\alpha ,\beta ,\gamma )$ of the machine tool is established, then the machining accuracy prediction model of the machine tool can be established as follows:

$$\begin{array}{c}\multicolumn{1}{c}{s(x,y,z)=\{\sum _{i=1}^n\left[{\mathit{M}}_{n+1}^{-1}{\mathit{C}}_{n+1}^{-1}{\mathit{K}}_i^{-1}\left(\underset{p=i+1}{\overset{n+1}{\Pi }}{\mathit{M}}_p{\mathit{C}}_p\right)\right]}\\ -{\mathit{M}}_{m+1}\sum _{j=1}^m\left({\mathit{M}}_{m+1}^{-1}{\mathit{C}}_{m+1}^{-1}{\mathit{K}}_j^{-1}\left(\underset{q=j+1}{\overset{m+1}{\Pi }}{\mathit{M}}_q{\mathit{C}}_q\right)\right)\}\mathit{F}\hfill \end{array}$$

(11)

In the formula: ${\mathrm{\Delta }}_{S(x,y,z)}$ – Machining error deformation; F– Equivalent cutting load under dynamic machining trajectory of the machine tool.

Equivalent cutting load analysis

Can be based on the “Cutting Dosage Concise Manual” ²⁸ test specimen milling process parameters, combined with the milling process cutting load equivalent conversion of the relevant theoretical formulas, can be solved to calculate the equivalent load.

When the milling process type and process parameters are determined, because the action point of the cutting load is at the tangent point of the tool circumference and the machined surface, this will cause the size and direction of the cutting load to change with the dynamic change of the tool path. On the other hand, due to the unique spatial geometric characteristics of the ruled surface of the S-shaped test specimen, the opening and closing angle between the ruled surface and the bottom surface of the test specimen is not a constant right angle, and the value of the opening and closing angle will change in real time along the direction of the guide line. During the cutting process, the tool circumference is always tangent to the machined surface, and the change of the opening and closing angle will lead to the change of the spatial angle between the tool coordinate system x′oy′ and the machine coordinate system xoy. Therefore, these two factors should be considered comprehensively, and the cutting load should be converted into the equivalent force (F _x , F _y , F _z ) in the three directions of X, Y, and Z and the equivalent torque (M _a , M _b , M _c ) in the three directions of A, B, and C in the machine tool system.

As shown in Figure 10, according to the characteristics of cutting load, the load action point can be selected as the coordinate origin of the tool coordinate system. The direction from the action point to the tool axis is taken as the x′-axis positive direction of the tool coordinate system. The reverse direction of the tool feed direction is taken as the y′-axis positive direction of the tool coordinate system. The direction parallel to the tool axis and deviating from the spindle is the z′-axis positive direction, and the tool coordinate system is established. At this time, there are spatial angles α and γ between the tool coordinate system and the machine tool coordinate system, and the coordinate r _i of each point coordinate r _d in the tool coordinate system in the machine tool coordinate system can be obtained by the following pose transformation:

$${\mathit{r}}_{\mathrm{j}}={\mathit{R}}_{\mathrm{jd}}{\mathit{r}}_{\mathrm{d}}=\left(\begin{array}{ccc}\cos \gamma \hfill & -\sin \gamma \hfill & 0\hfill \\ \multicolumn{1}{c}{\cos \alpha \sin \gamma }& -\sin \alpha \cos \gamma \hfill & -\sin \alpha \hfill \\ \multicolumn{1}{c}{\sin \alpha \sin \gamma }& \sin \alpha \cos \gamma \hfill & \cos \alpha \hfill \end{array}\right){(x_{\mathrm{d}},y_{\mathrm{d}},z_{\mathrm{d}})}^{\mathrm{T}}$$

(12)

Figure 10.

Diagrammatic sketch of flank milling of the S-shaped test piece.

In the formula: R _jg– Rotation transformation matrix from tool coordinate system to machine coordinate system.

The equivalent conversion theoretical formula of cutting load in side milling process in machine system is as follows:

$$\begin{array}{c}F(F_x,F_y,F_z,M_a,M_b,M_c)=\hfill \\ \{\begin{array}{c}{F}_x=F_{x\prime }\cos \alpha \sin \gamma -F_{y\prime }\sin \alpha \sin \gamma \hfill \\ \multicolumn{1}{c}{F_y=F_{x\prime }\sin \alpha +F_{y\prime }\cos \alpha }\\ \multicolumn{1}{c}{F_z=F_{x\prime }\cos \alpha \cos \gamma +F_{y\prime }\sin \alpha \cos \gamma }\\ \multicolumn{1}{c}{M_a=F_{x}D}\\ \multicolumn{1}{c}{M_b=F_{y}d}\\ \multicolumn{1}{c}{M_c=F_{z}d}\end{array}\hfill \end{array}$$

(13)

In the formula: D– The distance between the top surface of the S-shaped test piece and the end face of the spindle box is (mm). d– Milling cutter radius (mm).

In the side milling of straight grain surface, a cylindrical milling cutter should be selected as the processing tool, tool material selection of high-speed steel, tool diameter selection of d = 40 mm, the number of tool teeth selection of z = 6. At this time, if the cutting speed of v = 20 m/min, the feed is used f = 0.06 mm/r, the depth of cut is used a_p = 2 mm, then according to the “Cutting Dosage Concise Handbook,” under the parameters of this process. The selection of correlation coefficients in the cutting load exponential type empirical formulas are listed in Tables 1 to 3:

Table 1.

Correlation coefficient of empirical formula of the main cutting force in flank milling straight surface.

Processing forms	Workpiece materials	Cutting tool materials	CFc	xFc	yFc	nFc	KFc
Side milling	Structural steel	High speed steel	967	1.0	0.75	0.88	1

Table 2.

Correlation coefficient of empirical formula of the radial thrust force in flank milling straight surface.

Processing forms	Workpiece materials	Cutting tool materials	CFp	xFp	yFp	nFp	KFp
Side milling	Structural steel	High speed steel	1100	0.9	0.75	0	1

Table 3.

Correlation coefficient of empirical formula of the axial thrust force in flank milling straight surface.

Processing forms	Workpiece materials	Cutting tool materials	CFf	xFf	yFf	nFf	KFf
Side milling	Structural steel	High speed steel	590	1.2	0.65	0	1

By substituting the coefficient values listed in Tables 1 to 3 into the formula (13), the equivalent load formula (14) of the ruled surface process of the side milling S-shaped test specimen in the machine tool system can be obtained by calculation.

$$\begin{array}{c}F(F_x,F_y,F_z,M_a,M_b,M_c)=\hfill \\ \{\begin{array}{c}{F}_x=0.27\cos \alpha \sin \gamma -0.30\sin \alpha \sin \gamma \hfill \\ \multicolumn{1}{c}{F_y=0.27\sin \alpha +0.30\cos \alpha }\\ \multicolumn{1}{c}{F_z=0.27\cos \alpha \cos \gamma +0.30\sin \alpha \cos \gamma }\\ \multicolumn{1}{c}{M_a=42.92\cos \alpha \sin \gamma -48.39\sin \alpha \sin \gamma }\\ \multicolumn{1}{c}{M_b=10.73\sin \alpha +12.10\cos \alpha }\\ \multicolumn{1}{c}{M_c=10.73\cos \alpha \cos \gamma +12.10\sin \alpha \cos \gamma }\end{array}\hfill \end{array}$$

(14)

Figure 11 shows the process of determining the size error of the workpiece by using the empirical formula and the stiffness model.

Figure 11.

Flow chart.

Machining error prediction and evaluation analysis

According to the clamping position of the S-shaped test specimen in the machine tool coordinate system, the process characteristics of the ruled surface side milling process, and the structural characteristics of the processing equipment, it can be seen that the deformation of the machine tool in the x and y directions has a great influence on the machining accuracy of the ruled surface, and the influence in the z direction is much smaller than that in the x and y directions, which can be ignored.

Therefore, taking the side milling process of the straight-grained face A of the S-shaped test specimen as an example, and taking the clamping angle α = 0° around the X-axis, the machining error of this face is solved, and the specific distribution obtained from the prediction can be shown in Figure 12:

Figure 12.

Distribution of prediction errors for the straight-lined surface in side milling process: (a) x-direction and (b) y-direction.

From the above results, it can be concluded that the maximum machining error in the x-direction under this machining process is −204.97 μm, and the maximum machining error in the y-direction is 332.02 μm, and the advantages and disadvantages of the machining accuracy of the machine tool are further evaluated based on the face contour error. As shown in Figure 13, where the color in the figure represents the coordinates of the straight-grained surface A in the z-direction. Along the height direction of the straight surface of the S-shaped test specimen, three inspection intercepts are selected at equal spacing, the intercepts are numbered 1#, 2#, 3# from top to bottom, the spacing of the intercepts is 15 mm, and the distance to the bottom surface of the bottom 3# intercepts is 5 mm. Twenty-five measurement points are selected at equal spacing on each of the intercepts, and the points are numbered 1, 2, … , 24, 25, along the positive direction of the Y-axis, respectively, so that the 75 points selected for the inspection are selected at equal distances. The contouring errors of the 75 selected measurement points are then detected and evaluated.

Figure 13.

Selection of measuring points for contour error evaluation of straight surface.

Extract the actual contour coordinates of each measurement point, as well as the machining error of each measurement point in the x and y directions, respectively, as shown in Figure 14:

Figure 14.

Comparative analysis of machining error of measuring point of straight surface A: (a) x-direction and (b) y-direction.

The above results are analyzed. In the process of side milling ruled surface A of machine tool, the maximum machining errors in x and y directions will occur near the No.25 measuring point of each section, that is, they all occur at the end of the specimen away from the workbench in the length direction.

The farther away from the worktable clamping position on the ruled surface (near the No.25 measuring point), the greater the contour error is, and the closer to the worktable area (near the No.1 ~5 measuring point), the smaller the error is. However, the error of the part with large curvature of the ruled surface will also increase (near the No.8–11 and No.18–21 measuring points).

There are two reasons for this error distribution: on the one hand, the area near the 8–11 and 18–21 measuring points is the area with the largest curvature of the ruled surface, and the processing difficulty is high; on the other hand, due to the uneven distribution of the spatial stiffness field of the machine tool for side milling, the support stability of the area near the No.25 measuring point far from the clamping position of the worktable is poor.

Experiment

In order to verify the accuracy, feasibility, and validity of the proposed spatial stiffness distribution model of the machine tool, a high-precision machine tool is used as the research object, and the side milling process of the straight grain surface A of the S-shaped test specimen is taken as an example to detect the machining error of the S-shaped specimen by utilizing the coordinate measuring machine.

The material type of the S-shaped test piece is structural steel, and the dimensional parameters are shown in Table 4. In the side milling process of straight grain surface, the cylindrical milling cutter is selected as the processing tool, the tool material is selected high-speed steel, the tool diameter is selected as d = 40 mm, and the number of teeth of the tool is selected as z = 6. At this time, if the cutting speed is adopted as v = 20 m/min, the feed is adopted as f = 0.06 mm/r, and the depth of cut is adopted as a_p = 2 mm.

Table 4.

Size parameters of S-shaped specimen.

Parameter types	Parameter value
Overall external dimensions: length (mm) × width (mm)	400 × 200
S-shaped specimen external dimensions: length (mm) × width (mm)	250 × 180
Square base height (mm)	30
S-shaped edge height (mm)	40

After calibrating the plane and zero reference of the specimen, according to the pre-written measurement program, the measuring machine automatically completes the measurement of the specimen’s straight surface error, measuring 75 test points. The process of measuring the test piece is shown in Figure 15.

Figure 15.

Three coordinate measuring machine measuring specimen process.

The X and Y direction deformation errors of the S-shaped specimen measured by the three-coordinate measuring instrument are compared with the predicted X and Y direction deformation errors. As shown in Figure 16 and Tables 5 to 7, it can be seen that the deformation error measured by the three coordinate measuring instrument at 75 test points has the same trend as the deformation error predicted by the machining accuracy, and the larger deformation error is the two-segment curve part of the S-shaped specimen. The reason is that the bending radian of the ruled surface in these two regions is larger.

Figure 16.

Comparison of deformation error between prediction and test.

Table 5.

Numerical deviation of measuring points in 1# X direction and Y direction.

Station serial number	X relative error (%)	Y relative error (%)	Station serial number	X relative error (%)	Y relative error (%)	Station serial number	X relative error (%)	Y relative error (%)
1	0.047	0.027	10	0.024	0.014	19	0.043	0.014
2	0.020	0.064	11	0.036	0.015	20	0.009	0.014
3	0.065	0.167	12	0.023	0.043	21	0.030	0.014
4	0.400	0.000	13	0.016	0.015	22	0.020	0.034
5	0.063	0.034	14	0.023	0.017	23	0.012	0.013
6	0.000	0.012	15	0.030	0.017	24	0.011	0.000
7	0.012	0.009	16	0.029	0.050	25	0.020	0.000
8	0.019	0.011	17	0.021	0.017
9	0.016	0.007	18	0.058	0.008

Table 6.

Numerical deviation of measuring points in 2# X direction and Y direction.

Station serial number	X relative error (%)	Y relative error (%)	Station serial number	X relative error (%)	Y relative error (%)	Station serial number	X relative error (%)	Y relative error (%)
1	0.054	0.039	10	0.041	0.014	19	0.055	0.007
2	0.073	0.125	11	0.019	0.016	20	0.045	0.010
3	0.130	0.500	12	0.023	0.022	21	0.046	0.000
4	2.000	0.057	13	0.046	0.040	22	0.020	0.007
5	0.095	0.069	14	0.064	0.011	23	0.018	0.013
6	0.056	0.012	15	0.079	0.016	24	0.016	0.000
7	0.048	0.000	16	0.051	0.029	25	0.010	0.006
8	0.047	0.012	17	0.059	0.013
9	0.034	0.011	18	0.058	0.015

Table 7.

Numerical deviation of measuring points in 3# X direction and Y direction.

Station serial number	X relative error (%)	Y relative error (%)	Station serial number	X relative error (%)	Y relative error (%)	Station serial number	X relative error (%)	Y relative error (%)
1	0.063	0.200	10	0.034	0.015	19	0.034	0.007
2	0.091	1.000	11	0.038	0.024	20	0.028	0.000
3	0.188	0.069	12	0.023	0.026	21	0.031	0.024
4	0.750	0.058	13	0.030	0.029	22	0.020	0.000
5	0.077	0.044	14	0.040	0.011	23	0.006	0.010
6	0.036	0.041	15	0.095	0.016	24	0.005	0.006
7	0.048	0.027	16	0.093	0.029	25	0.010	0.012
8	0.059	0.024	17	0.057	0.022
9	0.035	0.015	18	0.058	0.016

Conclusion

In this paper, based on the deformation transfer relationship of multi-body system, a spatial stiffness distribution model of machine tool end machining is established. The stiffness condensation theory and finite element method are used to solve the stiffness distribution information of machine tool efficiently. The established spatial stiffness distribution model of machine tool is used to predict the machining accuracy of the component surface of the S-shaped test specimen of machine tool milling. The following conclusions are drawn:

1. Based on the theory of multi-body system, taking a high-precision CNC machine tool as an example, according to the spatial assembly relationship of the machine tool, the topological structure of the multi-body system of the machine tool is described by the low-order body array method. The spatial stiffness of the machine tool is further characterized, and the spatial stiffness distribution model of the end processing of the machine tool is established;

2. The stiffness polycondensation theory and finite element analysis method are used to obtain the equivalent stiffness matrix of the components of the machine tool to characterize its stiffness characteristics, so as to efficiently solve the stiffness distribution characteristics of the machine tool processing space, and intuitively show the stiffness distribution characteristics of the machine tool processing space;

3. Based on the characterization model of the spatial stiffness distribution of the machine tool and the mathematical characterization model of the geometric characteristics of the S-shaped specimen surface, combined with the equivalent conversion theoretical formula of the cutting load under the machining trajectory of the machine tool, the prediction model of the machining accuracy of the machine tool is established, and the profile error of the side milling ruled surface of the specimen is predicted;

4. Taking a high-precision machine tool as the research object, the contour error detection test of S-shaped specimen is carried out. The measurement results are consistent with the prediction results of machining accuracy in Chapter 4. The accuracy, feasibility, and effectiveness of the spatial stiffness distribution model of the machine tool are verified by experiments.