Free vibration and chatter stability of a rotating thin-walled composite bar

Abstract

The free vibration and chatter stability of a rotating thin-walled composite bar under the action of regenerative milling force were investigated in this article. The free vibration equations of the rotating thin-walled composite bar were presented and solved by the Galerkin method according to the derived Lagrange equation for the rotating thin-walled composite bar, and the influences of the taper ratio, ply angle, and mode number on the first natural frequency were analyzed. In order to study the chatter stability of the rotating thin-walled composite bar, the time-delay motion equation of the system was derived with consideration of the effects of regenerative milling force and internal/external damping. The two-mode Galerkin method was used to discretize the time-delay free vibration equations, and the semi-discrete method in the time domain was employed to predict the stability lobes. Finally, the correctness of the method was validated, and the stability analysis of the rotating thin-walled composite bar was conducted. Emphatically, the influences of the ply angle, taper ratio, and internal/external damping on the chatter stability of the rotating thin-walled composite bar were analyzed.

Keywords

Chatter cutter bar composite thin-walled variable cross-section

The most common method of mechanical cutting is to install the cutter head on the cutter bar. The rotating cutter bar can drive the cutter head to machine the fixed workpiece and thus complete the cutting process. In certain specific machining environments, such as boring deep holes and milling deep grooves, the cutter bar shall have a slender structure. Consequently, its stiffness is low and the cutter bar is prone to chatter under the action of periodically changing cutting force.^1–3 In practice, the techniques to improve the stability of the slender cutter bar include active damping,^4–6 passive damping,^7–9 and spindle speed variation.^10,11 However, active and passive damping control technologies both need to add devices and structures to the cutting system, which increases the volume of the processing system and makes the structure complicated. The advantage of the spindle speed variable method lies in its simple structure and wide suitability for frequency band adjustment. Unfortunately, the precision of the machined surface by this method is significantly affected by the changing speed.

In order to improve the stability of machine tool structures including slender cutter bars, some researchers have investigated the possibility of replacing metal materials with composite materials in the design of boring bars because of low weight, high stiffness, and high damping of composite materials.^12–16 Lee and Suh¹² designed a composite boring bar and proved that the composite boring bar had better stability than traditional metal boring bars. Nagano et al.¹³ studied the dynamic characteristics of composite boring bars with different shapes of metal core. Kim et al.¹⁶ derived the time-delay vibration differential equation of variable cross-section milling cutter bars based on the Timoshenko beam theory and drew the diagram of stability lobes by the frequency-domain analysis method. Especially, the influences of the taper ratio (TR), damping, and material properties on the stability of the bars were also analyzed. However, due to the high rotating speed of boring bars, it is also necessary to consider the gyroscopic effects.¹⁷

From the existing literature, little work has been done on the modeling of rotating composite bars and their chatter stability. Thus, it is of great importance to investigate the influences of various properties on the dynamic behaviors of composite bars in order to gain further understanding of composite bars in the metal cutting process.

In this article, the composite thin-walled beam theory proposed by Berdichevsky et al.¹⁸ and Badir¹⁹ was used to model composite boring bars. Based on the displacement function of Berdichevsky and the Lagrange equation of continuous system, free vibration differential equations of a rotating thin-walled composite bar (RTWCB) with bending–bending coupling were presented. Furthermore, dynamic equations of the RTWCB were developed by taking internal damping, external damping, and the cutting force with time-delay effect into account. Then, on the basis of analyzing the free vibration, the influences of the ply angle, TR, and internal/external damping on the chatter stability were analyzed. Numerical results showed that the ply angle and TR had obvious influences on the chatter stability, while the internal damping affected the stability only at high speeds. The correctness of the stability analysis method could be verified by comparison with the frequency-domain method as well as the calculation results of the time-domain displacement response.

The rest of this article is organized as follows. In section “Dynamic mechanical models of the RTWCB,” the dynamic mechanical models of the RTWCB are established. In section “Governing equations,” the Lagrange equation for the RTWCB is derived to deduce the governing equation of the RTWCB. Then, in section “Solution,” the Galerkin method is used to transform the vibration equation into an ordinary differential equation with respect to time, and the semi-discrete method in the time domain is used to predict the stability lobes. In section “Numerical results,” the free vibration characteristics of the RTWCB rods are discussed, and the influences of the ply angle, TR, and internal/external damping on the chatter stability are analyzed based on the correctness of the stability analysis method. Finally, a summary and some conclusions are provided in section “Summary and conclusion.”

Dynamic mechanical models of the RTWCB

Figure 1 shows the cutter bar model and coordinate systems. The RTWCB rotates around its longitudinal axis with a constant angular speed W. The coordinate system (X, Y, Z) stands for a fixed inertial coordinate system, and its unit vector is ( I , J , K ); (x, y, z) stands for a moving coordinate system, which is rotating with the cutter bar, and its unit vector is ( i , j , k ); (x, s, x) represents the local coordinate system, where s is along the midline of the contour tangent direction and x is along the normal direction; h represents the wall thickness; r₁ and r₂ are the contour line radii of the root and the free cross-section, respectively; L is the axial length; and q is the ply angle of the composite fibers. Here, the TR of the cutter bar is defined as

TR = 100 (\frac{r_{2} - r_{1}}{L})

(1)

Figure 1.

Cutter bar model and coordinate systems.

It is assumed that (1) the cutter bar is a slender thin-walled elastic cylinder, namely, d ≪ L, h ≪ d, h ≪ r₁, and h ≪ r₂, where d is the maximum outer diameter of the cutter bar cross-section and (2) the transverse shear effect of the thin-walled structure can be neglected.

Governing equations

The displacement field of an anisotropic closed-section thin-walled structure derived by the variational approach is very applicable to the dynamic analysis of the RTWCB. It can be expressed as follows^18,19

\begin{matrix} u_{1} (x, s) = U_{1} (x) - y {U'}_{2} (x) - z {U'}_{3} (x) + g (s, x) \\ u_{2} (x, s) = U_{2} (x) - z φ (x) \\ u_{3} (x, s) = U_{3} (x) + y φ (x) \end{matrix}

(2)

where u₁, u₂, and u₃ denote the displacements along the X-, Y-, and Z-axes of any point on the RTWCB, respectively; U₁(x), U₂(x), and U₃(x) denote the average displacements along the X-, Y-, and Z-axes of the x cross-section, respectively; $φ (x)$ is the twist angle of the same section; the function g(s, x) added to the classical displacement field represents the out-of-plane warping of the cross-section, which can be determined from the continuity of the displacement field and the vanishing of the resultant hoop stress and constant shear flow; and ${U'}_{2} (x)$ and ${U'}_{3} (x)$ denote the differential forms of U₂(x) and U₃(x) with respect to x, respectively.

Based on the displacement functions in equation (2), the strain energy density (the strain energy per unit length) of the RTWCB can be calculated as follows^20,21

\begin{matrix} U = \frac{1}{2} C_{11} {U^{'}}_{1}^{2} + C_{12} ϕ^{'} {U^{'}}_{1} + C_{13} {U^{″}}_{3} {U^{'}}_{1} + C_{14} {U^{″}}_{2} {U^{'}}_{1} + \frac{1}{2} C_{22} {φ^{'}}^{2} + C_{23} {U^{″}}_{3} φ^{'} \\ + C_{24} {U^{″}}_{2} φ^{'} + \frac{1}{2} C_{33} {U^{″}}_{3}^{2} + C_{43} {U^{″}}_{2} {U^{‴}}_{3} + \frac{1}{2} C_{44} {U^{″}}_{2}^{2} \end{matrix}

(3)

where $C_{ij} (i, j = 1, 4)$ can be expressed in terms of the cross-section geometry and material properties.²¹

The kinetic energy on the unit length of the cutter bar can be expressed as²²

\begin{matrix} T = \frac{1}{2} [I_{m} ({\overset{\cdot}{U}}_{1}^{2} + {\overset{\cdot}{U}}_{2}^{2} + {\overset{\cdot}{U}}_{3}^{2}) + I_{d} ({\overset{\cdot}{U}}^{'}_{2}^{2} + {\overset{\cdot}{U}}^{'}_{3}^{2}) - 2 Ω I_{P} {U^{'}}_{2} {\overset{\cdot}{U}}^{'}_{3} + I_{P} {\overset{\cdot}{φ}}^{2} \\ + 2 Ω I_{P} \overset{\cdot}{φ} + Ω^{2} I_{P} + Ω^{2} I_{d} ({U^{'}}_{2}^{2} + {U^{'}}_{3}^{2})] \end{matrix}

(4)

where

\begin{matrix} I_{m} = π \sum_{i = 1}^{n} ρ_{i} (r_{i + 1}^{2} - {r_{i}}^{2}) \\ I_{d} = \frac{π}{4} \sum_{i = 1}^{n} ρ_{i} (r_{i + 1}^{4} - {r_{i}}^{4}) \\ I_{P} = \frac{π}{2} \sum_{i = 1}^{n} ρ_{i} (r_{i + 1}^{4} - {r_{i}}^{4}) \end{matrix}

(5)

where n is the number of layers in the laminate and $r_{i}$ and $r_{i + 1}$ are the inner and outer radii of the ith layer, respectively.

Because $Ω^{2} I_{P} >> Ω^{2} I_{d} ({U'}_{2}^{2} + {U'}_{3}^{2})$ , the last term in equation (4) can be ignored. Then

\begin{matrix} T = \frac{1}{2} [I_{m} ({\overset{\cdot}{U}}_{1}^{2} + {\overset{\cdot}{U}}_{2}^{2} + {\overset{\cdot}{U}}_{3}^{2}) + I_{d} ({\overset{\cdot}{U}'_{2}}^{2} + {\overset{\cdot}{U}'_{3}}^{2}) \\ - 2 Ω I_{P} {U'}_{2} \overset{\cdot}{U}'_{3} + I_{P} {\overset{\cdot}{φ}}^{2} + 2 Ω I_{P} \overset{\cdot}{φ} + Ω^{2} I_{P}] \end{matrix}

(6)

Based on the expressions of the kinetic energy density (equation (6)) and the strain energy density (equation (3)) of the cutter bar, the Lagrange equation of the continuum can be established as²³

\begin{matrix} - \frac{d}{d t} \frac{\partial ζ}{\partial {\overset{\cdot}{U}}_{1}} - \frac{d}{d x} \frac{\partial ζ}{\partial {U'}_{1}} = 0 \\ - \frac{d}{d t} \frac{\partial ζ}{\partial {\overset{\cdot}{U}}_{2}} + \frac{d^{2}}{d x^{2}} \frac{\partial ζ}{\partial {U ″}_{2}} - \frac{d}{d x} \frac{\partial ζ}{\partial {U'}_{2}} + \frac{d}{d x} \frac{d}{d t} \frac{\partial ζ}{\partial \overset{\cdot}{U}'_{2}} = 0 \\ - \frac{d}{d t} \frac{\partial ζ}{\partial {\overset{\cdot}{U}}_{3}} + \frac{d^{2}}{d x^{2}} \frac{\partial ζ}{\partial {U ″}_{3}} + \frac{d}{d x} \frac{d}{d t} \frac{\partial ζ}{\partial \overset{\cdot}{U}'_{3}} = 0 \\ - \frac{d}{d t} \frac{\partial ζ}{\partial \overset{\cdot}{φ}} - \frac{d}{d x} \frac{\partial ζ}{\partial φ'} = 0 \end{matrix}

(7)

where

ζ = T - U

(8)

The free vibration equations can be expressed in terms of kinematic variables by substituting equation (8) into equation (7) as follows

\begin{matrix} I_{m} {\overset{\cdot\cdot}{U}}_{1} - C'_{11} {U'}_{1} - C_{11} {U ″}_{1} - C'_{12} φ' - C_{12} φ ″ \\ - C'_{13} {U ″}_{3} - C_{13} U ″'_{3} - C'_{14} {U ″}_{2} - C_{14} U ″'_{2} = 0 \end{matrix}

(9)

\begin{matrix} I_{P} \overset{\cdot\cdot}{φ} - C'_{12} {U'}_{1} - C_{12} {U ″}_{1} - C'_{22} φ' - C_{22} φ ″ \\ - C'_{23} {U ″}_{3} - C_{23} U ″'_{3} - C'_{24} {U ″}_{2} - C_{24} U ″'_{2} = 0 \end{matrix}

(10)

\begin{matrix} m_{c} {\overset{\cdot\cdot}{U}}_{3} - I_{d} {\overset{\cdot\cdot}{U} ″}_{3} - I'_{d} \overset{\cdot}{U}'_{3} - Ω I'_{P} \overset{\cdot}{U}'_{2} - Ω I_{P} {\overset{\cdot}{U} ″}_{2} \\ + {C ″}_{13} {U'}_{1} + 2 C'_{13} {U ″}_{1} + C_{13} U ″'_{1} + {C ″}_{23} φ' \\ + 2 C'_{23} φ ″ + C_{23} φ ″' + {C ″}_{33} {U ″}_{3} + 2 {C ″}_{33} {U ″}_{3} \\ + C_{33} U' ″'_{3} + {C ″}_{43} {U ″}_{2} + 2 {C ″}_{43} {U ″}_{2} + C_{43} U' ″'_{2} = 0 \end{matrix}

(11)

\begin{matrix} m_{c} {\overset{\cdot\cdot}{U}}_{2} - I_{d} {\overset{\cdot\cdot}{U} ″}_{2} - I'_{d} \overset{\cdot}{U}'_{2} + Ω I'_{P} \overset{\cdot}{U}'_{3} + Ω I_{P} {\overset{\cdot}{U} ″}_{3} \\ + {C ″}_{14} {U'}_{1} + 2 C'_{14} {U ″}_{1} + C_{14} U ″'_{1} + {C ″}_{24} φ' + \\ 2 C'_{24} φ ″ + C_{24} φ ″' + {C ″}_{43} {U ″}_{3} + 2 {C ″}_{43} {U ″}_{3} \\ + C_{43} U' ″'_{3} + {C ″}_{44} {U ″}_{2} + 2 {C ″}_{44} {U ″}_{2} + C_{44} U' ″'_{2} = 0 \end{matrix}

(12)

Since the RTWCB has a circular cross-section and circumferential symmetric stiffness ( $θ (y) = θ (- y)$ , $θ (z) = θ (- z)$ ), it is available that C₁₃ = C₁₄ = C₂₃ =C₂₄ = C₃₄ = 0. Then equations (9)–(12) can be simplified as

I_{m} {\overset{\cdot\cdot}{U}}_{1} - C'_{11} {U'}_{1} - C_{11} {U ″}_{1} - C'_{12} φ' - C_{12} φ ″ = 0

(13)

I_{P} \overset{\cdot\cdot}{φ} - C'_{12} {U'}_{1} - C_{12} {U ″}_{1} - C'_{22} φ' - C_{22} {φ ″}_{2} = 0

(14)

\begin{matrix} m_{c} {\overset{\cdot\cdot}{U}}_{3} - I_{d} {\overset{\cdot\cdot}{U} ″}_{3} - I'_{d} \overset{\cdot}{U}'_{3} - Ω I'_{P} \overset{\cdot}{U}'_{2} - Ω I_{P} {\overset{\cdot}{U} ″}_{2} \\ + {C ″}_{33} {U ″}_{3} + 2 {C ″}_{33} {U ″}_{3} + C_{33} {U ″ ″}_{3} = 0 \end{matrix}

(15)

\begin{matrix} m_{c} {\overset{\cdot\cdot}{U}}_{2} - I_{d} {\overset{\cdot\cdot}{U} ″}_{2} - I'_{d} \overset{\cdot}{U}'_{2} + Ω I'_{P} \overset{\cdot}{U}'_{3} + Ω I_{P} {\overset{\cdot}{U} ″}_{3} + {C ″}_{44} {U ″}_{2} \\ + 2 {C ″}_{44} {U ″}_{2} + C_{44} {U ″ ″}_{2} = 0 \end{matrix}

(16)

Considering equations (12)–(15), for an RTWCB with circumferentially uniform stiffness (CUS) configuration, the bending motion can be decoupled from the axial and torsional motions. Here, the bending motion is analyzed separately.

The damping terms play an important role in the suppression/enhancement of instability, but they are not included in equations (15) and (16). When the RTWCB works as a milling cutter bar, it suffers periodic milling force. Taking these terms into account, equations (15) and (16) can be transformed as

\begin{matrix} m_{c} {\overset{\cdot\cdot}{U}}_{2} - I_{d} {\overset{\cdot\cdot}{U} ″}_{2} - I'_{d} \overset{\cdot}{U}'_{2} + Ω I'_{P} \overset{\cdot}{U}'_{3} + Ω I_{P} {\overset{\cdot}{U} ″}_{3} + (c_{e} + c_{i}) {\overset{\cdot}{U}}_{2} \\ + c_{i} Ω U_{3} + {C ″}_{44} {U ″}_{2} + 2 {C ″}_{44} {U ″}_{2} + C_{44} {U ″ ″}_{2} = f_{y} \end{matrix}

(17)

\begin{matrix} m_{c} {\overset{\cdot\cdot}{U}}_{3} - I_{d} {\overset{\cdot\cdot}{U} ″}_{3} - I'_{d} \overset{\cdot}{U}'_{3} - Ω I'_{P} \overset{\cdot}{U}'_{2} - Ω I_{P} {\overset{\cdot}{U} ″}_{2} \\ + (c_{e} + c_{i}) {\overset{\cdot}{U}}_{3} - c_{i} Ω U_{2} + {C ″}_{33} {U ″}_{3} + 2 {C ″}_{33} {U ″}_{3} \\ + C_{33} {U ″ ″}_{3} = f_{z} \end{matrix}

(18)

where $f_{y}$ and $f_{z}$ represent the applied forces per unit length in the y- and z-directions, respectively; $c_{e}$ is an external damping coefficient per unit length of the cutter bar and is proportional to the speed in the fixed coordinate system; and $c_{i}$ is an internal damping coefficient per unit length of the cutter bar and is proportional to the speed in the rotating coordinate system.²⁴ Here, the damping values are assigned based on the concepts from an oscillator with a single degree of freedom, namely

c_{e} + c_{i} = \frac{2 ς ω_{n} M_{33}}{\int_{0}^{L} {(ψ_{1})}^{2} d z}

(19)

where $ς$ is the damping ratio, which is assumed to be 0.015 for the problem at hand; $ω_{n}$ is the lowest natural bending frequency at $Ω = 0$ ; and $ψ_{1}$ and $M_{33}$ are the one-mode bending Galerkin function and the approximate bending mass, respectively, and their detailed expressions will be given in the following section.

Solution

Galerkin approximation

The Galerkin method is an effective and inexpensive method for converting a continuous differential equation to a discrete problem. It can provide a powerful numerical solution to the dynamic analysis of a continuous system. In this study, the Galerkin method is used to find the approximate solution of equations (17) and (18). By employing this method, the governing equations (17) and (18) are approximated by a system of ordinary differential equations, namely

U_{2} = \sum_{j = 1}^{N} W_{j} (t) ψ_{j} (x), U_{3} = \sum_{j = 1}^{N} V_{j} (t) ψ_{j} (x)

(20)

The Galerkin function $ψ_{j} (x)$ for the variables U₂ and U₃ are taken to be the mode shapes of a non-rotating, non-coupling, and uniform cantilever Euler–Bernoulli beam as follows

\begin{matrix} \cos (β_{j}) \cosh (β_{j}) = - 1, λ_{j} = - \frac{\cos β_{j} + \cosh β_{j}}{\sin β_{j} + \sinh β_{j}}, \\ j = 1, 2, \dots, N \\ ψ_{j} = \cos (\frac{β_{j} x}{L}) - \cosh (\frac{β_{j} x}{L}) \\ + λ_{j} (\sin (\frac{β_{j} x}{L}) - \sinh (\frac{β_{j} x}{L})) \end{matrix}

(21)

The Galerkin method can yield a set of equations in the form of

[M] \overset{\cdot\cdot}{X} + [G] \overset{\cdot}{X} + KX = {F}

(22)

where

\begin{matrix} M = {\begin{matrix} M_{33} & 0 \\ 0 & M_{44} \end{matrix}}, G = {\begin{matrix} G_{33} & G_{34} \\ G_{43} & G_{44} \end{matrix}}, \\ K = {\begin{matrix} K_{33} & K_{34} \\ K_{43} & K_{44} \end{matrix}}, {F} = {\begin{matrix} {f_{yj}} \\ {f_{zj}} \end{matrix}} \end{matrix}

(23)

X = {(\begin{matrix} \sum_{j = 1}^{N} V_{j} (t) & \sum_{j = 1}^{N} W_{j} (t) \end{matrix})}^{T}

The elements in the above matrix are listed as follows

\begin{matrix} K_{33 ij} = \int_{0}^{L} ({C ″}_{33} {ψ ″}_{j} ψ_{i} + 2 C'_{33} ψ ″'_{j} ψ_{i} + C_{33} ψ' ″'_{j} ψ_{i}) d x, \\ K_{44 ij} = \int_{0}^{L} ({C ″}_{44} {ψ ″}_{j} ψ_{i} + 2 C'_{44} ψ ″'_{j} ψ_{i} + C_{44} ψ' ″'_{j} ψ_{i}) d x \\ K_{34 ij} = - K_{43 ij} = \int_{0}^{L} - c_{i} Ω ψ_{i} ψ_{j} d x, M_{33 ij} = M_{44 ij} \\ = \int_{0}^{L} (I_{m} ψ_{i} ψ_{j} - I_{d} ψ_{i} {ψ ″}_{j} + I'_{d} ψ_{i} ψ'_{j}) d x \end{matrix}

(24)

\begin{matrix} G_{34 ij} = - G_{43 ij} \\ = \int_{0}^{L} (- I_{P} Ω ψ_{i} {ψ ″}_{j} - I'_{P} Ω ψ_{i} ψ'_{j}) d x, G_{34 ij} = G_{43 ij} \\ = \int_{0}^{L} (c_{e} + c_{i}) ψ_{i} ψ_{j} d x \\ f_{yj} = \int_{0}^{L} f_{y} ψ_{j} d x, f_{zj} = \int_{0}^{L} f_{z} ψ_{j} d x \end{matrix}

When the RTWCB works as a milling cutter bar, it suffers milling force. The milling forces are treated as concentrated loads at its free end. Then $f_{y}$ and $f_{z}$ in equation (24) can be expressed as

f_{y} = F_{y} (t) δ (z - L), f_{z} = F_{z} (t) δ (z - L)

(25)

where $δ$ denotes a Dirac delta function. Considering a half immersion, the up-milling operation with a four-fluted end-milling operation is adopted, and the deflection-dependent cutting force model is employed. The tangential and radial forces based on the work of Smith and Tlusty²⁵ are given by

\begin{matrix} F_{y} = F_{y 0} - F_{y 1} {U_{2} (t) - U_{2} (t - T)} \\ - F_{y 2} {U_{3} (t) - U_{3} (t - T)} \end{matrix}

(26)

\begin{matrix} F_{z} = F_{z 0} - F_{z 1} {U_{2} (t) - U_{2} (t - T)} \\ - F_{z 2} {U_{3} (t) - U_{3} (t - T)} \end{matrix}

(27)

where

\begin{matrix} F_{y 0} (t) = \frac{K_{s} d f_{t}}{2} {- 1 + \frac{1}{\sin β_{0}} \sin (2 θ_{t} + β_{0})}, \\ F_{y 1} (t) = \frac{K_{s} d}{2} {\frac{1}{\tan β_{0}} + \frac{1}{\sin β_{0}} \cos (2 θ_{t} + β_{0})} \\ F_{y 2} (t) = \frac{K_{s} d}{2} {- 1 + \frac{1}{\sin β_{0}} \sin (2 θ_{t} + β_{0})}, \\ F_{z 0} (t) = \frac{K_{s} d f_{t}}{2} {\frac{1}{\tan β_{0}} - \frac{1}{\sin β_{0}} \cos (2 θ_{t} + β_{0})} \\ F_{z 1} (t) = \frac{K_{s} d}{2} {1 + \frac{1}{\sin β_{0}} \sin (2 θ_{t} + β_{0})}, \\ F_{z 2} (t) = \frac{K_{s} d}{2} {\frac{1}{\tan β_{0}} - \frac{1}{\sin β_{0}} \cos (2 θ_{t} + β_{0})} \\ θ_{t} = remainder of (\frac{Ω t}{0.5 π}) \\ = {\begin{matrix} Ω t if 0 \leq Ω t < 0.5 π \\ Ω t - 0.5 π if 0.5 π \leq Ω t < π \\ ⋮ \end{matrix} \\ T = \frac{2 π}{N_{t} Ω} \end{matrix}

(28)

where $K_{s}$ is the specific cutting force per unit chip area and its value is $K_{s} = 551.6 MPa$ ;¹ $β_{0}$ is the constant angle between the cut surface normal and the direction of cutting force, and its value is 73.3°;² $θ_{t}$ represents the azimuth at which the cutting edge is located; d represents the cutting depth; $f_{t}$ is the feed per tooth; and $T = 2 π / (N_{t} Ω)$ denotes the tooth period, where the number of teeth is $N_{t} = 4$ .

By combining equations (22), (26), and (27), it can be obtained that

\begin{matrix} [M] \overset{\cdot\cdot}{X} (t) + [G] \overset{\cdot}{X} (t) + [K] X (t) + ([K_{f}] \\ + [K_{P} (t)]) {X (t) - X (t - T)} = {F} \end{matrix}

(29)

where

\begin{matrix} [K_{f}] = \frac{K_{s} d}{2} [\begin{matrix} \frac{1}{\tan β_{0}} [z_{ij}] & - [z_{ij}] \\ [z_{ij}] & \frac{1}{\tan β_{0}} [z_{ij}] \end{matrix}], [K_{p} (t)] = \frac{K_{s} d}{2 \sin β_{0}} \\ [\begin{matrix} \cos (2 θ_{t} + β_{0}) [z_{ij}] & \sin (2 θ_{t} + β_{0}) [z_{ij}] \\ \sin (2 θ_{t} + β_{0}) [z_{ij}] & - \cos (2 θ_{t} + β_{0}) [z_{ij}] \end{matrix}] \\ {f_{0 j}^{y}} = F_{y 0} {ψ_{j} (L)}, {f_{0 j}^{z}} = F_{z 0} {ψ_{j} (L)}, \\ z_{ij} = ψ_{i} (L) ψ_{j} (L) \end{matrix}

(30)

Considering the effect of the cutting force, equation (29) contains one traditional stiffness matrix K , as well as two additional matrices $K_{f}$ and $K_{P}$ caused by time delay.

Semi-discrete stability solution

Ignoring the constant force and transporting the time delay to the right-hand side, equation (29) can be expressed as

[M] \overset{\cdot\cdot}{X} (t) + [G] \overset{\cdot}{X} (t) + [K + K_{F}] X (t) = [K_{F}] X (t - T)

(31)

where

\begin{matrix} [K_{F}] = [\begin{matrix} K_{F 11} d & K_{F 12} d \\ K_{F 21} d & K_{F 22} d \end{matrix}] \\ K_{F 11} = \frac{K_{s}}{2} z (\frac{1}{\tan β_{0}} + \frac{\cos (2 θ_{t} + β_{0})}{\sin β_{0}}), \\ K_{F 12} = \frac{K_{s}}{2} z (- 1 + \frac{\sin (2 θ_{t} + β_{0})}{\sin β_{0}}) \end{matrix}

(32)

\begin{matrix} K_{F 21} = \frac{K_{s}}{2} z (1 + \frac{\sin (2 θ_{t} + β_{0})}{\sin β_{0}}), \\ K_{F 22} = \frac{K_{s}}{2} z (\frac{1}{\tan β_{0}} - \frac{\cos (2 θ_{t} + β_{0})}{\sin β_{0}}) \end{matrix}

From equation (31), it can be obtained that

\begin{matrix} \overset{\cdot\cdot}{X} (t) + [M]^{- 1} [G] \overset{\cdot}{X} (t) + [M]^{- 1} [K + K_{F}] X (t) \\ = [M]^{- 1} [K_{F}] X (t - T) \end{matrix}

(33)

Using the semi-discrete method, in the ith time interval, equation (33) can be expressed as²⁶

\begin{matrix} \overset{\cdot\cdot}{X} (t) + [M]^{- 1} [G] \overset{\cdot}{X} (t) + [M]^{- 1} [K + K_{F}] X (t) \\ = [M]^{- 1} [K_{F}] X_{τ, i} \end{matrix}

(34)

By Cauchy transformation, equation (38) can be written in the form of

\overset{\cdot}{q} (t) = A_{i} q (t) + w_{a} B_{i} q_{i - m + 1} + w_{b} B q_{i - m + 1}

(35)

where

\begin{matrix} A_{i} = (\begin{matrix} 0 & I \\ - M^{- 1} (K + K_{F}) & - M^{- 1} G \end{matrix}) \\ B_{i} = (\begin{matrix} 0 & 0 \\ M^{- 1} K_{F} & 0 \end{matrix}), q (t) = (\begin{matrix} X (t) \\ \overset{\cdot}{X} (t) \end{matrix}), \\ q_{j} = q (t_{j}) = (\begin{matrix} X (t_{j}) \\ \overset{\cdot}{X} (t_{j}) \end{matrix}) = (\begin{matrix} X_{j} \\ {\overset{\cdot}{X}}_{j} \end{matrix}) \end{matrix}

(36)

and $w_{a}$ and $w_{b}$ are the weighting factors. For the milling system, the time delay is equal to the cutting cycle, so $w_{a} = w_{b} = 1 / 2$ . If the initial condition is $q (t_{i}) = q_{i}$ , $q_{i + 1}$ can be expressed as

q_{i + 1} = P_{i} q_{i} + w_{a} R_{i} q_{i - m + 1} + w_{b} R_{i} q_{i - m}

(37)

where

\begin{matrix} P = \exp (A_{i} Δ t) \\ R = (\exp (A_{i} Δ t) - I) A_{i}^{- 1} B_{i} \end{matrix}

(38)

and m is an approximation parameter regarding the time delay.

It can be found that $X_{i}$ and ${\overset{\cdot}{X}}_{i}$ depend on $X_{i - m + 1}$ and $X_{i - m}$ instead of ${\overset{\cdot}{X}}_{i - m + 1}$ and ${\overset{\cdot}{X}}_{i - m}$ . Therefore, the elements in the column n_j of the matrices $B_{i}$ and $R_{i}$ are both 0. Thus, let

z_{i} = col (\begin{matrix} X_{i} & {\overset{\cdot}{X}}_{i} & X_{i - 1} & \dots & X_{i - m} \end{matrix})

(39)

After a series of discretization, the final value of the system is

z_{i + 1} = D_{i} z_{i}

(40)

where

D_{i} = (\begin{matrix} P_{i, 11} & P_{i, 12} & P_{i, 13} & P_{i, 14} & 0 & \dots & 0 & w_{a} R_{i, 11} & w_{a} R_{i, 12} & w_{b} R_{i, 11} & w_{b} R_{i, 12} \\ P_{i, 21} & P_{i, 22} & P_{i, 23} & P_{i, 24} & 0 & \dots & 0 & w_{a} R_{i, 21} & w_{a} R_{i, 22} & w_{b} R_{i, 21} & w_{b} R_{i, 22} \\ P_{i, 31} & P_{i, 32} & P_{i, 33} & P_{i, 34} & 0 & \dots & 0 & w_{a} R_{i, 31} & w_{a} R_{i, 32} & w_{b} R_{i, 31} & w_{b} R_{i, 32} \\ P_{i, 41} & P_{i, 42} & P_{i, 43} & P_{i, 44} & 0 & \dots & 0 & w_{a} R_{i, 41} & w_{a} R_{i, 42} & w_{b} R_{i, 41} & w_{b} R_{i, 42} \\ I & 0 & 0 & 0 & 0 & \dots & 0 & 0 & 0 & 0 & 0 \\ 0 & I & 0 & 0 & 0 & \dots & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & I & \dots & 0 & 0 & 0 & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & 0 & 0 & \dots & I & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \dots & 0 & I & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \dots & 0 & 0 & I & 0 & 0 \end{matrix})

(41)

where $P_{i, hj}$ and $R_{i, hj}$ are the h row and j column elements of the matrices $P_{i}$ and $R_{i}$ , respectively. Thus, the transfer matrix $Φ$ of (4m + 8) dimensions can be obtained as

Φ = [\begin{matrix} D_{k} & D_{k - 1} & \dots & D_{1} D_{0} \end{matrix}]

(42)

If the eigenvalues of the transfer matrix $Φ$ are less than 1, the system is stable; otherwise the system is unstable. When the eigenvalues are equal to 1, the system is in the critical stable state.

Numerical results

Free vibration characteristics

According to equation (22) and ignoring the cutting force, the first-mode natural frequency of the free vibration of the RTWCB can be obtained using MATLAB software.

Figure 2 shows the variation of the first-mode natural frequency of the RTWCB with ply angle (TR = 2, L = 0.4 m, $r_{1} + r_{2} = 0.03175 m$ , ply thickness = 0. 386806 mm, and layer number = 6). The mechanical properties of the composite materials are listed in Table 1. The analysis results from equation (21) are compared with the simulation results obtained using the ANSYS finite element method. As illustrated in Figure 2, the analysis results agree well with the simulation results, and the ply angle has a significant influence on frequencies up to 40°.

Figure 2.

Variation of the first-mode natural frequency of the RTWCB with ply angle.

Table 1.

Properties of the boron–epoxy composite material and steel.

	Boron–epoxy	Steel
E₁₁ (GPa)	211	207
E₂₂ (GPa)	24.1	207
G₁₂ = G₁₃ (GPa)	6.9	80
G₂₃ (GPa)	6.9	80
v ₁₂	0.36	0.2975
ρ(kg/m³)	1967	7700

Figure 3 illustrates the variation of the first-mode natural frequency of the RTWCB with TR when the ply angle is 20°. As can be seen, with the increase of the TR, the natural frequency gradually increases and the analysis results agree well with the simulation results obtained using the ANSYS finite element method.

Figure 3.

Variation of the first-mode natural frequency of the RTWCB with the taper ratio (ply angle = 20°).

When the Galerkin method is employed to solve the problem, the number of modes has a certain influence on the accuracy of the results. Table 2 shows the influence of the number of modes on the analysis results. As can be seen, when the number of modes is larger than 2, the first- and second-mode natural frequencies both change slightly. Consequently, the mode number of 2 is used in the subsequent analysis by the Galerkin method.

Table 2.

Influence of the number of modes on the analysis results ([20°]₆, TR = 2).

Mode number	First-mode natural frequency (Hz)	Second-mode natural frequency (Hz)
1	262.816323933815	262.816323933815
2	256.937085040506	256.937085040506
3	256.921689345760	256.921689345760
4	256.910620352933	256.910620352933

TR: taper ratio.

Stability analysis

Verification of the solution method

In this section, the lobes of the RTWCB computed by the time-domain method in this study and by the frequency-domain method in Movahhedy and Mosaddegh¹⁵ will be compared with those of the steel milling bar shown in Figure 7 in Movahhedy and Mosaddegh.¹⁵

Figure 4 shows the time-domain stability results from equation (42) and the frequency-domain stability results from Movahhedy and Mosaddegh.¹⁵ It is obvious that the curve obtained from the two-mode time-domain method is in good agreement with that from the first-mode frequency-domain method within the most speed range.

Figure 4.

Comparison between the time-domain and frequency-domain methods.

In order to calculate the displacement response curve, the Runge–Kutta method is used to reduce equation (33) to a one-order differential equation, and the dde23 function provided by MATLAB can be used to solve the one-order delay differential equation. When the speed is 3000 rad/s, the cutting depth is 0.9 and 0.7 mm corresponding to Point A and Point B in Figure 4, respectively, and the displacement response curves of Point A and Point B are illustrated in Figure 5. According to the displacement response curves, Point A is divergent, indicating that there is a chatter phenomenon, while Point B is convergent, indicating that the cutting process is stable. These results are consistent with the analysis results shown in Figure 4. Point A locates in the unstable zone based on two and three modes, while Point B locates in the stable zone. The results demonstrate that the two-mode Galerkin method is necessary to obtain high-precision results. Figure 6 shows the convergence situations of different algorithms. As can be seen, the result of the two modes coincides with that of the three modes, and they show a certain difference from the result of one mode. This confirms the convergence and high precision of the proposed algorithm.

Figure 5.

Bending displacement response curves of Point A and Point B.

Figure 6.

Convergence situations of different algorithms.

Stability lobes of the RTWCB

In order to obtain the stability lobes of the RTWCB with different parameters, the program is written based on equation (42) using the two-mode Galerkin method in MATLAB software.

Figures 7 –9 illustrate the influence of the ply angle on the cutting stability of the RTWCB with different TRs when $c_{i} / (c_{e} + c_{i}) = 0.5$ . As can be seen, with the decrease of the ply angle, the scope of an unconditionally stable zone extends gradually and the lobes move to the right-hand side, which leads to the scope extension of the conditionally stable zone in the high rotating speed range. These changes are helpful to improve the milling efficiency.

Figure 7.

Influence of the ply angle on the cutting stability of the RTWCB when TR = 0 $(c_{i} / (c_{e} + c_{i}) = 0.5)$ .

Figure 8.

Influence of the ply angle on the cutting stability of the RTWCB when TR = 2 $(c_{i} / (c_{e} + c_{i}) = 0.5)$ .

Figure 9.

Influence of the ply angle on the cutting stability of the RTWCB when TR = 4 $(c_{i} / (c_{e} + c_{i}) = 0.5)$ .

Figures 10 –12 illustrate the influence of TR on the cutting stability of the milling cutter with different ply angles when $c_{i} / (c_{e} + c_{i}) = 0.5$ . From these figures, it is easy to find out that the scope of an unconditionally stable zone extends gradually with the increase of TR, and that the lobes also move to the right-hand side, which leads to the scope extension of the conditionally stable zone in the high rotating speed range. These changes are helpful to provide a larger spindle speed.

Figure 10.

Influence of the TR on the cutting stability of the milling cutter when θ = 0° $(c_{i} / (c_{e} + c_{i}) = 0.5)$ .

Figure 11.

Influence of the TR on the cutting stability of the milling cutter when θ = 20° $(c_{i} / (c_{e} + c_{i}) = 0.5)$ .

Figure 12.

Influence of the TR on the cutting stability of the milling cutter when θ = 90° $(c_{i} / (c_{e} + c_{i}) = 0.5)$ .

Figure 13 shows the influence of internal damping on the cutting stability. As can be seen, the internal damping has a weaker influence on the cutting stability than the ply angle or the TR. As the rotating speed increases, the influence of internal damping is enhanced gradually. Since the cutting speed is conventionally not very high, the influence of internal damping on the cutting stability can generally be ignored.

Figure 13.

Influence of internal damping on the cutting stability of the milling cutter when TR = 2 and θ = 20°.

Summary and conclusion

In summary, free vibration equations were derived and analyzed to investigate the free vibration characteristics of an RTWCB with bending–bending coupling, and chatter stability analysis was performed by incorporating regenerative milling force and internal/external damping into the free vibration equations. It was found that the ply angle and TR of the composite milling bar had a significant influence on the natural frequency and that a smaller ply angle and a larger TR could increase the natural frequency. Also, the two-mode Galerkin method could obtain accurate results for chatter stability analysis of the RTWCB. Based on the semi-discrete method and the two-mode Galerkin method, the stability lobes of the RTWCB for different parameters could be predicted, and a smaller ply angle and a larger TR could extend the stable cutting zone and improve the cutting efficiency. The internal damping had an increasingly significant influence on the chatter stability with the increase of the cutting speed, and its influence on the cutting stability was generally not obvious considering the relatively low cutting speed in the conventional machining process.

Footnotes

Handling Editor: Nima Mahmoodi

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 was supported by the National Natural Science Foundation of China (Grant Nos 11672166 and 51605264).

ORCID iD

Jingmin Ma

References

Tlusty

Dynamics of high-speed milling. J Eng Ind 1986; 108: 59–67.

Tlusty

Macneil

Dynamics of cutting forces in end milling. Ann CIRP 1975; 24: 21–25.

Munoa

Beudaert

Dombovari

et al . Chatter suppression techniques in metal cutting. CIRP Ann: Manuf Tech 2016; 65: 785–808.

Tanaka

Obata

Matsubara

et al . Active chatter suppression of slender boring bar using piezoelectric actuators. JSME Int J 1994; 37: 601–606.

Matsubara

Maeda

Yamaji

Vibration suppression of boring bar by piezoelectric actuators and LR circuit. CIRP Ann: Manuf Tech 2014; 63: 373–376.

Chen

Altintas

Magnetic actuator for active damping of boring bars. CIRP Ann: Manuf Tech 2014; 63: 369–372.

Madoliat

Hayati

Ghalebahman

AG.

Modeling and analysis of frictional damper effect on chatter suppression in a slender endmill tool. J Adv Mech Des Syst 2011; 5: 115–128.

Yiqing

Qiang

Min

Optimization of the tuned mass damper for chatter suppression in turning. J Vib Eng Technol 2010; 23: 717–724.

Ema

Marui

Suppression of chatter vibration of boring tools using impact dampers. Int J Mach Tool Manu 2000; 40: 1141–1156.

10.

Bediaga

Hernández

Munoa

et al . Comparative analysis of spindle speed variation techniques in milling. J Membrane Sci 2004; 443: 69–82.

11.

Al-Regib

Lee

SH.

Programming spindle speed variation for machine tool chatter suppression. Int J Mach Tool Manu 2003; 43: 1229–1240.

12.

Lee

Suh

NP.

Manufacturing and testing of chatter free boring bars. CIRP Ann: Manuf Tech 1988; 37: 365–368.

13.

Nagano

Koizumi

Fujii

et al . Development of a composite boring bar. Compos Struct 1997; 38: 531–539.

14.

Dai Gil

Hui Yun

Jin Kook

. Design and manufacture of a carbon fiber epoxy rotating boring bar. Compos Struct 2003; 60: 115–124.

15.

Movahhedy

Mosaddegh

Prediction of chatter in high speed milling including gyroscopic effects. Int J Mach Tool Manu 2006; 46: 996–1001.

16.

Kim

Argento

Scott

RA.

Forced vibration and dynamic stability of a rotating tapered composite Timoshenko shaft: bending motions in end-milling operations. J Sound Vib 2001; 246: 583–600.

17.

Tian

Hutton

SG.

Chatter instability in milling systems with flexible rotating spindles—a new theoretical approach. J Manuf Sci Eng 2001; 123: 1–9.

18.

Berdichevsky

Armanios

Badir

Theory of anisotropic thin-walled closed-cross-section beams. Compos Eng 1992; 2: 411–432.

19.

Badir

AM.

Analysis of advanced thin-walled composite structures. PhD Thesis, Georgia Institute of Technology, Atlanta, GA, 1992.

20.

Armanios

Badir

AM.

Free vibration analysis of anisotropic thin-walled closed-section beams. AIAA J 1994; 33: 1905–1910.

21.

Dancila

Armanios

EA.

The influence of coupling on the free vibration of anisotropic thin-walled closed-section beams. Int J Solids Struct 1998; 35: 3105–3119.

22.

Min-young

Jeng-Keag

Chih-Yung

A simple spinning laminated composite shaft model. Int J Solids Struct 2004; 41: 637–662.

23.

Jingmin

Yongsheng

Wenli

Free vibration analysis of composite thin-walled cantilever beams with variable closed cross-section. J Vib Shock 2012; 31: 117–123 (in Chinese).

24.

Chandra

Singh

Gupta

Damping studies in fiber-reinforced composites—a review. Compos Struct 1999; 46: 41–45.

25.

Smith

Tlusty

An overview of modelling and simulation of the milling process. J Eng Ind 1991; 113: 169–175.

26.

Insperger

Stépán

Updated semi-discretization method of periodic delay-differential equations with discrete delay. Int J Numer Meth Eng 2004; 61: 117–141.