A three-dimensional model for the dynamics of micro-endmills including bending, torsional and axial vibrations

Precision Engineering 35 (2011) 24–37

Contents lists available at ScienceDirect

Precision Engineering journal homepage: www.elsevier.com/locate/precision

A three-dimensional model for the dynamics of micro-endmills including bending, torsional and axial vibrations Sinan Filiz, O. Burak Ozdoganlar ∗ Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh, PA 15213, United States

a r t i c l e

i n f o

Article history: Received 14 September 2009 Received in revised form 9 May 2010 Accepted 17 May 2010 Available online 27 May 2010 Keywords: Micro-endmill dynamics model Micromilling dynamic testing Micromachining dynamics

a b s t r a c t Attainable geometric accuracy and surface finish in a micromilling operation depends on predicting and controlling the vibrations of micro-endmills. The specific multi-scale geometry of micro-endmills results in complexities in dynamic behavior, including three-dimensional vibrations, which cannot be accurately captured using one-dimensional (1D) beam models. This paper presents an analytically based three-dimensional (3D) model for micro-endmill dynamics, including actual cross-section and fluted (pretwisted) geometry. The 3D model includes not only bending, but also coupled axial/torsional vibrations. The numerical efficiency is enhanced by modeling the circular cross-sectioned shank and taper sections using 1D beam equations without compromising in model accuracy, while modeling the complex cross-sectioned and pretwisted fluted section using 3D linear elasticity equations. The boundary-value problem for both 1D and 3D models are derived using a variational approach, and the numerical solution for each section is obtained using the spectral-Tchebychev (ST) technique. Subsequently, component mode synthesis is used for joining the individual sections to obtain the dynamic model for the entire tool. The 3D model is validated through modal experimentation, by comparing natural frequencies and modeshapes, for two-fluted and four-fluted micro-endmills with different geometries. The natural frequencies from the model was seen to be within 2% to those from the experiments for up to 90 kHz frequency. Comparison to numerically intensive, solid-element finite-elements models indicated that the 3D and FE models agree with less than 1% difference in natural frequencies. The 3D-ST model is then used to analyze the effect of geometric parameters on the dynamics of micro-endmills. © 2010 Published by Elsevier Inc.

1. Introduction The precision of machining processes depends greatly on the dynamic behavior of the machining system as reflected at cutting points, including the dynamics of the cutting tools [1,2]. The forced vibrations can cause surface-location errors, leading to reduced dimensional accuracies and increased surface roughness [3]. Furthermore, self-excited vibrations (chatter), again arising from the dynamic behavior of the machining system, significantly limits the productivity of a machining operation [4–6]. The magnitude of these errors depends on the dynamic response of the tool [7], and can be minimized if the dynamic behavior of cutting tools are well-characterized. A comprehensive dynamic model of the tool facilitates predicting the system response, and minimizing the errors by determining the most favorable cutting conditions (speed, feed, and depth of cut) and cutting-tool geometry.

∗ Corresponding author at: Department of Mechanical Engineering, Carnegie Mellon University, 303 Scaife Hall, 5000 Forbes Ave., Pittsburgh, PA 15213, United States. Tel.: +1 412 268 9890; fax: +1 412 268 3348. E-mail address: [email protected] (O.B. Ozdoganlar). 0141-6359/$ – see front matter © 2010 Published by Elsevier Inc. doi:10.1016/j.precisioneng.2010.05.003

In micromilling, which uses micro-scale milling tools to fabricate three-dimensional features on a wide range of materials, the accuracy and surface-finish requirements are exacerbated due to the micro-scale feature dimensions [8]. In addition to precise kinematic (machine) motions, spindle rotations and tool geometry, attaining the required level of precision necessitates controlling the vibrations during micromilling. Forced and self-excited vibrations not only deteriorate the part accuracy and surface roughness, but also limit the overall productivity of the system by causing cutting instability, tool breakage, and machine-tool damage. Being one of the most flexible components in a micromachining system, micro-scale cutting tools commonly dominate the dynamic behavior during micromilling. The specific (sectioned) geometry of the micro-endmills, including a macro-scale shank, a taper, and the micro-scale fluted section (see Fig. 1), results in a large number of closely coupled dynamic modes (natural frequencies) to be present within the frequency range of interest. In addition to bending modes, coupled axial/torsional modes become important. Therefore, successful prediction and control of vibrations of micro-endmills require derivation of accurate, yet numerically efficient, dynamic models that will capture bending, axial and torsional modes.

S. Filiz, O.B. Ozdoganlar / Precision Engineering 35 (2011) 24–37

Fig. 1. Micro-endmill geometry, and parameters for geometric description.

The dynamics of macro-scale endmills have been investigated in the literature, with the aim of minimizing the undesired effects of forced and self-excited vibrations during milling, including poor dimensional accuracy, tool failure, and reduced productivity, e.g. [2,9]. Earlier approaches described the tool (and the spindle/machine-tool) dynamics using multi-degree-of-freedom (MDOF) lumped parameter models [6,10,11]. Increased prediction accuracy was later obtained by using distributed-parameter models (e.g., Euler–Bernoulli beam model) while considering an equivalent circular cross-section for the fluted section [12–15]. The important effect of pretwisted (fluted) geometry with actual crosssection was recognized in modeling drill dynamics [16–18]. It was shown that the actual non-axisymmetric cross-section, pretwist, and number of flutes significantly affect the dynamics of the macroscale cutting tools [16,17]. Recent literature has recognized the complexities in obtaining dynamic models for micro-scale cutting tools, largely due to their sectioned multi-scale geometry that induces three-dimensional vibrations. Jun et al. [19,20] analyzed the micro-endmill dynamics by considering an equivalent circular cross-section for the fluted region, and used Timoshenko beam equations to derive the boundary value problem. The effect of mass unbalance and setup errors was included in the model. The numerical solution of the problem was obtained using the finite-elements (FE) technique. Huang [21] modeled micro-drill vibrations by considering the twisted tip geometry and axial drilling force. The micro-drills are modeled as sectioned beams (shank, taper, and the fluted region) using the Euler–Bernoulli beam model. The numerical solution was obtained using the Galerkin’s method with cantilever-beam mode shapes as the basis functions. It was concluded that the non-axisymmetric cross-section and twisted geometry of the fluted region has a considerable effect on the frequency response of micro-drills. Gong et al. [18] modeled the micro-drills as twisted rotating Timoshenko beams and solved the boundary-value problem by using the finite-elements method. Their detailed study on the geometric parameters and boundary conditions showed the significance of these parameters on the dynamics of micro-drills, specifically on the critical speeds and buckling loads. Filiz et al. [22–24] modeled the dynamics of micro-endmills as sectioned Timoshenko beams, where the actual geometry of fluted section is considered. The effect of attachment errors and gyroscopic effects due to rotation were included in the models. The models are solved using the spectralTchebychev technique, thereby resulting a significantly simpler system description as compared to that from the FE technique [25]. Although one-dimensional beam models are widely used for modeling the bending vibrations of cutting-tools, their accuracy decreases when the pretwist rate is increased. Filiz et al. [24] have shown that increased pretwist rate induces large errors (>10%) to natural frequency predictions, especially for the bending modes that are dominated by the fluted tip-section. Additionally, the torsional and axial vibrations cannot be accurately modeled using one-dimensional beam models when the actual cross-section and the pretwisted geometry of the fluted section of the micro-endmills are considered. For instance, finding the torsion constant requires calculation of the warping function over the cross-section, which,

25

in turn, requires solution of a two-dimensional boundary value problem [26,27]. The pre-twist also causes coupled torsional-axial vibrations [28,29]. This paper presents an analytically based three-dimensional (3D) model for micro-endmill dynamics. The model includes the actual cross-sectional geometry and pretwist of the fluted section, and captures the bending, torsional, and axial vibrations of microendmills. To attain increased numerical efficiency, portions of the micro-endmill with circular cross-section (shank and tapered sections) are modeled using 1D beam models without compromising the model accuracy. Three-dimensional linear elasticity equations are used to describe the fluted section, including the pretwisted shape and actual cross-sectional geometry. The boundary-value problem (BVP) for the fluted section is derived using a variational approach based on the extended Hamilton’s principle. The spectral-Tchebychev technique is subsequently used to spatially discretize the BVPs, and component mode synthesis is used to join the sections to obtain the model for the entire micro-endmill. The derivation of receptances (frequency response functions) from the model is then described. The complete 3D model is validated by comparing natural frequencies and mode shapes from 3D model to those from the modal tests on micro-endmills with different geometries. The 3D model is also compared to a solid-element FE model. The validated model is then used to analyze the effect of geometric parameters on the dynamics of the micro-endmills. The original contributions of this work include (1) a new 3D model and associated spectral-Tchebychev solution that not only capture the bending motions of micro-endmills more accurately than previous 1D solutions (e.g., in [22–24]), but also incorporate the coupled torsional-axial motions, actual cross-sectional geometry and the pretwist of the tools, (2) combination of 1D (for the shank and taper regions) and 3D (for the fluted region) models through component mode synthesis to attain a computationally efficient model, (3) a new experimental approach that enables modal testing for both bending and torsional vibrations of micro-tools, and (4) utilization of the model to determine the effect of tool geometry (diameter, aspect ratio, and helix angle) on the natural frequencies of the micro-endmills. 2. Derivation of the boundary-value problem The overall geometry of micro-endmills include a macro-scale shank section, a micro-scale fluted section, and a taper section that connects the shank and fluted sections (see Fig. 1). Having circular cross-sections, the dynamic behavior of the shank and taper sections is simplified, resulting in independent (uncoupled) axial, torsional, and bending deflections. For those sections, one-dimensional beam equations accurately capture the dynamic behavior [30]. However, complex non-axisymmetric cross-section and the pretwisted geometry of the fluted section induces complicated dynamic behavior, including coupled axial, torsional, and bending motions, which cannot be captured with one-dimensional models. Being lower order, one-dimensional equations can be solved with significantly higher numerical efficiency than three-dimensional equations. Therefore, to simultaneously attain numerical efficiency and modeling accuracy, the shank and taper sections are modeled using one-dimensional models, and the fluted section is modeled using a three-dimensional model. The complete dynamic model for the micro-endmill is then obtained by joining individual sections using a component mode synthesis approach. 2.1. Derivation of the one-dimensional model In this work, one-dimensional Timoshenko beam equations derived in [22,23] are used to model the bending motions of the

26

S. Filiz, O.B. Ozdoganlar / Precision Engineering 35 (2011) 24–37

shank and taper sections. To model the axial and torsional motions, the extended Hamilton’s principle was also used as in [23] by writing the strain and kinetic energies for axial motions as 1 2

ESax =





L

EA(z) 0

∂w ∂z

2



1 2

EKax =

dz,



L

A(z) 0

∂w ∂t

2

1 2

EStr =



L

GJs (z) 0

∂ z ∂z

2 dz,

1 2

EKtr =





L

Ip (z) 0

∂ z ∂t

dz,

(1)

2 dz, (2)

where z is the torsional deflection, Ip (z) is the polar moment of inertia, G is the shear modulus, and Js (z) is the torsion constant [30]. In [22,23], the BVP for the 1D Timoshenko beam was derived in the differential equation form. Here, we use an alternative approach, which lends itself to a simpler representation and becomes critically effective for derivation of the BVP for 3D cases. At this point, rather than deriving the differential BVP, the extended Hamilton’s principle is applied to obtain the integral form of the BVP for bending, axial and torsional motions of shank and taper sections as

 L 0

∂2 y ˆ ∂2 ∂2 u ∂2 v ∂2 w ˆ A 2 uˆ + Ix vˆ + Iy 2 x ˆ x + A 2 w y + A 2 2 ∂t ∂t ∂t ∂t ∂t

+Ip

∂2 ∂t 2

z

ˆ z + EIx ∂ y + ks GA ∂z ∂z

∂ x ∂ˆx +EIy + ks GA ∂z ∂z

 +Pa (z)

+EA

∂ˆy



∂u ∂uˆ ∂v ∂ˆv + ∂z ∂z ∂z ∂z

∂v + ∂z



 x

∂u − ∂z

∂ˆv + ˆx ∂z

 y

∂uˆ − ˆy ∂z



t2

ı(EK − ES + Wnc ) dt = 0,

ıqi (x, y, z, t) = 0

at t = t1 , t2 ,

ES =

 L

1 2

{}T {} dx dy dz,

ES =

 L

1 2

{q}T [B]T [BC ]{q} dx dy dz,

where [BC ] = [C][B]. For the three-dimensional structure, the kinetic energy and the work done by non-conservative forces is expressed as [31] EK =

 L

1 2

T

˙ {q} ˙ dx dy dz,  {q} 0

 L

{Fq }T {q} dx dy dz,

t2

 L

¨ T {ıq} − {q}T [B] [BC ]{ıq} + {Fq }T {ıq}] dx dy dz dt = 0. [−A{q} L

(3)

0 T

T ˆ ˆ x ˆ y ˆ z } are the functions representing the where {ˆq} = {uˆ vˆ w variation of the terms. The boundary value problem in Eq. (3) is then discretized using the spectral-Tchebychev technique [25,23]. The discretized system of equations are in the form of

qˆ (M1 q¨ + K1 q − F1 ) = 0

(4)

where q and qˆ are the discretized vectors of the deflections and variation of the deflections, respectively. Since this equation must be valid for an arbitrary qˆ , the equations of motion can be written as M1 q¨ + K1 q = F1 ,

(10)

Area

respectively, where {Fq } = {Fx Fy Fz }T . The variational approach will be applied by first substituting Eqs. (8)–(10) into Eq. (6). Upon integration by parts, we obtain the variational form as

t1

{Fq }T {ˆq} dz,

dz =

(9)

Area

T



(8)

Area

0

0



(7)

Area

0

Expressing {q} = {u v w}T and writing the strain-displacement relation as {} = [B]{q}, and the constitutive equations as {} = [C]{}, where [C] is the constitutive matrix [26], the strain energy can be rewritten in terms of the displacement vector {q} as





(6)

where EK , ES and Wnc represent the kinetic energy, strain (potential) energy, and the work done by nonconservative forces, respectively, t1 and t2 are two instants of time, and z is the spatial variable. Here, q1 = u, q2 = v, and q3 = w, and ı is the variational operator. Using the 3D linear elasticity theory, the strain energy for this case can be written as [31]

Wnc =



ˆ ∂w ∂w ∂ z ∂ˆz + GJs ∂z ∂z ∂z ∂z



t1

where z is the axial coordinate, w is the deflection along the tool axis, L is the beam length,  is the density, E is the Young’s modulus, and A(z) is the cross-sectional area. Similarly, the strain and kinetic energies arising from the torsional motions are written as



expressing the extended Hamilton’s principle as

(5)

where M1 , K1 and F1 are mass, stiffness, and forcing matrices, respectively, which are expressed in terms of geometric and material properties, and Tchebychev matrices [25,23].

(11)

Area

0

To impose the variational condition [30], the deflection will be expressed as ¯ {q(x, y, z, t)} = {q(x, y, z, t)} + [C0 ]{ˆq(x, y, z, t)},

(12)

where q(x, y, z, t) is the solution of the boundary value problem, ¯ q(x, y, z, t) represents a family of admissible functions, [C0 ] is a (diagonal) matrix of arbitrary scalars, and qˆ (x, y, z, t) is a timeinvariant arbitrary function [30,32]. Eq. (12) states that when an ¯ is used for approximating the arbitrary admissible function ({q}) solution, there will be a variation [C0 ]{ˆq} from the true solution {q}. Using Eq. (12), the variations (ı terms) can expressed as {ıq} = [C0 ]{ˆq}. Inserting this expression into Eq. (11), and considering the validity of the resulting expression for any t1 and t2 , we obtain the integral form of the 3D boundary value problem as

 L

¨ T {ˆq} + {q}T [B]T [BC ]{ˆq}] dx dy dz [A{q} 0

Area

 L {Fq }T {ˆq} dx dy dz = 0.

= 0

(13)

Area

2.2. Derivation of the three-dimensional model 2.3. Application to the pretwisted complex cross-section Unlike the one-dimensional model that defines deflections of and rotations about the tool axis with a total of six degrees-offreedom, the three-dimensional model defines three translational degrees of freedom (uk , vk , wk ) for each point k within the volumetric domain. The derivation of the boundary value problem begins by

Application of the volume integrals obtained in Eq. (13) to complex cross-sectional geometries and pretwisted shapes poses significant challenges. However, the domain of the problem can be simplified by using coordinate transformations. For this purpose,

S. Filiz, O.B. Ozdoganlar / Precision Engineering 35 (2011) 24–37

27

Fig. 2. (a) The domain in the physical coordinates, (b) the domain after the twist is removed using a local reference frame, and (c) the domain after mapping.

two coordinate transformation are applied here. First, to remove ¯ z¯ ) that rotates with the the pretwist, a local reference frame (¯x, y, pretwist is defined (Fig. 2(b)). This transformation can be achieved by using a rotational transformation matrix that imposes a righthand rotation of ˛ z about the z axis, where ˛ is the twist rate. ¯ and z¯ Next, the complex cross-section is mapped from x¯ , y, coordinates onto a rectangular cross-section defined by ,  and  coordinates (Fig. 2(c)). The details of the polynomial mapping used for this transformation are described in Appendix A. For an effective transformation, a one-to-one mapping is required, which is guaranteed if the Jacobian is positive for all sampling points. For either of the transformation, the derivative operators (in [B] and [BC ]) in the physical (x, y, z) domain are related to the final mapped domain (, , ) using Jacobian matrices as ∂  x¯ ∂ = Jix¯jx Jij ∂xj ∂i

and

∂ ∂ = Ji−1 , j ∂ ∂xj i

i, j = 1, 2, 3,

(14)

 x¯

where Jix¯jx and Jij ’s are the elements of the Jacobian matrices.

Using the differential operators in Eq. (14) along with dx dy dz = |J| d d d, the functions in Eq. (13) can be expressed in the mapping coordinates ,  and . It should be noted that, although simplifying the domain of the problem, the coordinate transformations increase the complexity of the differential equations. 2.4. Spectral-Tchebychev solution of the three-dimensional model To apply the spectral-Tchebychev technique to the threedimensional model in Eq. (13), the displacements in the (,,¯z ) domain is expressed as a triple Tchebychev expansion as N

qi (, , , t) 

 N N  

I

O

O



q = Iu q, v =



q qˆ dx dy dz = qT V xyz qˆ

and

i

0

(17)

Area

where qni,x is the n th spatial derivative of q with respect to the i

variable xi , Qnxi is the n th extended derivative matrix with respect to coordinate xi , and V xyz is the inner product matrix with dimension (N N N × N N N ) (see Appendix B). Substituting Eq. (17) into Eq. (13), and applying the derivative and inner-product expressions in the Tchebychev domain (Eq. (17)), we obtain T

qˆ (M q¨ + Kq − F) = 0,

(18)

where M, K and F are mass, stiffness, and forcing matrices, respectively. Since this equation must be true for any admissible function qˆ , the term in the parenthesis must be zero, yielding M q¨ + Kq = F,

(19)

where M = (ITu V xyz Iu + ITv V xyz Iv + ITw V xyz Iw ), K = BT Vxyz BC , F = V xyz F q ,

(20)

and Vxyz is a matrix assembled with the inner product matrix V xyz on its diagonal, such that the inner product operation can be performed on the discretized derivative operator matrices B and BC . The system matrices expressed in Eq. (20) only capture the natural (force) boundary conditions directly, and the procedure described below (with projection matrices) is used to incorporate the essential and time-dependent boundary conditions.

(15)

where T’s are the scaled (orthogonal) Tchebychev polynomials, aqi are the time-dependent expansion coefficients, and Nj is the number of polynomials associated with the coordinate j = , , . To express the equations in a spatially sampled form, Gauss–Lobatto sampling is used in three-dimensional (,,¯z ) domain. Although the sampling produces a third-rank tensor, a tensor-vector mapping is used here to express the sampled functions as vectors (see Appendix B). The sampled (discretized) displacements are now written as



 L

qni,x = Qnxi Ii q

3. System matrices for the entire micro-endmill aqi (t)Tl−1 ()Tm−1 ()Tn−1 (),

l=1 m=1 n=1

u=

[25]

O

I

O



q = Iv q, w =



O

O

I



q = Iw q,

(16)

where I and O are (N N N × N N N ) identity and zero matrices, respectively, and q = {u v w}T . Underline indicates the spatially sampled version of a variable. The derivative (in [B] and [BC ]) and the inner product (integration) operations are imposed in the sampled domain by writing

To obtain a complete model for the micro-endmill with the shank, taper, and fluted sections, a component mode synthesis technique will be used. The technique is based on equating the deflections and slopes along the shared boundaries between the sections, and assembling the individual mass, stiffness, and force matrices for each section into global mass, stiffness, and forcing matrices. While the application of component mode synthesis is straightforward between the shank and taper sections that are modeled using one-dimensional equations [33,23,24], connecting the one-dimensional taper section with the three-dimensional fluted section is more complicated. Although the natural (force) boundary conditions are accounted for within the integral form of the boundary value problem in Eq. (13), essential (displacement) and time-dependent boundary conditions must be imposed to complete the system equations. For applying the essential and time-dependent boundary conditions, the spectral-Tchebychev technique uses basis recombination as {q} = P{qd } + R{qb },

(21)

28

S. Filiz, O.B. Ozdoganlar / Precision Engineering 35 (2011) 24–37

where P and R are the projection matrices [25]. In this expression, P{qd } satisfies the essential boundary conditions. Furthermore, since the admissible functions (functions for the variation of the deflection terms) are required to satisfy the essential boundary conditions, they can be expressed as {ˆqa } = P{ˆq}. The boundary conditions are automatically satisfied at the free ends. The displacement and slope compatibility conditions have to be satisfied at the boundaries combining the micro-endmill sections. The one-dimensional solution is expressed as three displacements and three rotations about the axis of the cross-section, whereas the three-dimensional solution is expressed as three displacements for every sampled point distributed within the cross-section. For a sampling point located at (x, y) within the crosssection, the three displacement degrees of freedom (u, v, w) of the three-dimensional model can be obtained in terms of the six ¯ v¯ , w, ¯ degrees-of-freedom (u, x , x , x ) of the one-dimensional model as {u, v, w} = {u¯ − y

v¯ + x

z,

z,

¯ +x w

y

+y

x }.

(22)

It should be noted that this expression assumes that the rotational transformations are independent from the order of rotations, which is valid for small rotational deformations seen during vibrations of the micro-endmills. Furthermore, small-angle assumption has been imposed in arriving at Eq. (22). Projection matrices are calculated by applying singular value decomposition on the boundary conditions matrix ˇ, which is formed from the boundary conditions as



ˇ=

ˇ1 0

0 ˇ2



,

(23)

where the compatibility conditions (the deflection and slope equality) on the boundary between the shank and taper sections are satisfied through the matrix ˇ1 , and the compatibility conditions (Eq. (22)) on the boundary between the taper and fluted section are satisfied through the ˇ2 matrix. The construction of ˇ1 and ˇ2 matrices are explained in Appendix C. The mass and stiffness matrices, and the forcing vector, for the entire of the micro-endmill are obtained by applying the projection matrices on the global mass and stiffness matrices and on the global forcing vector. Thus, the equations of motion for the entire endmill can be written as Mq¨ d + Kqd = F.

(24)

The global mass M and stiffness K matrices are formed using the mass and stiffness matrices computed for the shank (1) and taper (2) (see Eq. (5)) and fluted (3) section (see Eq. (19)) under free-free boundary conditions such that M = P T Mg P,

K = P T Kg P,

where



Mg =

M1 0 0

0 M2 0

0 0 M3

(25)



,

Kg =

K1 0 0

0 K2 0

0 0 K3



T

I q¨ m + [ωn2 ] qm = [ M ] F,

(29)

where orthogonality and mass-normalization of modal vectors T are use to obtain the identity matrix I = [ M ] M[ M ], the diagonal eigenvalue (squared natural frequency) matrix [ωn2 ] = T T [ M ] K[ M ], and the modal force vector [ M ] F. Eq. (29) represents the system in the modal coordinates by a decoupled (independent) set of single-degree-of-freedom (SDOF) systems. Therefore, a model reduction can be realized by simply ignoring those SDOF equations that correspond to natural frequencies above the frequency range of interest. For instance, if the system equations is ordered from low-to-high eigenvalues, and only 20 modes are deemed to be sufficient for the model, only the first 20 equations can be considered. It should be noted that it is not required to find all the modes for applying the model reduction, since numeric eigenvalue solvers can provide modal matrices and natural frequencies up to a specific number of modes. These decoupled system representation also allows introduction of modal damping as T

I q¨ m + 2[d ] [ωn ] q˙ m + [ωn2 ]qm = [ M ] F,

(30)

where [n ] is a diagonal matrix with modal damping ratios (obtained from modal testing) associated with mode n in its diagonal, and [ωn ] is the diagonal natural frequency matrix. One way to obtain the complex receptance matrix for the system is to substitute q˙ m = Qm ejωt into Eq. (30) as T

[−I ω2 + 2j ω [n ] [ωn ] + [ωn2 ]] Qm = [ M ] F. .

(26)

The forcing vector is computed as F = P T Fg − P T Mg Rq¨ b − P T Kg Rqb , T {F1T F2T F3T }

(31)

Considering the coordinate transformation (to modal coordinates) Qd = [ M ]Qm , after manipulation, we obtain Qd = [G(jω)] F,

(27)

where Fg =

lating the vibration response of the system. However, the model presented in Eq. (24) must be appended by adding the damping characteristics of the tool-holder-spindle system. Furthermore, as presented, the model will yield a large number of natural frequencies, many of which will be outside the frequency range of interest. In this section, a model reduction technique, as well as addition of modal damping to system equations, are described. Moreover, the reduced and damped model is used to derive the system receptance (i.e., frequency response function between the applied force and displacement response) in the frequency domain, since the receptance is the preferred representation of system dynamics for stability analysis. The spectral-Tchebychev technique described above results in symmetric system matrices M and K. Furthermore, the mass matrix M is positive definite, and the stiffness matrix K is positive semi-definite. Considering these properties of system matrices, a coordinate transformation to modal coordinates can be defined as qd = [ M ]qm , where qm is the vector of modal coordinates, and [ M ] is the modal matrix that includes mass-normalized modal vectors obtained from the algebraic eigenvalue problem associated with Eq. (24). Applying this coordinate transformation to Eq. (24), and pre-multiplying with the transpose (i.e., inverse) of the modal matrix, we obtain

where [G(jω)] = [ M ] [−Iω2 + 2j ω [n ] [ωn ] + [ωn2 ]]

(28)

4. Obtaining frequency response function in the form of receptance As discussed in the introduction section, our main interest in deriving dynamic solution for the micro-endmills lies in simu-

(32)

−1

[ M ]

T

(33)

is the complex receptance matrix for the system. To demonstrate this approach, the receptance for the F2 tool (see Table 1) is calculated. A modal reduction is realized by considering only the first 20 modes of the system. Subsequently, a 0.5 percent modal damping ratio (i.e., n = 0.005) is applied to each mode of the system (clearly, the formulation allows adding different damping ratios to each mode.) Using these parameters, the complex 20 x 20

S. Filiz, O.B. Ozdoganlar / Precision Engineering 35 (2011) 24–37

29

Table 1 Geometric parameters of the micro-endmills used during the validation experiments. Tool

Tool diameter (␮m)

Number of flutes

Shank length (mm)

Taper length (mm)

Fluted length (Lt ) (mm)

Helix angle (◦ )

T1 T2 F1 F2

245 811 245 747

2 2 4 4

29.80 30.83 30.80 29.28

7.22 5.40 6.05 5.74

0.94 2.70 1.05 2.95

27.5 46.8 27.9 31.0

ces with their scaled forms is written as K = K∗ Lr3 ,

Fig. 3. Simulated receptance of the direct response of the F2 tool along the x direction at the tool tip. The endmill is fixed at the shank end.

receptance matrix [G(jω)] is obtained from Eq. (33). Fig. 3 presents the amplitude and phase of the receptance considering the direct response at the tip along the x direction (i.e., only considering the Glk (jω), where l corresponds to tip displacement on the tool axis in the x direction (see Fig. 2(a)), and k corresponds to the force input along the same direction). 5. Efficiency and convergence of the numerical solutions For attaining the optimum numerical performance, a scaling and non-dimensionalization will be applied to both the one-dimensional and three-dimensional equations. For one dimensional equations, a reference length Lr is selected. The displacement/length, area, and moment of inertia terms are then non-dimensionalized dividing them by Lr , Lr2 , and Lr4 , respectively. Furthermore, the moment terms are divided by Lr , and force and angular deflections are not manipulated. For the three-dimensional solution, all the displacement and length terms are divided by the reference length Lr , including the derivative terms such as ∂x. The force terms for the three-dimensional case are divided by Lr2 . Using the scaled parameters, the relationships between the system matri-

M = M∗ Lr5 ,

F = F∗ Lr2 ,

(34)

where the starred quantities represent the scaled forms. As such, the numerical solution produces the scaled natural frequencies, and the physical (actual) natural frequencies can be obtained by dividing the scaled natural frequencies by Lr . The number of polynomials to be used in Tchebychev expansion for both one-dimensional and three-dimensional (Eq. (15)) cases must be determined in a way to insure both the accuracy and the efficiency of the numerical solution. For the one-dimensional case, the convergence was studied in [25], and 12 and 14 terms were deemed sufficient to obtain convergence for shank and taper sections, respectively. Furthermore, it was shown that the convergence, as expected by the nature of the spectral-Tchebychev solution, is exponential [25]. The numerical convergence of the three-dimensional solution is studied in the mapped (, , ) domain (which is the domain used for numerical solution). During the convergence study, the entire micro-endmill is considered. A preliminary analysis has indicated that a larger number of polynomials are required along the the axial () direction than along the two lateral directions. Furthermore, since the cross-section is mapped onto a square shape, the same number of polynomials are selected for the  and  directions. To analyze the convergence, the difference between the natural frequencies of the selected number of polynomials and those for (9, 9, 22) polynomials (corresponding to , , ) is calculated as CRk = |ωi,k − ωi,r |/ ωi,r , where CRk is the absolute non-dimensional difference between the ith natural frequency when using current set of polynomials (k) with that when using the reference (r) set of polynomials. The convergence is studied by creating contour plots of the CRk for increased number of polynomials. Fig. 4 provides sample convergence plots for the third bending mode and the first torsional-axial mode of an 811 ␮m microendmill. As seen from the contour plots, a region (gray areas) can be identified for obtaining the optimum number of polynomials. Based on this convergence study, (7, 7, 16) polynomials are used to examine the fluted section of the micro-endmills used in this work. For the 811 ␮m micro-endmill, the difference between the natural frequencies obtained when using (7, 7, 16) and (9, 9, 22) polynomials are found to be 1.9 × 10−3 %, 2.5 × 10−3 %, −1.3 × 10−3 % and

Fig. 4. Numerical convergence for the T2-micro-endmill. Base 10 logarithm of the differences are plotted. (a) First bending-mode pair, first direction (B11), (b) first torsion-axial mode (TA1).

30

S. Filiz, O.B. Ozdoganlar / Precision Engineering 35 (2011) 24–37

Fig. 5. SEM images of micro-endmills. (a)–(e) T1, (b)–(f) T2, (c)–(g) F1, and (d)–(h) F2.

−7.5 × 10−3 % for the first bending, second bending, third bending, and first torsion modes, respectively. 6. Experimental validation of the model A set of modal test are conducted to validate the developed model and the associated numerical solution. Due to their small size, validation of micro-endmill dynamics brings forth various challenges [34–36], and renders the modal testing approaches used for macro-scale cutting tools (e.g., [20,34,35,37–44]) inapplicable. The excitation forces (or displacement) must provide sufficient energy to a wide range of frequencies (generally up to 100 kHz), which cannot be obtained through traditional impact hammers or conventional shakers. Non-contact response measurement is required to ensure that dynamics of the micro-endmills are not altered due to relatively large measurement sensors. Furthermore, well-characterized boundary conditions must be selected for the experimentation. Ultimately, the dynamic behavior is represented in the form of a frequency response function (FRF). Only a few attempts have been made in the literature to measure the dynamic response of micro-endmills. Jun et al. [20] used acoustic emission sensors for measurement of the micro-endmill vibrations during the operation of the machine tool. They compared the experimental and simulated results within the frequency range of the machining forces, which were used as the excitation. Filiz et al. [24,33] tested the dynamics of micro-endmills and micro-drills by using piezoelectric elements for excitation and a microscopecoupled laser Doppler vibrometer for measurement. The tools were suspended with elastic bands such that free-free boundary conditions are achieved. They also measured the operating deflection shapes (mode-shapes) by scanning the vibrations along the tool. In this work, the validation experiments are conducted in a manner similar to those in [24,33]. Both two-fluted and four-fluted micro-endmills with different geometries (cross-sectional shape, length of each sections, and twist rate; see Fig. 1) given in Table 1 are tested. To facilitate the interpretation of the results, and to compare the three-dimensional spectral-Tchebychev (3D-ST) model with that of a finite-element (FE) model, a three-dimensional FE model (using ANSYS® Version 11) is also developed for each microendmill geometry considered in Table 1. 6.1. Capturing the fluted geometry of the micro-endmill In order to accurately model the micro-endmill dynamics, the actual cross-section and pretwist of the fluted section must be

measured. To measure the cross-sectional geometry of the microendmills used for the validation experiments, a set of scanning electron microscope (SEM) images of the fluted cross-sections are obtained from the bottom (tip) of each tool (see Fig. 5). An image-processing procedure is then used in Matlab® to select the boundary of the cross-section, and to identify a number of sample points along the boundary. Similarly, SEM images from the side are used to determine the pretwist angle, section diameters, and section lengths. Table 1 gives the measured geometric parameters for the tools. Using the above-mentioned mapping technique, two-fluted, four fluted, and circular cross-sections are mapped (see Appendix A for details). To obtain a continuous description of the cross-section for the fluted section, as well as the transition from the fluted to the taper section (with circular cross-section), the coefficients for the cross-sectional mapping are obtained at seven locations of the micro-endmill fluted section. The first location is chosen as the intersection of fluted and taper sections, where the cross-section is circular, whereas the second location is selected within the fluted region at the end of the transition from the taper to fluted section. From the second to the seventh sampling location, which is at the tip of the micro-endmill, the cross-sectional geometry remains constant, and those locations are selected with a uniform spacing. After obtaining the mapping coefficients at seven locations, a sixthorder polynomial is fitted to describe the axial variation of each coefficient continuously along the micro-endmill axis. This process yields a three-dimensional continuous function that describes the geometry of the fluted section, including a smooth transition region from the fluted cross-section to the circular taper crosssection. The geometry of the F2 micro-endmill captured using this approach is given in Fig. 6. This geometric description is used in both the spectral-Tchebychev and the FE modeling of the microendmills. 6.2. Experimental setup Fig. 7 shows the setup used for validation experiments. To obtain free-free boundary conditions, micro-endmills are suspended using elastic bands, dynamic effects of which were previously determined to be negligible [24]. The excitation is obtained by using piezoelectric elements glued to the shank of each micro-endmill as shown in Fig. 7(b). The attachment is arranged such that the contraction and expansion of the piezoelectric element generate both force and torque at the center of the tool shank (see Fig. 7(c)) exciting bending and torsional

S. Filiz, O.B. Ozdoganlar / Precision Engineering 35 (2011) 24–37

31

along the tool axis enabled capturing operating deflection shapes (mode shapes). The collected data is processed to obtain voltage-todisplacement FRFs between the voltage to the piezoelectric element and the displacement response for each measurement point. Compiling the FRFs for all points enabled obtaining the mode shapes. To determine the natural frequencies from the experimental data, a single-degree-of-freedom curve-fitting procedure is used in the vicinity of each (well-separated) natural frequency. 6.3. Model assessment

Fig. 6. (a) The geometry of the F2 micro-endmill captured in the model and (b) zoomed view of the transition region from the fluted to taper section.

motions. A pseudo-random voltage is supplied to the piezoelectric elements within the frequency range of 1 Hz to 100 kHz. One hundred averages are used to minimize the effect of unbiased noise on the data. The vibration response is measured in a non-contact manner using a laser Doppler vibrometer. To reduce the spot size of the laser, the fiber-optic lasers are fed through a microscope. To differentiate bending and torsional/axial modes, measurements were conducted on the shim attached to the shank-end of microendmills (see Fig. 7(b)). For a bending mode, the shim experiences a translation motion, and for a torsional mode, the shim experiences a rotary motion. The displacement of the shim is scanned along the radial direction for each resonance peak to identify whether a mode is bending or torsional type. The measured displacements from the shim are given in Figs. 8 and 10 for each resonance peak. A one-axis computer-controlled slide with 1 ␮m resolution is used to control the measurement location. Obtaining data from multiple locations

When obtaining the natural frequencies and mode shapes from the model for the tungsten carbide micro-endmills, the material properties of  = 14450 kg/m3 , E = 580 GPa,  = 0.2 are used [20]. Based on the convergence study above, the shank and tapered sections are expanded with 12 and 14 polynomials, respectively, and the fluted section is expanded with 7, 7, and 16 polynomials in the two lateral and axial directions, respectively. To compare the 3D-ST model predictions, solid-element FE models of micro-endmills were developed in a commercial FEM solver (ANSYS® V11, Solid187 3D 10-node tetrahedral structural solid-type elements). The convergence of the FE model is studied, and for the F2 micro-endmill 69 thousand elements were seen to be sufficient (for up to the 8th natural frequency). The magnitude of experimentally obtained voltage-todisplacement FRF for the T1 micro-endmill is given in Fig. 8. Using the small shims attached to the shank, the third mode was identified as the only torsional mode up to 90 kHz for the T1 micro-endmill. For the T2 micro-endmill, the torsional mode was seen to be the fourth mode. The non-axisymmetric geometry of the 2-fluted micro-endmills causes each bending mode to split into two modes, similar to bending behavior of a rectangular beam (natural frequencies in two principal directions for the same mode shape would differ). For this reason, each mode is associated with a pair of closely coupled natural frequencies (e.g.,

Fig. 7. Experimental setup for micro-endmill experiments. (a) Laser Doppler vibrometer coupled with a microscope, (b) measurement area, and (c) the piezoelectric element placement on the shank for exciting both bending and torsional modes.

32

S. Filiz, O.B. Ozdoganlar / Precision Engineering 35 (2011) 24–37

Fig. 8. Voltage-to-displacement FRF for the T1 micro-endmill. Measured deflection shapes are indicated for each mode.

Table 2 Comparison of the experimentally obtained and predicted natural frequencies for the two-fluted micro-endmills. Tool

Mode

T1

B11 B12 B21 B22 TA1 B31 B32 B11 B12 B21 B22 B31 B32 TA1 B41 B42

T2

Piezo/LDV Nat. Freq. (Hz)

Error (%) 3D-ST vs. Piezo/LDV

Difference (%) 3D-ST vs. FEM

16,322 16,322 42,196 42,196 65,562 75,915 75,915

−1.0 −1.0 −0.9 −0.9 −0.5 −0.8 −0.8

−0.07 −0.07 −0.17 −0.17 −0.03 −0.29 −0.24

15,430 15,430 38,950 38,466 50,721 56,153 63,271 77,056 78,269

−0.9 −0.8 −0.8 −1.0 0.5 1.1 0.1 −1.6 −1.5

−0.04 −0.03 −0.01 −0.05 0.86 0.87 0.00 −0.09 −0.03

B11 and B12 are the first pair of frequencies for the first bending mode. The comparison of natural frequencies between the 3D-ST model and experiments, and between the 3D-ST model and the FE model are given in Table 2 for both T1 and T2 micro-endmills. The 3D-ST model captured the natural frequencies with less than 1% error for T1, and with less than 1.6% error for T2. In contrast, the 1D pretwisted model presented in [24] produced a 7.7% error for the third bending mode (B31), which arises from the motions of the fluted section (the 3D-ST model has less than 1% error for the same mode.) In addition, the torsional/axial mode was not predicted in [24]. Compared to FE model (which is used as the reference), the differences between the 3D-ST and FE models are less than 0.3% and 0.9% for T1 and T2, respectively. Fig. 9 provides the mode-shape comparison for the T2 microendmill. As seen in Fig. 9, the mode shapes obtained from the experiments and those calculated using the 3D-ST model match well for the bending modes. Since the torsional mode-shapes could not be measured accurately from the micro-scale fluted sections of the micro-endmill, the torsional mode shape from the 3D-ST model is compared to that from the FE model, showing a good match. Fig. 10 gives the magnitude of experimentally obtained voltageto-displacement FRF for the F1 micro-endmill. For an ideal cross-sectional geometry with four-fold symmetry, the mode splitting does not occur (any two mutually perpendicular directions within the cross-section are the principal directions). In reality, however, grinding errors may disrupt the four-fold symmetry of the four-fluted micro-endmills. In Fig. 10 (zoomed images in Fig. 10(b) and (c)), slight mode-splitting on two bending modes are seen. Since the cross-section is considered to be (four-fold) symmetric for the four-fluted tools during modeling, the model only produce

Fig. 9. Mode-shapes of the T2 micro-endmill. Solid lines represent the modelprediction, triangles represent the experimental mode-shapes, and dotted line represent the mode-shape from the FEM solution. (a) First bending mode, (b) second bending mode, (c) third bending mode, (d) first torsional/axial mode, and (e) fourth bending mode. Table 3 Comparison of the experimentally obtained and predicted natural frequencies for the four-fluted micro-endmills. Tool

Mode

Piezo/LDV Nat. Freq. (Hz)

Error (%) 3D-ST vs. Piezo/LDV

Difference (%) 3D-ST vs. FEM

F1

B11 B12 B21 B22 TA1 B31 B32

15,444 15,534 40,653 40,653 62,998 74,637 74,724

2.0 1.4 1.4 1.4 1.1 1.2 1.0

−0.03 −0.03 −0.10 −0.10 −0.01 −0.20 −0.20

F2

B11 B12 B21 B22 B31 B32 TA1 B41 B42

16,353 16,426 40,307 40,307 52,293 52,293 65,774 81,184 81,512

0.9 0.4 −1.6 −1.3 −1.8 −1.2 1.0 0.5 0.2

−0.05 −0.05 0.13 0.11 0.58 0.57 0.00 −0.10 −0.10

a single natural frequency for each bending mode. It should be noted, however, that if the actual geometry is precisely provided, the model will capture the mode splitting behavior (as it does for the two-fluted tools). The natural frequencies from the experiments and the 3D-ST model are presented in Table 3. The 3D-ST model predicts the natural frequencies of the F1-micro-endmill with better than 2% accuracy. The differences between the 3D-ST predictions and the FE results agree with less than 0.2% difference. For F2, the largest difference between the 3D-ST model and the experiments is −1.8% on the third mode. The differences between the 3D-ST and FE results are less than 0.6% for F2. The mode-shapes of the F2-micro-endmill are obtained from the 3D-ST model and compared with the mode-shapes from the experiments (for the bending modes) and from the FE model (for the torsional/axial mode). The mode-shape comparison presented in Fig. 11 shows that the model successfully captures the modeshapes for the four-fluted micro-endmill.

S. Filiz, O.B. Ozdoganlar / Precision Engineering 35 (2011) 24–37

33

Fig. 10. Voltage-to-displacement FRF for the F1 micro-endmill.

7. Model application In this section, the validated 3D-ST model is used to analyze the effect of tool diameter, aspect ratio, and helix angle of the fluted section on the natural frequencies of micro-endmills. The tool material is considered as tungsten carbide, using the material properties mentioned above. 7.1. The effect of micro-endmill diameter A parametric study is performed to analyze the effect of tool diameter (i.e., the diameter of the fluted section) on the natural frequencies of the micro-endmills. The shank diameter and the total length of the micro-endmill (Ls + Lp + Lt ) are kept constant at the T2 values given in Table 1. To resemble commercially available microtooling, which have the same total tool length, the total length of

the micro-endmill is kept constant. The taper angle ( ) is set to be 12◦ (see Fig. 1). The taper and shank lengths are specified based on the selected diameter. The aspect ratio of the tip (Lt /dt ) and helix angle () are kept constant at 4 and 30 degrees, respectively. The fluted section cross-sectional geometry of T2 is used in the analysis (see Fig. 5). The variation of the bending and torsional/axial natural frequencies for varying tool diameter is plotted in Fig. 12. Fig. 12(a)–(c) indicates that the mode-splitting becomes more pronounced when the tool diameter is increased. The reason for this is the increased influence of the fluted section (which causes the mode splitting) on the dynamic behavior of the entire micro-endmill with increased tool diameter. Since the aspect ratio of the tip section is kept constant, the length of the fluted section increases for increasing tool diameter. The fluted-section dynamics dominate the dynamic behavior on the first, second and third mode for diameters larger than 800, 600 and 400 ␮m, respectively. It is observed from Fig. 12(a) and (b) that the bending natural frequencies first increase, reach a peak value at a certain diameter, and then decrease. The decrease in bending natural frequencies is coincident with fluted-section deformations becoming dominant. When the deformation of a mode is not dominated by the fluted-section deformation, the increase in the diameter causes a reduction in the total modal mass of the micro-endmill (due to the constant aspect ratio and length) leading to increasing natural frequencies. When the fluted-section deformation becomes dominant, the increase in the aspect ratio decreases the modal stiffness, leading to a decrease in natural frequencies. As seen in Fig. 12(d), the torsional/axial natural frequency monotonously increases with increasing tool diameter. Increasing diameter (and length) of the fluted section leads to a reduction in the modal mass of the first torsional/axial mode, which arises from the shank of the micro-endmill, while not affecting the stiffness of the mode. Consequently, the natural frequency for this mode increases with increasing aspect ratio. 7.2. The effect of aspect ratio

Fig. 11. Mode-shapes for the F2 micro-endmill. Solid lines represent the modelprediction, triangles represent the experimental mode-shapes, and dotted line represent the mode-shape from the FEM solution. (a) First bending mode, (b) second bending mode, (c) third bending mode, (d) first torsional/axial mode, and (e) fourth bending mode.

In this analysis, the tip section aspect ratio (Lt /dt ) is varied from 2 to 5 while keeping the diameter constant at 500 ␮m. The shank diameter and the total length of the micro-endmill (Ls + Lp + Lt ) are kept constant at the T2 values. The taper angle ( ) and helix angle () are kept at 12 and 30 degrees, respectively.

34

S. Filiz, O.B. Ozdoganlar / Precision Engineering 35 (2011) 24–37

Fig. 12. The effect of varying tip diameter on the natural frequencies of the micro-endmill. (a) First bending mode, (b) second bending mode, (c) third bending mode, and (d) first torsional/axial mode. Solid lines represent the first direction in the bending pair, dashed lines represent the second direction in the bending pair.

The fluted section cross-sectional geometry of T2 is used in the analysis. The mode-shapes in Fig. 9 show that the first bending mode is arising from the shank of the micro-endmill (as observed from the change in slope along the shank section). Since the total length and diameter of the micro-endmill are kept constant, increasing aspect ratio causes the fluted-section length to increase and the shank length to decrease. Increasing the aspect ratio of the fluted section (and hence its length) does not significantly reduce the modal stiffness of the first mode since the stiffness of the first mode is dominated by the shank stiffness. However, increasing tip aspect

ratio leads to a reduction in the first mode’s modal mass since it reduces the shank length. Consequently, the first bending natural frequency monotonously increase with increasing tip aspect ratio as observed from Fig. 13(a). Fig. 13(b) and (c) shows that the natural frequencies first increase then decrease. The mode-splitting of the second and third modes shows that these modes are dominated by the tip deformations for tip aspects ratios larger than 4 and 2.5, respectively. When the tip deformation is not dominant, the increase in the aspect ratio causes a reduction in the modal mass, leading to an increase in natural frequencies. When the tip deformation is dominant, the

Fig. 13. The effect of varying tip aspect ratio on the natural frequencies of the micro-endmill. (a) First bending mode, (b) second bending mode, (c) third bending mode, and (d) first torsional/axial mode. Solid lines represent the first direction in the bending pair, dashed lines represent the second direction in the bending pair.

S. Filiz, O.B. Ozdoganlar / Precision Engineering 35 (2011) 24–37

35

Fig. 14. The effect of varying tip helix angle on the natural frequencies of the micro-endmill. (a) First bending mode, (b) second bending mode, (c) third bending mode, and (d) first torsional/axial mode. Solid lines represent the first direction in the bending pair, dashed lines represent the second direction in the bending pair.

level of mode-splitting. The effect of helix angle on the torsionalaxial mode is found to be negligible as this mode is dominated by the deformation of the shank section (and the modal mass of the fluted section does not vary with helix angle).

8. Summary and conclusions

Fig. A.1. A fourth order polynomial mapping of the actual cross-section to a rectangular domain.

increase in the aspect ratio decreases the modal stiffness, leading to a decrease in natural frequencies. Similar to the tool-diameter analysis, the torsional/axial mode monotonously increase with increasing aspect ratio since increasing aspect ratio causes a bigger reduction in the modal mass than in the modal stiffness. 7.3. The effect of helix angle To study the effect of helix angle on the micro-endmill dynamics, the helix angle is varied from 0◦ to 50◦ . The tip diameter (dt ) and fluted section aspect ratio (Lt /dt ) of the micro-endmill are selected to be 750 ␮m and 4, respectively. The shank diameter and total length are kept constant at T2 values. Taper angles and helix angle are selected as 12◦ and 30◦ , respectively. The cross-sectional and fluted section geometry of T2 is used in this analysis. Fig. 14 shows that the pairs of bending modes are seen to show a larger split at low helix angles, with a maximum split at the helix angle of zero. For the helix angle of zero, the non-axisymmetry is maximum leading to the maximum mode-splitting. As the helix angle increases, the total angle of twist also increase, and the effective stiffness along the two principal directions become increasingly uniform, leading to increased symmetry and a reduced

This paper presented a combined one-dimensional/threedimensional model for bending, torsional, and axial vibrations of micro-endmills, while considering the actual cross-section and pretwisted geometry. Shank and taper sections of micro-endmills are modeled using one-dimensional beam models, and the flutedsection is modeled using the three-dimensional model, which is based on three-dimensional linear elasticity. A variational (the extended Hamilton’s principle) approach is used in deriving the integral boundary value problem. The model for the entire microendmill is then obtained using a component mode synthesis technique. For the numerical solution of the models, the spectralTchebychev technique is used to discretize the system equations. The model is the validated against modal tests conducted on two- and four-fluted micro-endmills with different geometries, by comparing the natural frequencies and mode shapes. The natural frequencies are also compared to those from finite-elements models. The validated model is then used to analyze the effects of tool diameter, aspect ratio, and helix angle on micro-endmill dynamics. The following specific conclusions are drawn from this work; • The presented model can accurately capture the bending and torsional/axial natural frequencies of micro-endmills. The experimental validation study showed that the natural frequencies can be obtained with better than 2% accuracy for the range of microendmill geometries tested (both two- and four-fluted), and for all the modes up to 90 kHz. • As seen from the comparison between the experimental and model results, the model also captures the mode shapes successfully. This capability is critical for accurate micromilling simulations.

36

S. Filiz, O.B. Ozdoganlar / Precision Engineering 35 (2011) 24–37

• In modeling the dynamics of micro-endmills with high accuracy, capturing the effects of actual cross-sectional geometry and the pre-twist is critical. As compared to a (complete) one-dimensional model, the accuracy of predictions of the threedimensional model is found to be significantly better. • The analytically based 3D-ST model can capture the dynamic behavior of micro-endmills as accurately as an FE model, while the latter requires much larger computational effort. The natural frequencies predicted from the 3D-ST model and those from the FE model differ less than 0.9% for the analyzed cases. • The effect of geometric and material properties on the dynamics of micro-endmills can be analyzed using the 3D-ST model. The tool diameter and aspect ratio have significant affect in bending modes. Increased tool diameter and aspect ratio also increase the axial/torsional natural frequency. The change in helix angle affects the bending modes, especially the mode splitting in nonaxisymmetric geometries, but does not affect the axial/torsional mode significantly.

The constant A0 is found by setting the polynomial value to 1 for the mapping coordinates of point i (A0 = 1/( i (i , i ))). Appendix B. Physical frame derivative and inner product matrices for the three-dimensional beam First, the mapping domain derivative matrices are extended to three-dimensional sampled domain as i Qce = Qli l 1 c2

1 2

l1 = 1, . . . , Ni ,

l2 = 1, . . . , Ni ,

m = 1, . . . , Nj ,

c1 = (l1 − 1)Nj Nk + (m − 1)Nk + n,

n = 1, . . . , Nk ,

(37)

c2 = (l2 − 1)Nj Nk + (m − 1)Nk + n,



where Qi is (Ni × Ni ) is extended to the derivative matrix Q ei with a size of (Ni Nj Nk × Ni Nj Nk ). In Eq. (37), first subscript represents the row number and second subscript represent the column number such that Qli l is the element of the Q i matrix on the l1th 1 2

row and l2th column. The Jacobian is extended as Jce1 c2 = Jl m n , l = 1, . . . , N , m = 1, . . . , N , n = 1, . . . , N , c1 = c2 = (l − 1)N N + (m − 1)N + n.

Acknowledgments The authors wish to acknowledge Dr. Louis Romero of Sandia National Laboratories and Prof. Jacobo Bielak of Carnegie Mellon University for their invaluable assistance in developing the spectral-Tchebychev solution for the three-dimensional model. This work was supported in part by the National Science Foundation CAREER award CMMI-0547534 (Ozdoganlar). e





Vc1 c2 = Vl

1 l2

The components of the inverse of the Jacobian matrix is extended in a similar way. The physical frame derivative matrices are obtained by multiplying the inverse Jacobian matrix components with the mapping domain derivative matrices (see Eq. (14)). Next, the extended inner product matrix in the mapping frame is computed as



Vm1 m2 Vn1 n2

l1 = 1, . . . , N , l2 = 1, . . . , N , m1 = 1, . . . , N , m2 = 1, . . . , N , n1 = 1, . . . , N , c1 = (l1 − 1)N N + (m1 − 1)N + n1 , c2 = (l2 − 1)N N + (m2 − 1)N + n2 .

n2 = 1, . . . , N ,

(38)

The inner product matrix in the physical reference frame is obtained as V xyz = J e V e . Appendix A. Polynomial mapping

Appendix C. Matrices for the compatibility equations

The physical coordinates can be written in terms of a set of mapping points and mapping polynomials as

The ˇ1 matrix can be expressed as

⎡ eN 1



(m+1)(n+1)

x¯ =



(m+1)(n+1)

x¯ i i ,

y¯ =

i

y¯ i i

(35)

i

where x¯ i and y¯ i are mapping points, m and n are the polynomial order in  and  directions, respectively. The accuracy of the mapping depends on the selection of the mapping points and the order of mapping polynomials [45]. A fourth order mapping (with m = n = 4 evenly spaced lines in mapping coordinates) is used in this analysis (see Fig. A.1), which leads to 25 mapping polynomials. A minimum number of mapping points ((m + 1)(n + 1) = 25) are used in this analysis as shown in Fig. A.1 to capture the cross-section of the micro-endmills. There are m + 1 = 5 and n + 1 = 5 lines in  and  coordinates, respectively. ¯ coordinate frame correspond to the The mapping points in the (¯x, y) points at the intersection of the lines in the mapping coordinates. The mapping points are located such that a one-to-one mapping is obtained. The value of the mapping polynomial is 1 at its corresponding point and 0 at other points. The method of constructing the polynomials is recursive, and explained in Eq. (36) for a point i located at the intersection of lines mi and ni .

⎛

i (, ) = A0 ⎝



m+1

k=1,k = / mi

⎞⎛

( − k )⎠ ⎝

n+1  l=1,l = / ni

⎞

( − l )⎠ .

(36)

Z ⎢ Z11 ˇ1 = ⎢ ⎣ Z1 Z1 Z1

Z1 e1N Z1 Z1 Z1 Z1

Z1 Z1 e1N Z1 Z1 Z1

Z1 Z1 Z1 e1N Z1 Z1

Z1 Z1 Z1 Z1 e1N Z1

Z1 Z1 Z1 Z1 Z1 e1N

−e21 Z2 Z2 Z2 Z2 Z2

Z2 −e21 Z2 Z2 Z2 Z2

Z2 Z2 −e21 Z2 Z2 Z2

Z2 Z2 Z2 −e21 Z2 Z2

Z2 Z2 Z2 Z2 −e21 Z2

Z2 Z2 Z2 Z2 Z2 −e21

⎤ ⎥ ⎥. ⎦

(39)

j

In the ˇ1 matrix, ei represent a vector for the i th section (the first section is the shank and the second section is the taper) whose j th element is unity and all other elements are zeros (N being the number of Tchebychev polynomials used in the expansion), and Zi is a vector for the i th section with all elements zero. The construction of the ˇ2 matrix is illustrated here for the first compatibility equation u − u¯ + y z = 0. The first term of the constraint equation represents the fluted section and the second and third terms represent the taper section. The −u¯ term is inserted into the matrix as ˇ2 c1 c2 = − e2N c

3

d1 =1, d2 =1, l1 =1, . . . , N , m1 =1, . . . , N , k1 =1, . . . , Nz2 , c1 =(d1 −1)N N +(l1 − 1)N + m1 , c2 = (d2 − 1)Nz2 + k1 , c3 = k1 , (40)

where c1 and c2 represents the row and column number of the ˇ2 matrix, respectively, c3 represents element number of the e2N vector, d1 represents the dimensions of the three-dimensional deflections of the fluted section (u, v, w), d2 represents the degree of freedom for the taper section (u¯ is the first degree of freedom), l1 is for the sampling on the fluted section in  direction, m1 is for the sampling

S. Filiz, O.B. Ozdoganlar / Precision Engineering 35 (2011) 24–37

on the fluted section in  direction and k1 is for the sampling on the fluted section in z direction of the fluted section, and Ni represent the number of sampling points along the ith direction. The term y z is inserted into the ˇ2 as ˇ2 c1 c2 =yl

1 m1

e2N c

3

d1 =1, d2 =6, l1 =1, . . . , N , m1 =1, . . . , N , k1 =1, . . . , Nz2 , c1 =(d1 − 1)N N +(l1 − 1)N +m1 , c2 =(d2 − 1)Nz2 + k1 , c3 =k1 , (41) where y is a matrix with the y position of the sampling points on the fluted section. The u term is inserted as ˇ2 c1 c2 = e1

c3

d1 =1, l1 =1, . . . , N , m1 =1, . . . , N , n1 =1, . . . , N , c1 =(d1 − 1)N N +(l1 − 1)N +m1 , c2 =6Nz2 +(d1 − 1)N N N +(l1 − 1)N N +(m1 − 1)N +n1 ,

c3 =n1 , (42)

where n1 is for the sampling on the fluted section in  direction. The second and third compatibility equations are expressed similarly. The dimension (d1 ) is 2 and 3 for the second and third equations, respectively. The degree of freedom parameters (d2 ) are also updated accordingly. The rest of the ˇ2 matrix elements are zeros. References [1] Sutherland J, DeVor R. An improved method for cutting force and surface error prediction in flexible end milling systems. Journal of Engineering for Industry, Transactions of the ASME 1986;108:269–79. [2] Altintas Y, Montgomery D, Budak E. Dynamic peripheral milling of flexible structures. Journal of Engineering for Industry, Transactions of the ASME 1992;114:137–45. [3] Montgomery D, Altintas Y. Mechanism of cutting force and surface generation in dynamic milling. Journal of Engineering for Industry, Transactions of the ASME 1991;113:160–8. [4] Tobias S, Fishwick W. The chatter of lathe tools under orthogonal cutting conditions. Transactions of the ASME 1958;80:1079–85. [5] Tlusty J, Polacek M, p. 465–74 Stability of machine tool against self-excited vibration in machining. In: International Production Engineering Research Conference – Proceedings. 1963. [6] Smith S, Tlusty J. An overview of modeling and simulation of the milling process. Journal of Engineering for Industry, Transactions of the ASME 1991;113:169–75. [7] Schmitz T, Ziegert J, Canning J, Zapata R. Case study: a comparison of error sources in milling. Precision Engineering 2008;32(2):126–33. [8] Filiz S, Conley C, Wasserman M, Ozdoganlar O. An experimental investigation of micro-machinability of copper 101 using tungsten carbide micro-endmills. International Journal of Machine Tools and Manufacture 2007;47(7–8):1088–100. [9] Schmitz TL, Duncan GS. Receptance coupling for dynamics prediction of assemblies with coincident neutral axes. Journal of Sound and Vibration 2006;289(4–5):1045–65. [10] Ema S, Fujii H, Marui E. Whirling vibration in drilling. Part 3. Vibration analysis in drilling workpiece with a pilot hole. Journal of Engineering for Industry, Transactions of the ASME 1988;110(4):315–21. [11] Liu X, Jun MBG, Devor RE, Kapoor SG, p. 583–92 Cutting mechanisms and their influence on dynamic forces, vibrations and stability in micro-endmilling. In: Proceedings of IMECE. 2004. [12] Sutherland JW. A dynamic model of the cutting force system in the end milling process, American Society of Mechanical Engineers. Production Engineering Division (Publication) PED 1988;33:53–62. [13] Jun MB, Kapoor SG, Devor RE. The effects of end mill alignment errors on vibrations at high spindle speeds. Transactions of NAMRI/SME 2004;32:9–16. [14] Gupta K, Ozdoganlar OB, Kapoor SG, DeVor RE. Modeling and prediction of hole profile in drilling. Part I. Modeling drill dynamics in the presence of drill alignment errors. Journal of Manufacturing Science and Engineering, Transactions of the ASME 2003;125:6–13. [15] Kivanc EB, Budak E. Structural modeling of end mills for form error and stability analysis. International Journal of Machine Tools and Manufacture 2004;44(11):1151–61.

37

[16] Magrab EB, Gilsinn DE. Buckling loads and natural frequencies of drill bits and fluted cutters. Journal of Engineering for Industry, Transactions of the ASME 1984;106(3):196–204. [17] Rincon DM, Ulsoy AG. Complex geometry, rotary inertia and gyroscopic moment effects on drill vibrations. Journal of Sound and Vibration 1995;188(5):701–15. [18] Gong Y, Ehmann KF, Lin C. Analysis of dynamic characteristics of micro-drills. Journal of Materials Processing Technology 2003;141:16–28. [19] Jun MBG, Liu X, Devor RE, Kapoor SG. Investigation of the dynamics of microend milling. Part I. Model development. Journal of Manufacturing Science and Engineering, Transactions of the ASME 2006;128:893–900. [20] Jun MBG, Liu X, Devor RE, Kapoor SG. Investigation of the dynamics of microend milling. Part II. Model validation and interpretation. Journal of Manufacturing Science and Engineering, Transactions of the ASME 2006;128:901– 12. [21] Huang B. The drilling vibration behavior of a twisted microdrill. Journal of Manufacturing Science and Engineering, Transactions of the ASME 2004;126(4):719–26. [22] Filiz S, Ozdoganlar OB, Romero LA. An analytical model of microendmill dynamics. Journal of Vibration and Controls 2008;14(8):1125– 50. [23] Filiz S, Ozdoganlar OB. Micro-endmill dynamics including the actual fluted geometry and setup errors. Part I. Model development and numerical solution. Journal of Manufacturing Science and Engineering, Transactions of the ASME 2008;130, 031119–1–10. [24] Filiz S, Ozdoganlar OB. Micro-endmill dynamics including the actual fluted geometry and setup errors. Part II. Model validation and application. Journal of Manufacturing Science and Engineering, Transactions of the ASME 2008;130, 031120–1-13. [25] Yagci B, Filiz S, Romero LA, Ozdoganlar OB. A spectral-Tchebychev technique for solving linear and non-linear beam equations. Journal of Sound and Vibration 2009;321(1–2):375–404. [26] Timoshenko SP, Goodier JN. Theory of elasticity. New York: McGraw-Hill; 1970. [27] Leipholz H. Theory of elasticity. Leyden: Noordhoff International Publishing; 1974. [28] Rosen A. Theoretical and experimental investigation of the nonlinear torsion and extension of initially twisted bars. Journal of Applied Mechanics, Transactions of the ASME 1983;50:321–6. [29] Hodges DG. Torsion of pretwisted beams due to axial loading. Journal of Applied Mechanics, Transactions of the ASME 1980;47:393–7. [30] Meirovitch L. Principles and techniques of vibrations. New Jersey: Prentice-Hall Inc.; 1997. [31] Rao SS. Vibration of continuous systems. Hoboken, New Jersey: John Wiley and Sons Inc.; 2007. [32] Zienkiewicz OC, Taylor RL, Zhu JZ. The finite element method: its basis and fundamentals. Burlington: Elsevier Butterworth-Heinemann; 2005. [33] Filiz S, Ozdoganlar OB. Experimental modal analysis of micro-drills. Transactions of NAMRI/SME 2008:185–92. [34] Lee SW, Mayor R, Ni J. Dynamic analysis of a mesoscale machine tool. Journal of Manufacturing Science and Engineering, Transactions of the ASME 2006;128:194–203. [35] Mascardelli BA, Park SS, Freiheit T, p. 145–50 Substructure coupling of microend mills. In: Proceedings of IMECE2006, Chicago, IL, USA (13129). 2006. [36] Ozdoganlar OB, Hansche BD, Carne TG. Experimental modal analysis for microelectromechanical systems. Experimental Mechanics 2005;45(6):498– 506. [37] Erturk A, Ozguven HN, Budak E. Analytical modeling of spindle-tool dynamics on machine tools using Timoshenko beam model and receptance coupling for the prediction of tool point FRF. International Journal of Machine Tools and Manufacture 2006;46:1901–12. [38] Jackson ME, Hyde LJ, Robinson GM, Ahmet W. Dynamic response of a tetrahedral nanomachining machine tool structure. International Journal of Nanomanufacturing 2006;1(1):26–46. [39] Schmitz TL, Duncan GS. Three-component receptance coupling substructure analysis for tool point dynamics prediction. Journal of Manufacturing Science and Engineering, Transactions of the ASME 2005;127:781–90. [40] Filiz S, Cheng C, Powell KB, Schmitz TL, Ozdoganlar OB. An improved tool-holder model for RCSA tool-point frequency response prediction. Precision Engineering 2009;33:26–36. [41] Chou YF, Wang LC. On the modal testing of microstructures: its theoretical approach and experimental setup. Journal of Vibration and Acoustics 2001;123:104–9. [42] Montalvao JM, Silva E, Araujo Gomes AJM. Experimental dynamic analysis of cracked free-free beams. Experimental Mechanics 1990;30(1):20–5. [43] Duncan GS, Tummond MF, Schmitz TL. An investigation of the dynamic absorber effect in high-speed machining. International Journal of Machine Tools and Manufacture 2005;45:497–507. [44] Schmitz TL, Davies MA, Medicus K, Snyder J. Improving high-speed machining material removal rates by rapid dynamic analysis. CIRP Annals—Manufacturing Technology 2001;50(1):263–8. [45] Becker EB, Carey GF, Oden JT. Finite elements. vol. 1. An introduction. New Jersey: Prentice-Hall; 1986.