Chatter-free cutter postures in five-axis machining

Abstract

In five-axis machining, stability of the machining process is determined not only by the combination of depth of cut and spindle speed but also the cutter postures (orientations). In this article, we propose to construct posture stability graphs to guide the selection of cutter postures during tool-path generation. Posture stability graph provides a partitioning of the allowable range of cutter postures, dividing it into stable and chatter zones. Using postures only from the stable zone, we can achieve chatter avoidance without the trouble of tuning other process parameters like spindle speed. Posture stability graphs are constructed by identifying the cutter orientations that would render the machining system borderline stable under the given machining conditions. These postures would make the boundaries between stable and unstable zones on the posture stability graph. Such a process is achieved by modelling the machining system as a 2-degrees-of-freedom spring-mass-damper system and applying the full-discretization method for chatter identification. Many experiments have been carried out, based on which the effectiveness of the posture stability graph approach has been verified.

Introduction

Possible existence of machining chatter is a key adverse factor that affects the surface finish quality for all machining processes. Chatter is a kind of self-excited and regenerative vibration, usually caused by lack of stiffness of the machining system, made of the cutter, holder, workpiece, and the machine tool itself.^1–3 Much research^4–11 has been done regarding the mechanism of chatter, which shows that a bad combination of spindle speed and depth of cut would trigger chatter in three-axis machining. Stability lobe diagram (SLD) has since been developed as a practical tool to guide machinists to select the proper machining condition for chatter avoidance.^12–14

Meanwhile, for five-axis machining, cutter orientation or posture is probably the most important process parameter. ¹⁵ A cutter posture is usually described by its lead and tilt angles, defined against the cross-feed and feed directions during machining (see Figure 1(a)). Intuitively, the introduction of cutter posture as an extra variable would make chatter prediction more complicated. An example is provided in Figure 2, showing a ball-nose cuter machining, a straight ‘ditch’ on a flat surface with different cutter postures. The machining parameters are kept exactly the same in both cases (depth of cut: 0.5 mm, spindle speed: 4800 r/min, material: AL5075) except the cutter posture. As shown in Figure 2, completely different surface finish qualities can be observed from the microscopic images of the experiments, machined with two cutter postures that have only 20° difference in the lead angle. A clear distinction can be made between stable and chatter machining. To address the influence of cutter postures on machining stability, we propose the concept of posture stability graph (PSG), which will partition the allowable range of cutter postures into stable and unstable zones. The target is to use the knowledge of machining dynamics to guide the selection of cutter postures during process planning in an integrated manner.

Figure 1.

Cutting force prediction in five-axis machining: (a) CCS, (b) top view of CCS, (c) discretization of cutter geometry, and (d) calculation of UCT.

Figure 2.

Influence of cutter posture on chatter stability and surface finish quality: (a) slotting with a ball-nose cutter, (b) lead angle: 60°; tilt angle: 0°, and (c) lead angle: 40°; tilt angle: 0°.

Construction of PSG starts with the prediction of machining chatter in five-axis machining. Over the years, much research^5,10,16–19 has been done to predict chatter in three-axis machining, most of which focuses on determining maximum chatter-free depth of cut for a specified spindle speed by modelling the machining system dynamically. Sridhar et al. ¹⁰ modelled the milling operations using linear differential difference equations with periodic coefficients, the stability of which can be identified analytically. Later, Minis et al. ¹⁶ solved the stability problem for two-dimensional (2D) dynamic milling problems with an iterative method, using Nyquist stability criterion to identify the stability limits. Budak and Altintas ⁵ determined the chatter stability limit analytically by approximating the time-varying directional machining coefficients using their Fourier series components. Many other methods stem from this approach, aimed at better accuracy and lower computational cost. A more detailed discussion on them is provided in section ‘Stability analysis with FD method’ later in this article when we explain the reason behind our choice of method.

In the meantime, limited study is reported for chatter stability of five-axis milling. Generally speaking, the introduction of different cutter orientations other than the one normal to the machining surface affects chatter stability in two ways: (1) it changes the cutter–workpiece engagement condition, which effectively determines the dynamic machining force. (2) It affects the distribution of modal parameters regarding the feed direction, that is, the two orthogonal degrees-of-freedom (DOF) of the cutter will no longer be aligned with the feed and cross-feed directions as in three-axis machining. Some of the early works to identify these influences are those of Ozturk et al.,^20,21 in which the cutter–workpiece engagement region is identified analytically under given lead and tilt angles and the modal parameters are transformed into the cutter frame. Both single-frequency and multi-frequency methods have been applied to determine the stability of the chatter equation. ²² Similar work includes that of Shamoto and Akazawa, ²³ who used a variation of the zeroth order method (by predicting the mean of the cutting force) for chatter prediction. These works are significant in identifying the influence of cutter orientation on machining stability. Yet, such approaches can only deal with machining with a fixed posture, which has limited applications. Mousseigne et al. ²⁴ addressed a special machining process when a torus cutter is used in climb milling. The experiment setting makes sure that the tool–workpiece engagement as well as the cutter’s orientation remains constant. This enables the traditional SLD to be built for the case. Ahmadi and Ismail established the stability maps relating cutter location (CL) and spindle speeds for five-axis peripheral milling of ruled surfaces. With this machining set-up, the cutter will be tangent to the machining surface. Thus, a change in CL is effectively a change in cutter orientation. ²⁵ This piece of work is a great demonstration for the influence of cutter orientation on chatter stability. Yet, simplification by assuming the cutter–workpiece engagement condition remains constant makes it not directly applicable to five-axis end milling.

During five-axis machining of freeform surfaces, the cutter postures would be continuously adjusted. The depth of cut, as affected by residue material from previous pass, is also changing continuously between CLs. Following the approach of existing work, the only viable solution is to construct the SLD and adjust the spindle speed at each CL. However, there would be many adverse effects, such as cutter marks from the spindle speed adjustment, increased power consumption, and deviation from optimal cutter speed. A better option is to make sure that all the cutter postures are stable during process planning. This requires the partitioning of the posture solution space into chatter/stable zones, which is the idea behind PSG. An analogy can be made between the proposed method and the division of posture solution space into accessible/inaccessible regions based on interference avoidance requirements. ²⁶

In the remaining of this article, the model of the five-axis machining system is to be presented and the method to solve for the stability of the system using the full-discretization (FD) method is to be introduced. Based on the result of the stability analysis, PSG can be constructed for any combination of depth of cut and spindle speed. Several experiments, including simultaneous five-axis machining, have been designed and conducted to validate the proposed method.

Chatter prediction for five-axis machining

Chatter is predicted based on the stability of the machining system under the excitation of cyclic machining forces. In this section, the dynamic model of the machining system is to be established, followed by the introduction of the proposed method to determine the stability conditions of any given cutter posture.

Cutting force model in five-axis machining

Chatter prediction starts with the prediction of cutting force, which acts as the excitation input for the machining system. Cutting force prediction is conducted under the assumption that only ball-part of the cutter is engaged with the workpiece, which is generally true for five-axis machining as it is mainly intended for finish machining purposes.

For five-axis machining, the cutter’s orientation is described with a pair of lead and tilt angles with respect to the feed direction of machining. For convenience, a workpiece coordinate system (WCS) O_W−X_WY_WZ_W is constructed, with X_W coincident with the feed direction and Z_W coincident with the surface normal (see Figure 1). A cutter coordinate system (CCS) O_T−X_TY_TZ_T is also constructed, with Z_T along cutter axis. When the cutter is placed normal to the machining surface, WCS and CCS are assumed aligned. A lead or tilt angle means the cutter should rotate around Y_W (cross-feed direction) or X_W (feed direction), which gives us the transformation matrix between CCS and WCS as

$$\begin{array}{c}\multicolumn{1}{c}{T_{\mathrm{cTw}}=\left[\begin{array}{ccc}c(\mathrm{tilt})\hfill & 0\hfill & s(\mathrm{tilt})\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{-s(\mathrm{tilt})}& 0\hfill & c(\mathrm{tilt})\hfill \end{array}\right]}\\ ·\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& c(\mathrm{lead})\hfill & -s(\mathrm{lead})\hfill \\ \multicolumn{1}{c}{0}& s(\mathrm{lead})\hfill & c(\mathrm{lead})\hfill \end{array}\right]\hfill \end{array}$$

(1)

where c(tilt) is short for cos(tilt), and s(tilt) is short for sin(tilt), and this abbreviation will be used in the remaining of this article.

Cutting force prediction is conducted in CCS. As cutting edges are helical on ball-nose cutters, the cutter is discretized along the cutter axis into N layers.

A cutting edge segment on a layer can be assumed straight and undergoing oblique cutting conditions. ²⁷ The cutting force generated by a single elemental cutting edge is decomposed into tangential, radial, and axial cutting forces, denoted as dF_t, dF_r, and dF_a, respectively, and given as ¹

$$\begin{array}{c}\mathrm{d}{F}_{\mathrm{t}}=K_{\mathrm{te}}\mathrm{d}S+K_{\mathrm{tc}}{t}_{\mathrm{uct}}\mathrm{d}b\hfill \\ \mathrm{d}{F}_{\mathrm{r}}=K_{\mathrm{re}}\mathrm{d}S+K_{\mathrm{rc}}{t}_{\mathrm{uct}}\mathrm{d}b\hfill \\ \mathrm{d}{F}_a=K_{\mathrm{ae}}\mathrm{d}S+K_{\mathrm{ac}}{t}_{\mathrm{uct}}\mathrm{d}b\hfill \end{array}$$

(2)

where K_te, K_re, and K_ae are the edge force coefficients and K_tc, K_rc, and K_ac are the shear coefficients, which can be calibrated for any given combination of cutter and material. dS is the differential element length of the cutting edge, db is the uncut chip width, and t is the uncut chip thickness (UCT) at the elemental cutting edge. From a kinematic point of view, t_uct can be calculated as the displacement of the cutter projected along the surface normal at location of the elemental cutting edge given as

$$\begin{array}{c}{t}_{\mathrm{uct}}=\mathrm{ud}\hfill \\ u={\left[\begin{array}{c}s\left(\kappa \right)s\left(\psi \right)\hfill \\ \multicolumn{1}{c}{s\left(\kappa \right)c\left(\psi \right)}\\ \multicolumn{1}{c}{c\left(\kappa \right)}\end{array}\right]}^{\mathrm{T}}\hfill \\ d=\left[\begin{array}{c}x\left(t\right)-x(t-T)\hfill \\ \multicolumn{1}{c}{y\left(t\right)-y(t-T)}\\ \multicolumn{1}{c}{z\left(t\right)-z(t-T)}\end{array}\right]\hfill \end{array}$$

(3)

where κ and ψ are the height and angular position angles of the cutter element, respectively, (see Figure 1(b) and (c)), [x(t) y(t) z(t)]^T gives the location of the cutter centre at time t, and T is the tooth passing period of the cutter.

During machining, a cutting edge will enter and exit from the engagement region. When mapped onto each cutter disc, the engagement region is represented with one or several pairs of entry and exit angles (φ_st, φ_ex). Cutting force will only be generated by certain cutting edge elements whose angular positional angle satisfies φ_st < ψ < φ_ex. Hence, the following step function g(ψ) is proposed

$$\begin{array}{c}g(\psi )=1\mathrm{when}{\phi }_{\mathrm{st}}<\psi <{\phi }_{\mathrm{ex}}\hfill \\ g(\psi )=0\mathrm{when\; otherwise}\hfill \end{array}$$

(4)

Finally, the elemental cutting force can be transformed to the cutter frame as

$${\left[\begin{array}{c}\mathrm{d}{F}_{\mathrm{x}}^{\mathrm{j}}\hfill \\ \multicolumn{1}{c}{\mathrm{d}{F}_{\mathrm{y}}^{\mathrm{j}}}\\ \multicolumn{1}{c}{\mathrm{d}{F}_{\mathrm{z}}^{\mathrm{j}}}\end{array}\right]}_{\mathrm{c}}=g({\varphi }_{\mathrm{j}})T_{\mathrm{l}}\left[\begin{array}{c}\mathrm{d}{F}_{\mathrm{r}}\hfill \\ \multicolumn{1}{c}{\mathrm{d}{F}_{\mathrm{t}}}\\ \multicolumn{1}{c}{\mathrm{d}{F}_{\mathrm{a}}}\end{array}\right]$$

(5)

where T_l is the transformation matrix given as

$$T_{\mathrm{l}}=\left[\begin{array}{ccc}-s(\kappa )s(\psi )\hfill & -c(\psi )\hfill & -c(\kappa )s(\psi )\hfill \\ \multicolumn{1}{c}{-s(\kappa )c(\psi )}& s(\psi )\hfill & -c(\kappa )c(\psi )\hfill \\ \multicolumn{1}{c}{c(\kappa )}& 0\hfill & -s(\kappa )\hfill \end{array}\right]$$

(6)

A summation of the elemental cutting forces contributed by each cutter disc will give the total cutting force acting on the cutter. The engagement region can be obtained via geometric analysis. It is worth noting that when machining with a ball-nose cutter, the change in cutter postures will not affect the shape of the engagement region, but only the distribution of the region on the cutter surface. This means engagement information for any posture can be reconstructed from that of a single input posture with reduced computation load (see Figure 3). ²⁸

Figure 3.

Construction of engagement information for prediction of cutting force: (a) engagement information of three-axis posture and (b) engagement information of five-axis posture reconstructed from three-axis posture.

Dynamic model of machining system

In our model, the cutter is assumed to have two DOF in two orthogonal directions, both normal to the cutter axis. The modal parameters along these two directions can be obtained with two separate impact tests. A modal frame O_M−X_MY_MZ_M is attached to the cutter, where X_M and Y_M are along the two impact directions and Z_M is along the cutter axis (see Figure 4).

Figure 4.

Illustration of dynamic machining system.

During machining, the static part of the machining force, that is, the part proportional to dS, which is not affected by the UCT, will not contribute to the regenerative mechanism (the cause behind chatter) and will thus be dropped. Substituting db = dz/s(κ), where dz is the thickness of one cutter disc, the following equation can be obtained

$$\begin{array}{c}{\left[\begin{array}{c}\mathrm{d}{F}_{\mathrm{x}}^{\mathrm{j}}\hfill \\ \multicolumn{1}{c}{\mathrm{d}{F}_{\mathrm{y}}^{\mathrm{j}}}\\ \multicolumn{1}{c}{\mathrm{d}{F}_{\mathrm{z}}^{\mathrm{j}}}\end{array}\right]}_{\mathrm{C}}=T_{\mathrm{l}}\left[\begin{array}{c}\mathrm{d}{F}_{\mathrm{r}}\hfill \\ \multicolumn{1}{c}{\mathrm{d}{F}_{\mathrm{t}}}\\ \multicolumn{1}{c}{\mathrm{d}{F}_{\mathrm{a}}}\end{array}\right]=g\left(\psi \right)\frac{\mathrm{d}z}{s\left(\kappa \right)}{T}_{\mathrm{l}}\hfill \\ ·\left[\begin{array}{c}{K}_{\mathrm{rc}}\hfill \\ \multicolumn{1}{c}{K_{\mathrm{tc}}}\\ \multicolumn{1}{c}{K_{\mathrm{ac}}}\end{array}\right]u(\left[\begin{array}{c}{x}_{\mathrm{C}}(t)\hfill \\ \multicolumn{1}{c}{y_{\mathrm{C}}(t)}\\ \multicolumn{1}{c}{z_{\mathrm{C}}(t)}\end{array}\right]-\left[\begin{array}{c}{x}_{\mathrm{C}}(t-T)\hfill \\ \multicolumn{1}{c}{y_{\mathrm{C}}(t-T)}\\ \multicolumn{1}{c}{z_{\mathrm{C}}(t-T)}\end{array}\right])\hfill \end{array}$$

(7)

The additional subscript in x_C(t) indicates that the calculation takes place in CCS. Summing up the elemental cutting force from all the involved edges, the total dynamic forces are calculated as

$${\left[\begin{array}{c}{F}_{\mathrm{x}}\hfill \\ \multicolumn{1}{c}{F_{\mathrm{y}}}\\ \multicolumn{1}{c}{F_{\mathrm{z}}}\end{array}\right]}_C=\sum _{i=1}^n\sum _{j=1}^N\left[\begin{array}{c}\mathrm{d}{F}_{\mathrm{x}}^{\mathrm{j}}\hfill \\ \multicolumn{1}{c}{\mathrm{d}{F}_{\mathrm{y}}^{\mathrm{j}}}\\ \multicolumn{1}{c}{\mathrm{d}{F}_{\mathrm{z}}^{\mathrm{j}}}\end{array}\right]$$

(8)

Defining

$$B(t)=\sum _{i=1}^n\sum _{j=1}^N\left(\frac{g\left(\psi \right)}{s\left(\kappa \right)}{T}_{\mathrm{l}}\left[\begin{array}{c}{K}_{\mathrm{rc}}\hfill \\ \multicolumn{1}{c}{K_{\mathrm{tc}}}\\ \multicolumn{1}{c}{K_{\mathrm{ac}}}\end{array}\right]n\right)$$

Equation (8) can be rewritten as

$${\left[\begin{array}{c}{F}_{\mathrm{x}}\hfill \\ \multicolumn{1}{c}{F_{\mathrm{y}}}\\ \multicolumn{1}{c}{F_{\mathrm{z}}}\end{array}\right]}_{\mathrm{C}}=\Delta \mathrm{zB}(t)(\left[\begin{array}{c}{x}_{\mathrm{C}}(t)\hfill \\ \multicolumn{1}{c}{y_{\mathrm{C}}(t)}\\ \multicolumn{1}{c}{z_{\mathrm{C}}(t)}\end{array}\right]-\left[\begin{array}{c}{x}_{\mathrm{C}}(t-T)\hfill \\ \multicolumn{1}{c}{y_{\mathrm{C}}(t-T)}\\ \multicolumn{1}{c}{z_{\mathrm{C}}(t-T)}\end{array}\right])$$

(9)

So far, the cutter displacements and corresponding cutting forces are all calculated in the CCS. However, it should be noted that it is impossible in practice to align the modal frame with every CCS corresponding to each cutter posture when impact test is done on the machine. For the sake of convenience, the impacting directions are selected to be along the machine’s X and Y directions when every rotary axis is at zero position on the machine (see Figure 5). The modal frame will also rotate with the cutter when a certain posture is taken by the cutter. The transformation matrix between CCS and the modal frame will therefore be determined by the machine’s structure.

Figure 5.

Determination of impact testing directions: (a) five-axis machining centre and (b) illustration of machine structure.

For instance, in our experiments, the used machine (model: DMG 80P duo Block) has a non-orthogonal swivel head configuration. The transformation between the modal frame and the workpiece frame is given as

$$T_{\mathrm{mTw}}=T_{\mathrm{tTw}}{T}_{\mathrm{C}}{T}_{135}{T}_{\mathrm{B}}{T}_{45}$$

(10)

where T_tTw is the transformation matrix from the machine table to the workpiece, determined by the way the workpiece is set up on the table and the other matrices are given as

$$T_{\mathrm{C}}=\left[\begin{array}{ccc}c\left(C\right)\hfill & -s\left(C\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{s\left(C\right)}& c\left(C\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right]$$

$$T_{135}=\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& c\left({135}^{\circ }\right)\hfill & -s\left({135}^{\circ }\right)\hfill \\ \multicolumn{1}{c}{0}& s\left({135}^{\circ }\right)\hfill & c\left({135}^{\circ }\right)\hfill \end{array}\right]$$

$$T_{\mathrm{B}}=\left[\begin{array}{ccc}c\left(B\right)\hfill & s\left(B\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{-s\left(B\right)}& c\left(B\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right]$$

$$T_{45}=\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& c\left({45}^{\circ }\right)\hfill & -s\left({45}^{\circ }\right)\hfill \\ \multicolumn{1}{c}{0}& s\left({45}^{\circ }\right)\hfill & c\left({45}^{\circ }\right)\hfill \end{array}\right]$$

Following this, the transformation matrix between the cutter and modal coordinate systems can be obtained as

$$T_{\mathrm{cTm}}=T_{\mathrm{cTw}}{\left(T_{\mathrm{mTw}}\right)}^{-1}$$

(11)

Multiplying T_cTm on both sides of equation (9) gives us

$${\left[\begin{array}{c}{F}_{\mathrm{x}}\hfill \\ \multicolumn{1}{c}{F_{\mathrm{y}}}\\ \multicolumn{1}{c}{F_{\mathrm{z}}}\end{array}\right]}_{\mathrm{M}}=\Delta z{T}_{\mathrm{cTm}}B\left(t\right)(\left[\begin{array}{c}{x}_{\mathrm{C}}(t)\hfill \\ \multicolumn{1}{c}{y_{\mathrm{C}}(t)}\\ \multicolumn{1}{c}{z_{\mathrm{C}}(t)}\end{array}\right]-\left[\begin{array}{c}{x}_{\mathrm{C}}(t-T)\hfill \\ \multicolumn{1}{c}{y_{\mathrm{C}}(t-T)}\\ \multicolumn{1}{c}{z_{\mathrm{C}}(t-T)}\end{array}\right])$$

(12)

The dynamic equation for our model of machining system can be written as

$$\begin{array}{c}\left[\begin{array}{ccc}{m}_{\mathrm{x}}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& m_{\mathrm{y}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & m_{\mathrm{z}}\hfill \end{array}\right]\left[\begin{array}{c}{\stackrel{··}{x}}_{\mathrm{M}}\hfill \\ \multicolumn{1}{c}{{\stackrel{··}{y}}_{\mathrm{M}}}\\ \multicolumn{1}{c}{{\stackrel{··}{z}}_{\mathrm{M}}}\end{array}\right]+\left[\begin{array}{ccc}{c}_{\mathrm{x}}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& c_{\mathrm{y}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & c_{\mathrm{z}}\hfill \end{array}\right]\left[\begin{array}{c}{\stackrel{·}{x}}_{\mathrm{M}}\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{y}}_{\mathrm{M}}}\\ \multicolumn{1}{c}{{\stackrel{·}{z}}_{\mathrm{M}}}\end{array}\right]\hfill \\ +\left[\begin{array}{ccc}{k}_{\mathrm{x}}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& k_{\mathrm{y}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & k_{\mathrm{z}}\hfill \end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{M}}\hfill \\ \multicolumn{1}{c}{y_{\mathrm{M}}}\\ \multicolumn{1}{c}{z_{\mathrm{M}}}\end{array}\right]={\left[\begin{array}{c}{F}_{\mathrm{x}}\hfill \\ \multicolumn{1}{c}{F_{\mathrm{y}}}\\ \multicolumn{1}{c}{F_{\mathrm{z}}}\end{array}\right]}_{\mathrm{M}}\hfill \end{array}$$

(13)

which can be rewritten as

$$\begin{array}{c}M{\stackrel{··}{X}}_{\mathrm{M}}\left(t\right)+C{\stackrel{·}{X}}_{\mathrm{M}}\left(t\right)+K{X}_{\mathrm{M}}\left(t\right)\hfill \\ =\Delta z{T}_{\mathrm{cTm}}B\left(t\right)T_{\mathrm{cTm}}^{-1}(X_{\mathrm{M}}\left(t\right)-X_{\mathrm{M}}(t-T))\hfill \end{array}$$

(14)

where M, C, and K represent inertial, damping ratio, and stiffness matrices of the system, respectively, obtained through impact testing. X_M is a column vector representing the displacement of the cutter in the modal frame. Equation (14) is the general expression for five-axis machining. In reality, as the stiffness in the Z direction is very high, the terms in Z direction are usually dropped.

Stability analysis with FD method

It is obvious that equation (14) is a delay differential equation (DDE). At a specified spindle speed, which corresponds to a specific value of the tooth passing period T in equation (14), the stability of the equation determines whether chatter exists during machining. Over the years, many attempts have been made by researchers to solve for stability condition. The time domain simulation is proposed by Smith and Tlusty. ²⁹ Although it is the most accurate method to calculate the critical value, it is highly time-consuming. Altintas and Budak ³⁰ presented the zeroth order approximation (ZOA) method using the zeroth order Fourier term to approximate the time-varying cutting force coefficients in milling. However, this method cannot accurately predict the critical parameters at low radial immersion milling operations. ³¹ Subsequently, the multi-frequency method was proposed to deal with this problem by extending the ZOA method. ³² However, the order of the Fourier term is still difficult to choose. Semi-discretization (SD) method was proposed to convert the DDE into a series of autonomous ordinary differential equations (ODEs), and then the stability of the system is determined by calculating the eigenvalues of the transition matrix, constructed using the coefficients of the series ODEs.^31,33 Temporal finite element analysis (TFEA) was proposed to deal with highly interrupted milling process. ³⁴ However, this method is only suitable for very low radial immersion milling process. Ding et al. ³⁵ presented a FD method which is similar to the SD method in the sense that both methods approximate the original DDE by a series of ODEs but produces better computational efficiency. In this article, the FD method is used to determine the stability in five-axis milling in our work.

First, the state variable U is defined as

$$U={\left[\begin{array}{cc}{X}_{\mathrm{M}}\left(t\right)\hfill & M{\stackrel{·}{X}}_{\mathrm{M}}\left(t\right)+C{X}_{\mathrm{M}}\left(t\right)/2\hfill \end{array}\right]}^T.$$

With the state variables U, equation (14) can be rewritten as

$$U\left(t\right)=A_0^{\mathrm{u}}U\left(t\right)+A^{\mathrm{u}}\left(t\right)U\left(t\right)+B^{\mathrm{u}}\left(t\right)U(t-T)$$

(15)

where

$$A_0^{\mathrm{u}}=\left[\begin{array}{cc}-M^{-1}C/2\hfill & M^{-1}\hfill \\ \multicolumn{1}{c}{C{M}^{-1}C/4-K}& -C{M}^{-1}/2\hfill \end{array}\right]$$

$$A^{\mathrm{u}}\left(t\right)=\left[\begin{array}{cc}\begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\hfill & \begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\hfill \\ \multicolumn{1}{c}{-\Delta z{T}_{\mathrm{cTm}}B\left(t\right)T_{\mathrm{cTm}}^{-1}}& \begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\hfill \end{array}\right]$$

$$B^{\mathrm{u}}\left(t\right)=\left[\begin{array}{cc}\begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\hfill & \begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\hfill \\ \multicolumn{1}{c}{\Delta z{T}_{\mathrm{cTm}}B\left(t\right)T_{\mathrm{cTm}}^{-1}}& \begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\hfill \end{array}\right]$$

The main idea of FD method is to discrete the time period T into m small time interval such that T = µτ. Based on the direct integration scheme, on arbitrary time interval kτ ≤ t ≤ (k+1)τ, with the initial condition u_k = u(kτ), the response of equation (15) can be obtained via the direct integration scheme as

$$\begin{array}{c}\multicolumn{1}{c}{u\left(t\right)=e^{{\mathrm{A}}_0^{\mathrm{u}}(\mathrm{t}-\mathrm{k}\mathrm{\tau })}u\left(k\tau \right)}\\ +\underset{k\tau }{\overset{t}{\int }}\left\{e^{{\mathrm{A}}_0^{\mathrm{u}}(\mathrm{t}-\mathrm{\xi })}[A^{\mathrm{u}}\left(\xi \right)u\left(\xi \right)+B^{\mathrm{u}}\left(\xi \right)u(\xi -T)]\right\}d\xi \hfill \end{array}$$

(16)

Equivalently, it can be expressed as

$$\begin{array}{c}u(k\tau +t)=e^{{\mathrm{A}}_0^{\mathrm{u}}\mathrm{t}}u\left(k\tau \right)\hfill \\ +\underset{0}{\overset{t}{\int }}\{[A^{\mathrm{u}}(k\tau +t-\xi )u(k\tau +t-\xi )+B^{\mathrm{u}}(k\tau +t-\xi )u(k\tau +t-\xi -T)]\}d\xi \hfill \end{array}$$

(17)

where 0 ≤ t ≤ τ. Then, u_k+1, then u(kτ+τ) can be obtained using

$$\begin{array}{c}{u}_{\mathrm{k}+1}=e^{{\mathrm{A}}_0^{\mathrm{u}}\tau }{u}_{\mathrm{k}}\hfill \\ +\underset{0}{\overset{\tau }{\int }}\{e^{{\mathrm{A}}_0^{\mathrm{u}}\xi }[A^{\mathrm{u}}(k\tau +\tau -\xi )u(k\tau +\tau -\xi )+B^{\mathrm{u}}(k\tau +\tau -\xi )u(k\tau +\tau -\xi -T)]\}d\xi \hfill \end{array}$$

(18)

Equation (18) can be rewritten as

$$\begin{array}{c}\multicolumn{1}{c}{u_{\mathrm{k}+1}=(F_0+F_{0,1})u_{\mathrm{k}}+F_{\mathrm{k}+1}{u}_{\mathrm{k}+1}}\\ +F_{\mathrm{m}-1}{u}_{\mathrm{k}+1-\mathrm{m}}+F_{\mathrm{m}}{u}_{\mathrm{k}-\mathrm{m}}\hfill \end{array}$$

(19)

where

$$\begin{array}{c}{F}_0={\Phi }_0\hfill \\ F_{0,1}=({\Phi }_2/\tau ){\left(A_0^{\mathrm{u}}\right)}^{\left(\mathrm{k}\right)}+({\Phi }_3/\tau ){\left(A_1^{\mathrm{u}}\right)}^{\left(\mathrm{k}\right)}\hfill \\ F_{\mathrm{k}+1}=({\Phi }_1-{\Phi }_2/\tau ){\left(A_0^{\mathrm{u}}\right)}^{\left(\mathrm{k}\right)}+({\Phi }_2-{\Phi }_3/\tau ){\left(A_1^{\mathrm{u}}\right)}^{\left(\mathrm{k}\right)}\hfill \\ F_{\mathrm{m}-1}=({\Phi }_1-{\Phi }_2/\tau ){\left(B_0^{\mathrm{u}}\right)}^{\left(\mathrm{k}\right)}+({\Phi }_2-{\Phi }_3/\tau ){\left(B_1^{\mathrm{u}}\right)}^{\left(\mathrm{k}\right)}\hfill \\ F_{\mathrm{m}}=({\Phi }_2/\tau ){\left(B_0^{\mathrm{u}}\right)}^{\left(\mathrm{k}\right)}+({\Phi }_3/\tau ){\left(B_1^{\mathrm{u}}\right)}^{\left(\mathrm{k}\right)}\hfill \\ {\Phi }_0=e^{{\mathrm{A}}_0^{\mathrm{u}}\tau };{\Phi }_1={\left(A_0^{\mathrm{u}}\right)}^{-1}({\Phi }_0-I)\hfill \\ {\Phi }_2={\left(A_0^{\mathrm{u}}\right)}^{-1}(\tau {\Phi }_0-{\Phi }_1)\hfill \\ {\Phi }_3={\left(A_0^{\mathrm{u}}\right)}^{-1}({\tau }^2{\Phi }_0-2{\Phi }_2)\hfill \\ {\left(A_0^{\mathrm{u}}\right)}^{\left(\mathrm{k}\right)}=A_{\mathrm{k}+1}^{\mathrm{u}};{\left(A_1^{\mathrm{u}}\right)}^{\left(\mathrm{k}\right)}=(A_{\mathrm{k}}^{\mathrm{u}}-A_{\mathrm{k}+1}^{\mathrm{u}})/\tau \hfill \\ {\left(B_0^{\mathrm{u}}\right)}^{\left(\mathrm{k}\right)}=B_{\mathrm{k}+1}^{\mathrm{u}};{\left(B_1^{\mathrm{u}}\right)}^{\left(\mathrm{k}\right)}=(B_{\mathrm{k}}^{\mathrm{u}}-B_{\mathrm{k}+1}^{\mathrm{u}})/\tau \hfill \end{array}$$

If the matrix [I−F_k+1] is non-singular, u_k+1 can be expressed in the following explicit form

$$\begin{array}{c}\multicolumn{1}{c}{u_{\mathrm{k}+1}={[I-F_{\mathrm{k}+1}]}^{-1}(F_0+F_{0,1})u_{\mathrm{k}}}\\ +{[I-F_{\mathrm{k}+1}]}^{-1}{F}_{\mathrm{m}-1}{u}_{\mathrm{k}+1-\mathrm{m}}\hfill \\ +{[I-F_{\mathrm{k}+1}]}^{-1}{F}_{\mathrm{m}}{u}_{\mathrm{k}-\mathrm{m}}\hfill \end{array}$$

(20)

In the case [I−F_k+1] is singular, u_k+1 can be obtained using

$$u_{\mathrm{k}+1}=(F_0+F_{0,2})u_{\mathrm{k}}+F_{\mathrm{m}-1}{u}_{\mathrm{k}+1-\mathrm{m}}+F_{\mathrm{m}}{u}_{\mathrm{k}-\mathrm{m}}$$

(21)

where

$$F_{0,2}={\left({\Phi }_1{A}_0^{\mathrm{u}}\right)}^{(\mathrm{k})}+{\left({\Phi }_2{A}_1^{\mathrm{u}}\right)}^{(\mathrm{k})}$$

A discrete map can be defined as

$$y_{\mathrm{k}+1}=D_{\mathrm{k}}{y}_{\mathrm{k}}$$

(22)

where the vector y_k denotes

$$y_{\mathrm{k}}=\mathrm{col}(u_{\mathrm{k}}{u}_{\mathrm{k}-1}\dots u_{\mathrm{k}+1-\mathrm{m}}{u}_{\mathrm{k}-\mathrm{m}})$$

(23)

D_k is defined as

$$D_{\mathrm{k}}=\left[\begin{array}{ccccccc}{[I-F_{\mathrm{k}+1}]}^{-1}(F_0+F_{0,1})\hfill & 0\hfill & 0\hfill & \dots \hfill & 0\hfill & {[I-F_{\mathrm{k}+1}]}^{-1}{F}_{\mathrm{m}-1}\hfill & {[I-F_{\mathrm{k}+1}]}^{-1}{F}_{\mathrm{m}}\hfill \\ \multicolumn{1}{c}{I}& 0\hfill & 0\hfill & \dots \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& I\hfill & 0\hfill & \dots \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\dots }& \dots \hfill & \dots \hfill & \dots \hfill & \dots \hfill & \dots \hfill & \dots \hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & \dots \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & \dots \hfill & I\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & \dots \hfill & 0\hfill & I\hfill & 0\hfill \end{array}\right]$$

(24)

Now, the transition matrix Φ over one periodic time interval can be constructed using the sequence of discrete maps D_k, (k = 0, …, m−1), that is

$$y_{\mathrm{m}}=\Phi y_0$$

(25)

where Φ is defined by

$$\Phi =D_{\mathrm{m}-1}{D}_{\mathrm{m}-2}\dots D_1{D}_0$$

(26)

According to the Floquet theory, the stability of the system can be determined in the following way: if the modulus of all the eigenvalues λ_Φ of the transition matrix Φ is less than 1, the system is stable, otherwise, it is unstable.

Construction and application of PSGs

With the proposed methods, given a combination of spindle speed, depth of cut, and cutter posture, whether machining will be stable can be determined. Spindle speed will be fixed by recommended cutting speed, which leaves two important parameters during tool-path generation, that is, depth of cut and cutter orientation. For any depth of cut, the feasible range for cutter orientations can be divided into stable and unstable regions. Postures should only come from stable regions, which is the idea behind PSG.

Generation of PSG at a certain depth of cut relies on a sampling algorithm. The boundary between stable and unstable machining will be the set of cutter orientations for which the ∣λ_Φ∣ of the transition matrix Φ are equal to 1. For cutter postures with ∣λ_Φ∣ larger than 1, the postures will be unstable and for those with ∣λ_Φ∣ smaller than 1, the postures will be stable. ³⁵

The chattering behaviour of a certain machining process is heavily affected by the engagement area between the cutter and the workpiece. In this article, experiments are carried out at the simple scenario of slotting on a plane (see Figure 2(a)). The engagement information for such a scenario is easy to obtain. The tool of PSG construction could be applied to freeform machining scenarios as well, as long as the cutter–workpiece engagement condition is provided. The machine tool used is a DMG 80P Duo Block five-axis machining centre (see Figure 6). The material to be machined is AL7075 and cutter is a Mitsubishi carbide cutter with 10 mm diameter. Spindle speed is set at 4800 r/min. The cutting force coefficients and modal parameters are provided in Tables 1 and 2, respectively.

Figure 6.

Setup of cutting experiments and impact test: (a) for cutting experiments and (b) for impact test.

Table 1.

Calibrated cutting force coefficients.

Kte (N/mm)	Ktc (N/mm2)	Kre (N/mm)	Krc (N/mm2)	Kae (N/mm)	Kac (N/mm2)
11.4	1079.2	30.2	373.03	−2.3	293.3

Table 2.

Modal parameters.

Natural frequency (Hz)		Process damping (%)		Stiffness (N/m2)
X M	Y M	X M	Y M	X M	Y M
2213.75	1613.75	0.016713	0.013928	6012175.53	6292140.24

Following the sampling algorithm, PSG is first generated for a cutting depth of 0.5 mm as shown in Figure 7(a). The sampling is conducted at 1° for both lead and tilt angles. The range of sampling for both lead and tilt angles are set as 60°, large enough to cover the range of practical machining. The horizontal and vertical axes of in Figure 7 represent the tilt and lead angles, the combination of which determines a cutter posture. Thus, the rectangular area defined by these two axes gives the allowable range for cutter postures. The contours on the graph give the postures that share the same value of ∣λ_Φ∣. The contour with ∣λ_Φ∣ = 1, as explained before, is the boundary between stable and chatter machining, whereas contours with value of ∣λ_Φ∣ smaller or larger than 1 correspond to stable and chatter postures, respectively. The PSGs are also generated for the cutting depth of 0.6 and 1 mm as well, as shown in Figure 7(b) and (c). It can be seen clearly that the unstable region has enlarged in accordance with our expectations. In Figure 8, the PSGs at these three different depths of cut are summarized into a complete PSG. Each contour represents the stability boundary at the depth of cut indicated by the number on it, taken as the line with ∣λ_Φ∣ = 1 from the corresponding PSG. The benefit of such a summarized PSG is two-fold: first, it demonstrates clearly how the stable region enlarges with the reduction in depth of cut; second, for depth of cut other than the three indicated values, the stability region of the next larger depth of cut can be chosen as a conservative measure. The complete work flow for construction of PSGs is illustrated in Figure 9.

Figure 7.

PSGs for slotting on a flat surface at spindle speed 4800 r/min: (a) depth of cut = 0.5 mm, (b) depth of cut = 0.6 mm, and (c) depth of cut = 1.0 mm.

Figure 8.

Summarized PSG corresponding to Figure 7(a)–(c).

Figure 9.

Flowchart for PSG generation.

To test the validity of generated PSGs, a series of machining experiments have been carried out under indexed machining mode, with which the cutter’s posture is pre-set and will remain the same during one pass of machining, leaving a consistent machined surface. The first set of the experiments was carried out with a depth of cut at 0.5 mm. The results are shown in Figure 10, where the cutter postures experimented with have been overlaid over the PSG. It is obvious that this is the same PSG as shown in Figure 7(a) but showing only the boundary of the chatter zone, that is, the contour with ∣λ_Φ∣ = 1. Microscopic images for some of the finished machining surfaces are also provided in Figure 10. It can be seen that a high level of conformity has been achieved between our prediction and the experimental results.

Figure 10.

Stability comparisons between predicted and experimental measurements.

Another experiment is carried out under the simultaneous machining mode as shown in Figure 11 when the cutter posture is continuously adjusted from (L: 0°, T: 0°) to (L: 60°, T: 0°) smoothly with a depth of cut 0.5 mm over 21 CLs with an interval of 5 mm. The interpolation from the starting to the ending cutter posture is conducted by commercial CAM software (Unigraphics NX 8.5). At the intermediate CLs, the lead and tilt angles have been retrieved as given in Figure 11.

Figure 11.

Machined surface under varying cutter posture with depth of cut 0.5 mm.

According to the PSG, the machining will enter chatter zone after the 10th CL (corresponding posture: L: 26.7°, T: 0°). This is generally in accordance with the microscopic image of the machined surface, in which chatter emerges at around 22° as shown in Figure 11, which demonstrates the effectiveness of our algorithm. In addition, according to the PSG, the machining will enter stable area after the 19th CL (corresponding posture: L: 54.8°, T: 0°). This is also in accordance with the microscopic image of the machined surface, in which chatter disappear at around 58° as shown in Figure 11. The minor difference between the predicted key posture and experimental key posture may be caused by the measurement deviation for modal parameters and ignored influence of cutting force along the cutter spindle.

Conclusion

In this article, PSG is proposed to guide the determination of chatter-free cutter postures during five-axis machining. A procedure to obtain the PSG has also been established. At first, a dynamic model of the machining system under cyclic cutting force is established, in which the influences of the cutter posture, depth of cut, and spindle speed have all been accounted for. The stability of the machining process is then derived from the transition matrix of the system. By identifying cutter postures that will lead to borderline stability under a certain machining condition, the feasible range of cutter orientations will be divided into stable and unstable zones, which makes the PSG. Machining experiments have been conducted and the results verified the validity of PSG under both indexed machining and the simultaneous machining. The development of PSG is meaningful for the integration of process planning and machining dynamics. Our future work will focus on applying our algorithm for freeform surface machining with varying cutter–workpiece engagement conditions.