The finite element modeling and stability prediction of high-speed spindle system dynamics with spindle-holder-tool joints

2.1. FEM modeling principle of parallel rotors

The parallel rotors system consists of inner rotor and outer rotor. The unit nodes are numbered continuous when the rotor system is discrete. A transition section is added at the end of outer rotor, the information of which is given randomly. The assembly process includes giving the mass and stiffness matrixes of all the units firstly, assembling all the matrixes and removing the matrixes of transition section, finally, the mass matrix $\mathbf{M}$ and stiffness matrix $\mathbf{K}$ of parallel rotors could be got.

Assuming $p$ units in outer rotor and $q$ units in inner rotor, there are $p+q+1$ units in the parallel rotor system and $N$ nodes, $N=p+1+1+1+q+1$. The mass and stiffness matrices of each unit are 12×12 orders, which is consist of two nodes and each node owns six degree freedoms. The assembly of adjacent units is shown as Fig. 2. The mass and stiffness matrices of whole system are 6$N$×6$N$ orders, from which the transition section matrices will be removed. The final complete mass and stiffness matrices of parallel rotor are shown as Fig. 3.

Fig. 2Schemes of adjacent unit assembly

Fig. 3Schemes of parallel rotors assembly

2.2. Dynamic simulation of spindle-holder-tool joints

The spindle-holder-tool joints consist of double parallel rotor systems, one is spindle and holder, the other is holder and tool. The sectional drawing of complete spindle system is shown as Fig. 4. Dynamic simulation of joints is illustrated taken the spindle-holder joint as example. The distribution-spring model of spindle-holder joint is shown as Fig. 5 (Fig. 5(a) is the toper joint structure and Fig. 5(b) is the joint FEM).

Fig. 4Sectional drawing of spindle system

Assuming the spindle and holder joint is divided into $n$ stepped shafts, in the joint there are $n+1$ nodes, $i$, $i+1$,…, $i+n$ on the spindle and $j$, $j+1$,…, $j+n$ on the holder, respectively. The coupling forces on the spindle and holder are:

where, ${\mathbf{K}}_{coup}$ and ${\mathbf{C}}_{coup}$ are coupling stiffness and damping matrices at spindle-holder joint, ${\mathbf{X}}_p$ and ${\mathbf{X}}_q$ are the displacement vectors of spindle and holder at the joint, $\mathrm{\Delta }{\mathbf{F}}_i$ and $\mathrm{\Delta }{\mathbf{F}}_q$ are the coupling forces at spindle-holder joint, $p$ and $q$ are node numbers at spindle-holder joint.

The stiffness and damping matrices of spindle-holder joint expressed with ${\mathbf{K}}_c$ and ${\mathbf{C}}_c$ are:

Fig. 5The distribution-spring model of spindle-holder joint

a) Structure

b) FEM

The stiffness characteristic of distribution-spring is shown in Eq. (2) and the damping characteristic is shown in Eq. (3). The assembly method of stiffness on joint is shown in Fig. 6.

Fig. 6Stiffness assembly theory of spindle-holder joint

2.3. Modeling process of spindle-holder-tool joints

The spindle-holder-tool joints modeling should be combined with the system modeling, which is not independent exit. Based on the parallel rotors idea, the joints are built with distributed spring model. The parameters of spring numbers, density, stiffness and damping are set, the spring numbers and density of which could be amendment with experiment. The modeling process of spindle-holder-tool joints is shown as Fig. 7.

Step 1: building the FEM model of spindle, holder and tool respectively. Firstly, the different sections are divided according to inner and outer diameters, secondly, two transition sections are added at the end of spindle and holder respectively, finally, the nodes of units are numbered continuously including the transition sections.

Step 2: creating the spindle-holder-tool joints model. The spring density, number, stiffness and damping are chosen primary. Then the distribution-spring model is created with the divided units of joints using the spring stiffness and damping. The complete spindle system model will be got assembling the joints and spindle, holder, tool model.

Step 3: the modal test and dynamic simulation. The FRF at tool point of the built model is got though modal analysis, which is verified by modal test. The number and density of distribution-spring will be corrected through the test and simulation results.

Step 4: the model application. The dynamic and stability prediction will be done with the complete spindle model considering the spindle-holder-tool joints.

Fig. 7Modeling process of spindle-holder-tool joints

2.4. FEM of complete spindle system

The high-speed spindle system is established with isoparametric beam element, each unit of which exits 12-DOF including 3 translational DOFs and 3 rotational DOFs on each node. The unit mass, stiffness and damping matrixes are assembled into the dynamic equation of system considering bearing force and spindle-holder-tool joints as follow:

where, $\mathbf{m}$ is mass matrix; $\mathbf{c}$ is damping matrix; $\mathbf{k}$ is stiffness matrix; $\mathbf{D}$ is damping matrix; $\mathbf{J}$ is gyroscopic matrix; ${\mathbf{k}}_H$ is coupling stiffness matrix of joints; ${\mathbf{k}}_s$ is stiffness matrix of spindle; ${\mathbf{k}}_C$ is centrifugal stiffness matrix of spindle; ${\mathbf{F}}_R$ is supporting bearing force vector; ${\mathbf{F}}_H$ is hysteretic force vector of spindle-holder-tool; ${\mathbf{F}}_P$ is cutting force vector; $\mathbf{G}$ is gravity vector; $\mathbf{X}$ is displacement vector; ${\mathbf{q}}_i$ is displacement vector of node $i$; $x_i$, $y_i$, $z_i$, ${\theta }_{xi}$, ${\theta }_{yi}$, ${\theta }_{zi}$ are 3 translational and 3 rotational DOFs of node $i$.

The FEM model of high-speed spindle system in the free-state is shown as Fig. 8. There are 66 units and 66 nodes in the FEM including 2 transition sections. The structures of spindle system are shown in Table 1 and Table 2. The dynamic parameters of system are: elasticity modulus $E=$2.06×10¹¹ Pa, material density $\rho =$7850 kg/m³. The FEM of test object is built with 5-distribution-spring connection, 10-distribution-spring connection and rigid connection respectively.

Fig. 8The model of high speed spindle system

3.1. Test object

The modal experiment is done on VMC0540d by hammering to prove the effective of modeling method for spindle-holder-tool joints and the whole spindle system. The results obtained from the experiment, 5-distribution-spring model. 10-distribution-spring model and rigid model will be compared by the frequency response function of tool point.

The test system is set up and shown in Fig. 9. The test object is a spindle system which consists of HCS150 motorized spindle from GMN, BBT50-MEGA32D-165 holder and HLXX32 end milling cutter. The main test equipments include the data acquisition, hammer and vibration sensors in Table 3.

Fig. 9The test system

The FRFs at tool point of rigid connection joints, 5-distribution-spring joints and test results are show as Fig. 10(a), and the FRFs of 5-distribution-spring joints, 10-distribution-spring joints and test results are show as Fig. 10(b).

Fig. 10Comparison of different joints models and test results

a) Different model types

b) Different spring density

The inherence frequencies of joints different forms are shown in Table 4 including rigid connection, concentrated spring model, different distribution-spring models and test results. There is no rigid modal in test results. From the Table, the simulation errors of rigid model are largest than other connections from the test results, which is from 82.65 % to 158 %.

The Fig. 10 shows that the FRF of 5-distribution-spring joint model is closer to the test result than 10-distribution-spring model, and also closer than the rigid one. The Table 4 shows that the distribution-spring model results are much closer to the test than the concentrated. For the spring layout form, the 5-distribution-spring joint model are more perfect than the 10-distribution-spring joints for the size of joints which is closer to the other unit in the model, which illustrates the influence of FRF on distribution-spring form is obvious and the 5-distribution-spring joint model is verified effectively.

4.1. Influence factors of stability

The cutting stability is evaluated through lobe diagram generally, which is characterized by cutting limit depth, as Eq. (5) [14]:

where, $b_{lim}$ is the maximum cutting depth without chatter, $Z$ is cutter teeth, $K_f$ is cutting force coefficient in feed direction, $G\left({\omega }_c\right)$ is real part of FRF, ${\omega }_c$ is angular frequency of negative FRF.

From the Eq. (5), the zone under $b_{lim}$ is unconditionally stability zone affected by $Z$, $K_f$ and $G$. The system will chatter when the cutting depth is higher than the $b_{lim}$ under specific speed theoretically. However, the emergence of conditionally stability zone will make the system stability with higher cutting depth. The conditionally stability zone is affected by lobe distribution, which still is controlled by FRF. Therefore, the FRF is a significant factor of cutting stability.

The FRF of system could be obtained by direct method and inherent characteristics method except modal test. The direct method is using the concept of frequency response function, which is the Laplace transformation ratio of response (or output) and excitation (or input) when the initial condition is zero in a linear constant system as Eq. (6):

where, $y\left(s\right)$ is Laplace transformation ration of response, $F\left(s\right)$ is Laplace transformation ration of excitation, $\psi \left(s\right)$ is frequency response function of system.

The inherent characteristics method begins with the essence of FRF. When the system equation is $\mathbf{m}\ddot{\mathbf{x}}+\mathbf{c}\ddot{\mathbf{x}}+\mathbf{k}\mathbf{x}=\mathbf{F}\left(t\right)\text{,}$ Laplace transformation of two sides is $\left(\mathbf{m}{s}^2+\mathbf{c}s+\mathbf{k}\right)\stackrel{-}{\mathbf{x}}\left(s\right)=\stackrel{-}{\mathbf{F}}\left(s\right)$, then there is:

where, ${\mathbf{H}}_d\left(s\right)$ is transfer function matrix of system, $\mathbf{F}\left(s\right)$ is Laplace transformation of excitation vector, $\stackrel{-}{\mathbf{x}}\left(s\right)$ is Laplace transformation of response vector:

Take the Eq. (8) into Eq. (7), $s=j\omega$, the displacement FRF of excitation at $q$, and response at $p$ is as follow:

where, ${\phi }_{pi}$ and ${\phi }_{qi}$ are the $i$ order modal shape of $p$ and $q$, ${\xi }_i$ is damping ratio, ${\omega }_i$ is inherence angular frequency, $\omega$ is excitation angular frequency.

From the Eq. (7) to (9), there are many influence parameters on FRF including excitation $F$, response $y$, mass $m$, stiffness $k$, damping $c$, modal shape $\phi$, inherence frequency $\omega$, damping ratio $\xi$, even the spindle-holder-tool joints dynamic characteristics discussed in the paper above and rotational speed verified by [15]. This paper focuses on the influence of cutting stability on cutting force (excitation) and speed (process condition) with specific system.

4.2. Stability prediction process

Based on modeling the system with distribution-spring joint model and solving the FRF, the lobe diagram will be drawn as follow:

Step 1: the FRF is divided into two parts: the real and imaginary part as $\varphi \left(j{\omega }_c\right)=G+jH$. Solve the value of ${\omega }_c$ to satisfy $G\left({\omega }_c\right)<0$ and obtain the $b_{lim}$ according to the $K_f$ and Eq. (5).

Step 2: positive and negative of imaginary part $H\left({\omega }_c\right)$ is judged by ${\omega }_c$. If $H\left({\omega }_c\right)<0$ then the phase angle of system will be $\psi =-\pi +\text{arctan}(H/G)$ or $\psi =-\pi -\text{arctan}\left(\right|H/G\left|\right)$. The rotational speed $n$ will be got according to $\epsilon =3\pi +2\psi$, $T=(l+\epsilon /2\pi )/f_c$ and $n=60/T$.

Step 3: the $b_{lim}$ and $n$ will be obtained in each lobe. The lobe diagram will be formed through drawing several lobes in the picture.

5.1. Influence of cutting force

The cutting force during the process is simulated by harmonic force, and the influence rule of stability is analyzed by changing the force amplitude. The three-dimensional lobe diagram with changing forces is shown as Fig. 11(a) and stability lobe diagram under two different forces is shown as Fig. 11(b), both of which show the limit cutting depth reduces as the force increases and the conditional zone affected by the force.

The limit cutting depths and lobe intersections correspond to the critical speeds changing with forces are shown as Fig. 12. From the Fig. 12(a), the limit cutting depths are reduced as the force increasing, which are affected greatly when the cutting force is small or gently when it is bigger. The first intersection is the one of $l=$0 and $l=$1, then the second intersection is the one of $l=$1 and $l=$2 from the Fig. 12(b), which shows the intersections move to the left as the force increasing. Therefore, the conditional zones are changed with the forces.

Fig. 11Stability lobes of different cutting forces

a) Three-dimensional lobe

b) Lobe under different cutting forces

Fig. 12Stability influence rules of cutting forces

a) Minimum limit cutting depth

b) Lobe intersection

5.2. Influence of rotational speed

The three-dimensional lobe diagram of stability with changing rotational speeds is shown as Fig. 13(a) and stability lobe diagram under two different sweep speeds are shown as Fig. 13(b). As the rotational speed rise, the limit cutting depths increase, which illustrates the cutting stability strengthen.

Fig. 13Stability lobes of different sweep speeds

a) Three-dimensional lobe

b) Lobe under different sweep speeds

The limit cutting depths and lobe intersections correspond to the critical speeds changing with rotational speed are shown as Fig. 14. From the Fig. 14(a), the limit cutting depths are enlarged as the speed increasing, but the first order remarkably and the secondly order is slow. Fig. 14(b) shows the intersections move to the right as the speed increasing, which states the conditional zones change complicated.

Fig. 14Stability influence rules of sweep speeds

a) Minimum limit cutting depth

b) Lobe intersection