Chatter stability prediction of milling considering nonlinearities

Abstract

Nonlinearities have been evidenced during the chatter vibration of milling. Machinability of the thin-walled part is feed rate and position-dependent, and is subject to process damping at low cutting speed. Therefore, chatter stability prediction of milling considering nonlinear cutting force, nonlinear structural stiffness and process damping is investigated. The cutting force and stiffness are established based on the polynomial model and the process damping is investigated based on the dissipated energy. The dynamic cutting force and stability lobes are solved in the time domain with coefficients updated at each iteration. By formulating the displacement as an expanded form via the perturbation method, the time-consuming solution of delay differential equations is avoided. After formulating the identification of the nonlinear model via cutting tests and modal tests, numerical simulations considering nonlinearities are carried out and compared with the analytical method. The proposed method attains high accuracy of classic time-domain solution, but with an improved computational efficiency. Finally, cutting tests are conducted to verify the prediction of cutting force and stability lobes.

Keywords

Thin-walled part chatter stability milling nonlinear perturbation method

Introduction

Chatter vibration is apt to occur during the milling process of thin-walled workpiece, which reduces machining efficiency and deteriorates workpiece surface. Many strategies for chatter avoidance have been proposed, such as cutting process optimization, structural dynamics modification through passive and active control. Si and Wang¹ presented an iterative compensation strategy to reduce the tool/workpiece deformation-induced surface error during the five-axis flank milling of thin-walled workpiece. Liu et al.² proposed a tool inclination method based on the adaption of the cutting force component to guarantee the machined surface finish of the thin-walled workpiece. Especially, cutting parameter optimization through the stability lobes is an efficient method and has been widely employed in the industry.³

The chatter stability prediction based on linear model achieves high accuracy for most cases. However, nonlinearities have been evidenced during the milling process. The machinability is feed rate dependent while the linear cutting force model is unable to predict the effect of feed rate. Meanwhile, the machinability is changing with the tool path due to the variation of workpiece stiffness, especially for a thin-walled part. Besides, low cutting speed is usually selected during machining hard-to-cut materials, and the machining stability is quite affected due to the introduction of process damping.

A large amount of research work has been conducted to investigate the nonlinearities during the milling process. The cutting force models formulated in a nonlinear mathematical form are mostly investigated. Lamikiz et al.⁴ presented calculation of linear, quadratic and cubic polynomial shear cutting coefficients of a ball-end cutting tool. Muñoa et al.⁵ studied the stability prediction in frequency domain based on the exponential instantaneous force model and demonstrated the influence of feed rate on the stability. Guzeev and Pimenov⁶ proposed a mathematical model that includes the cutting forces due to processes in the shear zone and force component due to processes at the real surface of tool. Besides, investigations into the nonlinear dynamics of the machine tool/cutter/workpiece are reported. Moradi et al.⁷ investigated the nonlinear dynamics of regenerative chatter and presented an extended dynamic model of the peripheral milling process. Yang et al.⁸ proposed a dynamic model for tool/workpiece assembly, and the influence of dynamic behaviours and tool feed direction were considered. Karimi and Nategh⁹ developed optimal tool path and workpiece setup of a hexapod machine tool by considering the kinematics nonlinearity. Furthermore, process damping has been another important factor resulting in the nonlinear behaviours of chatter. Budak and Tunc¹⁰ presented an identification and modelling method of process damping coefficients from chatter tests, after determining the indentation force coefficient through energy analysis. Feng et al.¹¹ studied the generation mechanism of process damping in the milling of thin-walled part. The current work aims to establish a nonlinear chatter model by incorporating the nonlinear cutting force, nonlinear stiffness of workpiece and process damping together. Meanwhile, the identification of the nonlinear model is presented.

The mathematical solution of chatter stability lobes has been investigated quite a lot. Dun et al.¹² proposed a numerical difference method based on Adams-Bashforth scheme and multi-modal scheme for stability prediction. Yang et al.¹³ developed an efficient decomposition-condensation method to predict the in-process workpiece dynamics for chatter prediction of thin-walled structure. Li et al.¹⁴ extended dynamic model of milling system by considering regeneration, helix angle and process damping into the high-order time domain algorithm. Song et al.¹⁵ proposed a time-space discretization method for milling stability prediction of thin-walled component based on thin-plate theory and mode superposition principle. Sellmeier and Denkena¹⁶ proposed a semi-discretization method considering the time dependency of the system by combining Ackermann’s method and the method of the piecewise constant subsystems. The time domain methods solve the delay differential equations directly, while the semi-discretization method utilizes a finite-dimensional transformation matrix to approximate a single-valued operator of infinite dimensions. The above numerical methods have been demonstrated as high accuracy, but low efficiency compared to the frequency domain method. Current work has devoted to solving the nonlinear problems in the prediction of milling stability. However, the introduction of the nonlinearities to the chatter model presents many difficulties. The solution is much more complicated as more parameters are related to the chatter model, and sacrifice of high order terms representing nonlinearities for the purpose of improving computational efficiency causes the reduction of accuracy.

Methods to solve nonlinearity have been presented in literature. Afazov et al.¹⁷ solved the equation of motion by integrating numerically in the time domain using the Runge-Kutta fourth-order method, considering nonlinearities of the micro-milling cutting force. Zhou et al.¹⁸ analysed nonlinear behaviours and developed a semi-analytic method to predict the shearing, plowing and friction forces in primary and tertiary shear zones. Moradi et al.¹⁹ investigated the nonlinear effect of structural and cutting force, tool wear and process damping, and simulated the forced vibration in milling process via perturbation method. Olvera et al.²⁰ presented the stability analysis of a thin-walled workpiece by using the enhanced multistage homotopy perturbation. Urbikain et al.²¹ proposed multimode analysis to predict chatter in large horizontal lathes by using the modern collocation method based on the Chebyshev polynomials. As the perturbation method is beneficial for simplifying mathematical calculation, it is employed to formulate the dynamic differential equation in current work. Afterward, the iteration method is utilized to solve stability lobes due to coefficients updating of the chatter model.

As machinability of the thin-walled part is feed rate and position-dependent, and is subject to process damping at low cutting speed, a nonlinear chatter model considering nonlinear cutting force, nonlinear stiffness and process damping is established. By formulating the displacement as an expanded form via the perturbation method, the dynamic cutting force and stability lobes are solved in the time domain with an improved computational efficiency. The remainder of the paper is organized as follows: the nonlinear chatter stability model and its solution are presented in Section 2. In Section 3, the identification procedure of the nonlinear chatter model is proposed based on the least square method. In Section 4, experiments of identifying the cutting force coefficients, stiffness coefficients and process damping coefficients are presented. In Section 5, the cutting force, vibration and chatter stability lobes are predicted and compared with the analytical method. Cutting tests are conducted in Section 6, and the paper is concluded in the end.

Theoretical formulation

Dynamic model of the milling process

Regenerative chatter vibration is easy to occur during the milling process of thin-walled part and is harmful to the machining process. As fluctuation of chip thickness is induced due to the vibration of machine tool/cutter and workpiece assembly, dynamic model of the milling process is demonstrated in Figure 1.

Figure 1.

Dynamic model of the milling process of thin-walled part.

Dynamic differential equations of the milling process can be established in x and y directions as equation (1).

MQ ″ + CQ' + KQ = F

(1)

where M, C and K are the mass, damping and stiffness matrix. F and Q are the cutting force and displacement vector.

F can be written as equation (2), and its components in x and y directions are calculated by the projection of tangential force F_t and radial force F_r.

F = [\begin{matrix} F_{x} \\ F_{y} \end{matrix}] \begin{matrix} , \end{matrix} {\begin{matrix} F_{x} = - F_{t} \cos ϕ - F_{r} \sin ϕ \\ F_{y} = + F_{t} \sin ϕ - F_{r} \cos ϕ \end{matrix}

(2)

ϕ is the immersion angle (Figure 1). For a discretized cutting flute, it is

\begin{matrix} ϕ (i, j, k) = ϕ_{st} + (i - 1) d ϕ + (j - 1) ϕ_{p} \\ - 2 (k - 1) \tan (β) dz / D \end{matrix}

(3)

where ϕ_st is the entry angle, ϕ_p is the pitch angle, D is the diameter of cutter and β is the helix angle. i, j and k are increments denoting the time, tooth and axial depth.

In order to incorporate the nonlinear influence of feed rate on cutting force, F_t and F_r are expressed in a polynomial form as equation (4).

{\begin{matrix} F_{t} = \sum_{i = 1}^{m} K_{tc, i} a_{p} h^{i} + K_{te} a_{p} \\ F_{r} = \sum_{i = 1}^{m} K_{rc, i} a_{p} h^{i} + K_{re} a_{p} \end{matrix} (m = 1, 2, . . .)

(4)

where K_tc,i and K_rc,i are the tangential and radial cutting force coefficients, K_te and K_re are edge constants. a_p is the axial depth of cut. h is dynamic chip thickness and is expressed as

h (ϕ) = h_{0} + Δ υ = h_{0} + Δ x \sin ϕ + Δ y \cos ϕ

(5)

where h₀ = csinϕ corresponds to static part of the chip thickness, and c is the feed rate per tooth. Δυ is the regenerative chip thickness, and Δx and Δy are the fluctuations of chip thickness in x and y directions.

According to equations (2)–(5), F_x and F_y can be rewritten as equation (6) when m = 3 (equation (4)). A detailed formulation of coefficients μ₁₁, $μ_{11}'$ , μ₁₂,…, μ₇, ${μ'}_{7}$ are presented in Appendix 1.

{\begin{matrix} F_{x} = - \frac{N}{2 π} (μ_{11} Δ {x_{c}}^{3} + μ_{12} Δ y_{c}^{3} + μ_{21} Δ {x_{c}}^{2} \\ + μ_{22} Δ y_{c}^{2} + μ_{31} Δ x_{c} + μ_{32} Δ y_{c} + 3 μ_{4} Δ x_{c}^{2} Δ y_{c} \\ + 3 μ_{5} Δ x_{c} Δ y_{c}^{2} + 2 μ_{6} Δ x_{c} Δ y_{c} + μ_{7}) \\ F_{y} = \frac{N}{2 π} (μ_{11}' Δ x_{c}^{3} + μ_{12}' Δ y_{c}^{3} + μ_{21}' Δ x_{c}^{2} + μ_{22}' Δ y_{c}^{2} \\ + μ_{31}' Δ x_{c} + μ_{32}' Δ y_{c} + 3 μ_{4}' Δ x_{c}^{2} Δ y_{c} \\ + 3 μ_{5}' Δ x_{c} Δ y_{c}^{2} + 2 μ_{6}' Δ x_{c} Δ y_{c} + μ_{7}') \end{matrix}

(6)

where Δx_c = (c+Δx)sinϕ and Δy_c = Δycosϕ are the chip thickness components. Meanwhile, the stiffness of thin-walled workpiece is established based on the polynomials as it varies with the tool path. Then, the stiffness matrix K can be expanded as equation (7).

K = [\begin{matrix} k_{x} & 0 \\ 0 & k_{y} \end{matrix}], {\begin{matrix} k_{x} = \sum_{i = 0}^{m} b_{x, i} x^{i} \\ k_{y} = \sum_{i = 0}^{m} b_{y, i} y^{i} \end{matrix} (m = 0, 1, 2, . . .)

(7)

where b_x,i and b_y,i are stiffness coefficients in x and y directions respectively.

The process damping is generated due to the contact between the wavy surface finish and round edge of the cutter, and it is dominant at low spindle speed. The tangential and radial process damping force is written as equation (8)) according to Jin and Altintas.²²

{\begin{matrix} F_{pt} (t) = C_{t} a_{p} υ' (t) / V_{c} \\ F_{pr} (t) = C_{r} a_{p} υ' (t) / V_{c} \end{matrix}

(8)

where C_t and C_r are the process damping coefficients, υ′(t) is the vibration velocity along the direction of chip thickness and V_c is the cutting speed.

By moving the term of process damping to the left side of equation (1), the damping matrix C can be formulated as

C = [\begin{matrix} C_{0 x} + (C_{t} \sin ϕ \cos ϕ + C_{r} \sin^{2} ϕ) a_{p} / V_{c} & (C_{t} \cos^{2} ϕ + C_{r} \sin ϕ \cos ϕ) a_{p} / V_{c} \\ - (C_{t} \sin^{2} ϕ - C_{r} \sin ϕ \cos ϕ) a_{p} / V_{c} & C_{0 y} - (C_{t} \sin ϕ \cos ϕ - C_{r} \cos^{2} ϕ) a_{p} / V_{c} \end{matrix}]

(9)

where C_0x and C_0y are the inherent damping which can be obtained via modal tests.

Prediction of chatter stability

The dynamic cutting force, vibration and chatter stability lobes during the milling process of thin-walled part are predicted in time domain with coefficients updated at each iteration, enabling the simulation of entire tool path. As the introduction of high order polynomials to represent nonlinearities complicates the mathematical modelling and solution of chatter stability, the perturbation method is employed to formulate the displacement as an expanded form for the purpose of improving computational efficiency. According to the perturbation method, an approximate analytical solution of displacement x and y (equation (1)) can be formulated as equation (10).

{\begin{matrix} x = x_{0} + x_{1} ε + x_{2} ε^{2} + \dots; \\ Δ x = Δ x_{0} + Δ x_{1} ε + Δ x_{2} ε^{2} + \dots \\ y = y_{0} + y_{1} ε + y_{2} ε^{2} + \dots; \\ Δ y = Δ y_{0} + Δ y_{1} ε + Δ y_{2} ε^{2} + \dots \end{matrix}

(10)

where Δx_j = x_j(t)−x_j(t−T) and Δy_j = y_j(t)−y_j(t−T) (j = 0, 1, 2, …). ε is the perturbation coefficient, and T is the tooth period.

Then equation (1) is rewritten in Appendix 2 according to the order of ε. For the simplification of the solution, only the zero and first orders are remained. The solution of x₀, y₀, x₁, y₁ is given in equation (11) based on the Euler formula.

{\begin{matrix} x_{0} = a_{0} \cos (ω t) + A_{0}; y_{0} = b_{0} \cos (ω t) + B_{0} \\ x_{1} = a_{1} \exp [(μ_{x} + i ω) t]; y_{1} = b_{1} \exp [(μ_{y} + i ω) t] \end{matrix}

(11)

where ω is angular velocity of the spindle. According to the initial condition

x (0) = 0, y (0) = 0

(12)

And, it is derived as A₀+a₀+a₁ = 0 and B₀+b₀+b₁ = 0. The a₀, b₀, A₀, B₀, a₁, b₁ are solved in Appendix 3 according to the Taylor series expansion. Similarly, μ_x and μ_y can be solved. Finally, an approximate solution of x and y is shown as equation (13).

{\begin{matrix} x = a_{0} \cos (ω t) + A_{0} + ε a_{1} \exp [(μ_{x} + i ω) t] \\ y = b_{0} \cos (ω t) + B_{0} + ε b_{1} \exp [(μ_{y} + i ω) t] \end{matrix}

(13)

As a result, the cutting force and vibration can be solved, and stability lobes are plotted by determining whether the process is unstable (Figure 2).

Figure 2.

Flow chart of chatter stability prediction.

Formulation of coefficients identification

Cutting force coefficients

The cutting force coefficients are identified during a stable cutting when the fluctuation of the chip thickness has little influence on cutting force and can be ignored. The instantaneous cutting force F_x/F_y/F_z is expanded to mth order polynomial of static chip thickness h₀ as equation (14).

{\begin{cases} F_{x} (ϕ) = - [(\sum_{i = 0}^{m} K_{t c, i} \times h_{0}^{i + 1} + K_{t e}) \cos ϕ \\ + (\sum_{i = 0}^{m} K_{r c, i} \times h_{0}^{i + 1} + K_{r e}) \sin ϕ] a_{p} \\ F_{y} (ϕ) = + [(\sum_{i = 0}^{m} K_{t c, i} \times h_{0}^{i + 1} + K_{t e}) \sin ϕ \\ - (\sum_{i = 0}^{m} K_{r c, i} \times h_{0}^{i + 1} + K_{r e}) \cos ϕ] a_{p} \\ F_{z} (ϕ) = (\sum_{i = 0}^{m} K_{a c, i} \times h_{0}^{i + 1} + K_{a e}) a_{p} \end{cases}

(14)

The resultant force can be calculated by the superposition of F_x(ϕ), F_y(ϕ) and F_z(ϕ), and the average cutting force $\bar{F_{x}},$ $\bar{F_{y}}$ and $\bar{F_{z}}$ are obtained by dividing 2π in equation (15).

[\begin{matrix} \frac{\bar{F_{x}}}{N a_{p}} \\ \frac{\bar{F_{y}}}{N a_{p}} \\ \frac{\bar{F_{z}}}{N a_{p}} \end{matrix}] = [\begin{matrix} - \frac{K_{rc, m}}{2 π} \int_{0}^{π} {(\sin ϕ)}^{m + 2} d ϕ & \dots & - \frac{K_{rc, 0}}{2 π} \int_{0}^{π} {(\sin ϕ)}^{2} d ϕ & - \frac{K_{re}}{2 π} \int_{0}^{π} \sin ϕ d ϕ \\ \frac{K_{tc, m}}{2 π} \int_{0}^{π} {(\sin ϕ)}^{m + 2} d ϕ & \dots & \frac{K_{tc, 0}}{2 π} \int_{0}^{π} {(\sin ϕ)}^{2} d ϕ & \frac{K_{te}}{2 π} \int_{0}^{π} \sin ϕ d ϕ \\ \frac{K_{ac, m}}{2 π} \int_{0}^{π} {(\sin ϕ)}^{m + 1} d ϕ & \dots & \frac{K_{ac, 0}}{2 π} \int_{0}^{π} \sin ϕ d ϕ & \frac{1}{2 π} \int_{0}^{π} K_{ae} d ϕ \end{matrix}] [\begin{matrix} c^{m + 1} \\ \dots \\ c \\ 1 \end{matrix}]

(15)

The identification algorithm of cutting force coefficients is provided in Appendix 4, and the identified coefficients when m is even (equation (16)) and odd (equation (17)) are presented as follows.

{\begin{matrix} K_{re} = - π a_{x, 0} & K_{rc, 0} = - 4 a_{x, 1} & \dots & K_{rc, m} = - 2 a_{x, m + 1} Π_{i = 0}^{m / 2} \frac{2 i + 2}{2 i + 1} \\ K_{te} = π a_{y, 0} & K_{tc, 0} = 4 a_{y, 1} & \dots & K_{tc, m} = 2 a_{y, m + 1} Π_{i = 0}^{m / 2} \frac{2 i + 2}{2 i + 1} \\ K_{ae} = 2 a_{z, 0} & K_{ac, 0} = π a_{z, 1} & \dots & K_{ac, m} = π a_{z, m + 1} Π_{i = 1}^{m / 2} \frac{2 i + 2}{2 i + 1} \end{matrix} (m = 2, 4, . . .)

(16)

{\begin{matrix} K_{re} = - π a_{x, 0} & K_{rc, 0} = - 4 a_{x, 1} & \dots & K_{rc, m} = - π a_{x, m + 1} Π_{i = 1}^{(m + 1) / 2} \frac{2 i + 1}{2 i} \\ K_{te} = π a_{y, 0} & K_{tc, 0} = 4 a_{y, 1} & \dots & K_{tc, m} = π a_{y, m + 1} Π_{i = 1}^{(m + 1) / 2} \frac{2 i + 1}{2 i} \\ K_{ae} = 2 a_{z, 0} & K_{ac, 0} = π a_{z, 1} & \dots & K_{ac, m} = 2 a_{z, m + 1} Π_{i = 1}^{(m + 1) / 2} \frac{2 i}{2 i - 1} \end{matrix} (m = 1, 3, . . .)

(17)

where coefficients a_x,i, a_y,i and a_z,i (i=1, …, m) are identified with the least square method.

Stiffness coefficients

As the stiffness of thin-walled part is dependent on the position, the stiffness vector, that is in x direction, is

K_{x} = {[\begin{matrix} k_{x 0} & k_{x 1} & \dots & k_{xi} \end{matrix}]}^{T}

(18)

where k_xi denotes the stiffness at points x_i. The vector of stiffness coefficients is

B_{x} = {[\begin{matrix} b_{x, 0} & b_{x, 1} & \dots & b_{x, m} \end{matrix}]}^{T}

(19)

Then equation (20) is formulated based on the least square method.

X^{T} X B_{x} = X^{T} K_{x}

(20)

where $X = [\begin{matrix} 1 & x_{0} & {x_{0}}^{2} & \dots & {x_{0}}^{m} \\ 1 & x_{1} & {x_{1}}^{2} & \dots & {x_{1}}^{m} \\ 1 & x_{2} & {x_{2}}^{2} & \dots & {x_{2}}^{m} \\ \dots & \dots & \dots & \dots & \dots \\ 1 & x_{i} & {x_{i}}^{2} & \dots & {x_{i}}^{m} \end{matrix}]$ .

By substituting the stiffness identified via modal tests into above equation, the coefficients vector B_x are solved. Similarly, the stiffness coefficients in y direction can be identified.

Process damping coefficients

The process damping coefficients C_t and C_r are identified from the dissipated energy as²²

C_{t} = \frac{\int_{0}^{τ_{p}} F_{pt} (t) υ' (t) dt}{\frac{a_{p}}{V_{c}} \int_{0}^{τ_{p}} υ' {(t)}^{2} dt}, C_{r} = \frac{\int_{0}^{τ_{p}} F_{pr} (t) υ' (t) dt}{\frac{a_{p}}{V_{c}} \int_{0}^{τ_{p}} υ' {(t)}^{2} dt}

(21)

τ_p is the vibration period. The integration terms can be replaced by the summation of discrete sampling points via experiments. Among the components (i.e. static force, dynamic force due to regenerative effect and process damping force) of cutting force, the process damping force is effective at low cutting speed. Therefore, F_pt and F_pr can be estimated by the experimental force and theoretical force (equation (4)) calculated by cutting force coefficients identified at high cutting speed. The vibration velocity υ′ can be obtained by the vibration of thin-walled part in x and y directions according to equation (5).

Experimental identification of coefficients

Identification of cutting force coefficients

Milling tests are conducted to identify the cutting force coefficients on a machining centre of DMG DMC 635V. A carbide end mill with diameter D = 12 mm, number of teeth N = 3 and helix angle β = 30° is used. The workpiece is Aluminum alloy 7075. The spindle speed n is selected as 3000 rpm while the axial depth of cut a_p is 0.5 mm, 1.0 mm, 1.5 mm and feed rate per tooth c is 0.033 mm, 0.067 mm, 0.010 mm. The cutting force F_x, F_y and F_z are measured by the dynamometer Kistler 9257B. The cutting force coefficients based on the polynomials are identified in Table 1, according to equations (16) and (17).

Table 1.

Identified cutting force coefficients.

m = 1	K_tc,1 = 1054.8 N/mm²	K_rc,1 = 389.5 N/mm²	K_ac,1 = 239.9 N/mm²
	K_te = 12.2 N/mm	K_re = 10.1 N/mm	K_ae = 1.3 N/mm
m = 2	K_tc,2 = 2164.0 N/mm³	K_rc,2 = 5362.5 N/mm³	K_ac,2 = 1597.2 N/mm³
	K_tc,1 = 1298.8 N/mm²	K_rc,1 = 994.1 N/mm²	K_ac,1 = 406.5 N/mm²
	K_te = 6.9 N/mm	K_re = 2.9 N/mm	K_ae = 1.7 N/mm
m = 3	K_tc,3 = 163522.5 N/mm⁴	K_rc,3 = 57875.9 N/mm⁴	K_ac,3 = 25522.1 N/mm⁴
	K_tc,2 = 16612.6 N/mm³	K_rc,2 = 205.1 N/mm³	K_ac,2 = 750.0 N/mm³
	K_tc,1 = 597.4 N/mm²	K_rc,1 = 840.4 N/mm²	K_ac,1 = 315.2 N/mm²
	K_te = 4.8 N/mm	K_re = 6.1 N/mm	K_ae = 3.7 N/mm

Identification of stiffness coefficients

Modal tests are utilized to identify the changing stiffness of thin-walled workpiece at different positions. The impact damper PCB 086C03 with a sensitivity of 2.25 mV/N, accelerometer PCB 352A21 with a sensitivity of 106.17 mV/g and data acquisition card NI 9234 are used. Only the stiffness in x direction is identified in current work due to its flexibility compared to y direction. The stiffness at Points A, B, C and D (Figure 3) are 8.5e7 N/m, 6.3e7 N/m, 7.6e7 N/m and 9.2e7 N/m, respectively. By substituting them into equation (20), the identified stiffness coefficients (KN/m²) based on the second-order and third-order polynomials are presented in Table 2.

Figure 3.

Experimental setup.

Table 2.

Identified stiffness coefficients.

m = 2	b_x,0 = 9.9e8 N/m	b_x,1 = −1.9e10 N/m²	b_x,2 = 9.2e10 N/m³
m = 3	b_x,0 = 6.9e8 N/m	b_x,1 = −8.9e9 N/m²	b_x,2 = −8.4e9 N/m³	b_x,3 = 3.4e11 N/m⁴

Identification of process damping coefficients

Down milling tests are conducted to identify process damping coefficients. Three sets of parameters are utilized including axial depth of cut a_p = 0.5 mm and feed rate F = 300 mm/min, a_p=0.5 mm and F = 600 mm/min, a_p = 1 mm and F = 300 mm/min, while spindle speed n = 3000 rpm and width of cut a_e = 12 mm are kept constant. The cutting force and vibration velocity (Figure 4) in x and y directions during a single period are presented. The vibration velocity is measured by velocity sensor PCB 640B02 with a sensitivity of 0.31 mA/(mm/s). According to equation (22), the process damping coefficients are identified as C_t = 1.44e4 N/m and C_r = 3.98e4 N/m.

Figure 4.

Experimental cutting force and velocity: (a) a_p = 0.5 mm, F = 300 mm/min, (b) a_p = 0.5 mm, F = 600 mm/min and (c) a_p = 1 mm, F = 300 mm/min.

Numerical simulation

Based on the polynomial model of cutting force and stiffness, the milling force and stability lobes are predicted after identifying the coefficients, and then compared with analytical method and experiments. The cutter and workpiece are the same as above section.

Static simulation of the milling force

Three configurations of slot milling are simulated. The axial depth of cut a_p and feed rate per tooth c are changed while spindle speed n is 3000 rpm. It is seen that the simulated cutting force via proposed method is approximate to the analytical method¹³ and experiments with 5% difference on amplitude (Figure 5), but the computation time is much longer due to iterated solution. However, dynamic milling force and vibration are able to be simulated due to the coefficients updated during each iteration.

Figure 5.

Static simulation of milling force in x, y and z directions: (a) a_p = 0.5 mm, c = 0.05 mm, (b) a_p = 0.5 mm, c = 0.10 mm and (c) a_p = 1.0 mm, c = 0.05 mm.

Dynamic simulation of the milling force and vibration

Dynamic milling force and vibration of the thin-walled part (Figure 3) are simulated. Stiffness variation is considered to simulate actual milling process. For down milling of spindle speed n = 6000 rpm, feed rate per tooth c = 0.05 mm and width of cut a_e = 4 mm, the cutting force increases from 400 N, 800 N and 1000 N to 700 N, 1400 N and 1800 N from 0.1 s to 1.0 s when a_p = 4 mm, 6 mm and 8 mm (Figure 6). Meanwhile, the acceleration increases from 8 g, 10 g and 12 g to 21 g, 28 g and 42 g. The cutting force and workpiece vibration are maximum at t = 1.0 s when the tool position is located at the middle position and corresponding stiffness is the weakest. It is observed that there are 25%, 10% and 22% difference on the simulated workpiece vibration compared with experiments.

Figure 6.

Dynamic simulation of milling force and workpiece vibration considering stiffness variation: (a) a_p = 4 mm, (b) a_p = 6 mm and (c) a_p = 8 mm.

Prediction of stability lobes

Chatter stability simulation via the proposed method based on linear model is compared with analytical method¹³ and experiments. For slot milling of the thin-walled part, the critical depth of cut is 1.6 mm (Figure 7(a)) at x = 90 mm while it is 1.4 mm and 1.5 mm for analytical and semi-discretization methods respectively. Besides, the critical depth of cut for the three methods are 1.7 mm, 1.5 mm and 1.6 mm (Figure 7(b)) when down milling with a_e = 6 mm, and are 1.8 mm, 1.7 mm and 1.7 mm (Figure 7(c)) when down milling with a_e = 4 mm. Simulations of the three methods are quite close.

Figure 7.

Comparison of stability lobes: (a) slot milling, (b) down milling, a_e = 6 mm and (c) down milling, a_e = 4 mm.

Next, chatter stability simulations based on the third-order polynomial model of cutting force are demonstrated (Figure 8). For slot milling, the critical depth of cut increases from 1.3 mm to 2.5 mm as c increases from 0.04 mm to 0.08 mm. Meanwhile, the critical depth of cut increases from 1.4 mm and 1.7 mm to 2.7 mm and 2.8 mm as c increases for down milling, a_e = 6 mm and 4 mm. It is observed that c has an impact on the stability lobes, and the critical depth of cut increases as c increases.

Figure 8.

Chatter stability lobes when c = 0.04 mm, 0.06 mm and 0.08 mm: (a) slot milling, (b) down milling, a_e = 6 mm and (c) down milling, a_e = 4 mm.

Furthermore, the influence of nonlinear stiffness on chatter stability lobes is simulated at different tool positions. For down milling with a_e = 4 mm and c = 0.05 mm, three locations indicated in Figure 3 are simulated (Figure 9(a)). The critical depth of cut is 2.4 mm at x = 40 mm (Point A), 1.7 mm at x = 70 mm (Point B) and 2.1 mm at x = 90 mm (Point C). The minimum critical depth of cut is located at x = 70 mm. Besides, the influence of process damping on chatter stability lobes is simulated for down milling, a_e = 4 mm and c = 0.04 mm. For x = 90 mm, the critical depth of cut increases from 2.2 mm to 3.3 mm at n = 2200 rpm after process damping is considered (Figure 9(b)), corresponding to an increase of 60%. However, there is little difference after the spindle speed of 3800 rpm.

Figure 9.

Chatter stability lobes: (a) at different tool positions and (b) without and with process damping (x = 90 mm).

Stability analysis of thin-walled part

The stability lobes of thin-walled part considering nonlinear cutting force, nonlinear stiffness and process damping for down milling, a_e = 4 mm and c = 0.04 mm are illustrated in Figure 10(a), and the relation between the critical depth of cut and tool position is illustrated in Figure 10(b). Four spindle speeds are considered. The stability limit is maximum when n = 2000 rpm and is minimum when n = 8000 rpm. For n = 8000 rpm, the critical depth of cut decreases from 2.4 mm to 1.7 mm when the tool position moves from x = 40 mm to 70 mm, and increases to 2.1 mm at x = 90 mm. In order to verify the above simulation, the tool path is divided into three segments (Figure 10(b)) and corresponding machining signals are recorded. The Path I corresponds to x = 40 mm to 60 mm, the Path II corresponds to x = 60 mm to 80 mm, and the Path III corresponds to x = 80 mm to 100 mm.

Figure 10.

Stability lobes considering nonlinearities: (a) stability lobes at different tool positions and (b) relation between the critical depth of cut and tool positions.

Cutting tests

Down milling tests (Figure 3) are conducted to verify the above simulations. The width of cut a_e = 4 mm and feed rate F = 900 mm/min are kept constant, and four configurations of the axial depth of cut a_p and spindle speed n are utilized. The machining sound is collected by the microphone PCB 378B02.

When a_p = 3 mm and a_p = 5 mm, the maximum machining sound are 1.9 Pa (Figure 11(a)) and 2.8 Pa (Figure 11(b)) for n = 6000 rpm, and increase to 7.5 Pa (Figure 11(c)) and 9.8 Pa (Figure 11(d)) for n = 8000 rpm. The increase is quite large, and it is explained that chatter vibration occurs as can be observed from the simulation. According to the chatter stability lobes in Figure 10(b), the minimum critical depth of cut is 5.2 mm located at x = 70 mm when n= 6000 rpm, and critical depth of cut is 3.4 mm and 2.8 mm at x = 40 and 70 mm when n = 8000 rpm. Besides, it is noted that the maximum vibration occurs during path II as corresponding stiffness of thin-walled part is the weakest.

Figure 11.

The machining sound in time and frequency domain: (a) n = 6000 rpm, a_p = 3 mm, (b) n = 6000 rpm, a_p = 5 mm, (c) n = 8000 rpm, a_p = 3 mm and (d) n = 8000 rpm, a_p = 5 mm.

The time domain signals during tool path I, II and III (Figure 10(b)) are transformed to the frequency domain. For n = 6000 rpm, a_p = 3 mm (Figure 11(a)) and a_p = 5 mm (Figure 11(b)), the tooth passing frequency of 300 Hz is dominant, indicating a stable cutting process. When n = 8000 rpm and a_p = 3 mm (Figure 11(c)), the tooth passing frequency of 400 Hz is dominant for Path I. However, the chatter frequency of 1300 Hz increases by 60% for Path II. It indicates that the cutting process turns to be unstable from Path I to Path II due to the stiffness decrease. The machining process of Path III is evolved to be stable due to the stiffness increase. For n = 8000 rpm and a_p = 5 mm (Figure 11(d)), the chatter frequency of 1600 Hz and tooth passing frequency of 400 Hz are clear for entire tool path, and is treated as unstable by combining with the machining sound in the time domain.

Conclusion

The machining of thin-walled part with high quality and efficiency has been a challenge due to its low and position-dependent machinability. Chatter theory in the milling process of thin-walled part considering nonlinearities (i.e. cutting force, structural stiffness and process damping) is investigated, and stability lobes are solved in time domain with an improved computational efficiency. Cutting tests are carried out to verify the simulation. The paper is concluded as follows.

The nonlinearities are introduced to be able to predict real behaviours, and coefficients identification of nonlinear chatter model are presented. However, the introduction of high order polynomials to represent nonlinearities complicates the mathematical modelling and solution of chatter stability in time domain. By formulating the displacement as an expanded form via the perturbation method, the time-consuming solution of delay differential equations is avoided. It is demonstrated that the proposed method attains high accuracy of classic time-domain solution, while the computational efficiency is much improved.

The dynamic cutting force, vibration and stability lobes considering nonlinearities are solved in the time domain with coefficients updated at each iteration, enabling the simulation of entire tool path. The proposed method achieves approximated accuracy with the semi-discretization method, but with the sacrifice of efficiency compared to the analytical frequency domain method. Besides, the adoption of specific criterion to determine an unstable process could cause deviation on the accuracy.

After formulating the displacement as an expanded form, numerical terms of the expansion with minor magnitude are neglected during the simulation of milling process, and the effect on the accuracy needs further investigation. Besides, only the dynamics of thin-walled part are considered due to its flexibility compared to the machine tool/cutter assembly. The proposed model can be extended to incorporate the dynamics of machine tool/cutter/workpiece for future work.

Footnotes

Appendix 1

μ ₁₁, $μ_{11}'$ , μ₁₂, …, μ₇, $μ_{7}'$ are presented as

(22)

\begin{matrix} μ_{11} = a_{p} (\frac{1}{4} K_{tc 3} + \frac{3 π}{16} K_{rc 3}) \begin{matrix} , & μ_{12} = a_{p} (\frac{3 π}{16} K_{tc 3} + \frac{1}{4} K_{rc 3}) \begin{matrix} , & μ_{21} = a_{p} (\frac{1}{3} K_{tc 2} + \frac{2}{3} K_{rc 2}), \end{matrix} \end{matrix} \\ μ_{22} = a_{p} (\frac{2}{3} K_{tc 2} + \frac{1}{3} K_{rc 2}) \begin{matrix} , & μ_{31} = a_{p} (\frac{1}{2} K_{tc 1} + \frac{π}{4} K_{rc 1}) \end{matrix} \begin{matrix} , & μ_{32} = a_{p} (\frac{π}{4} K_{tc 1} + \frac{1}{2} K_{rc 1}) \end{matrix}, \\ μ_{4} = a_{p} (\frac{π}{16} K_{tc 3} + \frac{1}{4} K_{rc 3}) \begin{matrix} , & μ_{5} = a_{p} (\frac{1}{4} K_{tc 3} + \frac{π}{16} K_{rc 3}) \end{matrix} \begin{matrix} , & μ_{6} = a_{p} (\frac{1}{3} K_{tc 2} + \frac{1}{3} K_{rc 2}) \end{matrix}, \\ μ_{11}' = a_{p} (\frac{3 π}{16} K_{tc 3} - \frac{1}{4} K_{rc 3}) \begin{matrix} , & μ_{12}' = a_{p} (\frac{1}{4} K_{tc 3} - \frac{3 π}{16} K_{rc 3}) \begin{matrix} , & μ_{21}' = a_{p} (\frac{2}{3} K_{tc 2} - \frac{1}{3} K_{rc 2}), \end{matrix} \end{matrix} \\ μ_{22}' = a_{p} (\frac{1}{3} K_{tc 2} - \frac{2}{3} K_{rc 2}) \begin{matrix} , & μ_{31}' = a_{p} (\frac{π}{4} K_{tc 1} - \frac{1}{2} K_{rc 1}) \end{matrix} \begin{matrix} , & μ_{32}' = a_{p} (\frac{1}{2} K_{tc 1} - \frac{π}{4} K_{rc 1}) \end{matrix}, \\ μ_{4}' = a_{p} (\frac{1}{4} K_{tc 3} - \frac{π}{16} K_{rc 3}) \begin{matrix} , & μ_{5}' = a_{p} (\frac{π}{16} K_{tc 3} - \frac{1}{4} K_{rc 3}) \end{matrix} \begin{matrix} , & μ_{6}' = a_{p} (\frac{1}{3} K_{tc 2} - \frac{1}{3} K_{rc 2}) \end{matrix}, \\ μ_{7} = K_{te} + K_{re} \begin{matrix} , & μ_{7}' = \end{matrix} K_{te} - K_{re} \end{matrix}

Appendix 2

By writing the damping ratio ζ_x and ζ_y as εη_x and εη_y, equation (1) is formulated as equations (23)-(25) according to the order of ε.

(23)

ε^{0} : {\begin{matrix} x_{0}^{″} + ω_{x}^{2} x_{0} = - \frac{N}{2 π m_{x}} [μ_{11} (x_{0} - x_{0 T} + c)^{3} + μ_{12} (y_{0} - y_{0 T})^{3} + μ_{21} (x_{0} - x_{0 T} + c)^{2} + μ_{22} (y_{0} - y_{0 T})^{2} \\ + μ_{31} (x_{0} - x_{0 T} + c) + μ_{32} (y_{0} - y_{0 T}) + 3 μ_{4} (x_{0} - x_{0 T} + c)^{2} (y_{0} - y_{0 T}) + 3 μ_{5} (x_{0} - x_{0 T} + c) (y_{0} - y_{0 T})^{2} \\ + 2 μ_{6} (x_{0} - x_{0 T} + c) (y_{0} - y_{0 T}) + μ_{7}] \\ y_{0}^{″} + ω_{y}^{2} y_{0} = - \frac{N}{2 π m_{y}} [μ_{11}' (x_{0} - x_{0 T} + c)^{3} + μ_{12}' (y_{0} - y_{0 T})^{3} + μ_{21}' (x_{0} - x_{0 T} + c)^{2} + μ_{22}' (y_{0} - y_{0 T})^{2} \\ + μ_{31}' (x_{0} - x_{0 T} + c) + μ_{32}' (y_{0} - y_{0 T}) + 3 μ_{4}' (x_{0} - x_{0 T} + c)^{2} (y_{0} - y_{0 T}) + 3 μ_{5}' (x_{0} - x_{0 T} + c) (y_{0} - y_{0 T})^{2} \\ + 2 μ_{6}' (x_{0} - x_{0 T} + c) (y_{0} - y_{0 T}) + μ_{7}'] \end{matrix}

(24)

ε^{1} : {\begin{matrix} x_{1}^{″} + ω_{x}^{2} x_{1} = - η_{x} x_{0}' - \frac{N}{2 π m_{x}} [3 μ_{11} (x_{0} - x_{0 T} + c)^{2} (x_{1} - x_{1 T}) + 3 μ_{12} (y_{0} - y_{0 T})^{2} (y_{1} - y_{1 T}) \\ + 2 μ_{21} (x_{0} - x_{0 T} + c) (x_{1} - x_{1 T}) + 2 μ_{22} (y_{0} - y_{0 T}) (y_{1} - y_{1 T}) + μ_{31} (x_{1} - x_{1 T}) + μ_{32} (y_{1} - y_{1 T}) \\ + 3 μ_{4} (x_{0} - x_{0 T} + c)^{2} (y_{1} - y_{1 T}) + 6 μ_{4} (x_{0} - x_{0 T} + c) (x_{1} - x_{1 T}) (y_{0} - y_{0 T}) + 3 μ_{5} (x_{1} - x_{1 T}) (y_{0} - y_{0 T})^{2} \\ 6 μ_{5} (x_{0} - x_{0 T} + c) (y_{0} - y_{0 T}) (y_{1} - y_{1 T}) + 2 μ_{6} (x_{0} - x_{0 T} + c) (y_{1} - y_{1 T}) + 2 μ_{6} (x_{1} - x_{1 T}) (y_{0} - y_{0 T}) \\ y_{1}^{″} + ω_{y}^{2} y_{1} = - η_{y} y_{0}' - \frac{N}{2 π m_{y}} [3 μ_{11}' (x_{0} - x_{0 T} + c)^{2} (x_{1} - x_{1 T}) + 3 μ_{12}' (y_{0} - y_{0 T})^{2} (y_{1} - y_{1 T})] \\ + 2 μ_{21}' (x_{0} - x_{0 T} + c) (x_{1} - x_{1 T}) + 2 μ_{22}' (y_{0} - y_{0 T}) (y_{1} - y_{1 T}) + μ_{31}' (x_{1} - x_{1 T}) + μ_{32}' (y_{1} - y_{1 T}) \\ + 3 μ_{4}' (x_{0} - x_{0 T} + c)^{2} (y_{1} - y_{1 T}) + 6 μ_{4}' (x_{0} - x_{0 T} + c) (x_{1} - x_{1 T}) (y_{0} - y_{0 T}) + 3 μ_{5}' (x_{1} - x_{1 T}) (y_{0} - y_{0 T})^{2} \\ 6 μ_{5}' (x_{0} - x_{0 T} + c) (y_{0} - y_{0 T}) (y_{1} - y_{1 T}) + 2 μ_{6}' (x_{0} - x_{0 T} + c) (y_{1} - y_{1 T}) + 2 μ_{6}' (x_{1} - x_{1 T}) (y_{0} - y_{0 T})] \end{matrix}

(25)

ε^{2} : {\begin{matrix} x_{2}^{″} + ω_{x}^{2} x_{2} = - η_{x} x_{1}' - \frac{N}{2 π m_{x}} [3 μ_{11} (x_{0} - x_{0 T} + c) (x_{1} - x_{1 T})^{2} + 3 μ_{11} (x_{0} - x_{0 T} + c)^{2} (x_{2} - x_{2 T}) \\ + 3 μ_{12} (y_{0} - y_{0 T}) (y_{1} - y_{1 T})^{2} + 3 μ_{12} (y_{0} - y_{0 T})^{2} (y_{2} - y_{2 T}) + 2 μ_{21} (x_{0} - x_{0 T} + c) (x_{2} - x_{2 T}) \\ + μ_{21} (x_{1} - x_{1 T})^{2} + 2 μ_{22} (y_{0} - y_{0 T}) (y_{2} - y_{2 T}) + μ_{22} (y_{1} - y_{1 T})^{2} + μ_{31} (x_{2} - x_{2 T}) + μ_{32} (y_{2} - y_{2 T}) \\ + 3 μ_{4} (x_{0} - x_{0 T} + c)^{2} (y_{2} - y_{2 T}) + 3 μ_{4} (x_{1} - x_{1 T})^{2} (y_{0} - y_{0 T}) + 6 μ_{4} (x_{0} - x_{0 T} + c) (x_{1} - x_{1 T}) (y_{1} - y_{1 T}) \\ + 6 μ_{4} (x_{0} - x_{0 T} + c) (x_{2} - x_{2 T}) (y_{0} - y_{0 T}) + 3 μ_{5} (x_{0} - x_{0 T} + c) (y_{1} - y_{1 T})^{2} + 3 μ_{5} (x_{2} - x_{2 T}) (y_{0} - y_{0 T})^{2} \\ + 6 μ_{5} (x_{0} - x_{0 T} + c) (y_{0} - y_{0 T}) (y_{2} - y_{2 T}) + 6 μ_{5} (x_{1} - x_{1 T}) (y_{0} - y_{0 T}) (y_{1} - y_{1 T}) \\ + 2 μ_{6} (x_{0} - x_{0 T} + c) (y_{2} - y_{2 T}) + 2 μ_{6} (x_{2} - x_{2 T}) (y_{0} - y_{0 T}) \\ y_{2}^{″} + ω_{y}^{2} y_{2} = - η_{y} y_{1}' - \frac{N}{2 π m_{y}} [3 μ_{11}' (x_{0} - x_{0 T} + c) (x_{1} - x_{1 T})^{2} + 3 μ_{11}' (x_{0} - x_{0 T} + c)^{2} (x_{2} - x_{2 T}) \\ + 3 μ_{12}' (y_{0} - y_{0 T}) (y_{1} - y_{1 T})^{2} + 3 μ_{12}' (y_{0} - y_{0 T})^{2} (y_{2} - y_{2 T}) + 2 μ_{21}' (x_{0} - x_{0 T} + c) (x_{2} - x_{2 T}) \\ + μ_{21}' (x_{1} - x_{1 T})^{2} + 2 μ_{22}' (y_{0} - y_{0 T}) (y_{2} - y_{2 T}) + μ_{22}' (y_{1} - y_{1 T})^{2} + μ_{31}' (x_{2} - x_{2 T}) + μ_{32}' (y_{2} - y_{2 T}) \\ + 3 μ_{4}' (x_{0} - x_{0 T} + c)^{2} (y_{2} - y_{2 T}) + 3 μ_{4}' (x_{1} - x_{1 T})^{2} (y_{0} - y_{0 T}) + 6 μ_{4}' (x_{0} - x_{0 T} + c) (x_{1} - x_{1 T}) (y_{1} - y_{1 T}) \\ + 6 μ_{4}' (x_{0} - x_{0 T} + c) (x_{2} - x_{2 T}) (y_{0} - y_{0 T}) + 3 μ_{5}' (x_{0} - x_{0 T} + c) (y_{1} - y_{1 T})^{2} + 3 μ_{5}' (x_{2} - x_{2 T}) (y_{0} - y_{0 T})^{2} \\ + 6 μ_{5}' (x_{0} - x_{0 T} + c) (y_{0} - y_{0 T}) (y_{2} - y_{2 T}) + 6 μ_{5}' (x_{1} - x_{1 T}) (y_{0} - y_{0 T}) (y_{1} - y_{1 T}) \\ + 2 μ_{6}' (x_{0} - x_{0 T} + c) (y_{2} - y_{2 T}) + 2 μ_{6}' (x_{2} - x_{2 T}) (y_{0} - y_{0 T}) \end{matrix}

Appendix 3

According to the Taylor series, the coefficients of ε¹ are expanded as

(26)

{\begin{matrix} a_{0} (ω_{x}^{2} - ω^{2}) [1 - \frac{{(ω t)}^{2}}{2}] + A_{0} ω_{x}^{2} = - \frac{N}{2 π m_{x}} [\frac{μ_{21} a_{0}^{2} ω^{4}}{4} (T^{2} - 2 Tt)^{2} + \frac{μ_{21} a_{0} ω^{2} c}{2} (T^{2} - 2 Tt) \\ + \frac{μ_{22} b_{0}^{2} ω^{4}}{4} (T^{2} - 2 Tt)^{2} + \frac{μ_{31} a_{0} ω^{2}}{4} (T^{2} - 2 Tt) + μ_{31} c_{f} + \frac{μ_{32} b_{0} ω^{2}}{4} (T^{2} - 2 Tt) \\ + \frac{μ_{6} a_{0} b_{0} ω^{4}}{2} (T^{2} - 2 Tt)^{2} + \frac{μ_{6} b_{0} c ω^{2}}{2} (T^{2} - 2 Tt) + μ_{7}] \\ b_{0} (ω_{y}^{2} - ω^{2}) [1 - \frac{{(ω t)}^{2}}{2}] + B_{0} ω_{y}^{2} = - \frac{N}{2 π m_{y}} [\frac{μ_{21}' a_{0}^{2} ω^{4}}{4} (T^{2} - 2 Tt)^{2} + \frac{μ_{21}' a_{0} ω^{2} c}{2} (T^{2} - 2 Tt) \\ + \frac{μ_{22}' b_{0}^{2} ω^{4}}{4} (T^{2} - 2 Tt)^{2} + \frac{μ_{31}' a_{0} ω^{2}}{4} (T^{2} - 2 Tt) + μ_{31}' c + \frac{μ_{32}' b_{0} ω^{2}}{4} (T^{2} - 2 Tt) \\ + \frac{2 μ_{6}' a_{0} b_{0} ω^{4}}{4} (T^{2} - 2 Tt)^{2} + \frac{μ_{6}' b_{0} c ω^{2}}{2} (T^{2} - 2 Tt) + μ_{7}'] \end{matrix}

Then, a₀, b₀, A₀, B₀ are obtained by solving equation (27).

(27)

{\begin{matrix} - \frac{N ω^{2} T^{2}}{8 π m_{x}} [ω^{2} T^{2} (μ_{21} a_{0}^{2} + μ_{22} b_{0}^{2} + 2 μ_{6} a_{0} b_{0}) + 2 μ_{21} a_{0} c + μ_{31} a_{0} + μ_{32} b_{0}] - \frac{N μ_{7}}{2 π m_{x}} = a_{0} (ω_{x}^{2} - ω^{2}) + A_{0} ω_{x}^{2} \\ μ_{21} a_{0}^{2} ω^{4} T^{3} - μ_{21} a_{0} c ω^{2} T + μ_{22} b_{0}^{2} ω^{4} T^{3} + \frac{μ_{31} a_{0} ω^{2} T}{2} + \frac{μ_{32} b_{0} ω^{2} T}{2} + 2 μ_{6} a_{0} b_{0} ω^{4} T^{3} = 0 \\ - \frac{N ω^{4} T^{2}}{2 π m_{x}} (μ_{21} a_{0}^{2} + μ_{22} b_{0}^{2} + 2 μ_{6} a_{0} b_{0}) = \frac{a_{0} ω^{2} (ω_{x}^{2} - ω^{2})}{2} \\ - \frac{N ω^{2} T^{2}}{8 π m_{y}} [ω^{2} T^{2} (μ_{21}' a_{0}^{2} + μ_{22}' b_{0}^{2} + 2 μ_{6} a_{0} b_{0}) + 2 μ_{21}' a_{0} c + μ_{31}' a_{0} + μ_{32}' b_{0}] - \frac{N μ_{7}'}{2 π m_{y}} = b_{0} (ω_{y}^{2} - ω^{2}) + B_{0} ω_{y}^{2} \\ μ_{21}' a_{0}^{2} ω^{4} T^{3} - μ_{21}' a_{0} c ω^{2} T + μ_{22}' b_{0}^{2} ω^{4} T^{3} + \frac{μ_{31}' a_{0} ω^{2} T}{2} + \frac{μ_{32}' b_{0} ω^{2} T}{2} + 2 μ_{6}' a_{0} b_{0} ω^{4} T^{3} = 0 \\ - \frac{N ω^{4} T^{2}}{2 π m_{y}} (μ_{21}' a_{0}^{2} + μ_{22}' b_{0}^{2} + 2 μ_{6}' a_{0} b_{0}) = \frac{b_{0} ω^{2} (ω_{y}^{2} - ω^{2})}{2} \end{matrix}

Appendix 4

The integral part in equation (15) is solved as equation (28).

(28)

\begin{matrix} [\begin{matrix} \frac{\bar{F_{x}}}{N a_{p}} \\ \frac{\bar{F_{y}}}{N a_{p}} \\ \frac{\bar{F_{z}}}{N a_{p}} \end{matrix}] = [\begin{matrix} - \frac{1}{2} K_{rc, m} Π_{i = 0}^{m / 2} \frac{2 i + 1}{2 i + 2} & \dots & - \frac{1}{4} K_{rc, 0} & - \frac{1}{π} K_{re} \\ \frac{1}{2} K_{tc, m} Π_{i = 0}^{m / 2} \frac{2 i + 1}{2 i + 2} & \dots & \frac{1}{4} K_{tc, 0} & \frac{1}{π} K_{te} \\ \frac{1}{π} K_{ac, m} Π_{i = 0}^{m / 2} \frac{2 i + 2}{2 i + 3} & \dots & \frac{1}{π} K_{ac, 0} & \frac{1}{2} K_{ae} \end{matrix}] [\begin{matrix} c^{m + 1} \\ \dots \\ c \\ 1 \end{matrix}] (m = 2, 4, . . .) \\ [\begin{matrix} \frac{\bar{F_{x}}}{N a_{p}} \\ \frac{\bar{F_{y}}}{N a_{p}} \\ \frac{\bar{F_{z}}}{N a_{p}} \end{matrix}] = [\begin{matrix} - \frac{1}{π} K_{rc, m} Π_{i = 1}^{(m + 1) / 2} \frac{2 i}{2 i + 1} & \dots & - \frac{1}{4} K_{rc, 0} & - \frac{1}{π} K_{re} \\ \frac{1}{π} K_{tc, m} Π_{i = 1}^{(m + 1) / 2} \frac{2 i}{2 i + 1} & \dots & \frac{1}{4} K_{tc, 0} & \frac{1}{π} K_{te} \\ \frac{1}{2} K_{ac, m} Π_{i = 1}^{(m + 1) / 2} \frac{2 i - 1}{2 i} & \dots & \frac{1}{π} K_{ac, 0} & \frac{1}{2} K_{ae} \end{matrix}] [\begin{matrix} c^{m + 1} \\ \dots \\ c \\ 1 \end{matrix}] (m = 1, 3, . . .) \end{matrix}

The feed rates are selected as c₀, c₁, c₂, ……, c_l, and the least square method is utilized to solve the cutting force coefficients.

(29)

A^{T} A a_{x} = A^{T} y_{x}

where $A = [\begin{matrix} 1 & c_{0} & {c_{0}}^{2} & \dots & {c_{0}}^{m + 1} \\ 1 & c_{1} & {c_{1}}^{2} & \dots & {c_{1}}^{m + 1} \\ 1 & c_{2} & {c_{2}}^{2} & \dots & {c_{2}}^{m + 1} \\ \dots & \dots & \dots & \dots & \dots \\ 1 & c_{l} & {c_{l}}^{2} & \dots & {c_{l}}^{m + 1} \end{matrix}]$ , $y_{x} = [\begin{matrix} \frac{\bar{F_{x 0}}}{N a_{p}} & \frac{\bar{F_{x 1}}}{N a_{p}} & \frac{\bar{F_{x 2}}}{N a_{p}} & \dots & \frac{\bar{F_{xl}}}{N a_{p}} \end{matrix}]$ .

The coefficients vector a_x is

(30)

\begin{matrix} a_{x} = {[\begin{matrix} a_{x, 0} & a_{x, 1} & \dots & a_{x, m + 1} \end{matrix}]}^{T} \\ = {[\begin{matrix} - \frac{1}{π} K_{re} & - \frac{1}{4} K_{rc, 0} & \dots & a_{x, m + 1} \end{matrix}]}^{T} \end{matrix}

where

a_{x, m + 1} = {\begin{matrix} - \frac{1}{2} K_{rc, m} Π_{i = 0}^{m / 2} \frac{2 i + 1}{2 i + 2} \begin{matrix} , & m = 0, 2, . . . \end{matrix} \\ - \frac{1}{π} K_{rc, m} Π_{i = 1}^{(m + 1) / 2} \frac{2 i}{2 i + 1} \begin{matrix} , & m = 1, 3, . . . \end{matrix} \end{matrix}

The coefficients vector a_y and a_z can be formulated similarly.

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: This work is supported by the National Natural Science Foundation of China (grant no. 51675032 and 91960108) and Fundamental Research Funds for the Central Universities (grant no. YWF-19-BJ-J-72).

ORCID iD

Yiqing Yang

References

Wang

. Error compensation in the five-axis flank milling of thin-walled workpieces. Proc IMechE, Part B: J Engineering Manufacture 2019; 233(4): 1224–1234.

Liu

Dang

, et al. Parametric study on the tool inclination angle for side milling thin-walled workpiece edges based on FE simulation. Proc IMechE, Part B: J Engineering Manufacture 2020; 234(3): 439–452.

Zhang

Luo

, et al. Dynamic modeling and stability prediction in milling process of thin-walled workpiece with multiple structural modes. Proc IMechE, Part B: J Engineering Manufacture Epub ahead of print July 2020. DOI: 10.1177/0954405420933710.

Lamikiz

De Lacalle

Sanchez

, et al. Calculation of the specific cutting coefficients and geometrical aspects in sculptured surface machining. Mach Sci Technol 2005; 9(3): 411–436.

Muñoa

Zatarain

Bediaga

, et al. Stability study of the milling process using an exponential force model in frequency domain. In: CIRP-2nd International HPC Conference, June 2006, Vancouver, Canada.

Guzeev

Pimenov

. Cutting force in face milling with tool wear. Russ Engin Res 2011; 31(10): 989–993.

Moradi

Movahhedy

Vossoughi

. Dynamics of regenerative chatter and internal resonance in milling process with structural and cutting force nonlinearities. J Sound Vib 2012; 331(16): 3844–3865.

Yang

Zhang

, et al. Chatter prediction for the peripheral milling of thin-walled workpieces with curved surfaces. Int J Mach Tool Manu 2016; 109: 36–48.

Karimi

Nategh

. Contour maps for developing optimal toolpath and workpiece setup in hexapod machine tools by considering the kinematics nonlinearity. Proc IMechE, Part B: J Engineering Manufacture 2016; 230(9): 1572–1583.

10.

Budak

Tunc

. Identification and modeling of process damping in turning and milling using a new approach. CIRP Ann-Manuf Techn 2010; 59(1): 403–408.

11.

Feng

Wan

Gao

, et al. Mechanism of process damping in milling of thin-walled workpiece. Int J Mach Tool Manu 2018; 134: 1–19.

12.

Dun

Zhu

Wang

, et al. Multi-modal method for chatter stability prediction and control in milling of thin-walled workpiece. Appl Math Model 2020; 80: 602–624.

13.

Yang

Zhang

, et al. An efficient decomposition-condensation method for chatter prediction in milling large-scale thin-walled structures. Mech Syst Signal Pr 2019; 121: 58–76.

14.

Sun

Guo

. Chatter prediction utilizing stability lobes with process damping in finish milling of titanium alloy thin-walled workpiece. Int J Adv Manuf Tech 2017; 89(1–12): 2663–2674.

15.

Song

Shi

Liu

, et al. A time-space discretization method in milling stability prediction of thin-walled component. Int J Adv Manuf Tech 2017; 89(9–12): 2675–2689.

16.

Sellmeier

Denkena

. Stable islands in the stability chart of milling processes due to unequal tooth pitch. Int J Mach Tool Manu 2011; 51(2): 152–164.

17.

Afazov

Ratchev

Segal

, et al. Chatter modelling in micro-milling by considering process nonlinearities. Int J Mach Tool Manu 2012; 56: 28–38.

18.

Zhou

Deng

Peng

, et al. Semi-analytic modelling of cutting forces in micro ball-end milling of NAK80 steel with wear-varying cutting edge and associated nonlinear process characteristics. Int J Mech Sci 2020; 169: 105343.

19.

Moradi

Vossoughi

Movahhedy

, et al. Forced vibration analysis of the milling process with structural nonlinearity, internal resonance, tool wear and process damping effects. Int J Nonlin Mech 2013; 54: 22–34.

20.

Olvera

Urbikain

Eliaszuniga

, et al. Improving stability prediction in peripheral milling of Al7075T6. Appl Sci 2018; 8(8): 1316–2018.

21.

Urbikain

Campa

Zulaika

, et al. Preventing chatter vibrations in heavy-duty turning operations in large horizontal lathes. J Sound Vib 2015; 340: 317–330.

22.

Jin

Altintas

. Chatter stability model of micro-milling with process damping. J Manuf Sci E-T ASME 2013; 135(3): 031011.