Milling Mechanism and Chattering Stability of Nickel-Based Superalloy Inconel718

1. Introduction

Nickel-based Superalloy Inconel718 is widely used in the aerospace industry due to its extraordinary high-temperature strength, thermal stability, heat fatigue resistance, corrosion resistance, creep resistance, and oxidation resistance, especially in the field of manufacture of key components in the aerospace industry, such as turbine blades, combustion chamber components, and gas turbines[1,2,3,4]. Common processing methods for Inconel718 are mainly casting[5], forging[6], laser processing[7] and electrical discharge machining[8]. However, each of these machining methods suffers from flaws like low dimensional accuracy, inability to process complex geometry, high machining costs and containing heat-affected areas, slow machining speed, and poor surface quality respectively. Although utilizing the milling method can effectively avoid the above-mentioned problems, Inconel718 in the cutting process is prone to serious tool wear, high thermal stress, high cutting forces, vibration, and other phenomena[9]. In order to address these challenges and improve the machinability of Inconel718, the milling mechanisms and chattering stability of Inconel718 need to be investigated. Based on this, the optimized processing parameters could be determined by analyzing the influence of different processing parameters, such as spindle speed, feed per tooth, radial depth of cut, and axial depth of cut on cutting force, temperature, stress, and deformation, including machining quality during manufacturing.

Accurate modeling of milling force can help adjust cutting parameters to ensure high-quality machining results, therefore, it is especially necessary to establish an accurate milling force model. At present, the milling force models proposed by scholars mainly include empirical models, physical analytical models, and finite element models.

The empirical model of milling force is a model derived by fitting milling test data and is mainly used to predict the magnitude of milling force. Ding[10], Zhao[11], and Lei[12] used orthogonal tests to obtain the milling force coefficients and established an empirical model of milling force. Bergmann[13] developed a nonlinear empirical model of milling force and improved the prediction of stability limits achieved by parameterizing a milling force model for stability analysis. Zhao[14] established an empirical model of milling force for ultrasonic vibration milling of titanium alloy and explored the effects of different machining parameters on milling force. Liang[15] proposed an empirical model of cutting forces in micro-end-milling operations that takes into account the parameters of feed per tooth per tool radius, the true trochoidal nature of the tool edge path, and the chip thickness. Li[16] established an empirical model for the prediction of milling force and cutting temperature based on the multiple linear regression method. In fact, the empirical model of milling force relies on a large amount of experimental data and has obvious limitations, which cannot reveal the dynamic characteristics and mechanism of cutting, and can only be used to calculate the average milling force, rather than accurately predict the instantaneous cutting force during the milling process.

A physical analytical model of milling force is a model that describes the cutting forces generated during the milling process with the help of mathematical analytical methods. The physical analytical model of milling force is usually based on the basic theories of machining principles and material mechanics. Jiang[17] proposed that the projected area of the milling shear surface in three directions determines the milling force in three directions, and established an analytical model of milling force based on the shear surface. Madajewski[18] presented a computational method based on a combination of finite element analysis and classical analytical methods for predicting the cutting force components in face milling processes. Wang[19] discretized the end milling cutter along the circumferential direction and established an analytical cutting force model for the end milling cutter and simplified analytical indexes for predicting and estimating cutting force fluctuations. Kao[20] developed a mathematical model of the cutting forces of a ball-end milling cutter that includes tangential, radial, and axial cutting force components, and characterized the relationship between average cutting force and feed per groove as a linear function. Zhou[21], Kazuki[22], and Li[23] modeled the milling force of a ball-end milling cutter considering changes of microelement friction angle, shear angle, and shear stress constraints using microelement diagonal shear theory. Falta[24] and Dong[25] used an analytical modeling method for predicting the static milling force of a ball-end milling cutter by establishing a tool-workpiece engagement region. The removal mechanism of materials can be deeply analyzed through the physical analytical model; however, the physical analytical model is closely related to material constitutive parameters, which are difficult to obtain. Therefore, the physical analytical model is not convenient for practical application, and its modeling process is complex and cumbersome.

Finite element model of milling force is a mathematical model that predicts milling forces and contact forces during the cutting process by discretizing the material regions being milled into small finite element regions. Li[26] used ABAQUS software to analyze the effects produced by polycrystalline diamond tools for machining thin walls of SiCp/Al composite material. Jin[27] established an optimized finite element model based on pendulum motion, considered the performance parameters of the workpiece, and used ABAQUS software to study the rule of change of milling force during the milling process of aerospace thin-wall component. Hu[28] modeled the thermal-mechanical coupling property of ultrasonic torsional vibration-assisted micro-milling (UTVAM) using ABAQUS software. Pratap[29] used the Johnson-Cook constitutive equation to develop a finite element model for micro-channel micro-end milling of titanium alloys using ABAQUS that takes into account the strain, strain rate, and temperature effects on material properties. Li[30] imported the blade model created using UG NX into ANSYS software for simulation of milling force and calculated the elastic deformation of the blade using both MATLAB and ANSYS. The finite element model of milling force can simulate the real milling process, consider the contact and interaction between the tool and the workpiece, and is able to predict the mechanical properties and deformation behavior. However, the use of different finite element models will lead to different accuracy results. Therefore, it is recommended that a higher precision, more targeted finite element model should be chosen to accurately and intuitively simulate the milling process.

Chattering is an unstable state of motion in metal cutting processing, mainly caused by cutting force, cutting heat, and material deformation. Inconel718 is difficult-to-machine material, and its high hardness and heat resistance make it more susceptible to chattering during the milling process. At present, the theoretical analysis methods for milling stability mainly include frequency domain analysis and time domain analysis.

The frequency domain analysis is a method used to evaluate and solve the chattering problem in the milling process by converting the dynamic response of the milling system to the frequency domain. Frequency domain analysis method can be further divided into zero-order frequency domain method and multi-frequency method. Altintaş and Budak[31] proposed a zero-order frequency domain method of milling stability, in which the oriented dynamic milling force coefficient matrix is transformed into the frequency domain by Fourier transform to solve for the stability of the system without taking into account the cross-transfer function of the system. Based on that, the critical axial depth of cut without causing chattering is determined. Merdol[32] considered the higher-order expansion terms of the Fourier series of the oriented dynamic milling force coefficient matrix and gave a multi-frequency method that can be used in different radial depths of cut cases.

The time domain analysis is another method used to evaluate and solve the chattering problem in the milling process by observing and analyzing the time response. The time domain analysis method can be further divided into fully-discrete and semi-discrete methods. Ding[33,34] verified the accuracy of the method by single-degree-of-freedom and two-degree-of-freedom milling models and then proposed the second-order full-discrete method. Insperger[35,36,37] proposed a semi-discrete method to solve a system of time-delay differential equations containing a matrix of periodic coefficients. Subsequently, the method is improved so that it can not only solve the stability limit of a single-degree-of-freedom vibration system, but also a two-degree-of-freedom vibration system, and the first-order semi-discrete method with better convergence is proposed. Based on the semi-discrete method, Sun[38] established the stability lobe diagram of rotary ultrasonic milling of titanium alloy by defining the ultrasonic function angle, and the results showed that the rotary ultrasonic milling can suppress the machining vibration and improve the stability of milling of titanium alloy.

In this paper, the milling force model of a ball-end milling cutter is established by analytical method, additionally, both the experimental and finite element analysis is carried out. Since the material studied in this paper is Inconel718, which is prone to chattering and can easily lead to machining instability. Besides, the zero-order frequency domain method is computationally efficient and suitable for general processing situations, so this paper adopts the zero-order frequency domain method.

2. Modeling of Milling Forces

2.1 Modeling of Milling Forces of Ball-End Milling Cutter

The milling state of the helix cutting edge of the ball-end milling cutter varies with different cutting positions. Thus, it is essential to develop a geometric model of the ball-end milling cutter to investigate the milling force changes when the cutting edge interacts with the workpiece. The schematic of the proposed model is illustrated in Figure 1.

Figure 1 Milling force model of ball-end milling cutter.

According to the literature[20], neglecting the tip radius of the ball-end milling cutter results in the cutting force on any microelement of the cutting edge being described by equation (1).

$$\{\begin{array}{l}d{F}_{t,j}\left({\varphi }_j\left(z\right)\right)=K_{te}dS\left({\varphi }_j\left(z\right)\right)+K_{tc}{h}_j\left({\varphi }_j\left(z\right)\right)db\\ d{F}_{r,j}\left({\varphi }_j\left(z\right)\right)=K_{re}dS\left({\varphi }_j\left(z\right)\right)+K_{rc}{h}_j\left({\varphi }_j\left(z\right)\right)db\\ d{F}_{a,j}\left({\varphi }_j\left(z\right)\right)=K_{ae}dS\left({\varphi }_j\left(z\right)\right)+K_{ac}{h}_j\left({\varphi }_j\left(z\right)\right)db\end{array}$$

(1)

The instantaneous chip thickness at the immersion angle $h_j\left({\varphi }_j\left(z\right)\right)$ can be expressed as:

$$h_j\left({\varphi }_j\left(z\right)\right)=f_t\sin \left({\varphi }_j\left(z\right)\right)\sin \left(\kappa \left(z\right)\right)$$

(2)

The instantaneous edge length of the cutting segment $dS\left({\varphi }_j\left(z\right)\right)$ can be expressed as:

$$dS=dz·\sqrt{\frac{{\left(R_0-z\right)}^2}{R_0{}^2-{\left(R_0-z\right)}^2}+\frac{{\tan }^2\beta }{R_0{}^2}\left[R_0{}^2-{\left(R_0-z\right)}^2\right]+1}=M·dz$$

(3)

The instantaneous chip width $db$ can be expressed as:

$$db=\frac{dz}{\sin \left(\kappa \left(z\right)\right)}=\frac{R_0}{\sqrt{R_0{}^2-{\left(R_0-z\right)}^2}}dz$$

(4)

By transforming the cutting force microelement from a coordinate system $O-tra$ to $O-fnz$ and integrating it, the cutting force of the milling cutter along feed, normal and axial directions can be obtained as:

$$\left[\begin{array}{c}{F}_f\left(\varphi \right)\\ F_n\left(\varphi \right)\\ F_z\left(\varphi \right)\end{array}\right]=\left[\begin{array}{c}{\displaystyle \sum _{j=1}^N{F}_{f,j}\left({\varphi }_j\right)}\\ {\displaystyle \sum _{j=1}^N{F}_{n,j}\left({\varphi }_j\right)}\\ {\displaystyle \sum _{j=1}^N{F}_{z,j}\left({\varphi }_j\right)}\end{array}\right]$$

(5)

Subsequently, the average milling force per rotation of the ball-end milling cutter can be calculated using the equation (6):

$$\{\begin{array}{l}\overline{F_f}=\frac{N_f}{2\pi }{\displaystyle {\int }_0^{2\pi }\left({\displaystyle {\int }_{z_1\left({\varphi }_j\right)}^{z_2\left({\varphi }_j\right)}d{F}_f\left({\varphi }_j\left(z\right)\right)}\right)d\varphi }\\ \overline{F_n}=\frac{N_f}{2\pi }{\displaystyle {\int }_0^{2\pi }\left({\displaystyle {\int }_{z_1\left({\varphi }_j\right)}^{z_2\left({\varphi }_j\right)}d{F}_n\left({\varphi }_j\left(z\right)\right)}\right)d\varphi }\\ \overline{F_z}=\frac{N_f}{2\pi }{\displaystyle {\int }_0^{2\pi }\left({\displaystyle {\int }_{z_1\left({\varphi }_j\right)}^{z_2\left({\varphi }_j\right)}d{F}_z\left({\varphi }_j\left(z\right)\right)}\right)d\varphi }\end{array}$$

(6)

Equation (7) is finally obtained by simplifying the equation of the microelement of milling force and substituting it into (6):

$$\{\begin{array}{l}{\overline{F}}_f=\left(C_1{K}_{tc}+C_3{K}_{rc}+C_5{K}_{ac}\right)f_t+C_2{K}_{te}+C_4{K}_{re}+C_6{K}_{ae}\\ {\overline{F}}_n=\left(C_7{K}_{tc}+C_9{K}_{rc}+C_{11}{K}_{ac}\right)f_t+C_8{K}_{te}+C_{10}{K}_{re}+C_{12}{K}_{ae}\\ {\overline{F}}_z=\left(C_{13}{K}_{rc}+C_{15}{K}_{ac}\right)f_t+C_{14}{K}_{re}+C_{16}{K}_{ae}\end{array}$$

(7)

The parameters $C_1-C_{16}$ in equation (7) have been given by reference[20].

2.2 Recognition Experiment of Milling Force Coefficients

In order to be able to accurately find out the milling force generated during milling processing, it is necessary to determine the milling force coefficients $K_{te}$, $K_{re}$, $K_{ae}$, $K_{tc}$, $K_{rc}$, $K_{ac}$ in equation (7), therefore, recognition experiments of milling force coefficient are required.

For this experiment, the Sandvik 1B240-0800-XA 1630 carbide ball-end milling cutter is chosen, the diameter of the cutter is 8mm, the overall length of the cutter is 100mm, the length of the cutting section is 25mm, the depth of cut is set to 0.5mm and the rotating speed of spindle is set to 1000rpm. In order to facilitate the calculation of milling force coefficients, slot milling is used as the milling method, and the selected machining parameters and the measured average milling force data are recorded in Table 1.

Based on data in Table 1, the scatter plot is plotted and the milling force coefficient is obtained by calculating the slope and intercept using the coefficient recognition equation of milling force (8).

$$\{\begin{array}{l}{K}_{ac}=\frac{{\overline{F}}_{zc}\left(C_3{C}_7{C}_{15}-C_1{C}_9{C}_{15}\right)-C_7{C}_{13}{C}_{15}{\overline{F}}_{fc}+C_1{C}_{13}{C}_{15}{\overline{F}}_{nc}}{C_3{C}_7{C}^2{}_{15}-C_5{C}_7{C}_{13}{C}_{15}-C_1{C}_9{C}^2{}_{15}+C_1{C}_{11}{C}_{13}{C}_{15}}\\ K_{ae}=\frac{{\overline{F}}_{ze}\left(C_4{C}_8{C}_{16}-C_2{C}_{10}{C}_{16}\right)-C_{14}{C}_8{C}_{16}{\overline{F}}_{fe}+C_2{C}_{14}{C}_{16}{\overline{F}}_{ne}}{C_4{C}_8{C}^2{}_{16}-C_6{C}_8{C}_{14}{C}_{16}-C_2{C}_{10}{C}^2{}_{16}+C_2{C}_{12}{C}_{14}{C}_{16}}\\ K_{rc}=\frac{C_7{C}_{15}{\overline{F}}_{fc}-C_1{C}_{15}{\overline{F}}_{nc}+\left(C_1{C}_{11}-C_5{C}_7\right){\overline{F}}_{zc}}{C_3{C}_7{C}_{15}-C_5{C}_7{C}_{13}-C_1{C}_9{C}_{15}+C_1{C}_{11}{C}_{13}}\\ K_{re}=\frac{C_8{C}_{16}{\overline{F}}_{fe}-C_2{C}_{16}{\overline{F}}_{ne}+\left(C_2{C}_{12}-C_8{C}_6\right){\overline{F}}_{ze}}{C_4{C}_8{C}_{16}-C_6{C}_8{C}_{14}-C_2{C}_{10}{C}_{16}+C_2{C}_{12}{C}_{14}}\\ K_{tc}=\frac{{\overline{F}}_{fc}-C_3{K}_{rc}-C_5{K}_{zc}}{C_1}\\ K_{te}=\frac{{\overline{F}}_{fe}-C_4{K}_{re}-C_6{K}_{ze}}{C_2}\end{array}$$

(8)

Afterwards, a linear regression equation is formulated, taking the feed per tooth as the independent variable and the measured average milling force as the dependent variable. Figure 2 illustrates the milling force variation with changes in the feed per tooth along each axis, and Table 2 presents the coefficients of the milling force.

Figure 2 The relationship between average milling force and feed per tooth.

2.3 Validation of Milling Force Model for Ball-End Milling Cutter and Results Analysis

Based on the above-proposed milling force model of the ball-end milling cutter, the Matlab program, which is illustrated in the Appendix, is written to calculate the theoretical milling force. The results of the theoretical analysis are compared with the experimentally measured results to verify the correctness of the proposed milling force model. The same manufacturing parameters are set for both theoretical and experimental analysis, namely, the spindle speed is 1000rpm, the axial depth of cut is 0.4mm and the feed per tooth is 0.04mm/z. The results of the comparative analysis of the theoretical and experimental values of milling forces of the X, Y, and Z-axis are shown in Figure 3, observing the milling force values in each axis, it is evident that both the theoretical and experimental results exhibit a similar change trend, validating the accuracy of the proposed theoretical model for calculating the ball-end milling cutter's milling force. Although a slight degree of error is present in the figures, they remain within an acceptable range.

Figure 3 Experimental and theoretical comparison of milling forces.

The feed per tooth plays a pivotal role in determining both the magnitude of the milling force and the quality of the machined surface. Consequently, to investigate the influence of the feed per tooth on the milling forces, the different feed rate per tooth is set separately as $f_z$=0.1mm/z, 0.2mm/z, 0.3mm/z, 0.4mm/z, while keeping the radial depth of cut and axial depth of cut unchanged, namely $a_e$=3mm, $a_p$=0.5mm. In addition to feed per tooth, the tool helix angle is another important parameter that affects the milling process. Change in helix angle in ball-end milling cutter affects the cutting angle and feed angle, which in turn affects the contact between tool and workpiece. The tool helix angle is set separately as $\beta =15°$, $\beta =25°$, $\beta =35°$, $\beta =45°$ in the case of ensuring other parameters remain unchanged. The results of above studies are shown in Figure 4 - Figure 6. When the feed per tooth of the ball-end milling cutter increases, there is a corresponding upward trend observed in the milling force for each axis. Similarly, an increase in the helix angle of the ball-end milling cutter results in a proportional rise in the cutting force. This can be attributed to the expansion of the contact area between the tool and the workpiece, leading to an overall mounting trend in cutting force along all axes. Due to the significant impact of the ball-end milling cutter's helix angle on the generated milling force in each axis, careful consideration and selection of the appropriate tool helix angle are imperative to meet specific machining requirements

3. Finite Element Analysis of Milling of Inconel718

In this study, firstly, a finite element model is established, as shown in Figure 7, including the workpiece and the tool. In the following, the temperature, stress and milling force in the milling process are calculated. Finally, the influence of different machining parameters on the machining condition of Inconel718 is evaluated through simulation results.

3.1. The Constitutive Model of Inconel718

The constitutive model is a critical prerequisite of accuracy and reliability of finite element analysis. The behavior of materials and results of engineering problem can be better provided by a right chosen constitutive model[39]. The constitutive models applicable to metal cutting are mainly Power Law and Johnson-Cook model[40]. The Power Law constitutive model is usually applicable to analysis of deformation behavior of materials at high temperatures and high strain rates, while the Johnson-Cook constitutive model is more suitable for analysis of deformation behavior of materials at low to medium temperatures and low to medium strain rates. Therefore, the deformation behavior of Inconel718 at high temperatures and high strain rates can be described using the Power Law constitutive model. The mathematical equation of Power Law model could be expressed as:

$$\sigma \left({\epsilon }^p,\dot{\epsilon },T\right)=g\left({\epsilon }^p\right)\Gamma (\dot{\epsilon })\Theta (T)$$

(9)

In equation (9), the strain enhancement function $g\left({\epsilon }^p\right)$ could be further expressed as:

$$\{\begin{array}{l}g\left({\epsilon }^p\right)={\sigma }_0{\left(1+\frac{{\epsilon }^p}{{\epsilon }_0^p}\right)}^{1/n}{\epsilon }^p<{\epsilon }_{cut}^p\\ g\left({\epsilon }^p\right)={\sigma }_0{\left(1+\frac{{\epsilon }_{cut}^p}{{\epsilon }_0^p}\right)}^{1/n}{\epsilon }^p\ge {\epsilon }_{cut}^p\end{array}$$

(10)

The strain enhancement parameters of Inconel718 are shown in Table 3.

In equation (9), the thermal softening function $\Theta (T)$ could be expressed as:

$$\{\begin{array}{l}\Theta (T)=c_0+c_{1}T+c_2{T}^2+c_3{T}^3+c_4{T}^4+c_5{T}^6 \mathrm{T}<T_{\mathrm{cut} }\\ \Theta (T)=\Theta \left(T_{\mathrm{cut} }\right)\left(1-\frac{T-T_{\mathrm{cut} }}{T_{\mathrm{melt} }-T_{\mathrm{cut} }}\right) \mathrm{T}\ge T_{\mathrm{cut} }\end{array}$$

(11)

The thermal softening parameters of Inconel718 are shown in Table 4.

In equation (9), the strain rate effect function $\Gamma (\dot{\epsilon })$ could be expressed as:

$$\{\begin{array}{l}\Gamma (\dot{\epsilon })={\left(1+\frac{\dot{\epsilon }}{{\dot{\epsilon }}_0}\right)}^{\frac{1}{m_1}}\dot{\epsilon }\le {\dot{\epsilon }}_t\\ \Gamma (\dot{\epsilon })={\left(1+\frac{\dot{\epsilon }}{{\dot{\epsilon }}_0}\right)}^{\frac{1}{m_2}}{\left(1+\frac{{\dot{\epsilon }}_t}{{\dot{\epsilon }}_0}\right)}^{\left(\frac{1}{m_1}-\frac{1}{m_2}\right)}\dot{\epsilon }>{\dot{\epsilon }}_t\end{array}$$

(12)

The strain rate parameters of Inconel718 are shown in Table 5.

The damage model for Inconel718 can be determined by the following equation:

$$D={\displaystyle \sum _i\frac{\Delta {\epsilon }_i^p}{{\epsilon }_{f_i}^p}}$$

(13)

In equation (13), the parameter $\Delta {\epsilon }_i^p$ could be expressed as:

$${\epsilon }_{f_0}^p=d_0+d_{1}T+d_2{T}^2+d_3{T}^3+d_4{T}^4+d_5{T}^6$$

(14)

In equation (14), the parameters $d_0\sim d_5$ are the strain value deduced by tensile test.

3.2 The Result of Finite Element Analysis

In order to study the effect of different processing parameters such as spindle speed, feed per tooth, radial and axial depth of cut on the temperature, stress and cutting force generated during the milling process of Inconel718, parameters as shown in Table 6 are taken for finite element analysis.

3.2.1 Temperature Field Analysis of Milling Process

Inconel718 belongs to the difficult-to-machine materials, which has low thermal conductivity. The phenomenon of thermal mechanical coupling could occur during machining process, which results in residual stress and reduces the surface quality. Therefore, in order to effectively control tool wear and improve the surface integrity of the workpiece, the milling temperatures must be systematically studied and analyzed.

The fitted curves of the maximum temperature at different spindle speeds, feed per tooth, radial and axial depth of cut are plotted, according to the results of finite element analysis, these curves are shown in Figure 8 (a) – (d).

From Figure 8, it can be seen that the temperature increases as the spindle speed increases, this is mainly due to the fact that the friction phenomenon between the front (or back) of tool and workpiece becomes more intense as the spindle speed increases. Moreover, with the increase of radial depth of cut, axial depth of cut and feed per tooth, the temperature also increases, because rising of all these parameters makes the tool cut more material per unit of time, so the heat continues to accumulate.

3.2.2 Stress Field Analysis of Milling Process

Due to the high hardness and strength properties of Inconel 718, the milling process of this material is susceptible to stress concentration and crack formation. Stress field simulation enables the visualization of stress distribution on both the tool and workpiece surfaces during machining, facilitating a deeper understanding of the physical phenomena involved. Consequently, the impact of various machining parameters on stress can be predicted to mitigate stress concentration and prevent crack formation.

To investigate the impact of processing parameters on stresses during machining, the stress values in the top 10% of elements with the highest stresses on the tool are averaged and plotted as a function of time. Fitted curves representing stress variations at different spindle speeds, feed per tooth, radial and axial depths of cut are presented in Figure 9 (a) – (d).

From Figure 9, it is evident that an increase in spindle speed does not significantly affect the stress, as the material volume cut by the tool remains unchanged despite the enlargement of the spindle speed. Conversely, with the increment of feed per tooth, radial depth of cut, and axial depth of cut, there is a simultaneous rise in stress. This observation can be attributed to the larger volume of material being removed due to changes in these parameters.

3.2.3 Force Analysis of Milling Process

Given that the magnitude of the milling force directly influences the wear and fracture condition of the tool, a comprehensive analysis of the milling force generated during the Inconel 718 milling process becomes crucial.

The milling forces measured through milling experiment are presented in Figure 10 (a), and the milling forces calculated by using the theoretical milling force model are presented in Figure 10 (b). The milling forces of each axis obtained from the theoretical model align closely with the experimental results, validating the accuracy and correctness of the proposed model.

In order to better study the influence trend of different machining parameters on milling force, spindle speed is set to 1000rpm, 1300rpm and 1600rpm respectively, other parameters are kept unchanged, the milling force data of each axis are shown in Figure 11.

Figure 11 illustrates that the effect of different spindle speeds on the magnitude of the milling force for each axis is not significant. This can be attributed to the nearly constant volume of material being removed, regardless of the spindle speed variation.

The feed per tooth is varied at values of 0.2 mm, 0.3 mm, 0.4 mm, and 0.5 mm, while keeping the remaining parameters unchanged. The corresponding milling force data for each axis is presented in Figure 12.

Figure 12 demonstrates a clear and consistent increasing trend in milling forces for each axis with the rise in feed per tooth. This observation can be attributed to the larger volume of workpiece material being removed by the cutter as the feed per tooth increases.

The radial depth of cut is varied at values of 2 mm, 3 mm, 4 mm, and 5 mm, while keeping the remaining parameters unchanged. The corresponding milling force data for each axis is presented in Figure 13.

Figure 13 reveals that the milling forces in the X-Axis, Y-Axis, and Z-Axis exhibit a notable increase with the rise in radial depth of cut. This observation can be attributed to the corresponding increase in the material removal rate, which occurs as the radial depth of cut is enlarged.

Similarly, the value of axial depth of cut is set to 2mm, 3mm, 4mm and 5mm respectively, and the rest of parameters are kept unchanged, the milling force data of each axis are shown in Figure 14.

4. Analysis of Milling Stability

Inconel718 belongs to difficult-to-machine materials, and is susceptible to the phenomenon of chattering during milling process. Chattering in the milling process will lead to a positive feedback effect between the cutting force and vibration, which will result in degradation of milling quality, cutting edge wear, tool fracture and other problems. Through the stability analysis of chattering, the stability condition of the milling process can be predicted, and whether chattering will occur during the milling process could also be determined. Therefore, it is necessary to carry out chattering stability analysis of Inconel718, so as to improve the milling efficiency and machining quality as well as reduce the production cost.

As shown in Figure 15, the machine-tool system is simplified as a vibration system with two degrees of freedom, where the vibrations in the X and Y directions can be described by the kinetic differential equations:

$$\{\begin{array}{l}{m}_x\ddot{x}+c_x\dot{x}+k_{x}x={\displaystyle \sum _{j=1}^N{F}_{xj}}=F_x(t)\\ m_y\ddot{y}+c_y\dot{y}+k_{y}y={\displaystyle \sum _{j=1}^N{F}_{yj}}=F_y(t)\end{array}$$

(15)

In the milling process, as shown in Figure 16, the current cutting process will be affected by the previous cutting process, since the previous machining process leaves vibration on the machined surface, and when the tool cut the surface again can lead to changes in the cutting thickness. Such changes in the cutting thickness can cause fluctuations in cutting force, leading to further vibration in the system. The total cutting thickness generated by the superposition of the static cutting thickness and the dynamic cutting thickness caused by tool vibration can be expressed as:

$$h\left({\varphi }_j\right)=\left[f_t\sin {\varphi }_j+\left(v_{j,0}-v_j\right)\right]g\left({\varphi }_j\right)$$

(16)

In equation (16), the parameter $g\left({\varphi }_j\right)$ is a unit step function, which is expressed below:

$$g\left({\varphi }_j\right)=\{\begin{array}{l}1 , \left({\varphi }_{st}<{\varphi }_j<{\varphi }_{ex}\right)\\ 0, \left({\varphi }_j<{\varphi }_{st}，{\varphi }_j>{\varphi }_{ex}\right)\end{array}$$

(17)

Therefore, the total cutting force on the cutter can be expressed as:

$$\left\{\begin{array}{c}{F}_x\\ F_y\end{array}\right\}=\frac{1}{2}{a}_p{K}_{tc}\left[\begin{array}{cc}{\alpha }_{xx}& {\alpha }_{xy}\\ {\alpha }_{yx}& {\alpha }_{yy}\end{array}\right]\left\{\begin{array}{c}\Delta x\\ \Delta y\end{array}\right\}$$

(18)

In equation (18) the parameters ${\alpha }_{xx}$, ${\alpha }_{xy}$, ${\alpha }_{yx}$ and ${\alpha }_{yy}$ are average directional force coefficient:

$$\{\begin{array}{l}{\alpha }_{xx}={\displaystyle \sum _{j=0}^{N-1}-g_j\left[\sin 2{\varphi }_j+K_r\left(1-\cos 2{\varphi }_j\right)\right]}\\ {\alpha }_{xy}={\displaystyle \sum _{j=0}^{N-1}-g_j\left[\left(1+\cos 2{\varphi }_j\right)+K_r\sin 2{\varphi }_j\right]}\\ {\alpha }_{yx}={\displaystyle \sum _{j=0}^{N-1}{g}_j\left[\left(1-\cos 2{\varphi }_j\right)+K_r\sin 2{\varphi }_j\right]}\\ {\alpha }_{yy}={\displaystyle \sum _{j=0}^{N-1}{g}_j\left[\sin 2{\varphi }_j-K_r\left(1+\cos 2{\varphi }_j\right)\right]}\end{array}$$

(19)

Thus, the equation of dynamic milling force can be expressed as:

$$\left\{F\left(t\right)\right\}=\frac{1}{2}{a}_p{K}_{tc}\left[A\left(t\right)\right]\left\{\Delta \left(t\right)\right\}$$

(20)

Equation (20) can be further simplified as:

$$\left\{F\left(t\right)\right\}=\frac{1}{2}{a}_p{K}_{tc}\left[A_0\right]\left\{\Delta \left(t\right)\right\}$$

(21)

In equation (21), the parameter $\left[A_0\right]$ can be expressed as:

$$\left[A_0\right]=\frac{1}{T}{\displaystyle {\int }_0^T\left[A\left(t\right)\right]dt}=\frac{1}{{\varphi }_p}{\displaystyle {\int }_{{\varphi }_{st}}^{{\varphi }_{ex}}\left[A\left(\varphi \right)\right]d\varphi }=\frac{N}{2\pi }\left[\begin{array}{cc}{\alpha }_{xx}& {\alpha }_{xy}\\ {\alpha }_{yx}& {\alpha }_{yy}\end{array}\right]$$

(22)

The stability of the system is solved using the zero-order frequency domain method, such that the frequency response transfer function matrix of the tool-workpiece contact zone is:

$$\left[\Phi \left(i\omega \right)\right]=\left[\begin{array}{cc}{\Phi }_{xx}\left(i\omega \right)& {\Phi }_{xy}\left(i\omega \right)\\ {\Phi }_{yx}\left(i\omega \right)& {\Phi }_{yy}\left(i\omega \right)\end{array}\right]$$

(23)

In order to solve the dynamic equation, it is necessary to obtain the modal parameters of the tool, including the natural frequency, damping ratio, and vibration mode. The approach of force hammer excitation is used in this experiment to obtain accurate modal parameters.

The schematic hardware installation of system in the modal experiment is shown in Figure 17, and the photo of connection of experimental equipment is shown in Figure 18. The experimental equipment includes a carbide ball-end milling cutter, a three-way acceleration sensor, a TST5928 Dynamic Strain Gauge, and a computer equipped with the TST5928 Distributed Dynamic Signal Test and Analysis System.

Figure 19 shows the distribution of measurement points of the tool, a force hammer is used to excite the tool at different measurement points and the response of the tool is measured. Since the stiffness of tool in Z direction is much greater than in X and Y directions, only modal experiments in X and Y directions are performed.

The results of modal experiments are recorded in Table 7.

Figure 20 and Figure 21 show the tool modal of X and Y direction respectively.

Combining results of Figure 20 and Figure 21 as well as Table 7, it is possible to derive that the modal shapes corresponding to first-order natural frequency in X and Y directions of the tool are both bent and deformed in middle part of the tool. The modal shapes corresponding to second-order natural frequency in X and Y directions of the tool are both bent and deformed in upper part of the tool, and the degree of deformation in the Y direction is slightly larger than that in the X direction. The bending deformation of the modal shapes corresponding to third-order natural frequency in X-direction occurs in middle and upper part of the tool and the deformation of the bottom of the tool is small, while the bending deformation of the modal shapes corresponding to third-order natural frequency in Y-direction occurs in the lower part of the tool.

As it is necessary to select suitable processing parameters to avoid chattering during the milling process, therefore, there is a need to investigate the effect of tool radius, radial depth of cut, natural frequency ${\omega }_n$ and damping ratio $\zeta$ on cutting stability. The stability lobes are plotted for different tool radius, radial depths of cut, natural frequencies ${\omega }_n$ and damping ratios $\zeta$, as shown in Figure 22, Figure 23, Figure 24 and Figure 25.

5. Milling Experiment of Inconel718 and Parameter Optimization

In order to verify the correctness of the established milling force model and finite element analysis, a series of experiments including single factor experiments and orthogonal experiments are designed. After that, the milling force measured from the experiments with that calculated from the theoretical model as well as the milling forces obtained from the finite element analysis are compared.

A DMC635V vertical machining center, 1B240-0800-XA 1630 carbide ball-end milling cutter, Kistler 9527B force gauge, 5073A charge amplifier and acquisition card as well as Dynoware computer analysis software are selected to build the experimental platform. The Schematic diagram of connection of experimental equipment is shown in Figure 26, and the device connection diagram is shown in Figure 27.

The single factor experiments and orthogonal experiments are designed as shown in Table 8 and Table 9.

5.1 Results Analysis of Single Factor Experiment

To investigate the influence of spindle speed on the milling of Inconel 718, the parameter of spindle speed is varied at 1000 rpm, 1200 rpm, 1400 rpm, 1600 rpm, and 1800 rpm, while keeping the other machining parameters unchanged. The axial depth of cut is set to 0.4 mm, and the feed per tooth is set to 0.02 mm/z. The experimental results of the milling force obtained under these conditions are presented in Figure 28.

Figure 28 illustrates that the average milling forces in the X-Axis, Y-Axis, and Z-Axis show minimal changes with increasing spindle speed, which aligns with the conclusions drawn from the finite element analysis.

To unveil the influence of axial depth of cut on the milling of Inconel 718, the parameter of axial depth of cut is varied at 0.2 mm, 0.3 mm, 0.4 mm, 0.5 mm, and 0.6 mm, while keeping the other machining parameters constant. The spindle speed is set to 1000 rpm, and the feed per tooth is set to 0.02 mm/z. The experimental results of the milling force obtained under these conditions are presented in Figure 29.

Figure 29 reveals that the average milling force in the X, Y, and Z-Axis all demonstrate an increasing trend with the augmentation of axial depth of cut. The Z-Axis exhibits the most significant increase, while the average milling forces in the X and Y-Axis show relatively smaller increments. These observations are consistent with the conclusions drawn from the finite element analysis.

To investigate the influence of feed per tooth on the milling of Inconel 718, the parameter of feed per tooth is varied at 0.02 mm/z, 0.04 mm/z, 0.06 mm/z, 0.08 mm/z, and 0.10 mm/z, while keeping the other machining parameters unchanged. The spindle speed is set to 1000 rpm, and the axial depth of cut is set to 0.4 mm. The experimental results of the milling force obtained under these conditions are presented in Figure 30.

It can be seen from Figure 30 that the average milling force in the X, Y, and Z-Axis all reveal an increasing trend with the rising of feed per tooth, in which the increase in the Z-Axis is the largest, and the increase in X and Y is relatively small, which is compatible with the conclusions drawn from the finite element analysis.

5.2 Results Analysis of Orthogonal Experiments

The orthogonal experiments are analyzed based on the parameters shown in Table 9 and the orthogonal experimental results obtained are shown in Table 10.

An extreme variance analysis of the results shown in Table 10 leads to the following conclusions:

• Axial depth of cut is the most influential factor on the average milling force in X-axis, followed by feed per tooth and spindle speed.
• Feed per tooth is the most influential factor on the average milling force in Y-axis, followed by axial depth of cut and spindle speed.
• Axial depth of cut is the most influential factor on the average milling force in Z-axis, followed by the feed per tooth and spindle speed.

Considering that the effect of feed per tooth and axial depth of cut on average milling force in the previous analysis is large, while the effect of spindle speed is relatively small, different parameters of axial depth of cut and feed per tooth are selected for verification experiments. Groups 1-4 shown in Table 9 are selected, the theoretical and finite element analysis value of average milling force of these groups are calculated in order to compare with the value of experiments shown in Table 10, the results are recorded in Table 11 and Table 12.

Table 11 and Table 12 reveal that the minimum relative error between the milling force calculated by the theoretical model and the experimental data is 6.48%, while the minimum relative error between the milling force obtained from finite element analysis and the experimental data is 4.24%. These results further demonstrate the reliability and validity of both the theoretical model and the finite element analysis.

5.3 Establishment of Optimization Objective Function

To enhance the machining efficiency and tool life expectancy of Inconel 718, optimizing the processing parameters becomes essential. In this study, surface roughness, material removal rate, and their combination are considered as the optimization objectives. The fmincon algorithm, a nonlinear minimization method employing the interior point approach, is utilized to determine the optimal processing parameters. Known for its high accuracy in solving nonlinear optimization problems, the fmincon algorithm is well-suited for this task. The optimization procedure employing the fmincon algorithm is depicted in Figure 31.

Surface roughness is an important evaluation index in the machining process of Inconel718, therefore, a model with surface roughness as the optimization objective function is developed, and the optimization objective function can be expressed as:

$$R_a=C\cdot n^{c_1}\cdot f_z^{c_2}\cdot a_p^{c_3}$$

(24)

In equation (24), the parameter $C$ is influenced by both geometry of the milling cutter and material of workpiece. Parameters $c_1$, $c_2$ and $c_3$ are the influencing coefficients of the milling parameters.

Taking logarithms on both sides of equation (24) yields:

$$\ln R_a=\ln C+c_1\ln n+c_2\ln f_z+c_3\ln a_p$$

(25)

Taking $y=\ln R_a$, $c_0=\ln C$, $x_1=\ln n$, $x_2=\ln f_z$, $x_3=\ln a_p$, and substituting these parameters into equation (25) yields:

$$y=c_0+c_1{x}_1+c_2{x}_2+c_3{x}_3$$

(26)

Based on that, multiple linear regression equations are built sequentially:

$$\{\begin{array}{l}{y}_1={\beta }_0+{\beta }_1{x}_{11}+{\beta }_2{x}_{12}+{\beta }_3{x}_{13}\\ y_2={\beta }_0+{\beta }_1{x}_{21}+{\beta }_2{x}_{22}+{\beta }_3{x}_{23}\\ \hspace{1em} ⋮\\ y_n={\beta }_0+{\beta }_1{x}_{n1}+{\beta }_2{x}_{n2}+{\beta }_3{x}_{n3}\end{array}$$

(27)

Equation (27) is changed to a matrix representation:

$$Y=X\beta +\epsilon$$

(28)

The parameters in equation (28) can be expressed as:

$$Y=\left[\begin{array}{l}{y}_1\\ y_2\\ ⋮\\ y_4\end{array}\right],X=\left[\begin{array}{l}1 x_{11} x_{12} x_{13} x_{14}\\ 1 x_{21} x_{22} x_{23} x_{24}\\ ⋮\\ 1 x_{161} x_{162} x_{163} x_{164}\end{array}\right],\beta =\left[\begin{array}{l}{\beta }_0\\ {\beta }_1\\ ⋮\\ {\beta }_4\end{array}\right],\epsilon =\left[\begin{array}{l}{\epsilon }_1\\ {\epsilon }_2\\ ⋮\\ {\epsilon }_4\end{array}\right]$$

(29)

Using experimental results in reference[41], the optimization objective function (24) can be changed to:

$$f_1\left(x\right)=R_a=\min \left(R_a\left(n,f_z,a_p\right)\right)$$

(30)

In metal milling, improving productivity is critical, and material removal rate is usually used to measure the productivity. Therefore, a model with material removal rate as the optimization objective function is developed:

$$f_2\left(x\right)=\max \left(MRR\right)=\max \left(n\cdot f_z\cdot a_p\cdot N\right)$$

(31)

The parameter $n$ in equation (31) can be expressed as:

$$n=\frac{1000v_c}{\pi D}$$

(32)

In order to ensure the machining stability, the processing parameters of Inconel718 need to be selected within the stability region of stability lobes as shown in Figure 22 - Figure 25.

5.4 Optimization Results Analysis

When minimum surface roughness is used as the optimization objective, it can be observed from Figure 32 that the solution of optimization objective function gradually converges as the number of iterations increases. According to the optimal milling parameters in Table 13, it can be concluded that the surface roughness reaches the optimal value at the 19^th generation, and the optimal surface roughness is 0.43$\mu m$.

When the maximum material removal rate is used as the optimization objective, it can be seen from Figure 33 that the solution of the optimization objective function gradually converges as the number of iterations increases. According to the optimal milling parameters in Table 14, it can be concluded that the material removal rate reaches the optimal value at the 11^th generation, the optimal material removal rate is 58788.32(mm³/min).

If the maximum material removal rate and minimum surface roughness are taken as the research objective of multi-objective optimization, the mathematical model of multi-objective optimization is as follow:

$$\begin{array}{l}\min f(x)=\min ({\omega }_1{f}_1(x)+{\omega }_2{f}_2(x))\\ s.t.\{\begin{array}{l}{n}_{\min }\le n\le n_{\max }\\ f_z\le \frac{\pi D{v}_{f\max }}{1000vz}\\ a_{p\min }\le a_p\le a_{p\max }\end{array}\end{array}$$

(33)

In equation (33) the parameter ${\omega }_1$ and ${\omega }_2$ are weighting values of machining efficiency and surface roughness respectively.

It can be seen from Figure 34 that the solution of the optimization objective function gradually converges as the number of iterations increases. According to the optimal milling parameters in Table 15, it can be concluded that the optimal material removal rate reaches 4199.2 mm³/min, while the optimal surface roughness reaches 3.5$\mu m$.

$\mathrm{Surface} \mathrm{roughness} (\mu m$

6. Conclusion

This paper takes Inconel718 as the research object, and carries out an in-depth discussion on the key problems in the process of machining Inconel718, specifically including milling force, milling temperature, milling stability, and optimization of processing parameters. The main conclusions are as follow:

• An increase in feed and tool helix angle leads to an increase in milling force, and the results of finite element analysis indicate that with the increase of the spindle speed, radial depth of cut, axial depth of cut and feed per tooth, the increase of temperature in the machining area becomes more conspicuous. However, with the increase of spindle speed, the change of tool stress is not obvious, but with the increase of feed per tooth, radial depth of cut and axial depth of cut, the tool stress increases significantly. The spindle speed does not have a significant effect on the milling force in each axis. With the increase of feed per tooth and radial depth of cut, the milling force in X-, Y- and Z-Axis increases significantly, herein, the axial depth of cut has a greater effect on milling force of X-Axis and a smaller effect on milling forces of both Y- and Z-Axis.
• An increase in tool radius improves stability of milling system, however, an increase in radial depth of cut is detrimental to the stability. With increase of natural frequency, the critical depth of cut of milling remains unchanged, however, increase in damping ratio improves the system stability.
• The order of influence on the average milling force in X-axis is: axial depth of cut > feed per tooth > spindle speed, and the order of influence on the average milling force in Y-axis is: spindle speed > feed per tooth > axial depth of cut.
• Using the fmincon algorithm, if the minimum surface roughness is considered as the optimization objective, the optimal machining parameters are: spindle speed 3999.63rpm, feed rate 80.01mm/min, axial depth of cut 0.25mm, at this time, the optimal value of surface roughness is 0.43 μm. If the maximum material removal rate is considered as the optimization objective, the optimal machining parameters are: spindle speed 4000rpm, feed speed 700mm/min, axial depth of cut 2.54mm, at this time, the optimal value of material removal rate is 58788.32mm³/min. If both the minimum surface roughness and maximum material removal rate are considered as the optimization objective, the optimal processing parameters are: spindle speed 3199.2rpm, feed speed 80mm/min, axial depth of cut 0.25mm. At this time, the value of surface roughness is 3.5 μm, and the value of material removal rate is 4199.2mm³/min.