Mechanistic model and probability characteristics of micro-milling force with a new parameter identification method

Abstract

As a critical part of the analysis of the micro-milling process, mechanical modeling directly affects the prediction of the surface finish of machining materials and the cutting stability. In this paper, a flexible force model of micro-milling is established, which accounts for the influence of the actual cycloid tool path, tool runout, elastic recovery, tool deformation, and chip separation state. The Harris Hawks Error Optimization parameter identification method is proposed to efficiently obtain the cutting force coefficient/parameter. Based on the uncertainty of system parameters caused by the manufacturing and assembly process, probability distribution characteristics of the micro-milling force are predicted and analyzed. Additionally, the Adaptive Kriging method is used to reconstruct the complex implicit relationship between the cutting parameters and milling forces, improving the computational efficiency of the probability analysis. Finally, a series of micro-milling experiments and computational results verified the accuracy of the proposed method.

Introduction

With the characteristics of small machining size, a wide range of machinable materials, complex surface machining, and low cost, micro-milling has become an essential technology in the manufacturing fields such as aerospace, biomedical, automotive, electronic, and others, which is the current research hotspot.^1–3 However, micro-milling is easily affected by minor variations such as tool deformation and instantaneous cutting thickness. ⁴ Consequently, understanding and modeling the micro-milling process are critical when aiming to improve the quality of machined parts and production efficiency.

Compared to macro-machining, micro-machining considers problems at a small scale (micro problems), making it difficult to study the underlying micro-milling mechanic mechanisms. Wang et al. ⁵ redefined and calculated the cutting radius and instantaneous uncut chip thickness in the machining process to predict the cutting force in the micro end milling process. Further, Guo et al. ⁶ proposed an instantaneous undeformed thickness model for micro-milling considering tool runout, variable pitch, and helix angle. Based on multi-scale modeling and analysis, Niu and Cheng ⁷ proposed an innovative expression of micro machining cutting force. Next, Liu et al. ⁸ investigated a micro-milling mathematical model of the thin-walled copper structure and verified the proposed model using finite element simulation. Considering the tool size, tool runout, and spatial tool position, Jing et al. ⁹ and Zhang et al. ¹⁰ established more accurate micro-milling tool runout and instantaneous force models. All of the above-listed models described the micromachining processes under some influencing factors, providing a basis for establishing a more comprehensive micro-milling force model, which will consider multiple related factors.

To create a more accurate cutting force model, the factors such as chip separation state, instantaneous cutting thickness, tool runout, workpiece elastic recovery, effectively included angle, and the interaction between tool deformation and cutting forces are integrated into the modeling. In addition, the accurate and efficient measurement or calculation of the cutting coefficient/parameter in the mechanical model is essential. Gonzalo et al. ¹¹ calculated the specific cutting force coefficient by substituting the identified simulated force and measured force for reverse calculation into the linear mechanics model. Jin and Altintas ¹² and Tsai et al. ¹³ predicted the cutting force coefficient based on the finite element analysis and milling experimental results. By combining experimental measurement with numerical calculation, ¹⁴ Li et al. ¹⁵ analyzed the ratio of total force corresponding to each cutting edge according to the spatial geometric relationship to determine the tool runout parameters. However, using finite element analysis, cyclic iteration, and experimental methods for calculating cutting the coefficient/parameters is relatively complex and somewhat time-consuming.

Hence, the Harris Hawkes Error Optimization (HEO) method is developed to identify cutting coefficients/parameters both concisely and efficiently. The square difference between the predicted and experimental forces is used to construct the optimization objective function. The process parameters are usually regarded as deterministic variables in the prediction methods of micro-milling force, but these parameters are variable and affect the machining process. ¹⁶ Consequently, in this paper, the probability characteristics of the micro-milling force are studied. The Adaptive-Kriging (AK) method is used to reconstruct the functional relationship between the machining parameters and the milling force to improve the calculation efficiency. Comparing the predicted results of the AK method and directly proposed model calculations verify the accuracy of the probabilistic characteristic analysis conclusions.

The paper is structured as follows: a novel comprehensive cutting force model, including the effects above, is proposed in Section 2. In Section 3, the HEO algorithm is presented and used to identify the cutting coefficients/parameters. Section 4 illustrates the AK method used in probability analysis. The experimental device is designed to verify the proposed prediction model and AK model, and the probability distribution characteristics of the micro-milling force are obtained in Section 5. Section 6 summarizes the conclusions of this work.

Mechanical modeling of micro-milling

Chip thickness analysis

The radial variation of the actual cutting point results in periodical changes in the force acting on the cutting edge, which can be defined as tool runout. As shown in Figure 1(a), the existence of runout causes the geometric central axis of the tool to no longer be collinear with the central rotation axis, resulting in the change of the actual tool position and rotation radius. The actual tool rotation radius can be calculated as:

$${R_{\left(k\right)}}^2={(R_0\cos ({\gamma }_{\mathrm{o}}+2\pi (k-1)/K)+r_{\mathrm{o}})}^2+{\left(R_0\sin {\gamma }_{\mathrm{o}}\right)}^2$$

(1)

Figure 1.

Micro geometric relationship: (a) cutting tool center offset, (b) the geometry calculation of IUCT, and (c) chip formation mechanism.

where R₀ is the tool radius, k is the number of the current flute, K is the total number of tools, γ_o and r_o are the runout angle and value. The trajectories O of the tool center are formulated as:

$$(x_O,y_O)=(f{t}_k/60+r_o\sin (\omega t_k+{\gamma }_o),r_o\cos (\omega t_k+{\gamma }_o))$$

(2)

where f is the feed rate, ω is the angular velocity. The trajectory of the kth edge can be written as:

$$\{\begin{array}{c}{x}_k=f{t}_k/60+R\sin (\omega t_k-2\pi k/K)+r_{\mathrm{o}}\sin (\omega t_k+{\gamma }_{\mathrm{o}})\hfill \\ \multicolumn{1}{c}{y_k=R\cos (\omega t_k-2\pi k/K)+r_{\mathrm{o}}\cos (\omega t_k+{\gamma }_{\mathrm{o}})}\end{array}$$

(3)

According to the geometric relationship in Figure 1, the kth tooltip can be denoted as M (x, y) with the tool center located at O _k . The line MO_k can be expressed as follows:

$$\begin{array}{c}\multicolumn{1}{c}{Y=\frac{y_k-y_{O_k}}{x_k-x_{O_k}}X=\frac{X-f{t}_k/60-r_o\sin (\omega t_k+{\gamma }_{\mathrm{o}})}{\tan (\omega t_k-2\pi k/K)}}\\ \multicolumn{1}{c}{+R\cos (\omega t_k+{\gamma }_{\mathrm{o}})}\end{array}$$

(4)

The coordinates of the k−1th cutting edge M_k−1(x_k−1, y_k−1) are obtained based on the above kth tooth trajectory analysis presented earlier:

$$\{\begin{array}{c}{x}_{k-1}=f{t}_{k-1}/60+R\sin (\omega t_{k-1}-2\pi (k-1)/K)\hfill \\ \multicolumn{1}{c}{+r_{\mathrm{o}}\sin (\omega t_{k-1}+{\gamma }_{\mathrm{o}})}\\ \multicolumn{1}{c}{y_{k-1}=R\cos (\omega t_{k-1}-2\pi (k-1)/K)}\\ \multicolumn{1}{c}{+r_{\mathrm{o}}\cos (\omega t_{k-1}+{\gamma }_{\mathrm{o}})}\end{array}$$

(5)

It can be seen from Figure 1(b) that the instantaneous uncut chip thickness (IUCT) is the min-distance of the intersection between the cutting trajectory of the current and the previous tooth. Let M_k−1 be the intersection of the MO_k and the k−1th cutting trajectory. Substituting equation (5) into the linear MO_k equation of equation (4):

$$\begin{array}{c}\multicolumn{1}{c}{G(t\prime )=R\tan (\omega t_k-2\pi k/K)\cos (\omega t_{k-1}-2\pi (k-1)/K)}\\ \multicolumn{1}{c}{+f(t_k-t_{k-1})/60+r_o\tan (\omega t_k-2\pi k/K)}\\ \multicolumn{1}{c}{[\cos (\omega t_{k-1}+{\gamma }_o)-\cos (\omega t_k+{\gamma }_o)]}\\ \multicolumn{1}{c}{-R\sin (\omega t_{k-1}-2\pi (k-1)/K)}\\ \multicolumn{1}{c}{-r_o\sin (\omega t_{k-1}+{\gamma }_o)+r_o\sin (\omega t_k+{\gamma }_o)=0}\end{array}$$

(6)

The nonlinear equations related to time t′ are solved by the Newton Raphson iteration method, ¹⁷ and the initial time t_i iteration update can be expressed as:

$$t_{i+1}\prime =t_i\prime -G\left(t_i\prime \right)/G\left(t_i\prime \right)\prime =t-2\pi /\omega K-G\left(t_i\prime \right)/G\left(t_i\prime \right)\prime$$

(7)

where G(t_i′)′ is the first derivative of G(t_i′). The flow of solving the IUCT can be derived as follows:

$$\begin{array}{c}\multicolumn{1}{c}{h_u=H-U}\\ \multicolumn{1}{c}{=(R_{\left(k\right)}+\sqrt{{(x_k-x_{k-1})}^2+{(y_k-y_{k-1})}^2}\sin (\omega t_k-2\pi k/K+\alpha ))}\\ \multicolumn{1}{c}{-\sqrt{{R_{\left(k\right)}}^2-L^2\cos {(\omega t_k-2\pi k/K+\arctan [(y_{O_k}-{y_O}_{k-1})/(x_{O_k}-x_{O_{k-1}})])}^2}}\end{array}$$

(8)

As shown in Figure 1(c), the chips are not formed when the chip thickness is less than the minimum chip thickness (MCT), and part of material elastic recovery except for the chip. The actual IUCT requires superposition calculation, which can be expressed as:

$$h_u\left(k\right)=\{\begin{array}{c}{h}_u\left(k\right)+P_e\times h_u(k-1)h_u\left(k\right)\le h_{min}=\lambda r_e\hfill \\ \multicolumn{1}{c}{h_u\left(k\right)h_u\left(k\right)>h_{min}=\lambda r_e}\end{array}$$

(9)

where P_e is the elastic recovery rate of the workpiece, h_min is the MCT, and λ is the experimentally-obtained coefficient. ¹⁸ Figure 2 shows the instantaneous cutting thickness curves of two tooth cutters considering different work conditions.

Figure 2.

IUCT with different work conditions (ψ_k is the cutting edge position angle).

Due to the existence of the MCT, the IUCT curve suddenly increases at tooth 1, resulting in a superposition effect, which will gradually disappear as the runout or feed per tooth increases. Moreover, the IUCT increases with the growth of elastic recovery rate (P_e), runout length (r_o), and feed per tooth (f_s) and decreases with the rise of runout angle (γ_o). The change of the tool helix angle (β) will change the cutting tool working area, not the size of the IUCT. In addition, the effect of the runout on the instantaneous cutting thickness gradually decreases with the f_s increases.

Tool deformation analysis

As shown in Figure 3, the cutting is divided into the shear/plough dominant region (SDR/PDR). The chip flow direction is the same as that of the effective rake angle. The force acting on AC can be expressed as:

$$\{\begin{array}{c}d{F}_{\mathrm{ns}}=K_{\mathrm{ns}}h\left({\psi }_{\mathrm{k}}\right)\mathrm{dz}\cos \left({\alpha }_e\right)+K_{\mathrm{fs}}h\left({\psi }_{\mathrm{k}}\right)\sin \left({\alpha }_e\right)\hfill \\ \multicolumn{1}{c}{d{F}_{\mathrm{fs}}=-K_{\mathrm{ns}}h\left({\psi }_{\mathrm{k}}\right)\mathrm{dz}\sin \left({\alpha }_e\right)+K_{\mathrm{fs}}h\left({\psi }_{\mathrm{k}}\right)\cos \left({\alpha }_e\right)}\end{array}$$

(10)

Figure 3.

The diagram of the SDR-(a) and RDR-(b).

where K_ns and K_fs are the shear cutting coefficients, dz is the height of the cutting edge discrete elements. Through geometric analysis and integration simplification, the instantaneous effective rake angle α_e can be written as:

$${\alpha }_e=\{\begin{array}{c}\alpha \frac{h\left({\psi }_{\mathrm{k}}\right)+r_e[(\frac{\alpha }{2}-\frac{{\pi }^2}{8\alpha }-\frac{(1/\cos (1-h_{min}/r_e)+\pi )}{2\alpha \cos (1-h_{min}/r_e)})\cos \left(\alpha \right)-\sin \left(\alpha \right)-1]}{h\left({\psi }_{\mathrm{k}}\right)+r_e[(\alpha +\pi /2-1/\cos (1-h_{min}/r_e))\cos \left(\alpha \right)-\sin \left(\alpha \right)-1]}h\left({\psi }_{\mathrm{k}}\right)\ge r_e(1+\sin \left(\alpha \right))\hfill \\ \multicolumn{1}{c}{(1/\cos (1-h\left({\psi }_{\mathrm{k}}\right)/r_e)+1/\cos (1-h_{min}/r_e)-\pi )/2h\left({\psi }_{\mathrm{k}}\right)\ge r_e(1+\sin \left(\alpha \right))}\end{array}$$

(11)

The plough force acting on AB is proportional to the MCT. Its components in the direction of the cutting speed and the radial direction are given as:

$$\{\begin{array}{c}d{F}_{\mathrm{np}}=K_{\mathrm{np}}{h}_{min}\mathrm{dz}\sin \left(\zeta \right)+K_{\mathrm{fp}}{h}_{min}\mathrm{dz}\cos \left(\zeta \right)\hfill \\ \multicolumn{1}{c}{d{F}_{\mathrm{fp}}=K_{\mathrm{np}}{h}_{min}\mathrm{dz}\cos \left(\zeta \right)-K_{\mathrm{fp}}{h}_{min}\mathrm{dz}\sin \left(\zeta \right)}\\ \multicolumn{1}{c}{\zeta =0.5/\cos ((r_e-h_{min})/r_e)}\end{array}$$

(12)

where K_np and K_fp are plough coefficients and ζ is the angle between AB and the horizontal plane. As shown in Figure 3(b), the cutting state does not form chips, and the material accumulates at the height h_ac in front of the blade. The cutting force components can be calculated as:

$$\{\begin{array}{l}d{F}_{np}\prime =K_{np}\left(h\left({\psi }_k\right)+h_{ac}\right)dz\cos \left({\alpha }_p\right)+K_{fp}\left(h\left({\psi }_k\right)+h_{ac}\right)dz\sin \left({\alpha }_p\right)\\ d{F}_{fp}\prime =K_{np}\left(h\left({\psi }_k\right)+h_{ac}\right)dz\sin \left({\alpha }_p\right)-K_{fp}\left(h\left({\psi }_{\text{k}}\right)+h_{ac}\right)dz\cos \left({\alpha }_p\right)\\ {\alpha }_p=1/\cos \left(\sqrt{\left(h\left({\psi }_k\right)+h_{bf}\right)/\left(2r_e\right)}\right)\end{array}$$

(13)

where α_p is the instantaneous rake angle. There is a linear relationship between the height of the accumulation under plowing and the rheological factor ¹⁹ :

$$\begin{array}{c}\multicolumn{1}{c}{h_{\mathrm{ac}}=(0.0101806E_y\sqrt{h_{min}/(2r_e-h_{min})}/{\sigma }_0-0.040468)}\\ \multicolumn{1}{c}{h\left({\psi }_{\mathrm{k}}\right)}\end{array}$$

(14)

where E_y and σ₀ are the Young’s modulus and yield stress of the workpiece. The friction can be decomposed in radial and thrust directions:

$$\{\begin{array}{c}d{F}_{\mathrm{ne}}=-K_{\mathrm{ne}}{P}_{e}h\left({\psi }_{\mathrm{k}}\right)\mathrm{dz}\sin \left({\zeta }_e\right)+K_{\mathrm{ne}}{P}_{e}h\left({\psi }_{\mathrm{k}}\right)\mathrm{dz}\cos \left({\zeta }_e\right)\hfill \\ \multicolumn{1}{c}{d{F}_{\mathrm{fe}}=K_{\mathrm{ne}}{P}_{e}h\left({\psi }_{\mathrm{k}}\right)\mathrm{dz}\cos \left({\zeta }_e\right)+K_{\mathrm{ne}}{P}_{e}h\left({\psi }_{\mathrm{k}}\right)\mathrm{dz}\sin \left({\zeta }_e\right)}\end{array}$$

(15)

where K_ne and K_fe are coefficients of friction. The clearance angle of the blade γ, and the angle ζ_e can be determined from:

$${\zeta }_e=\{\begin{array}{c}\arctan \left(\sqrt{h_{\mathrm{er}}/(2r_e-h_{\mathrm{er}})}\right)h_{\mathrm{er}}\le r_e(1-\cos \left(\gamma \right))\hfill \\ \multicolumn{1}{c}{\arctan [h_{\mathrm{er}}\sin (\gamma )/(h_{\mathrm{er}}\cos }\\ \multicolumn{1}{c}{(\gamma )+r_e(1-\cos (\gamma )))]h_{\mathrm{er}}>r_e(1-\cos \left(\gamma \right))}\end{array}$$

(16)

where h_er is the elastic recovery height. Based on the calculations presented above and after dividing the cutting edge into M sections, the cutting forces in the x and y directions can be written as:

$$\begin{array}{c}\multicolumn{1}{c}{\sum _{i=1}^M\left[\begin{array}{c}d{F}_x\hfill \\ \multicolumn{1}{c}{d{F}_y}\end{array}\right]=\left[\begin{array}{c}-\cos \left({\psi }_k\right)-\sin \left({\psi }_k\right)\hfill \\ \multicolumn{1}{c}{-\sin \left({\psi }_k\right)\cos \left({\psi }_k\right)}\end{array}\right]\sum _{i=1}^M}\\ \multicolumn{1}{c}{\left[\begin{array}{c}d{F}_{\mathrm{ns}}+d{F}_{\mathrm{np}}+d{F_{\mathrm{np}}}^{\prime }+d{F}_{\mathrm{ne}}\hfill \\ \multicolumn{1}{c}{d{F}_{\mathrm{fs}}+d{F}_{\mathrm{fp}}+d{F_{\mathrm{fp}}}^{\prime }+d{F}_{\mathrm{fe}}}\end{array}\right]}\end{array}$$

(17)

Due to the small tool diameter, the effect of tool deflection on the IUCT in micro-milling should be considered. As shown in Figure 4(a), the tool flute forces are transmitted to the tool shank, which can be regarded as an arc with radius ρ with the action of force/moment ( F / M ). The bending moment and axial curvature can be expressed by numerical calculation:

$$\begin{array}{c}\multicolumn{1}{c}{{\rho }^{-1}=\mathit{M}\left(z\right)/\mathrm{EI}=\mathit{M}\left(z\right)/(E\pi r^4/4)=({\mathrm{d}}^2\eta \left(z\right)/\mathrm{d}{z}^2)}\\ \multicolumn{1}{c}{/{(1+{(\mathrm{d}\eta \left(z\right)/\mathrm{d}z)}^2)}^{3/2}}\end{array}$$

(18)

Figure 4.

Simplified model of the micro-milling tool: (a) Micro-milling cutting process and (b) Deformation process of micro-milling tool.

where η(z) is the deflection value of the tool, and I is the moment of inertia. The simplified beam deflection and rotation angle equations can be obtained:

$${\rho }^{-1}=\mathit{M}\left(z\right)/\mathrm{EI}={\mathrm{d}}^2\eta \left(z\right)/\mathrm{d}{z}^2={\eta }^{\prime \prime }$$

(19)

The overall force and deformation analysis of the tool are shown in Figure 4(b). The deformation at a point z on the tool flute under the effect of cutting force can be expressed as:

$$\{\begin{array}{c}{\eta }_{x(y)}={\eta }_1+{{\eta }_1}^{\prime }(L_m+L_t-z)+{\eta }_2+{{\eta }_2}^{\prime }(L_t-z)+{\eta }_3\hfill \\ \multicolumn{1}{c}{{\eta }_{x(y)}^{\prime }={{\eta }_1}^{\prime }+{{\eta }_2}^{\prime }+{{\eta }_3}^{\prime }}\end{array}$$

(20)

where η₁, η₁^′, η₂, η₂^′, η₃, and η₃^′are the deflection value and angle of the tool shank/neck/flute, which can be calculated by:

$$\{\begin{array}{c}{\eta }_1=F{\left(E\pi {R_r}^4\right)}^{-1}[(L_m+L_t-z)Z^2/2+F{(Z-L_m-L_t)}^3/6+A_{1}Z+A_2]L_m+L_t<Z\le L\hfill \\ \multicolumn{1}{c}{{\eta }_2=4F{\left(E\pi \right)}^{-1}\left\{\begin{array}{c}[R_r{[R_r-(Z-L_t)\tan \left(\tau \right)]}^{-2}/6-{[R_r-(Z-L_t)\tan \left(\tau \right)]}^{-1}/2]/\tan \left(\tau \right)\hfill \\ \multicolumn{1}{c}{+(L_t-z){[R_r-(Z-L_t)\tan \left(\tau \right)]}^{-2}/6+B_1[R_r-(Z-L_t)\tan \left(\tau \right)]+B_2}\end{array}\right\}{L}_t<Z\le L_m+L_t}\\ \multicolumn{1}{c}{{\eta }_3=4F{\left(E\pi {R_t}^4\right)}^{-1}{(<Z-z>}^3/6+C_{1}Z+C_2)0\le Z\le L_t}\end{array}$$

(21)

$$\{\begin{array}{c}{{\eta }_1}^{\prime }=4F{\left(\pi {R_r}^4\right)}^{-1}[(L_m+L_t-z)Z+F{(Z-L_m-L_t)}^2/2+A_1]L_m+L_t<Z\le L\hfill \\ \multicolumn{1}{c}{{{\eta }_2}^{\prime }=4F{\left(E\pi \right)}^{-1}\left\{\begin{array}{c}-[-R_r{[R_r-(Z-L_t)\tan \left(\tau \right)]}^{-3}/3+{[R_r-(Z-L_t)\tan \left(\tau \right)]}^{-2}/2]/\tan \left(\tau \right)\hfill \\ \multicolumn{1}{c}{+(L_t-z){[R_r-(Z-L_t)\tan \left(\tau \right)]}^{-3}/3-B_1}\end{array}\right\}{L}_t<Z\le L_m+L_t}\\ \multicolumn{1}{c}{{{\eta }_3}^{\prime }=4F{\left(E\pi {r_t}^4\right)}^{-1}({<Z-z>}^2/2+C_1)0\le Z\le L_t}\end{array}$$

(22)

where R_r is the shank radius, τ is the inclination angle, Z is the actual force point on the tool, and z (z_x/y(ψ_k)) is the equivalent force point. By substituting the boundary conditions into equations (21) and (22), the constant value parameters A₁, A₂, B₁, B₂, C₁, and C₂ can be solved. According to the equivalence relationship, the concentrated force moment can be expressed as the sum of the moments of all microelements:

$$\begin{array}{c}\multicolumn{1}{c}{z_{x/y}\left({\psi }_{\mathrm{k}}\right)=L-\left\{\sum _{i=1}^N\sum _{j=1}^{S}d{F}_{x/y}\left({\psi }_{\mathrm{k}}\right)[L-z_{x/y}\left({\psi }_{\mathrm{k}}\right)]\right\}/}\\ \multicolumn{1}{c}{F_{x/y}\left({\psi }_{\mathrm{k}}\right)}\end{array}$$

(23)

where dF_x/y(φ) is the microelement force in the x/y-axis. N and S are the numbers of flute and microelements.

Flexible force model

Due to the tool deformation-induced changes in the tooth trajectory and IUCT, the mechanical model must be updated. The current tool deflection value is used to update the position of the next tool center to obtain an accurate tool tip trajectory. As shown in Figure 5, the red and green dotted lines indicate the adjacent tooth tool path before the offset, while the solid lines show its path following the offset. Tooltip paths of the preceding and the current tooth before and after offset are (x_k^j+1, y_k^j+1) and (x_k−1^j+1, y_k−1^j+1), respectively, while the actual IUCT is designated as h_nk^j+1. Then the tool center trajectory in actual processing can be expressed as:

$$\{\begin{array}{c}({x_{\mathrm{ok}}}^{j+1},{y_{\mathrm{ok}}}^{j+1})=(x+{{\eta }_{\mathrm{kx}}}^{j+1},y+{{\eta }_{\mathrm{ky}}}^{j+1})\hfill \\ \multicolumn{1}{c}{({x_{\mathrm{ok}-1}}^{j+1},{y_{\mathrm{ok}-1}}^{j+1})}\\ \multicolumn{1}{c}{=(x^{\prime }+{{\eta }_{(k-1)x}}^{j+1},y^{\prime }+{{\eta }_{{(k-1)}_y}}^{j+1})}\end{array}$$

(24)

Figure 5.

The coordinate change caused by tool deflection.

where η_kx^j+1, η_ky^j+1, η_(k−1)x^j+1, and η_(k−1)y^j+1 are the distance of the tool center offset distance, which can be calculated in Section 2.2. The mean square error (MSE) is the difference in the IUCT between consecutive tooth channels. The MSE is measured to evaluate the cutting steady-state and convergence. The iterative discriminant can be expressed as:

$$\begin{array}{c}\multicolumn{1}{c}{E_{\mathrm{MSE}}={[\sum {[h_u\left(n\right)-h_u(n-1)]}^2/N]}^{1/2}}\\ \multicolumn{1}{c}{\le \Delta n\ge 2,\Delta \to 0}\end{array}$$

(25)

where N is the sample size, n is the number of tooth paths, h_u(n) and h_u(n−1) are the IUCT removed between the n and n−1 tooth path at the position angle ψ_k. Once the iteration reaches convergence, the IUCT in the stable-state can be obtained, and the steps of the iterative algorithm are given in Algorithm 1.

Algorithm 1: Flexible force iterative algorithm

Inputs: Cutting and cutter geometry parameters R, Lt, Lm, Lr, Pe, αe, β, fs, Kns, Knf, Knp, Kfp, Kne, Kfe, ro, and γ0. The tool flute is divided into N points. Repeat: The EMSE of the instantaneous cutting thickness between the current tooth and the previous tooth is used to determine whether the iteration is terminated or not. Determine the most stable force.While (stopping condition is not met) do. Calculate the instantaneous cutting thickness h0 and tool center (x0, y0) without considering the effect of deflection, and set the initial mean square error value to 1.Updatehu, xo, and yo.• Fori = 1: N. • UpdateFx, Fy using equation (17). • Updateηx, ηy using equation (20). • Updatexo, yo using equation (3), equation (5), and equation (24). • UpdateEMSE using equation (25). • Return: hu(i), ηx(i), and ηy(i). • Until (Iteration termination).Outputs: Cutting force in steady-state Fx, Fy.

Coefficient identification and fitting method

The effective rake and clearance angles are variable parameters for calculating IUCT-h(ψ_k). Further, the traditional linear fitting and experimental measurement ²⁰ cannot be directly used to obtain parameters. To overcome these shortcomings, the Harris hawks Error Optimization (HEO) analysis algorithm is developed to efficiently identify the cutting force coefficients/parameters.

Harris-hawk optimization (HHO)

The Harris hawks can show various chase ways based on the dynamic environment and the prey escape pattern. ²¹ The main hunting behavior of the Harris hawk can be shown in Figure 6.

Figure 6.

Description of each phase of HHO.

The Harris hawk has two strategies for waiting for finding its prey, which can be modeled as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\mathit{X}(p+1)=}\\ \multicolumn{1}{c}{\{\begin{array}{c}{\mathit{X}}_{\mathrm{rand}}\left(p\right)-r_1|{\mathit{X}}_{\mathrm{rand}}\left(p\right)-2r_2\mathit{X}\left(p\right)|q\ge 0.5\hfill \\ \multicolumn{1}{c}{{\mathit{X}}_{\mathrm{prey}}\left(p\right)-{\mathit{X}}_m\left(p\right)-r_3[L_b+r_4(U_b-L_b)]q<0.5}\end{array}}\end{array}$$

(26)

where X (*) represents the position vectors of hawks in *-iteration, X _prey (p) are prey position, r_1–4 and q are random numbers, U_b and L_b are the upper and lower variables boundaries, X _rand (p) is any hawk from the population, and X_m (p) is the average population position:

$${\mathit{X}}_m\left(p\right)=\sum _{i=1}^N{X}_i\left(p\right)/N_p$$

(27)

where X_i (p) is the position of each hawk in p iterations, N_p is the hawks number. Based on the relationship between prey escape and its energy, the energy of the prey E is simulated as:

$$E=2E_0(P-p)/P$$

(28)

Where P is the maximum iterate, and E₀ is the initial energy of prey (−1, 1). The Harris hawks attacking strategies process is as follows:

• (1) |E| ≥ 0 and r ≥ 0.5 (r is the probability of the prey escaping)

The prey tries to escape, and the hawk positions can be expressed as:

$$\mathit{X}(p+1)=[{\mathit{X}}_{\mathrm{prey}}\left(p\right)-\mathit{X}\left(p\right)]-E|\mathit{J}{\mathit{X}}_{\mathrm{prey}}\left(p\right)-\mathit{X}\left(p\right)|$$

(29)

where J = 2(1−r₅) represents the prey movement pace during the escape process, and r₅∈(0,1).

• (2) |E| ≥ 0.5 and r < 0.5

The hawks rapidly dive around the prey several times and determine their position and direction according to the movement of the prey. The position of the hawks can be simulated as follows:

$$X(p+1)=\{\begin{array}{c}Y\mathrm{if}F\left(Y\right)<F\left(X\left(p\right)\right)\hfill \\ \multicolumn{1}{c}{Z\mathrm{if}F\left(Z\right)<F\left(X\left(p\right)\right)}\end{array}$$

(30)

$$\mathit{Y}={\mathit{X}}_{\mathrm{rabbit}}\left(p\right)-E|\mathit{J}{\mathit{X}}_{\mathrm{rabbit}}\left(p\right)-\mathit{X}\left(p\right)|$$

(31)

$$\mathit{Z}=\mathit{Y}+{\mathit{S}}_{\mathit{p}}\times \mathrm{Levy}\left(D_p\right)$$

(32)

where F(*) is the fitness function, D_p is the search agent space, S_p is a random vector, and Levy can be expressed as:

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{Levy}\left(x\right)=0.01u\sqrt[-\kappa ]{v}}\\ \multicolumn{1}{c}{{(\Gamma (1+\kappa )\sin (0.5\pi \kappa )/\Gamma (0.5+0.5\kappa )2^{0.5(\kappa -1)}\kappa )}^{1/\kappa }}\end{array}$$

(33)

where u, v∈(0,1), and κ = 1.5.

• (3) |E| < 0.5 and r >= 0.5

The prey is unable to escape, and the hawk position is updated as follows:

$$\mathit{X}(p+1)={\mathit{X}}_{\mathrm{rabbit}}\left(p\right)-E\left|\Delta \mathit{X}\left(p\right)\right|$$

(34)

• (4) |E| < 0.5 and r < 0.5

Hawks attempt to reduce the distance from the average position of the fleeing prey according to the following rules:

$$X(p+1)=\{\begin{array}{c}\mathrm{W}\mathrm{if}F\left(W\right)<F\left(X\left(p\right)\right)\hfill \\ \multicolumn{1}{c}{V\mathrm{if}F\left(V\right)<F\left(X\left(p\right)\right)}\end{array}$$

(35)

$$\mathit{W}={\mathit{X}}_{\mathrm{prey}}\left(p\right)-E|\mathit{J}{\mathit{X}}_{\mathrm{prey}}\left(p\right)-{\mathit{X}}_m\left(p\right)|$$

(36)

$$\mathit{V}=\mathit{W}+\mathit{S}\times \mathrm{Levy}\left(D\right)$$

(37)

The HHO algorithm removes the derivative part of the limit state function to solve high-dimensional calculation problems more efficiently.

Parameter identification

The HEO model is established by combining the experimental data and the HHO method:

$$\{\begin{array}{c}\mathrm{Min}:{\left[\begin{array}{c}{f}_1\left(\mathit{x}\right)\hfill \\ \multicolumn{1}{c}{f_2\left(\mathit{x}\right)}\\ \multicolumn{1}{c}{f_3\left(\mathit{x}\right)}\end{array}\right]}^T=[\sum _{i=1}^n\sum _{j=1}^m{({\mathit{F}}_{\exp }^{\mathrm{avg}}-{\mathit{F}}_{\mathrm{sim}}^{\mathrm{avg}})}^2,\sum _{i=1}^n\sum _{j=1}^m{[({\mathit{F}}_{\exp }^{\max }-{\mathit{F}}_{\exp }^{\mathrm{avg}})-({\mathit{F}}_{\mathrm{sim}}^{\max }-{\mathit{F}}_{\mathrm{sim}}^{\mathrm{avg}})]}^2,\hfill \\ \multicolumn{1}{c}{\sum _{i=1}^n\sum _{j=1}^m{[({\mathit{F}}_{\exp }^{\min }-{\mathit{F}}_{\exp }^{\mathrm{avg}})-({\mathit{F}}_{\mathrm{sim}}^{\min }-{\mathit{F}}_{\mathrm{sim}}^{\mathrm{avg}})]}^2]}\\ \multicolumn{1}{c}{\mathrm{Subject}\mathrm{to}:[g_1\left(\mathit{x}\right),g_2\left(\mathit{x}\right)]=[K_{\mathrm{ns}}{K}_{\mathrm{fs}}{K}_{\mathrm{np}}{K}_{\mathrm{fp}}{K}_{\mathrm{ne}}{K}_{\mathrm{fe}}{r}_o\ge 0,\left|{\gamma }_0\right|\le 2\pi ]}\end{array}$$

(38)

where x  = [K_ns, K_fs, K_np, K_ne, K_fe, r_o, γ_o], m and n are the numbers of selected f_s and experiments, f_*( x ) and g_*( x ) are the objective and constraint function, respectively, F _sim and F _exp are the prediction results under the proposed model with the iteration of the optimization error algorithm and the experimental force data in x and y directions, ${\mathit{F}}_{\exp }^{\mathrm{avg}},{\mathit{F}}_{\exp }^{\max },{\mathit{F}}_{\exp }^{\min },{\mathit{F}}_{\mathrm{sim}}^{\mathrm{avg}},{\mathit{F}}_{\mathrm{sim}}^{\max }$ , and ${\mathit{F}}_{\mathrm{sim}}^{\min }$ are the mean, maximum, and minimum experimental and simulation forces, respectively.

The decrease in the objective function is advantageous to find the solution vector conforming to constraints rapidly. ²² This process can be modeled by introducing a penalty function:

$$\begin{array}{c}\multicolumn{1}{c}{P\left(\mathit{x}\right)=f\left(\mathit{x}\right)+\sum _{i=1}^I{o}_i{G}_i\left(\mathit{x}\right)+\sum _{j=1}^J{s}_j{H}_j\left(\mathit{x}\right),}\\ \multicolumn{1}{c}{(G_i\left(\mathit{x}\right)=max{[0,g_i\left(\mathit{x}\right)]}^a,H_j\left(\mathit{x}\right)={\left|h_j\left(\mathit{x}\right)\right|}^{\beta })}\end{array}$$

(39)

where o_i and s_j are penalty function coefficients, g_i ( x ) and h_j ( x ) are the inequality/quality constraint, α and β are coefficients. The greater the degree of violation of constraints, the greater the corresponding G( x ) and H( x ), and the greater the value of the penalty function. ²³ The transformed objective function is expressed as:

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{Min}:P{\left(\mathit{x}\right)}_{\mathrm{new}}=u{f}_1\left(\mathit{x}\right)+v{f}_2\left(\mathit{x}\right)+w{f}_3\left(\mathit{x}\right)}\\ \multicolumn{1}{c}{+o_1{g}_1\left(\mathit{x}\right)+o_2{g}_2\left(\mathit{x}\right),(u+v+w=1)}\end{array}$$

(40)

Using the HEO algorithm to calculate the new error function value, the iterative input value is changed according to the feedback value. Once the iterative process is stopped, the parameters can be obtained with the calculating flow shown in Figure 7.

Figure 7.

Flow chart of micro-milling parameter identification.

Adaptive Kriging fitting

Adaptive Kriging (AK) is used to fit the established implicit model to save calculation time and facilitate the analysis of extensive sample data. In addition, Monte Carlo Simulation (MCS) ²⁴ method was used to verify the AK-results.

Initial Kriging

The Kriging method is an improved technology with linear regression analysis, including linear regression and nonparametric parts. ²⁵ The nonparametric part is the realization of the stochastic distribution process, and the Kriging model can be written as:

$$Y\left(\mathit{x}\right)=\sum _{j=1}^{k_0}{\mathit{q}}_{\mathit{x}}\mathbf{\beta }+z\left(\mathit{x}\right)$$

(41)

where q_x is the basis function of x , β is the coefficient, z( x ) is a random distribution, and its covariance matrix can be given as:

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{Cov}[z\left({\mathit{x}}^j\right),z\left({\mathit{x}}^h\right)]={\sigma }^2\underset{l=1}{\overset{m}{\Pi }}{R}_l({\theta }_l,{{\mathit{x}}_l}^j,{{\mathit{x}}_l}^h)}\\ \multicolumn{1}{c}{j,h=1,2,\dots ,m}\end{array}$$

(42)

where θ_l (l=1, 2, …, m) is the correlation parameter, R_l is the correlation function, which is expressed as follows:

$$R(\mathbf{\theta },{\mathit{x}}^j,{\mathit{x}}^h)=\exp (-\sum _{l=1}^m{\theta }_l{|{{\mathit{x}}_l}^j,{{\mathit{x}}_l}^h|}^2)$$

(43)

Based on the generalized least squares method, the estimated values of β and σ ² are:

$$\{\begin{array}{l}{\beta }^{\prime }={\left(Q^T{R}^{-1}Q\right)}^{-1}{Q}^T{Y}^{-1}G\\ {\sigma \prime }^2={\left(Y-Q^T{\beta }^{\prime }\right)}^T{R}^{-1}\left(Y-Q^T{\beta }^{\prime }\right)/m\end{array}$$

(44)

where S is a vector of the response values, and Q _m×k is the matrix composed of regression models. According to the maximum likelihood estimation method, the parameter θ can be solved as follows:

$$\theta =\arg \underset{\theta \ge 0}{\min }\left[-m\ln \left({\sigma \prime }^2\right)+\ln \left|R\right|\right]/2$$

(45)

Adaptive-Kriging (AK) method

AK method improves the prediction accuracy by analyzing the fault state of the output parameters obtained by the initial Kriging. ²⁶ The model is updated circularly until the adaptive convergence condition is satisfied. The learning process can be outlined as follows:

$$\mathrm{U}\left(\mathit{x}\right)=|{\mu }_Y\left(\mathit{x}\right)/{{\sigma }_Y}^2\left(\mathit{x}\right)|$$

(46)

The Adaptive Kriging update process is as follows:

1. A sample pool S_p composed of N samples is generated in the input variable sample space.

2. Samples of N₁ input variables from the sample pool S_p are substituted into the actual function to calculate their corresponding function values. The initial training set T is formed, and Kriging model G( X ) is established.

3. The learning function U(x) corresponding to each sample in the sample pool S_p is calculated and the calculation process of sample points is updated:

   $${\mathit{x}}^{\prime }=\arg \underset{\mathit{x}\in {\mathit{S}}_{\mathit{p}}}{min}U\left(\mathit{x}\right)$$

   (47)

4. When $\underset{\mathit{x}\in \mathit{S}}{min}\mathrm{U}\left(\mathit{x}\right)$  ≥ 2, the adaptive learning process is stopped. ²⁷ If $\underset{\mathit{x}\in \mathit{S}}{min}\mathrm{U}\left(\mathit{x}\right)$  < 2, [ x ′, Y( x ′)] is added to the current T. The algorithm returns to step (3) and updates the Kriging model.

The complex flexible mechanical model established in Section 2 is fitted by the AK method. Then, many input samples are substituted into the new fitting model for output, providing the basis for the next section – the rapid micro-milling force probability analysis (Figure 8).

Figure 8.

Flow chart of AK method.

Experimental and probability analysis

The cutting experiment is carried out on the experimental micro-milling machining center with a maximum electric spindle speed of 24,000 r/min. The workpiece material is Al6061 (90 × 10 × 120 mm), with an elastic recovery rate of 0.1. The micro-milling cutter is a tungsten carbon uncoated two-edged with a diameter of 1 mm. Table 1 outlines the parameters used in the calculation simulation and experimental processes.

Table 1.

Relevant parameters of the micro-milling experiment (Ω = 18,000 r/min, b = 0.1 mm).

Parameters	Value
Edge radius (re)	2 μm
Shank diameter (D)	4 mm
Tool length (L)	50 mm
Neck length (Lr)	10 mm
Flute length (Lt)	4 mm
The inclination angle of tool neck (τ)	10°
Helix angle (β)	45°
Rake angle (α)	5°
Clearance angle (ξ)	10°
Minimum chip thickness (hmin)	0.3re
Feed per tooth (fs)	3–10 μm/tooth

As shown in Figure 9, the cutting path is set in the computer control center. The machining coordinate of the machining system is found through the operating handle, and then the cutting command is executed. The force sensor, data collector, and charge amplifier record the force of the micro-milling tool. Finally, the force curve is visualized on the computer, and the data values are obtained.

Figure 9.

Experimental setup.

Model validation

The parameters obtained using the proposed parameter identification method are listed in Table 2. To verify the output accuracy, the average values of 50 groups of experimental milling force data were obtained and plotted in Figure 10(a).

Table 2.

Cutting coefficient and runout parameters with the minimum error.

Kns	Kfs	Knp	Kfp	Kne	Kfe	ro (mm)	γ 0 (°)
2926	2238	3203	2496	1597	1004	0.0012	50

Figure 10.

Comparison diagram of cutting force curve: (a) spindle speed 18,000 r/min, cutting depth 100 μm, (b) f_s = 4 μm/tooth, and (c) f_s = 8 μm/tooth.

As shown in Figure 10(a), the maximum, minimum, and average cutting force values obtained by the flexible milling force model are within the extreme ranges of experimentally-obtained values, with small error values. To further explore the accuracy of the milling force obtained by the established AK method, the cutting forces at the different positions during the milling process are compared and analyzed

The Ub-AK and Lb-AK in Figure 10(b) and (c) represent the extreme-confidence interval range of the superposition effect of the cutting force considering the 95% confidence interval obtained by using the central limit theorem and the big tree law and the error value of AK fitting. Mean-AK is the average value of AK fitting cutting force. Sim and Exp represent the force predicted by the flexible force model and measured experimentally. The results of the comparison are as follows:

1. The prediction curve agrees with the experimental results, verifying the accuracy of the proposed HEO parameter identification method.

2. The Mean-AK coincides with the predicted curve of the flexible force model (Sim), indicating that the AK fitting method applies to micro-milling force calculation.

3. The majority of the experimental data are found in the extreme-confidence interval. A small number of experimental data points are around the boundary, conforming to the nonideal working conditions inherent to the realistic machining process.

Probability characteristic prediction

Following the modeling and analyzing the micro-milling, the probability distribution of the milling force is discussed. In this section, the AK fitting method and the direct MCS method combined with the proposed flexible force model are used to conduct the comparative analysis and determine the accuracy of the result.

The cutting force at the position angle ψ_k = 80° was analyzed (100 sample points) to further discuss the differences between the AK and the established models. Based on the characteristic analysis of general machining parameters, ²⁸ the specific distribution of micro-milling parameters can be determined (Table 3). The milling forces are compared, and the corresponding scatters plot obtained using both methods are shown in Figure 11. The curve distribution in the figures shows that the output data obtained by the AK method is practically identical to that obtained by combining the MCS with the proposed flexible force model, which verifies the accuracy of the AK method.

Table 3.

Characteristics of random variables.

Random variable	Distribution type	Mean value	Coefficient of variation (Cov = Stand/Mean)
Cutting depth b	Normal	100 μm	0.05
Spindle speed Ω		18,000 r/min
Helix angle β		45°
Shank diameter D		4 mm
Flute diameter d		1 mm
Feed per tooth fs		4 μm/tooth
Runout length ro		1.2 μm/tooth
Runout angle γ0		50°
Total length L		50 mm
Flute length Lt		4 mm
Tool rake angle α		5°
Tool clearance angle ξ		10°

Figure 11.

Comparison of scatter plots (ψ_k = 80°).

Based on the distribution of random parameters in Table 3, each random variable is sampled by the MCS (10⁵) and substituted into the flexible force model and the trained AK model for prediction. Figure 12 shows the probability distribution characteristics of cutting force when position angle ψ_k = 80°, indicating that the milling force frequencies obtained under the two methods are basically the same. The cutting force distribution is concentrated near the average value with the greater probability value.

Figure 12.

Analysis diagram of cutting force characteristics.

The statistical calculation provided in Figures 10 to 12 shows that the predicted mean value of the cutting force is close to the actual value. Its Cov is greater than the Cov of a single parameter, indicating that the calculation result of the cutting force is unbiased and has strong dispersion. Moreover, the analysis time used by the AK method is 10 min. In comparison, the calculation time used by the flexible force model directly is 240 min, which proves that the analysis of probability problems AK is more efficient.

Conclusion

A new method for identifying the micro-milling parameters using the HEO method is proposed. Based on the identified parameters and cutting mechanism, an accurate flexible cutting force model for micro-milling is established. Combined with the proposed model and AK method, efficient calculation and analysis of cutting force probability characteristics in micro-milling are realized. The main contents of this study can be summarized as follows:

• A new micro-milling flexible force model framework is present. The model considers the influences of multiple factors, including the tool runout, IUCT, MCT, material elastic recovery rate, tool deformation and cutting force coupling, and chip separation state. After comparing experimental data and calculation results, it is observed that the proposed mechanical model has a high prediction accuracy.

• The HEO parameter identification method for solving the cutting coefficient/parameter of micro-milling is developed. This method integrates the experimental data and mechanical model into the objective optimization function. It directly obtains the coefficients/parameters by using the optimization operation, which overcomes the complexity of repeated experimental measurements.

• The AK method is used to reconstruct the functional relationship in the mechanical model, avoiding the time consumption stemming from repeatedly using the complex flexible force model directly. The probability characteristics of micro-milling force are analyzed via AK and direct MCS, respectively.

The curve shown in Figure 10 indicates that the experimental data are distributed within the extreme-confidence interval of cutting force obtained by the AK method, which illustrates that the proposed flexible force model is accurate. The curves in Figures 10 to 12 show that the average values of cutting force simulation is an unbiased estimate of cutting force, and the cutting force has a large dispersion (Cov >0.05), so it is necessary to consider the randomness of parameters when analyzing the cutting force state.

Compared to the traditional mechanical model, this study considers the influence of multiple factors and the uncertainty of processing parameters on the micro-milling force modeling. As such, it is more in line with the actual engineering situation. The actual distribution of parameters can be used as a research issue in the forthcoming studies. The above research results have positive guiding value for analyzing the precise micro-milling mechanism, improving machining quality, and enhancing manufacturing efficiency.