1. Introduction
Micromilling can create complex miniature 3D features with sub-micron accuracy in various engineering materials, such as metals, ceramics, polymers, and composites. The micromilling process is extensively used in many high-technology sectors like biomedical, defense, and aerospace. However, in micromilling, self-excited vibration or chatter may cause catastrophic tool failure and damage the spindle and workpiece. Chatter vibration is detrimental to the dimensional accuracy and surface finish of the workpiece.
The dynamic interaction between the tool and cutting process causes chatter in machining. Regenerative chatter occurs during the cutting process when a flexible tool leaves behind a wavy surface on the workpiece in the preceding cut. The phase difference between the waviness of the current and the preceding passes can result in a condition where the amplitude of the vibration grows rapidly. In the micromilling process, the key reason behind the occurrence of the regenerative chatter is the limited flexural stiffness of the small end mills. Very high spindle speeds can counter the limited flexural stiffness resulting in lower chip loads and cutting forces. However, the spindle dynamics can influence the chatter at high rotational speeds due to gyroscopic effects, especially for low stiffness tools (≤100 µm in diameter and high length to diameter ratio). Therefore, it is essential to investigate the impact of the gyroscopic couple on stability in high-speed machining.
The adverse effects caused by chatter vibration have inspired researchers to develop different models for predicting and avoiding the chatter in conventional milling and micromilling. Tobias [
1] and Tlusty [
2] presented the mechanism of chip regeneration and stability charts for single-point metal cutting that provides stable depths of cut as a function of spindle speed. Minis and Yanushevsky [
3] proposed a theoretical approach based on Fourier analysis and the transfer function. They have developed a comprehensive mathematical model for a milling process with two orthogonal degrees-of-freedom (DOF). Budak and Altintas [
4,
5] analytically predicted the stability limits in the frequency domain. They also proposed a multi-frequency solution to improve the prediction at small radial immersions where the process is highly intermittent. Insperger et al. [
6,
7] performed stability analysis in the discrete time domain, including the periodically varying parameters. It is important to note that these studies have considered single or two DOF systems to predict the stability limits.
Several researchers have studied the occurrence of bifurcation during the machining process. Bifurcation is the qualitative behavior of the solutions as the control parameters in a system are varied [
8,
9]. Davies et al. [
10] used the Poincare sectioning technique to characterize and identify chatter in high-speed milling. Zhao and Balachandran [
11] developed a time-domain model for both regenerative and loss of contact effects in a milling system. They studied the bifurcations diagrams for parameters, such as axial depth of cut. Schmitz and Honeycutt [
12] presented the extended milling bifurcation diagram using time-domain simulation and once-per-revolution sampling. They found both stability behavior such as Hopf and period-n bifurcations in bifurcation diagrams. Consequently, it is crucial to understand the system behavior with a control parameter of interest using bifurcation analysis.
It may be noted that the gyroscopic effect on dynamic stability becomes more pronounced in high-speed micromachining. Movahhedy and Mosaddegh [
13] used finite element analysis to determine the modal properties of the spindle system and its frequency response function (FRF). They assumed rigid connections between the tool, tool holder, and spindle. They predicted speed-dependent stability boundaries and observed that the gyroscopic effects become significant only at very high speeds. Bediz et al. [
14] used a custom-made impact excitation system to obtain speed-dependent dynamics of high-speed spindles through experimental modal analysis. They noticed that the spindle speed affected the dynamic response significantly, especially at higher spindle speeds. Cao et al. [
15] investigated the effect of the gyroscopic moment and centrifugal forces on spindle dynamics. They used the Nyquist stability criterion to generate the stability lobe diagram of high-speed milling. They observed a speed-dependent shift in stability lobes. Xiong et al. [
16] investigated the gyroscopic effect of the rotating spindle on the stability characteristics of the conventional milling system and found that it makes the stability prediction less conservative. Wan et al. [
17] investigated the effect of gyroscopic couple on the stability. They used speed-dependent FRFs of the tool. They considered the effect of runout on the stability. Mokhtari et al. [
18] investigated the effects of rotary inertial dynamics on prediction of chatter in the milling operation. They considered a 3D rotating cantilever Timoshenko beam for modeling the cutting tool. They found that without the rotary inertial dynamics, significant errors are in the prediction of chatter limits at high rotational speeds. Lu et al. [
19] obtained the FRF of the spindle system based on the rotating Timoshenko beam theory and the receptance coupling substructure analysis (RCSA). They investigated the centrifugal force and gyroscopic effect caused by the high-speed rotation of the micro-milling spindle. Jalili et al. [
20] modeled a cutting tool as a rotating clamped-free beam which is excited by cutting forces. They investigated various kinds of resonances in milling operation using a 3-D nonlinear dynamic model of the milling process including both structural and cutting force nonlinearities, gyroscopic moment, and rotary inertia. Totis [
21] introduced additional rotational degrees of freedom in the stability analysis of the conventional milling process. However, the effect of the gyroscopic couple was not investigated in this study. Very few studies have investigated the gyroscopic effect on stability limits in high-speed machining. These models have not explicitly considered the rotational degrees of freedom along with the translational degrees of freedom.
Hence, this study is focused on understanding and modeling the gyroscopic effect in a compliant micro-tool (at high rotational speeds) on the stability boundaries via a rotor dynamics-based analysis with the cutting forces at one end. A MDOF model includes rotational degrees of freedom with the translational degrees of freedom to investigate the effect of gyroscopic moment. Bifurcation analysis was performed to understand the system behavior for a control parameter in the MDOF system. The time-periodic delay differential equations governing the dynamics of the MDOF milling system were numerically solved using a linear multistep (LMS) method. The stability charts and bifurcation diagrams were constructed along with a comparison of displacement profiles at different depths of cut and spindle speeds. Poincare maps were also developed to observe phase space trajectories. The predicted stability limits were validated at different rotational speeds and depths of cut. This study will be beneficial to understand the long-term dynamic behavior of the micromilling system, which is affected by gyroscopic moment especially, at high rotational speeds.
2. Equations of Motion for MDOF Stability Model
For a low flexural stiffness end mill, the effect of the gyroscopic couple will be prominent at high rotational speeds. Both rotational and translational degrees of freedom have been considered in the stability model in the authors’ previous work [
22] to capture the effect of gyroscopic couple on stability. This model considers a micro end mill with mass ‘
m’ a rotor, as shown in
Figure 1. Flexible supports at the ends (clamping end and tooltip end) of the micro end mill are assumed. The center of mass is located at distance ‘
a’ and ‘
b’ from the clamping end and the tool tip-workpiece interface, respectively. The micro end mill has a fixed frame of reference ‘XYZ’ and a rotating frame of reference ‘xyz’. The micro end mill is rotating around the
z-axis with spindle speed ‘
ω’. The center of mass has orthogonal translational degrees of freedom in
and
, and rotational degrees of freedom in
and
rotations about the x and y axes, respectively.
,
and
,
are the displacements at the clamping end and tooltip end, respectively, which can be defined in terms of
,
and
by,
The dynamic equations of motion of the center of mass in the translational and rotational directions can be determined using a Lagrangian formulation described in Mittal et al. [
22] as,
where
where
and
are the transverse and the axial moments of inertia, respectively,
are the stiffness values of the flexible bearing supports in x and y directions,
are the damping coefficients of the flexible bearing supports in x and y directions, and
Fx and
Fy are the cutting forces in the x and y directions. It may be noted that in Equations (4) and (5),
and
are the gyroscopic moment acting about x and y axes, respectively.
The derived set of equations of motion for a higher-order multiple degrees of freedom stability model consists of translational and rotational degrees of freedom can be expressed in the following matrix form,
where
M, C, and
K are the mass matrix, damping matrix with the gyroscopic component, and stiffness matrix, respectively.
is the displacement vector and
is the force matrix, determined using dynamic force modeling.
M, C, K,
and
are defined as,
A mechanistic force was used to determine the dynamic cutting forces
Fx and
Fy in terms of dynamic chip thickness (as shown in
Figure 2) as described by Budak and Altintas [
4].
The dynamic chip thickness is given by,
where
and
are the difference between the dynamic displacement during the present (
j) and previous (
j − 1) cut.
The cumulative cutting forces in matrix form are expressed as,
where
,
,
and
are the directional factors.
By substituting
Fx and
Fy and
and
in terms of
,
, and
, the force matrix
F can be expressed as,
Equation (10) can be rewritten as
where
is given as
,
is cutting force coefficient in tangential direction and
A(
t) is the directional factor matrix.
By substituting force matrix in Equation (6), the dynamic equations of motion for a micro end mill system can be expressed in matrix form as,
3. Stability Analysis Using Linear Multistep Method (LMM)-Based Numerical Scheme
The multiple degrees of freedom model has been expressed in terms of the periodic delay-differential equation where
A(
t) is a time-periodic term (
. The frequency-domain approach will not predict the stability limits as accurately as the time-domain approach. Therefore, in this study, the periodic delay differential equations have been solved numerically using the fourth-order Explicit Adams–Bashforth method, a linear multistep method (LMM) [
23,
24]. Further, the stability limits have been identified using numerical solutions at different rotational speeds. The dynamic equations of motion for a micro end mill system can be expressed in matrix form as second-order delay differential equation (Equation (12)),
where
,
,
,
, and
are the matrices of mass, damping coefficients, stiffness, directional factors, and displacement variables, respectively.
By defining the state variables as,
The second-order delay differential equation for the periodic milling dynamics shown in Equation (12) is organized in the first-order delay differential equation using state variables as follows,
where
By multiplying the Equation (14) on both sides with the inverse of matrix
, Equation (14) becomes
where
The values of state variables in the previous cut are required to solve Equation (15). The first-order delay differential equations as shown in Equation (15) can be written as
where
is the initial function, T is the time period,
is called the delay argument, and
is the solution of the delay term.
The linear multistep methods are extensively used for the numerical solution of initial value problems without the delay term. The delay period T is divided into
r discrete time intervals of
. Assume the value to
and
at the current time
is
and
, respectively. According to the fourth-order explicit Adams–Bashforth method, the next value
from previous values of
can be computed as
While solving the delay equation (Equation (15)) using the Adams–Bashforth method, the solution for first interval will depend on the value of function . The values of at previous steps will be determined using . For subsequent time intervals, the values of in previous steps will be used from determined values of in previous intervals. Using Equation (17), the motion of the micro end mill system (see Equation (12)) can be numerically simulated for different spindle speeds and axial depths of cut.
A schematic to generate the stability plots is shown in
Figure 3. The stability limits are predicted by progressively increasing the axial depth of cut for a spindle speed given by ω = 2π/T until the numerical simulations for motion show instability. Root mean square (RMS) values of first half and second half of displacements data were calculated as RMS
1 and RMS
2, respectively. If the RMS
1 is greater than the RMS
2, then the system is assumed as stable, otherwise it will be unstable system. This procedure is repeated for the next spindle speeds to construct the stability lobe diagrams.
4. Results and Discussion
The numerical modeling is carried out for slot milling operation. The simulations were performed for different spindle speeds and depths of cut to construct the stability lobe diagrams. In this study, cutting force coefficients were assumed to be constant. To determine the tangential and radial cutting force coefficients, slot milling experiments were conducted on the high-speed micromachining center developed in the Machine Tools Lab, IIT Bombay, as shown in
Figure 4. The two-fluted micro end has a diameter of 100 µm. The length of flute is 200 µm. The overhang length of tool (length of tool outside the clamping) is 20 mm. The workpiece material is Ti
6Al
4V which finds application in the aerospace and biomedical industries. The cutting forces in x and y directions were measured using a 3-axis tabletop dynamometer (Minidyn 9256C1, Kistler, Pune, India). The experimental conditions as shown in
Table 1 were used to carry out experiments to determine the cutting force coefficients. Experiments were carried out at feed per rev of 2 µm/flute and depth of cut of 10 µm. The mechanistic force modeling was used to determine the radial and tangential cutting force coefficients from experimentally measured cutting forces. Average value of measured cutting force coefficients for different speeds was used. The value of radial and tangential cutting force coefficients were 9182 MPa and 7205 MPa
, respectively.
The modal parameters (natural frequencies, stiffnesses, and damping ratios) at the clamping end and tooltip end are required to predict the stability limits. Impact testing is performed in which excitation is done using a miniature impact hammer (Dytran Chatsworth, CA, USA, Model 5800SL) and displacement measured using a laser displacement sensor (Keyence, Pune, India, Model LK-H020). The stiffness and damping ratios of the machine tool structure at the ends of the rotor, which are assumed as flexible support bearings, are shown in
Table 2. The dynamic parameters in x and y directions are assumed to be isotropic because of symmetrical collet and tool.
4.1. Bifurcation Analysis
A bifurcation diagram is a tool to study the long-term qualitative behavior of the dynamic system with the variation of the control parameter. In the stability analysis of milling, the control parameters are axial depth of cut and spindle speed during the machining operation. To construct the bifurcation diagram, simulations were performed for the tooltip displacement at a given depth of cut. Initial values for state variables
are 10
−6, 10
−6, 10
−4, 10
−4, 0, 0, 0, and 0, respectively. The displacement values are collected for a particular time period and the process is repeated for the next depths of cut. The bifurcation diagrams are numerically constructed for a slot milling operation, as shown in
Figure 5, for a higher spindle speed of 120,000 rpm. The control parameter for these bifurcation diagrams is the axial depth of cut. It can be seen that bifurcation occurs at depth of cut
ap = 4.1 µm. The system is stable at depths of cut below 4.1 µm and the system shows unstable behavior at depths of cut greater than 4.1 µm.
To find the dynamic behavior of Equation (12), the numerical simulations are performed for different values of bifurcation parameter depth of cut (
ap).
Figure 6 and
Figure 7 show the waveform and Poincare maps for pre-bifurcation (depth of cut
ap = 4.1 µm) and post-bifurcation (depth of cut
ap = 4.6 µm) motions at a spindle speed of 120,000 rpm.
Figure 6 shows that if the axial depth of cut is less than the critical value, then the zero-equilibrium point is asymptotically stable. The amplitude of displacements (see
Figure 6a,b) and tool deflections (see
Figure 6c,d) approaches zero with time. It can also be seen from the orbital map (see
Figure 6e) and the Poincare maps (see
Figure 6f) that trajectory moves towards the equilibrium point, and the system is stable. When the axial depth of cut increases beyond 4.1 µm, a Hopf bifurcation occurs at
ap = 4.4 µm, as shown in
Figure 5. For post-bifurcation axial depth of cut
ap = 4.6 µm, the displacement profiles (see
Figure 7a,b) and tool deflection profiles (see
Figure 7c,d) diverge with time, indicating unstable behavior. The orbital map and Poincare maps as shown in
Figure 7e,f respectively, show that trajectory moves away from the equilibrium point, and the system is unstable.
Phase portrait for pre-bifurcation, bifurcation point, and post-bifurcation motions are shown in
Figure 8. It can be seen that at the bifurcation point (
ap = 4.4 µm), the vibration is not growing or decaying. A depth of cut below 4.4 µm, trajectory moves towards zero (shown in the center of the plot), and the system is stable, but depth of cut greater than 4.4 µm indicates unstable behavior.
Figure 9 shows a comparison of the orbital maps (
vs.
) for spindle speeds of 60,000 and 120,000 rpm at one fixed depth of cut
ap = 4.1 µm. Variation in
and
are significantly higher for spindle speed of 120,000 rpm as compared to spindle speed of 60,000 rpm. This can also be seen from the waveform plots as shown in
Figure 10a,b. The amplitude of deflection in x and y direction is significantly higher for a higher spindle speed of 120,000 rpm. Root mean square (RMS) of
and
increases by 37% and 45%, respectively, for spindle speed of 120,000 rpm as compared to 60,000 rpm. The higher deflection at higher speeds shows the presence of the gyroscopic effect in the machining process.
4.2. Stability Lobe Diagram and Validation
The stability lobe diagrams were constructed for slot milling operation (100% immersion) on Ti
6Al
4V using LMS based numerical scheme.
Table 2 shows the modal parameters at tool tip and clamping ends used in the numerical simulations. The feed per flute is fixed at 2 µm for the prediction of stability limits. The delay period T is divided in
r = 100, which is a parameter for the discretization of the tooth period.
Figure 11 shows the stability limits for the 2DOF and MDOF stability models. In the stability chart, the loss of stability is associated with the Neimark bifurcation of periodic motions, also known as Hopf bifurcation. The MDOF chatter model shows conservative stability limits, especially at high rotational speeds (>100,000 rpm). This conservative behavior of stability limits at high rotational speeds can be attributed to the inclusion of additional rotational degrees of freedom with gyroscopic components. The higher deflection at higher spindles (as shown in
Figure 10) shows the effect of gyroscopic couple. Due to effect of gyroscopic couple at higher rotational speeds, there is a reduction in stability limits at higher spindle speeds compared to the lower spindle speeds.
Experiments were carried out to verify the stability predictions at different depths of cut and spindle speeds selected from the stability chart shown in
Figure 11. The FFT spectrum of velocity and machined surface data were used to identify the stable and unstable depths of cut experimentally.
Figure 12 shows the FFT spectrum of the velocity data at spindle speeds of 50,000 and 120,000 and corresponding depths of cut of 20 µm and 5 µm. It can be seen that there is no high amplitude peak at a rotational speed of 50,000 rpm and a depth of cut of 20 µm which shows stable machining. At a spindle speed of 120,000 rpm and a depth of cut of 5 µm, the amplitude peaks of the velocity spectrum are higher at around 6000 Hz which is in the vicinity of one of the natural frequencies of the tool tip. The critical limits for stability obtained from numerical simulations and experiments are compared in
Figure 13. It can be seen that the predicted stability limits from the MDOF model are in better agreement with experimental limits at high rotational speeds (>100,000 rpm). Difference between predicted and experimental limits at high rotational speeds is less than 16%. Hence, using an MDOF model is imperative to predict stability limits at very rotational speeds as the effect of the gyroscopic couple becomes pronounced, especially for miniature tools (~100 µm in diameter).