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.

Keywords

Micro-milling parameter identification Harris hawks Error Optimization Adaptive-Kriging probability characteristics

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_{(k)}}^{2} = {(R_{0} \cos (γ_{o} + 2 π (k - 1) / K) + r_{o})}^{2} + {(R_{0} \sin γ_{o})}^{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 (ω t_{k} + γ_{o}), r_{o} \cos (ω t_{k} + γ_{o}))

(2)

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

{\begin{matrix} x_{k} = f t_{k} / 60 + R \sin (ω t_{k} - 2 π k / K) + r_{o} \sin (ω t_{k} + γ_{o}) \\ y_{k} = R \cos (ω t_{k} - 2 π k / K) + r_{o} \cos (ω t_{k} + γ_{o}) \end{matrix}

(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{matrix} Y = \frac{y_{k} - y_{O_{k}}}{x_{k} - x_{O_{k}}} X = \frac{X - f t_{k} / 60 - r_{o} \sin (ω t_{k} + γ_{o})}{\tan (ω t_{k} - 2 π k / K)} \\ + R \cos (ω t_{k} + γ_{o}) \end{matrix}

(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{matrix} x_{k - 1} = f t_{k - 1} / 60 + R \sin (ω t_{k - 1} - 2 π (k - 1) / K) \\ + r_{o} \sin (ω t_{k - 1} + γ_{o}) \\ y_{k - 1} = R \cos (ω t_{k - 1} - 2 π (k - 1) / K) \\ + r_{o} \cos (ω t_{k - 1} + γ_{o}) \end{matrix}

(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{matrix} G (t') = R \tan (ω t_{k} - 2 π k / K) \cos (ω t_{k - 1} - 2 π (k - 1) / K) \\ + f (t_{k} - t_{k - 1}) / 60 + r_{o} \tan (ω t_{k} - 2 π k / K) \\ [\cos (ω t_{k - 1} + γ_{o}) - \cos (ω t_{k} + γ_{o})] \\ - R \sin (ω t_{k - 1} - 2 π (k - 1) / K) \\ - r_{o} \sin (ω t_{k - 1} + γ_{o}) + r_{o} \sin (ω t_{k} + γ_{o}) = 0 \end{matrix}

(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}' = t_{i}' - G (t_{i}') / G (t_{i}')' = t - 2 π / ω K - G (t_{i}') / G (t_{i}')'

(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{matrix} h_{u} = H - U \\ = (R_{(k)} + \sqrt{{(x_{k} - x_{k - 1})}^{2} + {(y_{k} - y_{k - 1})}^{2}} \sin (ω t_{k} - 2 π k / K + α)) \\ - \sqrt{{R_{(k)}}^{2} - L^{2} \cos {(ω t_{k} - 2 π k / K + \arctan [(y_{O_{k}} - {y_{O}}_{k - 1}) / (x_{O_{k}} - x_{O_{k - 1}})])}^{2}} \end{matrix}

(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} (k) = {\begin{matrix} h_{u} (k) + P_{e} \times h_{u} (k - 1) h_{u} (k) \leq h_{min} = λ r_{e} \\ h_{u} (k) h_{u} (k) > h_{min} = λ r_{e} \end{matrix}

(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{matrix} d F_{ns} = K_{ns} h (ψ_{k}) dz \cos (α_{e}) + K_{fs} h (ψ_{k}) \sin (α_{e}) \\ d F_{fs} = - K_{ns} h (ψ_{k}) dz \sin (α_{e}) + K_{fs} h (ψ_{k}) \cos (α_{e}) \end{matrix}

(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:

α_{e} = {\begin{matrix} α \frac{h (ψ_{k}) + r_{e} [(\frac{α}{2} - \frac{π^{2}}{8 α} - \frac{(1 / \cos (1 - h_{min} / r_{e}) + π)}{2 α \cos (1 - h_{min} / r_{e})}) \cos (α) - \sin (α) - 1]}{h (ψ_{k}) + r_{e} [(α + π / 2 - 1 / \cos (1 - h_{min} / r_{e})) \cos (α) - \sin (α) - 1]} h (ψ_{k}) \geq r_{e} (1 + \sin (α)) \\ (1 / \cos (1 - h (ψ_{k}) / r_{e}) + 1 / \cos (1 - h_{min} / r_{e}) - π) / 2 h (ψ_{k}) \geq r_{e} (1 + \sin (α)) \end{matrix}

(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{matrix} d F_{np} = K_{np} h_{min} dz \sin (ζ) + K_{fp} h_{min} dz \cos (ζ) \\ d F_{fp} = K_{np} h_{min} dz \cos (ζ) - K_{fp} h_{min} dz \sin (ζ) \\ ζ = 0.5 / \cos ((r_{e} - h_{min}) / r_{e}) \end{matrix}

(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{cases} d F_{n p}' = K_{n p} (h (ψ_{k}) + h_{a c}) d z \cos (α_{p}) + K_{f p} (h (ψ_{k}) + h_{a c}) d z \sin (α_{p}) \\ d F_{f p}' = K_{n p} (h (ψ_{k}) + h_{a c}) d z \sin (α_{p}) - K_{f p} (h (ψ_{k}) + h_{a c}) d z \cos (α_{p}) \\ α_{p} = 1 / \cos (\sqrt{(h (ψ_{k}) + h_{b f}) / (2 r_{e})}) \end{cases}

(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{matrix} h_{ac} = (0.0101806 E_{y} \sqrt{h_{min} / (2 r_{e} - h_{min})} / σ_{0} - 0.040468) \\ h (ψ_{k}) \end{matrix}

(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{matrix} d F_{ne} = - K_{ne} P_{e} h (ψ_{k}) dz \sin (ζ_{e}) + K_{ne} P_{e} h (ψ_{k}) dz \cos (ζ_{e}) \\ d F_{fe} = K_{ne} P_{e} h (ψ_{k}) dz \cos (ζ_{e}) + K_{ne} P_{e} h (ψ_{k}) dz \sin (ζ_{e}) \end{matrix}

(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:

ζ_{e} = {\begin{matrix} \arctan (\sqrt{h_{er} / (2 r_{e} - h_{er})}) h_{er} \leq r_{e} (1 - \cos (γ)) \\ \arctan [h_{er} \sin (γ) / (h_{er} \cos \\ (γ) + r_{e} (1 - \cos (γ)))] h_{er} > r_{e} (1 - \cos (γ)) \end{matrix}

(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{matrix} \sum_{i = 1}^{M} [\begin{matrix} d F_{x} \\ d F_{y} \end{matrix}] = [\begin{matrix} - \cos (ψ_{k}) - \sin (ψ_{k}) \\ - \sin (ψ_{k}) \cos (ψ_{k}) \end{matrix}] \sum_{i = 1}^{M} \\ [\begin{matrix} d F_{ns} + d F_{np} + d {F_{np}}^{'} + d F_{ne} \\ d F_{fs} + d F_{fp} + d {F_{fp}}^{'} + d F_{fe} \end{matrix}] \end{matrix}

(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{matrix} ρ^{- 1} = M (z) / EI = M (z) / (E π r^{4} / 4) = (d^{2} η (z) / d z^{2}) \\ / {(1 + {(d η (z) / d z)}^{2})}^{3 / 2} \end{matrix}

(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:

ρ^{- 1} = M (z) / EI = d^{2} η (z) / d z^{2} = η^{''}

(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{matrix} η_{x (y)} = η_{1} + {η_{1}}^{'} (L_{m} + L_{t} - z) + η_{2} + {η_{2}}^{'} (L_{t} - z) + η_{3} \\ η_{x (y)}^{'} = {η_{1}}^{'} + {η_{2}}^{'} + {η_{3}}^{'} \end{matrix}

(20)

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

{\begin{matrix} η_{1} = F {(E π {R_{r}}^{4})}^{- 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 \leq L \\ η_{2} = 4 F {(E π)}^{- 1} {\begin{matrix} [R_{r} {[R_{r} - (Z - L_{t}) \tan (τ)]}^{- 2} / 6 - {[R_{r} - (Z - L_{t}) \tan (τ)]}^{- 1} / 2] / \tan (τ) \\ + (L_{t} - z) {[R_{r} - (Z - L_{t}) \tan (τ)]}^{- 2} / 6 + B_{1} [R_{r} - (Z - L_{t}) \tan (τ)] + B_{2} \end{matrix}} L_{t} < Z \leq L_{m} + L_{t} \\ η_{3} = 4 F {(E π {R_{t}}^{4})}^{- 1} {(< Z - z >}^{3} / 6 + C_{1} Z + C_{2}) 0 \leq Z \leq L_{t} \end{matrix}

(21)

{\begin{matrix} {η_{1}}^{'} = 4 F {(π {R_{r}}^{4})}^{- 1} [(L_{m} + L_{t} - z) Z + F {(Z - L_{m} - L_{t})}^{2} / 2 + A_{1}] L_{m} + L_{t} < Z \leq L \\ {η_{2}}^{'} = 4 F {(E π)}^{- 1} {\begin{matrix} - [- R_{r} {[R_{r} - (Z - L_{t}) \tan (τ)]}^{- 3} / 3 + {[R_{r} - (Z - L_{t}) \tan (τ)]}^{- 2} / 2] / \tan (τ) \\ + (L_{t} - z) {[R_{r} - (Z - L_{t}) \tan (τ)]}^{- 3} / 3 - B_{1} \end{matrix}} L_{t} < Z \leq L_{m} + L_{t} \\ {η_{3}}^{'} = 4 F {(E π {r_{t}}^{4})}^{- 1} ({< Z - z >}^{2} / 2 + C_{1}) 0 \leq Z \leq L_{t} \end{matrix}

(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{matrix} z_{x / y} (ψ_{k}) = L - {\sum_{i = 1}^{N} \sum_{j = 1}^{S} d F_{x / y} (ψ_{k}) [L - z_{x / y} (ψ_{k})]} / \\ F_{x / y} (ψ_{k}) \end{matrix}

(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−₁^j+1, y_k−₁^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{matrix} ({x_{ok}}^{j + 1}, {y_{ok}}^{j + 1}) = (x + {η_{kx}}^{j + 1}, y + {η_{ky}}^{j + 1}) \\ ({x_{ok - 1}}^{j + 1}, {y_{ok - 1}}^{j + 1}) \\ = (x^{'} + {η_{(k - 1) x}}^{j + 1}, y^{'} + {η_{{(k - 1)}_{y}}}^{j + 1}) \end{matrix}

(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{matrix} E_{MSE} = {[\sum {[h_{u} (n) - h_{u} (n - 1)]}^{2} / N]}^{1 / 2} \\ \leq Δ n \geq 2, Δ \to 0 \end{matrix}

(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, L_t, L_m, L_r, P_e, α_e, β, f_s, K_ns, K_nf, K_np, K_fp, K_ne, K_fe, r_o, and γ₀. The tool flute is divided into N points.
Repeat:
The E_MSE 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 h₀ and tool center (x₀, y₀) without considering the effect of deflection, and set the initial mean square error value to 1.Updateh_u, x_o, and y_o.• Fori = 1: N.
• UpdateF_x, F_y using equation (17).
• Updateη_x, η_y using equation (20).
• Updatex_o, y_o using equation (3), equation (5), and equation (24).
• UpdateE_MSE using equation (25).
• Return: h_u(i), η_x(i), and η_y(i).
• Until (Iteration termination).Outputs:
Cutting force in steady-state F_x, F_y.

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{matrix} X (p + 1) = \\ {\begin{matrix} X_{rand} (p) - r_{1} | X_{rand} (p) - 2 r_{2} X (p) | q \geq 0.5 \\ X_{prey} (p) - X_{m} (p) - r_{3} [L_{b} + r_{4} (U_{b} - L_{b})] q < 0.5 \end{matrix} \end{matrix}

(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:

X_{m} (p) = \sum_{i = 1}^{N} X_{i} (p) / 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 = 2 E_{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:

X (p + 1) = [X_{prey} (p) - X (p)] - E | J X_{prey} (p) - X (p) |

(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{matrix} Y if F (Y) < F (X (p)) \\ Z if F (Z) < F (X (p)) \end{matrix}

(30)

Y = X_{rabbit} (p) - E | J X_{rabbit} (p) - X (p) |

(31)

Z = Y + S_{p} \times Levy (D_{p})

(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{matrix} Levy (x) = 0.01 u \sqrt[- κ]{v} \\ {(Γ (1 + κ) \sin (0.5 π κ) / Γ (0.5 + 0.5 κ) 2^{0.5 (κ - 1)} κ)}^{1 / κ} \end{matrix}

(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:

X (p + 1) = X_{rabbit} (p) - E | Δ X (p) |

(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{matrix} W if F (W) < F (X (p)) \\ V if F (V) < F (X (p)) \end{matrix}

(35)

W = X_{prey} (p) - E | J X_{prey} (p) - X_{m} (p) |

(36)

V = W + S \times Levy (D)

(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{matrix} Min : {[\begin{matrix} f_{1} (x) \\ f_{2} (x) \\ f_{3} (x) \end{matrix}]}^{T} = [\sum_{i = 1}^{n} \sum_{j = 1}^{m} {(F_{\exp}^{avg} - F_{sim}^{avg})}^{2}, \sum_{i = 1}^{n} \sum_{j = 1}^{m} {[(F_{\exp}^{\max} - F_{\exp}^{avg}) - (F_{sim}^{\max} - F_{sim}^{avg})]}^{2}, \\ \sum_{i = 1}^{n} \sum_{j = 1}^{m} {[(F_{\exp}^{\min} - F_{\exp}^{avg}) - (F_{sim}^{\min} - F_{sim}^{avg})]}^{2}] \\ Subject to : [g_{1} (x), g_{2} (x)] = [K_{ns} K_{fs} K_{np} K_{fp} K_{ne} K_{fe} r_{o} \geq 0, | γ_{0} | \leq 2 π] \end{matrix}

(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, $F_{\exp}^{avg}, F_{\exp}^{\max}, F_{\exp}^{\min}, F_{sim}^{avg}, F_{sim}^{\max}$ , and $F_{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{matrix} P (x) = f (x) + \sum_{i = 1}^{I} o_{i} G_{i} (x) + \sum_{j = 1}^{J} s_{j} H_{j} (x), \\ (G_{i} (x) = max {[0, g_{i} (x)]}^{a}, H_{j} (x) = {| h_{j} (x) |}^{β}) \end{matrix}

(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{matrix} Min : P {(x)}_{new} = u f_{1} (x) + v f_{2} (x) + w f_{3} (x) \\ + o_{1} g_{1} (x) + o_{2} g_{2} (x), (u + v + w = 1) \end{matrix}

(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 (x) = \sum_{j = 1}^{k_{0}} q_{x} β + z (x)

(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{matrix} Cov [z (x^{j}), z (x^{h})] = σ^{2} Π_{l = 1}^{m} R_{l} (θ_{l}, {x_{l}}^{j}, {x_{l}}^{h}) \\ j, h = 1, 2, \dots, m \end{matrix}

(42)

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

R (θ, x^{j}, x^{h}) = \exp (- \sum_{l = 1}^{m} θ_{l} {| {x_{l}}^{j}, {x_{l}}^{h} |}^{2})

(43)

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

{\begin{cases} β^{'} = {(Q^{T} R^{- 1} Q)}^{- 1} Q^{T} Y^{- 1} G \\ {σ'}^{2} = {(Y - Q^{T} β^{'})}^{T} R^{- 1} (Y - Q^{T} β^{'}) / m \end{cases}

(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:

θ = \arg \min_{θ \geq 0} [- m \ln ({σ'}^{2}) + \ln | R |] / 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:

U (x) = | μ_{Y} (x) / {σ_{Y}}^{2} (x) |

(46)

The Adaptive Kriging update process is as follows:

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

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.

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:

x^{'} = \arg min_{x \in S_{p}} U (x)

(47)

When $min_{x \in S} U (x)$ ≥ 2, the adaptive learning process is stopped.²⁷ If $min_{x \in S} U (x)$ < 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 (r_e)	2 μm
Shank diameter (D)	4 mm
Tool length (L)	50 mm
Neck length (L_r)	10 mm
Flute length (L_t)	4 mm
The inclination angle of tool neck (τ)	10°
Helix angle (β)	45°
Rake angle (α)	5°
Clearance angle (ξ)	10°
Minimum chip thickness (h_min)	0.3r_e
Feed per tooth (f_s)	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.

K_ns	K_fs	K_np	K_fp	K_ne	K_fe	r_o (mm)	γ ₀ (°)
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:

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

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.

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 f_s		4 μm/tooth
Runout length r_o		1.2 μm/tooth
Runout angle γ₀		50°
Total length L		50 mm
Flute length L_t		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.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors gratefully acknowledge the support of the National Natural Science Foundation of China (51975110), Liaoning Revitalization Talents Program (XLYC1907171), and Fundamental Research Funds for the Central Universities (N2003005, N2203004), Key-Area Research and Development Program of Guangdong Province (2020B090928001).

ORCID iDs

Xianzhen Huang

YuXiong Li

References

Bhople

Mastud

Mittal

Metallurgical and machining performance aspects of cryotreated tungsten carbide micro-end mill cutters. Proc IMechE, Part B: J Engineering Manufacture. Epub ahead of print 22 June 2022. DOI: 10.1177/09544054221101765.

Zhang

Wang

Cutting forces modeling for micro flat end milling by considering tool runout and bottom edge cutting effect. Proc IMechE, Part B: J Engineering Manufacture 2019; 233(2): 470–485.

Cai

Yao

, et al. Identifying the transient milling force coefficient of a slender end-milling cutter with vibrations. J Manuf Process 2021; 67(2): 262–274.

Ding

Huang

Zhang

, et al. Reliability updating and parameter inversion of micro-milling. Mech Syst Signal Process 2022; 174: 109105.

Wang

Zhang

Huang

, et al. Cutting force modeling and experimental validation for micro end milling. Int J Adv Manuf Technol 2021; 117(3–4): 933–947.

Guo

Jiang

Yang

, et al. An accurate instantaneous uncut chip thickness model combining runout effect in micro-milling using cutters with the non-uniform helix and pitch angles. Proc IMechE, Part B: J Engineering Manufacture 2020; 234(5): 859–870.

Niu

Cheng

Improved dynamic cutting force modelling in micro milling of metal matrix composites. Part I: theoretical model and simulations. Proc IMechE, Part C: J Engineering Manufacture 2020; 234(9): 1733–1745.

Liu

, et al. Micro milling of copper thin wall structure. Int J Adv Manuf Technol 2017; 90(1–4): 405–412.

Jing

Song

, et al. A novel runout model based on spatial tool position for micro-milling force prediction. J Manuf Process 2021; 68: 739–749.

10.

Zhang

Zhu

Generic instantaneous force modeling and comprehensive real engagement identification in micro-milling. Int J Mech Sci 2020; 176(15): 105504.

11.

Gonzalo

Beristain

Jauregi

, et al. A method for the identification of the specific force coefficients for mechanistic milling simulation. Int J Mach Tools Manuf 2010; 50(9): 765–774.

12.

Jin

Altintas

Prediction of micro-milling forces with finite element method. J Mater Process Technol 2012; 212(3): 542–552.

13.

Tsai

Chang

Hung

, et al. Investigation of milling cutting forces and cutting coefficient for aluminum 6060-T6. Comput Electr Eng 2016; 51: 320–330.

14.

Jing

Tian

Yuan

, et al. A runout measuring method using modeling and simulation cutting force in micro end-milling. Int J Adv Manuf Technol 2017; 91(9–12): 4191–4201.

15.

Zhu

Micro-milling force modeling with tool wear and runout effect by spatial analytic geometry. Int J Adv Manuf Technol 2020; 107(1–2): 631–643.

16.

Chen

Ming

Cao

. Reliability-based NC milling parameters optimization using ensemble metamodel. Int J Adv Manuf Tech 2018; 97(9): 3359–3369.

17.

Yuan

Jing

Ehmann

, et al. Modeling of cutting forces in micro end-milling. J Manuf Process 2018; 31: 844–858.

18.

Aslantas

Alatrushi

Bedir

, et al. An experimental analysis of minimum chip thickness in micro-milling of two different titanium alloys. Proc IMechE, Part B: J Engineering Manufacture 2020; 234(12): 1486–1498.

19.

Jun

MBG

Goo

Malekian

, et al. A new mechanistic approach for micro end milling force modeling. J Manuf Sci Eng 2012; 134(1): 011006.

20.

Sahoo

Patra

Singh

, et al. Modeling dynamic stability and cutting forces in micro milling of Ti6Al4V using intermittent oblique cutting finite element method simulation-based force coefficients. J Manuf Sci Eng 2020; 142(9): 091005.

21.

Heidari

Mirjalili

Faris

, et al. Harris hawks optimization: algorithm and applications. Future Gener Comput Syst 2019; 97: 849–872.

22.

Horng

Lin

SS.

Merging crow search into ordinal optimization for solving equality constrained simulation optimization problems. J Comput Sci 2017; 23: 44–57.

23.

Andradottir

Constrained Search via penalization for continuous simulation optimization. IEEE Trans Automat Contr 2020; 65(11): 4741–4752.

24.

Veeramani

Muthuswamy

Reinforcement learning based path planning of multiple agents of SwarmItFIX robot for fixturing operation in sheetmetal milling process. Proc IMechE, Part B: J Engineering Manufacture 2022; 236: 1667–1680.

25.

Sambakhé

Rouan

Bacro

, et al. Conditional optimization of a noisy function using a Kriging metamodel. J Glob Optim 2019; 73(3): 615–636.

26.

Teixeira

Nogal

O’Connor

, et al. Reliability assessment with density scanned adaptive Kriging. Reliab Eng Syst Saf 2020; 199: 106908.

27.

Yun

Zhou

, et al. AK-SYSi: an improved adaptive Kriging model for system reliability analysis with multiple failure modes by a refined U learning function. Struct Multidiscipl Optim 2019; 59(1): 263–278.

28.

Huang

Zhang

, et al. Probabilistic analysis of chatter stability in turning. Int J Adv Manuf Technol 2016; 87(9–12): 3225–3232.