Prediction of Chatter Stability in Bull-Nose End Milling of Thin-Walled Cylindrical Parts Using Layered Cutting Force Coefficients

Abstract

During the bull-nose end milling operations of thin-walled structures, chatter usually occurs and adversely affects cutter performance, finished surface quality, and production efficiency. To accurately predict chatter stability, a suitable dynamic model with effective system parameters is required. In this article, a three-degree-of-freedom (3-DOF) dynamic model is developed to analyze the milling stability of the thin-walled cylinders, in which the dynamics of the bull-nose end mill along the x-axis and y-axis directions and the dynamic of the workpiece along the z-axis direction are taken into account. Then, the cutter-workpiece engagement (CWE) is extracted by employing a slice-intersection-based approach. And the layered cutting force coefficients are identified by considering the influences of varying cutter diameters on the cutting speed. Thereafter, the semi-discretization method (SDM) is adopted to compute the stability lobe diagram (SLD). In the end, a group of milling tests are carried out on a thin-walled cylinder to validate the accuracy and reliability of the proposed model, and the results show that the model predictions agree well with the experimental data.

1. Introduction

Thin-walled structures are commonly applied in the aerospace, marine, and automotive industries owing to their incomparable advantages, such as lightweight and high specific strength. Considering the high material removal rate and sleek finished surface, bull-nose end milling has become a popular operation for the manufacturing of thin-walled parts, such as aero-engine blades, impellers, and casings [1,2,3]. These parts have complex and thin-walled structures and often bring many machining difficulties and challenges. Chatter is essentially a self-excited vibration and usually appears in the milling process. It not only pulls down the material removal rate but also results in a coarse finished surface, premature tool wear, excessive cutting noise, and so on [4,5]. Chatter has long been one of the main problems during the milling operation. Thus, to avoid chatter, the milling operation must be performed in a steady state [6]. Furthermore, the calculation of SLD has become and still is an important research topic in high-quality machining.

Presently, two primary mathematical solutions are employed to construct the SLD, namely the analytical and numerical approaches. In terms of the analytical approaches, Altintas and Budak [7] proposed a single-frequency solution (SFS) on the basis of the truncated Fourier series of the cyclical cutting forces. The SFS is also known as the zeroth-order approximation (ZOA) approach, which can plot the SLD efficiently. However, it fails to compute the extra stability boundaries and the doubling bifurcations. For solving this question, Budak and Altintas [8,9] subsequently reported a multi-frequency solution (MFS), where the higher-order parts of the Fourier series are adopted. Merdol and Altintas [10] expanded the application of the MFS to machining conditions with small radial depths of cut. The disadvantage of the MFS is that it consumes more computing time. In comparison with the analytical approaches, more investigations have concentrated on numerical approaches [11]. Insperger and Stépán [12] reported an important SDM to evaluate milling stability. In this approach, the zeroth-order piecewise constant function is employed to approach the delayed part. To achieve better solution accuracy and quicker convergence rates, many high-order SDMs were developed. Insperger et al. [13] took into account the first-order approximation of the delayed term (namely the first SDM). Jiang et al. [14] investigated the second-order approximation of the delayed term (namely the second SDM). Yan et al. [15] suggested a third SDM and a fourth SDM. On the other hand, Ding et al. [16] pioneered the influential full-discretization method (FDM), which has a higher computational efficiency. In the FDM, in addition to the delayed term, the state term is also approximated by employing linear interpolation. Subsequently, many other approaches inspired by the FDM have been proposed, just like the 2nd FDM [17], the 3rd FDM [18], the updated FDM (UFDM) [19], the 2nd UFDM [20], and the 3rd UFDM [21].

During the bull-nose end finishing process of thin-walled or thin-floor parts, the influence of the workpiece dynamics cannot be neglected [22,23,24]. Campa et al. [25,26] modeled the dynamics of the thin-floor structure as a 1-DOF system and predicted the SLD in bull-nose end finish milling. Based on this study, the influence of the material removal on the bull-nose end milling stability is analyzed, and the three-dimensional (3D) SLD along the toolpath is calculated [27]. Later, a half-active magnetorheological damper is designed to dampen regenerative chatter for the bull-nose end finishing of thin floors [28]. Zhou et al. [29] simplified the dynamics of the aero-engine casings into a 1-DOF model along the workpiece’s radial direction. The model only takes into account the workpiece flexibility and is able to evaluate the SLD for different small lead angles. The above models only consider the dynamics of the workpiece, while the cutter is assumed to have a rigid body. It is essential to establish a multi-DOF milling system when neither the tool nor the workpiece dynamics can be ignored [30,31,32]. However, the 3-DOF stability prediction model for bull-nose end milling has been rarely explored. Especially for the cylindrical part, which has a ring-shaped and closed geometry structure. Yan et al. [33] modeled a 3-DOF milling system by taking into account the tool dynamics along the x and y directions and the workpiece dynamics along the z direction. Dang et al. [34] suggested a 3-DOF milling stability model for thin-floor components, in which the static and dynamic deflections of the tool and the workpiece are considered together.

In the chatter stability model, the cutting force coefficients are very important for the calculation of SLD [3,35]. Currently, most models ignore the effect of the cutter diameter variation and treat the cutting force coefficients as constants for the bull-nose or ball end milling [36], which can result in inaccurate stability predictions. Based on previous research [29], this article presents a 3-DOF chatter stability model that takes into account the effects of cutter diameter variation on the cutting speed. In Section 2, the 3-DOF milling system for the bull-nose end milling of a thin-walled cylinder is established. The dynamic cutting force model, CWE calculation, identification of layered cutting force coefficients, and SLD calculation with SDM are also conducted in this section. In Section 3, a group of milling tests are performed to validate the accuracy and reliability of the 3-DOF chatter stability model. Finally, some conclusions are proposed.

2. 3-DOF Chatter Stability Model

2.1. Milling System for a Thin-Walled Cylinder

During the bull-nose end milling process of thin-walled cylinders, both the cutter and workpiece are regarded as flexible. The dynamic milling system for a thin-walled cylinder is demonstrated in Figure 1, which can be delineated by a 3-DOF mass-spring-damper model. In the x-axis and y-axis directions, a 2-DOF dynamic model is employed to delineate the cutter, while an absolutely rigid body is applied to assume the workpiece. In the z-axis direction, a 1-DOF dynamic model is employed to delineate the workpiece, while an absolutely rigid body is applied to assume the cutter. O-xyz stands for the cutter coordinate system (CCS). The origin O is centered at the end of the cutter. The workpiece feed direction is specified as the y-axis. The cutter axis direction is specified as the z-axis. The right-hand rule is employed to identify the x-axis. Based on the previous study [29], the cutting contact point P is identified according to the eccentric amount from the cutter axis to the centerline of the workpiece.

The dynamic equation of the 3-DOF bull-nose end milling system is established as follows:

$$\left[\begin{array}{ccc}{m}_x& 0& 0\\ 0& m_y& 0\\ 0& 0& m_z\end{array}\right]\left[\begin{array}{c}{\ddot{x}}_c(t)\\ {\ddot{y}}_c(t)\\ {\ddot{z}}_w(t)\end{array}\right]+\left[\begin{array}{ccc}{c}_x& 0& 0\\ 0& c_y& 0\\ 0& 0& c_z\end{array}\right]\left[\begin{array}{c}{\dot{x}}_c(t)\\ {\dot{y}}_c(t)\\ {\dot{z}}_w(t)\end{array}\right]+\left[\begin{array}{ccc}{k}_x& 0& 0\\ 0& k_y& 0\\ 0& 0& k_z\end{array}\right]\left[\begin{array}{c}{x}_c(t)\\ y_c(t)\\ z_w(t)\end{array}\right]=\left[\begin{array}{c}{F}_x\\ F_y\\ -F_z\end{array}\right]$$

(1)

where m_x, c_x, and k_x denote the modal mass, damping, and stiffness of the cutter, respectively, along the x-axis direction. m_y, c_y, and k_y denote the modal mass, damping, and stiffness of the cutter, respectively, along the y-axis direction. m_z, c_z, and k_z denote the modal mass, damping, and stiffness of the workpiece, respectively, along the z-axis direction. x_c(t) and y_c(t) denote the cutter displacements along the x-axis and y-axis directions, respectively. z_w(t) denotes the workpiece displacement along the z-axis direction. F_x, F_y, and F_z denote the cutting forces along the x-axis, y-axis, and z-axis directions, respectively.

2.2. Dynamic Cutting Force Modeling for Bull-Nose End Milling

Figure 2 demonstrates a bull-nose end cutter, which is characterized by two main parameters: the fillet radius R_f and the cutter diameter D.

For the bull-nose end milling process, the chip thickness can be distributed between the static and dynamic components. The static component is disregarded because it has no influence on milling stability [37]. The dynamic component of cutting-edge j is determined as follows:

$$h_d({\phi }_j,\kappa )=\left[\begin{array}{c}\sin {\phi }_j\sin \kappa \hspace{1em}\cos {\phi }_j\sin \kappa \hspace{1em}-\cos \kappa \end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\end{array}\right]$$

(2)

where

$$\begin{array}{c}{\phi }_j(z)=\phi +(j-1){\phi }_p-\psi (z)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{2\pi n}{60}t+(j-1)\frac{2\pi }{N}-\frac{2z\tan \beta }{D}\end{array}$$

(3)

$$\begin{array}{c}\kappa (z)=\arccos \frac{R_f-z}{R_f},\hspace{1em}(z<R_f)\end{array}$$

(4)

$$\left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\end{array}\right]=\left[\begin{array}{c}{x}_c(t)-x_c(t-T)\\ y_c(t)-y_c(t-T)\\ z_w(t-T)-z_w(t)\end{array}\right]$$

(5)

where φ_j denotes the radial immersion angle of cutting-edge j. φ denotes the rotation angle of the cutting edge. φ_p denotes the pitch angle of the cutter. ψ(z) denotes the radial lag angle caused by the local helix angle. n denotes the spindle speed. t denotes the time parameter. N denotes the number of cutting edges. β denotes the constant helix angle of the cutter. κ denotes the axial immersion angle. Δx, Δy, and Δz denote the x-axis, y-axis, and z-axis dynamic displacements in the CCS. x_c(t) and y_c(t) denote the dynamic displacements of the cutter at the present time t. x_c(t − T) and y_c(t − T) denote the dynamic displacements of the cutter at previous tooth periods t − T. z_w(t) denotes the dynamic displacement of the workpiece at the present time. t. z_w(t − T) denotes the dynamic displacements of the workpiece at previous tooth periods t − T.

The dynamic cutting forces in the CCS can be established as follows:

$$\left[\begin{array}{c}{F}_x\\ F_y\\ F_z\end{array}\right]={\displaystyle \sum _{j=1}^N{\displaystyle \sum _{l=1}^{M}g({\phi }_j)T({\phi }_j,\kappa )\left[\begin{array}{c}{K}_{t,l}\\ K_{r,l}\\ K_{a,l}\end{array}\right]}}{h}_d({\phi }_j,\kappa )\frac{dz}{\sin \kappa }$$

(6)

where

$$g({\phi }_j)=\left\{\begin{array}{c}1\hspace{1em}{\phi }_{st,l}<{\phi }_j<{\phi }_{ex,l}\\ 0\hspace{1em}\hspace{1em}\mathrm{otherwise}\end{array}\right.$$

(7)

$$T({\phi }_j,\kappa )=\left[\begin{array}{ccc}-\cos {\phi }_j& -\sin {\phi }_j\sin \kappa & -\sin {\phi }_j\cos \kappa \\ \sin {\phi }_j& -\cos {\phi }_j\sin \kappa & -\cos {\phi }_j\cos \kappa \\ \hspace{1em}0\hfill & \cos \kappa & -\sin \kappa \end{array}\right]$$

(8)

where M denotes the number of axial cutting disk elements. g(φ_j) denotes a unit step function to decide whether the flute is in or out of the cutting state. φ_st,l and φ_ex,l denote the start and exit angles of the cutting disk element l, respectively. T(φ_j, κ) denotes the rotational transformation matrix. K_t,l, K_r,l, and K_a,l denote the layered tangential, radial, and axial cutting force coefficients on the cutting disk element l, respectively.

Equation (6) delineates a linear cutting force model where the dynamic chip thickness is not correlated with the feed rate [38]. It can be converted into the following matrix form:

$$\left[\begin{array}{c}{F}_x\\ F_y\\ F_z\end{array}\right]=H(t)\left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\end{array}\right]=\left[\begin{array}{ccc}{h}_{xx}& h_{xy}& h_{xz}\\ h_{yx}& h_{yy}& h_{yz}\\ h_{zx}& h_{zy}& h_{zz}\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\end{array}\right]$$

(9)

where

$$\left\{\begin{array}{l}{h}_{xx}={\displaystyle \sum _{j=1}^N{\displaystyle \sum _{l=1}^{M}g({\phi }_j)}}\left[-K_{t,l}\sin {\phi }_j\cos {\phi }_j\sin \kappa -K_{r,l}{\sin }^2{\phi }_j{\sin }^2\kappa -K_{a,l}{\sin }^2{\phi }_j\sin \kappa \cos \kappa \right]\frac{dz}{\sin \kappa }\\ h_{xy}={\displaystyle \sum _{j=1}^N{\displaystyle \sum _{l=1}^{M}g({\phi }_j)}}\left[-K_{t,l}{\cos }^2{\phi }_j\sin \kappa -K_{r,l}\sin {\phi }_j\cos {\phi }_j{\sin }^2\kappa -K_{a,l}\sin {\phi }_j\cos {\phi }_j\sin \kappa \cos \kappa \right]\frac{dz}{\sin \kappa }\\ h_{xz}={\displaystyle \sum _{j=1}^N{\displaystyle \sum _{l=1}^{M}g({\phi }_j)\left[K_{t,l}\cos {\phi }_j\cos \kappa +K_{r,l}\sin {\phi }_j\sin \kappa \cos \kappa +K_{a,l}\sin {\phi }_j{\cos }^2\kappa \right]}}\frac{dz}{\sin \kappa }\\ h_{yx}={\displaystyle \sum _{j=1}^N{\displaystyle \sum _{l=1}^{M}g({\phi }_j)}}\left[K_{t,l}{\sin }^2{\phi }_j\sin \kappa -K_{r,l}\sin {\phi }_j\cos {\phi }_j{\sin }^2\kappa -K_{a,l}\sin {\phi }_j\cos {\phi }_j\sin \kappa \cos \kappa \right]\frac{dz}{\sin \kappa }\\ h_{yy}={\displaystyle \sum _{j=1}^N{\displaystyle \sum _{l=1}^{M}g({\phi }_j)}}\left[K_{t,l}\sin {\phi }_j\cos {\phi }_j\sin \kappa -K_{r,l}{\cos }^2{\phi }_j{\sin }^2\kappa -K_{a,l}{\cos }^2{\phi }_j\sin \kappa \cos \kappa \right]\frac{dz}{\sin \kappa }\\ h_{yz}={\displaystyle \sum _{j=1}^N{\displaystyle \sum _{l=1}^{M}g({\phi }_j)}}\left[-K_{t,l}\sin {\phi }_j\cos \kappa +K_{r,l}\cos {\phi }_j\sin \kappa \cos \kappa +K_{a,l}\cos {\phi }_j{\cos }^2\kappa \right]\frac{dz}{\sin \kappa }\\ h_{zx}={\displaystyle \sum _{j=1}^N{\displaystyle \sum _{l=1}^{M}g({\phi }_j)}}\left[K_{r,l}\sin {\phi }_j\sin \kappa \cos \kappa -K_{a,l}\sin {\phi }_j{\sin }^2\kappa \right]\frac{dz}{\sin \kappa }\\ h_{zy}={\displaystyle \sum _{j=1}^N{\displaystyle \sum _{l=1}^{M}g({\phi }_j)}}\left[K_{r,l}\cos {\phi }_j\sin \kappa \cos \kappa -K_{a,l}\cos {\phi }_j{\sin }^2\kappa \right]\frac{dz}{\sin \kappa }\\ h_{zz}={\displaystyle \sum _{j=1}^N{\displaystyle \sum _{l=1}^{M}g({\phi }_j)}}\left[-K_{r,l}{\cos }^2\kappa +K_{a,l}\sin \kappa \cos \kappa \right]\frac{dz}{\sin \kappa }\end{array}\right.$$

(10)

where h_xx, h_xy, h_xz, h_yx, h_yy, h_yz, h_zx, h_zy, and h_zz are the matrix coefficients.

2.3. Calculation of CWE

CWE signifies the instantaneous contact interface while a bull-nose end mill is cutting a thin-walled cylinder. It is one of the foundations for modeling the dynamics of the milling system. During the semi-finish milling or finish milling of thin-walled cylinders, the cutting edges on the toroidal surface of a bull-nose end cutter are applied to end milling, giving rise to a smaller axial depth of cut and a larger machined strip width [39,40].

As shown in Figure 3, the CWE for the bull-nose end milling of thin-walled cylindrical parts consists of three boundaries. Boundary 1 is the cutter profile curve in the current feed direction view. Boundary 2 is the intersection curve between the uncut surface of the workpiece and the toroidal surface of the tool. Boundary 3 is the intersection curve between the workpiece transition surface left by the previous cutting path (CP) and the toroidal surface of the tool. Furthermore, the workpiece transition surface can be generated by a series of boundaries 1, which correspond to each cutting contact point along the previous cutting path. Therefore, the calculation of boundary 1 is the key to determining the CWE at an arbitrary cutting contact point.

As shown in Figure 4, a slice S_i can be obtained by offsetting the workpiece’s machined surface S_m outward. Furthermore, the intersection curve C_i between the slice S_i and the cutter’s toroidal surface S_c can also be calculated. The boundary points P_i and Q_i are determined by employing the closest lines from the two auxiliary planes P_r1 and P_r2 to the intersection curve C_i. The normal directions of the two reference planes P_r1 and P_r2 are simultaneously orthogonal to the workpiece feed direction and the cutter axis direction. The two boundary points are located on boundary 1 and represent the maximum view width of the intersection curve C_i in the workpiece feed direction view. Then, the CWE boundaries can be identified using the intersection curve C_a and a spline obtained by fitting a series of boundary points. The surface trimming tool can be used to extract the engagement geometry for the following cut condition.

The flowchart of the CWE calculation algorithm is shown in Figure 5. And the detailed steps of the CWE calculation algorithm are described in the following:

• Input the machined surface of the workpiece S_m and the toroidal surface of the tool S_c, and determine the cutting contact point P_c according to the cutter axis orientation and the milling line position. The machined surface of the workpiece S_m and the toroidal surface of the tool S_c are tangent to each other at the cutting contact point P_c;
• Input the total number of slices N_s and the machining allowance h, and set the initial value of an index variable i to 1;
• Offset the workpiece’s machined surface S_m towards the outside by the distance Δh = i∙h/N_s to obtain a slice S_i;
• Calculate the intersection curve C_i between the slice S_i and the toroidal surface of the tool S_c. The intersection curve C_i is essentially boundary 2 with the corresponding machining allowance;
• Calculate the minimum distance points P_i and Q_i from the intersection curve C_i to the reference planes P_r1 and P_r2, and store the minimum distance points into the set of boundary points {P};
• Go to Step (3) and repeat the process until the index variable i is greater than the total number of slices N_s;
• Fit the set of boundary points {P} of the CWE as a spline curve B;
• Trim the cutter’s toroidal surface S_c using the spline curve B and the intersection curve C_a to extract the CWE in bull-nose end slotting of the thin-walled cylindrical workpiece. The CWEs under the following cut conditions can be obtained by trimming operations.

2.4. Identification of Layered Cutting Force Coefficients

For the toroidal surface of a bull-nose end cutter, the diameter varies continuously along the cutter’s axis direction. The cutting force coefficients are identified by the fact that the measured average cutting forces per tooth are equal to the predicted ones [41,42]. The predicted average cutting forces acting on the lth cutting disk element are given as follows:

$$\begin{array}{c}{\overline{F}}_{q,l}=\frac{1}{{\phi }_p}{\displaystyle {\int }_{{\phi }_{st,l}}^{{\phi }_{ex,l}}{\displaystyle {\int }_{z_{1,l}}^{z_{2,l}}d{F}_{q,l}(\phi ,z)d\phi }},\hspace{1em}(q=x,y,z\end{array})$$

(11)

where z_1,l and z_2,l are the contact boundaries of the lth cutting disk element within the cut.

In the bull-nose end milling experiments, the measured cutting forces are obtained at various axial depths of the cut. The differences in average cutting forces between two adjacent axial depths of cut are calculated so as to acquire the average cutting forces corresponding to each cutting disk element:

$$\begin{array}{cc}{\overline{F}}_{q,l}={\overline{F}}_q(z_l)-{\overline{F}}_q(z_{l-1}),& (q=x,y,z; l=1,2,3,\dots ,N_l)\end{array}$$

(12)

where z_l denotes the lth axial depth of cut and z₀ = 0, representing non-cutting. N_l is the number of layers.

The instantaneous cutting forces generated on the lth axial cutting disk element at the radial immersion angle φ are expressed as follows:

$$\left[\begin{array}{c}{F}_{x,l}\\ F_{y,l}\\ F_{z,l}\end{array}\right]=\frac{f_t}{2}A\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}\right]$$

(13)

where

$$A=\left[\begin{array}{ccc}-K_{t,l}\sin 2{\phi }_j& -2K_{r,l}{\sin }^2{\phi }_j& -2K_{a,l}{\sin }^2{\phi }_j\\ 2K_{t,l}{\sin }^2{\phi }_j& -2K_{r,l}\sin 2{\phi }_j& -K_{a,l}\sin 2{\phi }_j\\ 0& -2K_{a,l}\sin {\phi }_j& 2K_{r,l}\sin {\phi }_j\end{array}\right]$$

(14)

$$u_1=z\left|\begin{array}{l}{z}_{2,l}\\ z_{1,l}\end{array}\right.$$

(15)

$$u_2=-\frac{R_f}{2}\left(E\sqrt{1-E^2}+\arcsin E\right)\left|\begin{array}{l}{z}_{2,l}\\ z_{1,l}\end{array}\right.$$

(16)

$$u_3=\frac{1}{R_f}\left(R_{f}z-\frac{z^2}{2}\right)\left|\begin{array}{l}{z}_{2,l}\\ z_{1,l}\end{array}\right.$$

(17)

$$E=\frac{R_f-z}{R_f}$$

(18)

where f_t is the normal feed per tooth.

The instantaneous cutting forces are averaged within the radial immersion interval and expressed as follows:

$$\left[\begin{array}{c}{\overline{F}}_{x,l}\\ {\overline{F}}_{y,l}\\ {\overline{F}}_{z,l}\end{array}\right]=\frac{f_t}{{\phi }_p}B\left[\begin{array}{c}{K}_{t,l}\\ K_{r,l}\\ K_{a,l}\end{array}\right]$$

(19)

where

$$B=\left[\begin{array}{ccc}{v}_3{u}_1& (v_2-v_1)u_2& (v_2-v_1)u_3\\ (v_1-v_2)u_1& v_3{u}_2& v_3{u}_3\\ 0& -v_4{u}_3& v_4{u}_2\end{array}\right]$$

(20)

$$v_1=\frac{1}{2}\phi \left|\begin{array}{l}{\phi }_{ex,l}\\ {\phi }_{st,l}\end{array}\right.$$

(21)

$$v_2=\frac{1}{4}\sin (2\phi )\left|\begin{array}{l}{\phi }_{ex,l}\\ {\phi }_{st,l}\end{array}\right.$$

(22)

$$v_3=\frac{1}{4}\cos (2\phi )\left|\begin{array}{l}{\phi }_{ex,l}\\ {\phi }_{st,l}\end{array}\right.$$

(23)

$$v_4=\cos \phi \left|\begin{array}{l}{\phi }_{ex,l}\\ {\phi }_{st,l}\end{array}\right.$$

(24)

Therefore, the cutting force coefficients of the lth axial cutting disk element are determined as follows:

$$\left\{\begin{array}{l}{K}_{t,l}=\frac{2\pi }{f_{t}N{u}_1}\cdot \frac{v_3{\overline{F}}_{x,l}-(v_2-v_1){\overline{F}}_{y,l}}{v_3^2+{(v_2-v_1)}^2}\\ K_{r,l}=\frac{2\pi }{f_{t}N(u_2^2+u_3^2)}\cdot \left\{\frac{u_2[(v_2-v_1){\overline{F}}_{x,l}+v_3{\overline{F}}_{y,l}]}{v_3^2+{(v_2-v_1)}^2}-\frac{u_3{\overline{F}}_{z,l}}{v_4}\right\}\\ K_{a,l}=\frac{2\pi }{f_{t}N(u_2^2+u_3^2)}\cdot \left\{\frac{u_3[(v_2-v_1){\overline{F}}_{x,l}+v_3{\overline{F}}_{y,l}]}{v_3^2+{(v_2-v_1)}^2}+\frac{u_2{\overline{F}}_{z,l}}{v_4}\right\}\end{array}\right.$$

(25)

2.5. Calculation of SLD with SDM

Equation (1) can be represented by the matrix equation as follows:

$$M\ddot{q}(t)+C\dot{q}(t)+Kq(t)=H(t)\left[q(t)-q(t-T)\right]$$

(26)

where

$$M=\left[\begin{array}{ccc}{m}_x& 0& 0\\ 0& m_y& 0\\ 0& 0& m_z\end{array}\right]$$

(27)

$$C=\left[\begin{array}{ccc}{c}_x& 0& 0\\ 0& c_y& 0\\ 0& 0& c_z\end{array}\right]$$

(28)

$$K=\left[\begin{array}{ccc}{k}_x& 0& 0\\ 0& k_y& 0\\ 0& 0& k_z\end{array}\right]$$

(29)

$$q(t)=\left[\begin{array}{c}{x}_c(t)\\ y_c(t)\\ z_w(t)\end{array}\right]$$

(30)

where M, C, and K denote the modal mass matrix, the modal damping matrix, and the modal stiffness matrix of the 3-DOF milling system, respectively. q(t) denotes the displacement vector.

At present, the SDM serves as an effective and widely applied method to calculate SLD. According to Insperger and Stépán’s theory [11], Equation (26) is rewritten in the following state-space form:

$$\dot{Q}(t)=R(t)Q(t)+S(t)Q(t-\tau )$$

(31)

with

$$R(t)=R(t+T)$$

(32)

$$S(t)=S(t+T)$$

(33)

where

$$Q(t)=\left[\begin{array}{l}q(t)\\ \dot{q}(t)\end{array}\right]$$

(34)

$$R(t)=\left[\begin{array}{cc}0& I\\ -M^{-1}K-M^{-1}H(t)& -M^{-1}C\end{array}\right]$$

(35)

$$S(t)=\left[\begin{array}{cc}0& 0\\ M^{-1}H(t)& 0\end{array}\right]$$

(36)

where Q(t) denotes the time-state vector. R(t) and S(t) denote the periodic coefficient matrixes describing the dynamics of the 3-DOF milling system. τ is the time delay.

The time delay τ is discretized into m numbers of equal time intervals Δt, namely, τ = mΔt, where m denotes an integer and is determined as follows:

$$m=\mathrm{int}\left(\frac{\tau +\Delta t/2}{\Delta t}\right)$$

(37)

where int is the function for rounding a number down to the closest integer.

When the initial condition satisfies Q(t_i) = Q_i, the solution of Equation (31) is stated as follows:

$$Q(t)=e^{R_i(t-t_i)}\left[Q_i+R_i^{-1}{S}_i{Q}_{\tau ,i}\right]-R_i^{-1}{S}_i{Q}_{\tau ,i}$$

(38)

Substituting Q_i+1 = Q(t_i+1) and Equation (38) into Equation (31) yields as follows:

$$Q_{i+1}=U_i{Q}_i+\alpha V_i{Q}_{i-m+1}+\beta V_i{Q}_{i-m}$$

(39)

where

$$U_i=e^{R_i\Delta t}$$

(40)

$$V_i=\left(e^{R_i\Delta t}-I\right)R_i^{-1}{S}_i$$

(41)

$$\alpha =\frac{m\Delta t+\frac{\Delta t}{2}-\tau }{\Delta t}$$

(42)

$$\beta =\frac{\tau +\frac{\Delta t}{2}-m\Delta t}{\Delta t}$$

(43)

According to Equation (39), a discrete map can be defined as follows:

$$w_{i+1}=D_i{w}_i$$

(44)

where

$$w_i=\left[\begin{array}{c}{Q}_i\\ Q_{i-1}\\ ⋮\\ Q_{i-m}\end{array}\right]$$

(45)

$$D_i=\left[\begin{array}{ccccc}{U}_i& 0& \cdots & \alpha V_i& \beta V_i\\ I& 0& \cdots & 0& 0\\ 0& I& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & I& 0\end{array}\right]$$

(46)

For evaluating the stability of the 3-DOF milling system, Equation (44) is written as follows:

$$w_k=\Phi w_0=D_{k-1}{D}_{k-2}\cdots D_1{D}_0{w}_0$$

(47)

According to Floquet’s theory, the eigenvalues of the matrix Φ are employed to evaluate the 3-DOF milling stability. When all modules of eigenvalue are smaller than 1, the 3-DOF milling system is in a stable state; when all modules of eigenvalue are equal to 1, the 3-DOF milling system is critically stable; when all modules of eigenvalue are greater than 1, the 3-DOF milling system is in an unstable state.

3. Predicted Results and Experimental Validation

In this section, the input parameters of the 3-DOF stability model, i.e., the CWE, the layered cutting force coefficients, and the modal parameters, are initially determined. Subsequently, these parameters are utilized to plot the SLD for the 3-DOF milling system using the SDM. In the end, a series of chatter verification tests are conducted to validate the predictions.

It should be pointed out that the tool wear effect is not considered in this article. In other words, the wear state of the cutter is assumed to be constant over a short machining time. This does not mean that tool wear can be ignored during machining. In the actual milling experiments, a new bull-nose end mill is used. Furthermore, the layered cutting force coefficients are also detected with this cutter.

The basic information about the cutter and the workpiece involved in the verification experiments is described here. A TiAlN-coated solid carbide bull-nose end mill is employed as the cutter, whose detailed specifications are presented in

Table 1. Cutter specifications.

Cutter Type	Cutter Material	Number of Flutes	Diameter (mm)	Fillet Radius (mm)	Helix Angle (°)	Overhang (mm)
Bull-nose end mill	3 μm TiAlN-coated carbide	4	12	2	38	75

. A thin-walled cylindrical part made of stainless steel is employed as the workpiece, whose geometric and physical parameters are presented in

Table 2. Workpiece parameters.

Workpiece Material	Geometries (mm)			Density (kg/m3)	Young’s Modulus (GPa)	Poisson’s Ratio
	Inner Diameter	Outer Diameter	Height
06Cr19Ni10	211	217	90	7930	194	0.3

.

3.1. CWE

In the bull-nose end milling of thin-walled cylindrical parts, the relative location relationship between the cutter and the workpiece is reflected by a lead angle [29]. The CWEs with six different lead angles are computed by applying the slice-intersection-based method in Section 2.3, as shown in Figure 6a. The results show that the engagement will shrink circumferentially and expand axially as the lead angle increases, which will cause the machined strip width to decrease. In other words, for the same scallop height, the smaller the lead angle, the greater the material removal rate. Therefore, a lead angle of 0 degrees is selected for practical machining in this paper. The corresponding calculated three-dimensional CWE in the CCS is shown in Figure 6b.

3.2. Layered Cutting Force Coefficients

For acquiring the layered cutting force coefficients, a group of slotting experiments is conducted on a VT850Z machining center, as demonstrated in Figure 7. The milling process parameters are listed in

Table 3. Machining parameters used in the slotting tests.

Spindle Speed (r/min)	Feed Rate per Tooth (mm/tooth)	Axial Depth of Cut (mm)
5000	0.008, 0.016, 0.024, 0.032	0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0

. According to the identification approach in Section 2.4, the layered cutting force coefficients for the selected combination of bull-nose end mill and workpiece material are calculated, as demonstrated in Figure 8.

The calculated specific coefficients F_t, F_r, and F_a are fitted to three cubic polynomials. The corresponding functions are expressed as follows:

$$\left\{\begin{array}{l}{K}_t(z)=-4.064z^3+14.767z^2-15.829z+8.955\hfill \\ K_r(z)=-6.887z^3+24.386z^2-25.779z+11.173\hfill \\ K_a(z)=1.279z^3-4.338z^2+4.347-0.364\hfill \end{array}\right.$$

(48)

Taking the height z as a variable, the layered cutting force coefficients associated with each cutting disk element can be obtained.

3.3. Dynamic Parameters

The modal parameters of the 3-DOF milling system along the x-axis, y-axis, and z-axis directions are important input parameters for the proposed model and need to be determined beforehand. The modal parameters of the bull-nose end mill respond to the dynamics of the cutter, the tool holder, and the machine spindle. Similarly, the modal parameters of the thin-walled cylinder respond to the dynamics of the workpiece, the fixture, and the machine table. In general, they are the integrated dynamic characteristics of the 3-DOF milling system. Therefore, the modal parameter results are dependent on the machine setup. A group of modal impact experiments is performed on the experimental setup, as demonstrated in Figure 9.

A Kistler 9722A impact hammer with a sensitivity of 10 mV/N is employed to produce the striking force. A Dytran 5225F1 accelerator sensor with a sensitivity of 10.04 mV/g is employed to collect the acceleration signals. The modal impact tests are performed 10 times along the x-axis, y-axis, and z-axis directions, respectively. The collected values are averaged to decrease random errors. The dynamic parameters of the 3-DOF milling system are determined by means of modal analysis.

Table 4. Modal parameters.

Cases	Frequency (Hz)	Damping (%)	Stiffness (N/m)
Cutter (x-axis)	1768	2.6	1.85 × 108
Cutter (y-axis)	1753	2.3	1.79 × 108
Workpiece (z-axis)	794	3.6	4.25 × 107

presents the identified results.

3.4. Chatter Stability Analysis and Experimental Validation

According to the SDM presented in Section 2.5, the SLD for bull-nose end milling of a thin-walled cylinder is calculated, as demonstrated in Figure 10.

The spindle speed varies within the interval from 2000 r/min to 8000 r/min. The axial depth of the cut varies within the interval from 0.0 mm to 2.0 mm. It can be found that there are three milling states: stable, critical, and unstable, which are distinguished by the stability boundary. The stable and unstable regions are below and above the stability boundary, respectively. The critical stable region is located near the stability boundary. For a depth of cut of 0.5 mm, the spindle speed within the intervals of 3204–4612 r/min and 5950–8000 r/min theoretically ensures stable milling. For a depth of cut of 1.0 mm, the spindle speed within the intervals of 3450–4372 r/min and 6369–8000 r/min theoretically ensures stable milling. However, the above ranges have to be contracted appropriately due to the critical stable region. In this paper, ±100 r/min is used to determine the critical stable region. Then, the spindle speed intervals for stable milling with a 0.5 mm axial depth of cut are 3304–4512 r/min and 6050–8000 r/min; the spindle speed intervals for stable milling with a 1.0 mm axial depth of cut are 3550–4272 r/min and 6469–8000 r/min.

For verifying the accuracy of the SLD, a group of bull-nose end milling tests are performed on a thin-walled cylinder, as demonstrated in Figure 9. The shared milling conditions are set as follows: the milling mode is down milling, the radial depth of cut is fixed at 6 mm, and the feed per tooth is fixed at 0.016 mm/tooth.

Table 5. Stability states with different milling conditions.

No.	Spindle Speed (r/min)	Feed Rate (mm/min)	Axial Depth of Cut (mm)	Stability State
1	3000	192	0.5	Critical
2	3000	192	1.0	Unstable
3	3500	224	0.5	Stable
4	3500	224	1.0	Critical
5	4000	256	0.5	Stable
6	4000	256	1.0	Stable
7	4500	288	0.5	Stable
8	4500	288	1.0	Unstable
9	5000	320	0.5	Unstable
10	5000	320	1.0	Unstable
11	5500	352	0.5	Unstable
12	5500	352	1.0	Unstable
13	6000	384	0.5	Critical
14	6000	384	1.0	Unstable
15	6500	416	0.5	Stable
16	6500	416	1.0	Stable
17	7000	448	0.5	Stable
18	7000	448	1.0	Stable

presents the remaining milling parameters.

The time domain and frequency domain acceleration signals corresponding to points A, B, and C are demonstrated in Figure 11. The spindle speeds at points A, B, and C are 4000 r/min, 5000 r/min, and 6000 r/min, respectively. The axial depths of the cut at points A, B, and C are all 0.5 mm. The bull-nose end milling state at point A is stable. The amplitude of the acceleration signal at point A is the smallest. And the spectral amplitude is primarily emphasized at the tooth passage frequency of 267 Hz and its multi-times-frequencies of 533 Hz and 1067 Hz. The bull-nose end milling state at point B is unstable. The amplitude of the acceleration signal at point B is the greatest. The corresponding spectrum has a chatter frequency of 1166 Hz, which is not an integral multiple of the tooth passage frequency of 333 Hz. The spectrum energy at the chatter frequency is greater than that at the tooth passage frequency. The bull-nose end milling state at point C is critical and stable. The amplitude of the acceleration signal is greater than that in the stable case but smaller than that in the unstable case. The chatter frequency of 1398 Hz occurs in the acceleration spectrum, whose energy is relatively low.

Figure 12 illustrates the finished surfaces of the thin-walled cylinder that correspond to the milling parameters at points A, B, and C in Figure 10. It can be clearly seen that significant chatter marks appear in the unstable case; slight chatter waves appear in the critical stable case; while in the stable case, the finished surface is sleek and its roughness (Ra) is very low.

4. Conclusions

Bull-nose end milling is an excellent finishing operation since it can provide greater material removal rates and smoother surface finishes. However, the regenerative chatter often occurs in the bull-nose end milling of thin-wall cylindrical parts, which have a ring-shaped and closed geometry structure. The dynamics of the cutter subsystem and the workpiece subsystem, the CWE, and the cutting force coefficients have a significant influence on the chatter stability. The objective of this article is to establish an accurate and applicable model to predict chatter stability during bull-nose end milling of thin-walled cylindrical parts. The main conclusions are summarized as follows:

(1)

A 3-DOF chatter stability model is developed, which is applicable to forecast the regenerative chatter in the bull-nose end finish milling of thin-walled cylindrical parts. It reflects the dynamics of the cutter subsystem in the x-axis and y-axis directions and the dynamic of the workpiece subsystem in the surface normal (namely, z-axis) direction. The SDM is employed to solve the dynamic equation describing the 3-DOF milling system;

(2)

A slice-intersection-based method for calculating the engagement boundary curve is established. The method is applicable to determining the CWE geometries with different milling conditions, such as slotting and following cuts. Furthermore, a lead angle of 0 degrees is recommended to maximize the machined strip width;

(3)

The approach to layered cutting force coefficient identification is presented considering the effect of varying cutter diameter. The specific coefficients for each cutting disk element can be determined by substituting the corresponding height into the cubic fitting polynomial;

(4)

The validation tests of the proposed model are performed on a four-axis CNC machine tool. The predicted SLD agrees well with the experimental data. The spindle speed and the axial depth of the cut can be optimally chosen using the SLD to effectively avoid regenerative chatter and achieve smooth surfaces.