MECHANICS AND DYNAMICS OF THIN WALL MACHINING by MUHITTIN CANER EKSIOGLU B.Sc., Istanbul Technical University, Turkey, 2007 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in The Faculty of Graduate Studies (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April 2011 © Muhittin Caner Eksioglu, 2011 ii Abstract Peripheral milling of thin walled aerospace components takes considerable amount of machining time as blank blocks are cut down to thin webs under excessive structural oscillations during the process. Unstable chatter vibrations and stable forced vibrations cause poor surface finish on the machined part. Predicting the process me- chanics in advance eliminates the time consuming trial and error approach in reducing the vibrations which are within the tolerance limits of the part. This thesis presents the mathematical modeling of the peripheral milling of the thin walls with helical end mills. The cutting forces, vibrations and dimensional form errors left on the finish surface are predicted under stable but forced vibration conditions. The chatter stability diagram of the operation is predicted by using both frequency and semi-discrete time domain models. The relative vibrations between the flexible part and slender end mill are consi- dered. The tool and the workpiece are discretized along the contact axis to include effect of varying cutting forces and structural dynamics. The differential milling forces are evaluated from the static chip loads contributed by the rigid body motion of the milling operation, and dynamic chip loads caused by the relative vibrations between the flexible tool and flexible thin part. The different cylindrical end mill geometries with regular and non-uniform pitch and helix angles, and low speed process damping effects are included in the dynamic force model. The dynamic properties of the flexible structures are represented by expe- rimentally evaluated modal model in order to reduce the number of linear, periodic, delayed differential equations solved in frequency and time domain computations. The periodic, delayed differential equations are solved by the semi discrete time domain method to predict the amplitude of vibrations and forces. The equations of motion are simplified to constant coefficient type ordinary differential equations, and surface location errors are calculated by frequency domain solver. iii Chatter stability lobes are calculated using semi discrete time domain and fre- quency domain methods. Chatter stability solvers are validated by conducting chatter tests for roughing and finishing stages of thin walled aluminum part at high cutting speeds, and low speed machining of rigid steel block. iv Table of Contents Abstract ............................................................................................................................ ii Table of Contents ............................................................................................................ iv List of Tables ................................................................................................................... vi List of Figures ................................................................................................................ vii Nomenclature ................................................................................................................... x Acknowledgements ....................................................................................................... xiv 1. Introduction................................................................................................................ 1 2. Literature Review ...................................................................................................... 3 2.1 Overview ............................................................................................................ 3 2.2 Static Forces in Milling ...................................................................................... 3 2.3 Surface Location Error ....................................................................................... 5 2.4 Chatter Stability Prediction Models for Intermittent Milling Operation ........... 9 2.5 Process Damping in Milling ............................................................................ 15 3. Dynamics of Peripheral Milling ............................................................................. 17 3.1 Introduction ...................................................................................................... 17 3.2 Discretized Representation of Process Dynamics ............................................ 18 3.2.1 Forces in Milling .................................................................................. 18 3.2.2 Equations of Motion in Laplace Domain ......................................... 32 3.3 Modal Space Formulation ................................................................................ 35 4. Time Domain Simulation ........................................................................................ 46 4.1 Introduction ...................................................................................................... 46 4.2 Semi Discrete Time Domain Solution ............................................................. 46 4.2.1 SD Time Domain Solution Example.................................................... 53 4.3 Surface Location Error ..................................................................................... 58 v 4.4 SLE Example ................................................................................................... 64 5. Chatter Stability in Milling ..................................................................................... 71 5.1 Introduction ...................................................................................................... 71 5.2 Semi Discrete Time Domain Chatter Stability................................................. 71 5.3 Frequency Domain Chatter Stability ................................................................ 73 5.4 Chatter Tests..................................................................................................... 77 5.4.1 Thin Wall Roughing Test for Chatter Analysis ................................... 77 5.4.2 Thin Wall Finishing Test for Chatter Analysis .................................... 88 5.4.3 Low Speed Cutting Test for Rigid Workpiece ..................................... 94 5.4.4 Summary of Experiments ..................................................................... 96 6. Conclusions ............................................................................................................... 97 Bibliography ................................................................................................................... 99 Appendix A: Single Point Formulation of Milling Equation ................................... 105 Appendix B: Fourier Coefficients for Time Dependent Force Terms .................... 110 Appendix C: Variable Helix Cutter Representation ................................................ 113 Appendix D: Accelerometer Mass Elimination ........................................................ 115 vi List of Tables Table 4-1: Parameters for regular cylindrical end mill .................................................... 54 Table 4-2: Modal parameters of the regular cylindrical end mill .................................... 54 Table 4-3: Least common period and frequency for different cutters ............................. 61 Table 4-4: Modal parameters of Al 7075 plate along its y direction ............................... 65 Table 4-5: Parameters of cylindrical cutter ..................................................................... 66 Table 5-1: Least common period and frequency for different cutters ............................. 72 Table 5-2: Modal Parameters of the 10.5 mm thick plate ............................................... 79 Table 5-3: Modal parameters for tool and workpiece ..................................................... 82 Table 5-4: Tested conditions and results ......................................................................... 83 Table 5-5: Modal Parameters of the 5 mm thick plate .................................................... 90 Table 5-6: Modal parameters of the tool ......................................................................... 91 Table 5-7: Tested conditions for finishing operation ...................................................... 92 Table 5-8: Modal parameters at the tool tip ..................................................................... 94 vii List of Figures Figure 2-1: Down milling and up milling operations ........................................................ 3 Figure 2-2: Cutter and workpiece deflection ..................................................................... 5 Figure 2-3: Schematic model of dynamic milling ............................................................. 7 Figure 2-4: An example on surface location error, a) chatter stability diagram for the 10% radial immersion down milling, b) the surface location errors in micrometers are calculated along the axial depth of cut (ADOC) for the fixed spindle speed indicated in a) ........................................................................................................................................ 9 Figure 2-5: Low immersion down milling operation with geometric parameters of an end mill ................................................................................................................................... 10 Figure 2-6: Cutting condition parameters that affect the intermittency of the engagement for down milling operation; a,c,e,g are cutting forces in x and y directions, and b,d,f,h are Fast Fourier Transforms of those forces .................................................................... 11 Figure 2-7: Regeneration mechanism .............................................................................. 12 Figure 2-8: Chatter stability lobe and movement of eigenvalues in discrete domain for points A (Hopf bifurcation) and B (Flip bifurcation) ...................................................... 14 Figure 2-9: Tool flank’s indentation into wavy surface in low cutting speeds ............... 15 Figure 3-1: a) Axially discretized general end mill geometry, b) milling forces and angle convention ........................................................................................................................ 19 Figure 3-2: Milling Forces Scheme ................................................................................. 21 Figure 3-3 : Demonstration of dynamic and static chip thicknesses ............................... 23 Figure 3-4 : Indentation of flank due to edge radius and undulated surface ................... 25 Figure 3-5: Indentation of flank at several positions of the wavy surface ....................... 28 Figure 3-6 : Discretized representation of the plate and the end mill .............................. 32 Figure 3-7: Relative displacement for axial point r ....................................................... 41 Figure 4-1: Approximation of the delayed term by discretization .................................. 49 Figure 4-2: Discretization of an element in dynamic modal force component matrix with different number of points ............................................................................................... 50 viii Figure 4-3: a) Stability chart is taken from [59] (with the permission of Merdol) b)-d) time domain simulation in x and y directions respectively for 7 mm ADOC and 6800 rpm c)-e) simulations in x and y directions respectively for 7 mm ADOC and 7500 rpm ......................................................................................................................................... 56 Figure 4-4: Resultant cutting forces on x-y plane a) unstable cutting b) stable cutting .. 57 Figure 4-5: FFT of force data a) unstable cutting b) stable cutting ................................. 58 Figure 4-6: Surface location error a) undercut example b) overcut example .................. 59 Figure 4-7: Flowchart of SLE method ............................................................................. 60 Figure 4-8: Front view of the Al 7075 plate, with four measurement points indicated along its right corner ........................................................................................................ 65 Figure 4-9: Chatter stability for low immersion down milling operation and the region for SLE calculation is indicated with dashed line. .......................................................... 67 Figure 4-10: a) Force in y direction for one spindle rotation, b) Fast Fourier Transform (FFT) of the force signal .................................................................................................. 68 Figure 4-11: Fast SLE calculation algorithm ................................................................... 69 Figure 4-12: Three dimensional SLE graph for low immersion down milling of thin plate and SLE at 16580 rpm ..................................................................................................... 70 Figure 5-1: Nyquist plot of the characteristic equation ................................................... 76 Figure 5-2: Side and front views of the plate with dimensions and measurement locations ......................................................................................................................................... 77 Figure 5-3: Frequency response functions of workpiece at different locations and tool tip a) first bending mode of the workpiece, 6 611,wp 12,wp1.93 10 N/m, 2.18 10 N/m,k k 6 6 13,wp 14,wp2.61 10 N/m, 2.74 10 N/mk k b) first torsional mode of the workpiece, 6 9 6 21,wp 22,wp 23,wp7.52 10 N/m, 2.34 10 N/m, 8.58 10 N/m,k k k 9 24,wp 1.92 10 N/mk c) second bending mode of the workpiece 831,wp 1.97 10 N/m,k 8 31,wp 1.97 10 N/m,k 7 32,wp 6.90 10 N/m,k 8 8 33,wp 34,wp2.40 10 N/m, 1.67 10 N/mk k d) tool’s FRFs at x and y directions .......................................................................................................................... 81 Figure 5-4: Setup for stability test ................................................................................... 82 Figure 5-5: Predicted and measured chatter stability results for the most flexible condition .......................................................................................................................... 84 ix Figure 5-6: Y direction measurement and simulation for point A, a) Force measurement b) FFT of force measurement c) simulated vibration displacement d) FFT of displacement .................................................................................................................... 85 Figure 5-7: Y direction measurement and simulation for point B, a) Force measurement b) FFT of force measurement c) simulated vibration displacement d) FFT of displacement .................................................................................................................... 86 Figure 5-8: Y direction measurement and simulation for point C, a) Force measurement b) FFT of force measurement c) simulated vibration displacement d) FFT of displacement .................................................................................................................... 88 Figure 5-9: Side and front views of the plate with designated modal test points ............ 89 Figure 5-10: Predicted stability lobes for the most flexible case and measured points are presented for 5% radial immersion thin wall machining ................................................. 93 Figure 5-11: Two orthogonal direction FRFs measured at the tool tip ........................... 95 Figure 5-12: Chatter stability prediction using different algorithms and test results ...... 96 Figure A-1: Resultant force lumped on point k ............................................................. 105 Figure C-1 : Crest cut end mill with varying helix angle along the flute ...................... 114 Figure D-1: Overlayed measurements indicating the torsional mode of the structure .. 115 Figure D-2: Curve fit for the flexible top right corner, the pink line is the fit curve and the blue line is the measurement. The discrepancy at torsional mode at 1720 Hz is apparent. ......................................................................................................................... 116 Figure D-3: The most of the frequencies are 1719 Hz and variation is around 2 Hz. ... 118 Figure D-4: The curve is fitted from the direct FRF and since there is no shift, the errors are less at the top corners of the plate. ........................................................................... 119 x Nomenclature a Axial depth of cut 1, ˆ jA Differential static cutting force vector of flute j 2, ˆ jA Differential dynamic cutting force vector of flute j 1, ˆ jB Differential static edge force vector of flute j 2, ˆ jB Differential process damping force matrix of flute j c Feed per tooth ,dA dV Indented area and volume of the flank face dS Length of differential cutting edge dz Differential axial depth of cut D Tool diameter f Force vector in modal force representation r t a, ,F F F Radial, tangential and axial milling forces x y z, ,F F F Milling forces in Cartesian coordinates G Displacement transfer function matrix h Fourier coefficient index jh Static chip thickness of flute j d,jh Dynamic chip thickness of flute j xi Η Truncated Toeplitz matrix ji Helix angle of flute j I Identity matrix j Tooth index rc tc ac, ,K K K Radial, tangential and axial milling cutting force coefficients re te ae, ,K K K Radial, tangential and axial milling edge force coefficients spK Material dependent process damping indentation coefficient RK Residual stiffness L State matrix for current time terms wL Average flank wear length t wp,m m Total number of modes of tool and workpiece respectively M Transfer function matrix in modal domain RM Residual mass n Spindle speed in rpm N Number of flutes q Differential physical displacement vector in time domain ( )sQ Physical displacement vector in Laplace domain containing all axial levels R State matrix for delayed terms xii S Input matrix(forces) in state space representation t Time step , jT T Time delay, time delay of flute j T Transformation matrix for Cartesian coordinates u Mass normalized mode shape U Mass normalized mode shape matrix cv Cutting velocity of the edge in cut V Indented volume of the flank face z Axial elevation from the tip of the tool Cutter rotation angle j Instantaneous immersion angle of flute j p Pitch angle Φ State transition matrix Modal coordinates vector T Delayed modal coordinates vector γ Mode index Axial immersion angle Coulomb friction coefficient Θ Motion vector containing modal displacements and velocities xiii Least common period c Chatter frequency n Natural frequency p Least common frequency s Spindle frequency in rad/s ψ j Lag angle of flute j Ψ Relative vibrations vector Damping ratio T[ ] Transpose operator det Determinant operator ADOC Axial depth of cut DDE Delayed differential equation DOF Degree of freedom MF Multi-frequency solution SD Semi discretization SLE Surface location error xiv Acknowledgements I would like to sincerely express my gratitude to my research supervisor Prof. Yusuf Altintas for the encouragement, support and invaluable guidance he has provided throughout my study at UBC. I would also like to thank my co-supervisor Dr. Doruk Merdol for his continuous support and help. I would like to express my appreciation to our sponsors NSERC-Pratt & Whitney Canada for the Industrial Research Chair Grant, Machine Tool Technologies Research Foundation (MTTRF) for letting us use the Mori Seiki NMV5000 CNC Machining Center, Sandvik Coromant and Data Flute for supplying the tools we used in the machin- ing tests. It was a great privilege to be a member of an active research group in Manufac- turing Automation Laboratory (MAL). The discussions I have made with my colleagues and their share of knowledge elaborated this research work. It was also a great fun to be with the brilliant people that came from different parts of the world. Finally, I would like to thank my parents, Nazli and Mahmut Eksioglu. I felt their presence all the time with the support, encouragement, friendship and love they have provided from thousands of miles away. I dedicate this work to them. Chapter 1. Introduction 1 1. Introduction Peripheral milling of thin walled components is a common process for aircraft manufacturers. Wing structures, fuselage sections, jet engine compressors and turbine blades are machined from blank blocks to thin webs with slender end mills to meet tight tolerances, achieve complex geometries, and obtain acceptable thermo-mechanical properties. Machining of these parts require large amount of material removal, hence affecting the overall production time and cost. The problems that emerge due to high flexibility of the cutting tool-workpiece system limits the productivity, especially when trial and error based approaches are employed to set cutting conditions. This has moti- vated researchers to understand the physics of the process, and optimize the parameters to attain higher production rates, and achieve good surface finish quality. The relative vibrations between flexible end mill and workpiece are the source of most of the problems that occur during manufacturing of thin webs. Intermittent en- gagement of cutter and workpiece excites a wide range of structural natural frequencies that result as unstable chatter vibrations and stable forced vibrations. Effect of cutting conditions, tool geometry, structure dynamics, material characteristics, high speed and low speed machining properties have been studied separately to minimize process vibrations in order to reduce operation cycle times. Most of the available models aim specific application, or general models consume sizable computation time which is not convenient for machining industry. This thesis investigates a general mathematical model for thin wall machining, which is computationally efficient, and has acceptable prediction accuracy within the physical assumptions made. The thesis is organized as follows; Chapter 2 provides the review of existing research on milling mechanics and ef- fect of vibrations on cutting. Brief overview on calculation of static milling forces, the models to predict errors on the finished surface, chatter stability calculation of intermit- tent milling operation, and process damping effect at low cutting speeds are presented in this chapter. Chapter 1. Introduction 2 The dynamics of peripheral milling is discussed in Chapter 3 using axially dis- cretized representation. Static milling forces created at each axial level regarding the general end mill geometry are summarized; furthermore the forces emerge from dynamic interactions, such as regeneration and process damping, are included to form general force model. Workpiece and cutter dynamics are presented in axially discretized form which is organized to retrieve the data from experimental modal analysis. Consequently the interactions between milling forces and structure dynamics are expressed via equa- tion of motion in modal domain. The modal equation of motion obtained in Chapter 3 has been solved by semi discrete time domain and frequency domain surface location error algorithms in Chapter 4. The semi discrete time domain method solves linear or linearized periodic coefficient delayed differential equation for given spindle speeds, axial and radial depths. This solution is then compared to time domain solution of CUTPRO to analyze its weak- nesses and strengths. For stable cutting conditions, a fast algorithm to obtain surface location errors utilizing the pace of frequency domain solvers is discussed, and a theoret- ical example is presented. Chapter 5 is dedicated to chatter stability prediction, and experimental validation of theoretical models. Semi discrete chatter stability solution is presented along with frequency domain solver that uses Nyquist stability criterion. High speed milling chatter tests are conducted for machining thin walled part to analyze change of structural dynamics, and to test solvers’ accuracy at different cutting conditions. Finally, a rigid steel block is machined to analyze frequencies dominant at low cutting speeds and verify the process damping model used. The thesis is concluded in Chapter 6, summarizing the contributions made and the future research directions. Chapter 2. Literature Review 3 2. Literature Review 2.1 Overview Milling of thin walled components is a common application for aerospace industry, and it is becoming more common for automotive, computer hardware and bioengineer- ing industries. The problems caused by vibrations and static deflections have been studied in details by previous researchers. In this chapter literature survey will be presented on the foundations such as modeling static milling forces, surface location errors, chatter stability prediction models for intermittent engagement, and process damping at low speed cutting. 2.2 Static Forces in Milling In milling, the cutting edge is discontinuously engaged with the workpiece along its periphery. Two different milling techniques are shown in Figure 2-1, which are down milling an up milling. In down milling, chip thickness has the maximum value at the start of the engagement and it reduces to zero as cutter rotates, and in up milling it is the opposite. The amount of material tool removes, hence the forces at each instant vary in time as indicated in Figure 2-1. This brings time dependent characteristic to milling operation unlike turning. Figure 2-1: Down milling and up milling operations Chapter 2. Literature Review 4 Predicting the forces created in milling is one of the most important tasks for re- searchers since problems as surface location errors and tool breakage are related on milling forces. Modeling the time dependent chip formation is the first step to calculate cutting forces in milling. The chip load is dependent on kinematics of milling, and fluctuations that are coming from structure dynamics. In this section, milling kinematics is surveyed, which is the foundations of static milling forces. Martelotti considered both rotating and translating motions of the tool and as- sumed a trochoidal tool path [1]. For planar milling, he expressed that when tool radius is comparatively larger than the feed value, circular tool path approximation has negligi- ble error, thus the formula for static chip thickness simplifies to: ( ) sin( )h c , (2.1) where h is the instantaneous chip thickness, c is feed rate per tooth, and is instantane- ous immersion angle of the cutting tooth. In order to eliminate excessive experimentation, analytical and semi analytical methods have been searched in the past [2]. Experimentally identified specific cutting pressures are presented in the work of Koeingsberger and Sabberwal [3]: F Kah , (2.2) where specific cutting pressure is expressed as xK Ch . Axial depth of cut (ADOC) is denoted by a . The empirical constants depending on workpiece and tool properties are denoted by C and x . In this thesis, due to convenience in analytical calculations, the linear edge force model [4] is used for milling force predictions: c eF K ah K a , (2.3) where cK is specific cutting pressure related with the shearing of material and eK is edge force coefficient which accounts the rubbing between tool flank face and work- piece. These constants can be identified experimentally by fitting a line for federate versus force graphs or analytically calculating by orthogonal cutting theory. Chapter 2. Literature Review 5 2.3 Surface Location Error The dimensional accuracy of the finished surface is one of the key objectives of thin wall machining. The static deflections caused by the forces, and relative vibrations between cutting tool and workpiece may result as surface location errors (SLE). Depend- ing on the error values and tolerance limits, another finishing operation can be required for undercuts, and finished part may be scrapped for severe overcuts, thus bringing additional cost to operation. The surface accuracy for milling of thin walled components investigated by Kline et al. [5] considering the relative static deflections between tool and flexible workpiece. Figure 2-2: Cutter and workpiece deflection Chapter 2. Literature Review 6 The tool is considered as a cantilever beam, and the deflection of the tool in y direction as demonstrated in Figure 2-2, is expressed as: y 3 3 2 y,t y ( ) [(CFY ) ( ) 3( ) ( CFY)] 6E F d z z L z L z L I , (2.4) where y,t ( )d z is the deflection of the tool in y direction at point z , yF is the static y cutting force, CFY is the y force center, yI is the moment of inertia of the tool in y direction, which is calculated from 4 y / 64I D , E is the modulus of elasticity and D is the diameter of the tool. In their work the deflection of the workpiece is calculated using finite element analysis software and the total deflection is expressed as the sum: y y,t y,w( ) ( ) ( )d z d z d z , (2.5) where y,w ( )d z and y ( )d z are workpiece and total displacements at point z . The cutter leaves the error marks when it enters to cut in up milling, and when it exits the cut in down milling at surface generation point. The force and force center are calculated for that instant, the total deflection is obtained by applying Eq. (2.5). Including the effect of helix, the authors plotted the deflections along the z axis. Sutherland and DeVor [6], introduced an improved model where effect of system deflections on the chip load were taken into account. They claimed that previous rigid force models, which calculated the cutting forces neglecting the system deflections, are over-estimating the errors for end milling thin walled cavities as much as 100 to 150 percent. Montgomery and Altintas [7] proposed a dynamic milling model where dynamics of the process is included. The workpiece’s and tool’s dynamics are represented with discrete equivalent masses ( tm , wm ), stiffnesses ( tk , wk ) and damping coefficients ( tc , wc ) as shown in Figure 2-3. The interaction between workpiece and the tool modeled through cutting stiffness cutk , and process damping cutc . Chapter 2. Literature Review 7 Figure 2-3: Schematic model of dynamic milling The following equations of motions of cutter and workpiece are solved for trochoidal motion of the tool through numerical integration iteratively by updating the chip load with vibrations: t t t t t t w w w w w w m q c q k q F m q c q k q F , (2.6) where sub-indices t and w indicate tool and workpiece, respectively. T(x,y)q is the displacement vector, T x y( , )F FF is the force vector, m,c,k are mass, damping and stiffness matrices, in order. After solving the differential equations in time domain for displacement, they obtained the surface errors similar to previous papers. Budak [8] [9] investigated the problem further for thin wall milling for static form errors and used two models which are static deflections and flexible milling model. The static deflections model, similar to Sutherland’s approach, considered change of chip load with static deflections and then plotted the form errors with respect to feed and axial depth. The flexible milling model on the other hand, also introduced change of cutter engagement boundaries; hence change of radial immersion with static deflections. He showed that flexible milling force model matches better with the test measurements, although bring- ing additional computational effort. Both models excluded effect of structural modes on the form errors. Frequency domain analysis of surface location errors has taken attention in recent studies due to fast calculation of steady state vibrations, and convenience at using measurement frequency response function without necessity to identify modal parame- Chapter 2. Literature Review 8 ters. Schmitz and Mann [10] presented a frequency domain SLE approach which only considers forced vibrations at stable cutting conditions without the regeneration effects. The Fourier coefficients of static cutting force are expressed as: 2 -i 0 1 ( )e 2 h hF F d , (2.7) where hF is the Fourier coefficient of cutting force F , ( , )h is the coefficient index, and is the instantaneous immersion angle of the cutting edge. The displacement in steady state is defined as: s i ( 2 / ) 1 ( ) ( ω) e N h t k N h k h t h q G F , (2.8) where displacement vector is T(x,y)q , the Fourier coefficient vector is T x, y,( , )h h hF FF , the number of flutes is N ,and s is the spindle speed in rad/s. The transfer function at the tool tip is expressed as a matrix with direct and cross terms: xx xy yx yy ( ω) ( ω) ( ω) ( ω) ( ω) G h G h h G h G h G . (2.9) At the instant when tool leaves mark on the surface, the location error is stored. The location errors at fixed spindle speed with respect to axial depth of cut (ADOC) are illustrated in Figure 2-6b. The authors [11] found a stable cutting point in a stability pocket as indicated with a dot in Figure 2-6a. Excluding the regeneration forces, cutter runout, effect of vibrations in chip thickness and radial immersion, they calculated the SLE values. Only tool tip transfer functions are taken into account, therefore effect of changing dynamics along the axial depth is neglected. A similar study has taken place helical end mills with semi discretization method [12]. Chapter 2. Literature Review 9 a) b) Figure 2-4: An example on surface location error, a) chatter stability diagram for the 10% radial immersion down milling, b) the surface location errors in micrometers are calculated along the axial depth of cut (ADOC) for the fixed spindle speed indicated in a) 2.4 Chatter Stability Prediction Models for Intermittent Milling Operation Intermittent cutting forces often occur in thin wall milling applications, especially in the finishing stage, where very low immersion rates cannot be compensated by increas- ing the flute number and/or helix angle. The definition of intermittent cutting is having pulse like cutting forces which excites the structure not only in its tooth passing frequen- cy but also at multiple harmonics of the tooth passing frequency. The intermittency of cutting condition is dependent on the immersion rate of the cutting tooth, number of flutes N , and helix angle i . Chapter 2. Literature Review 10 Figure 2-5: Low immersion down milling operation with geometric parameters of an end mill Immersion rate in milling is defined based on the ratio between radial depth of cut (RDOC) and cutter diameter D . As the ratio becomes smaller (i.e. less than 10%), the cutting forces instantly have changing values between peak to zero, which excites the structure in a broad frequency range as illustrated in Figure 2-6a. When the immersion rate is increased to 30%, the cutting forces become more continuous as demonstrated in Figure 2-6c. Similarly, helix angle and number of flutes in a cutter are inversely propor- tional to intermittency of the cutting condition. As helix angle increases, the forces become more continuous, the peak values drop and the number of multiple harmonics of tooth passing frequency diminish as in Figure 2-6f. Decreasing the flute number, creates impulsive forces, thus the number of frequencies swept becomes broader. Chapter 2. Literature Review 11 Parameters Forces FFT of Forces 10% immersion, 0 helix angle, 4N a) b) 30% immersion, 0 helix angle, 4N c) d) 10% immersion, 45 helix angle, 4N e) f) 10% immersion, 0 helix angle, 2N g) h) Figure 2-6: Cutting condition parameters that affect the intermittency of the engagement for down milling operation; a,c,e,g are cutting forces in x and y directions, and b,d,f,h are Fast Fourier Transforms of those forces The regenerative vibration, or chattering, is reported as one of the most signifi- cant problems in milling of thin wall components as high flexibility reduces the stability 0 0.01 0.02 0 200 400 600 F o rc e [ N ] Time [s] x direction y direction 0 1000 2000 3000 4000 5000 0 1 2 x 10 5 F o rc e [ N ] Frequency [Hz] 0 0.01 0.02 -200 0 200 400 600 F o rc e [ N ] Time [s] 0 1000 2000 3000 4000 5000 0 1 2 x 10 5 F o rc e [ N ] Frequency [Hz] 0 0.01 0.02 0 200 400 F o rc e [ N ] Time [s] 0 1000 2000 3000 4000 5000 0 1 2 x 10 5 F o rc e [ N ] Frequency [Hz] 0 0.01 0.02 0 200 400 F o rc e [ N ] Time [s] 0 1000 2000 3000 4000 5000 0 1 2 x 10 5 F o rc e [ N ] Frequency [Hz] Chapter 2. Literature Review 12 limits and intermittent cutting condition brings additional instability. The regenerative chatter vibrations are introduced by Tobias [13] and Tlusty [14] in late 1950s. The idea stated that the chip load fluctuates due to current and past vibrations as illustrated in Figure 2-7. When current and past vibrations are out of phase, then it leads to grow in chip thickness and force in a loop which drives the system to instability. Figure 2-7: Regeneration mechanism Most of the studies until 90s were concentrated on turning stability. Minis and Yanu- shevsky [15] proposed a comprehensive model where milling operation has been modeled as a two dimensional time dependent delayed differential equation. They solved the problem in frequency domain by truncating the infinite Fourier series, or namely using truncated Toeplitz representation. They obtained the following characteristic equation to identify the stability limit: c ω t c T 1 det[ (1 e ) (ω ω )] 2 T r kK a ik I H , (2.10) where I is the identity matrix, tK is the tangential cutting force coefficient, c is the chatter frequency that is searched in the algorithm,T is the time delay, ( , )r k are numbers indicating the harmonic index such as , 0, 1, 2,....,r k m , m is the number of Fourier terms considered in truncated Toeplitz representation, r kH is the oriented transfer function matrix that contains transfer functions, and directional milling coefficients, and Chapter 2. Literature Review 13 Tω is the tooth passing frequency. The characteristic equation subsequently solved using Nyquist stability criteria to determine instability with respect to axial depth of cut (ADOC) and spindle speed. Followed by Minis’s work, Altintas and Budak [16] intro- duced fast analytical way to plot stability lobes in frequency domain when only zeroth term of Fourier series considered. Eq. (2.10) has been simplified to: c ω t 1 det[ (1 e ) (ω )] 2 T cK a I H , (2.11) where H is the oriented transfer function of average milling directional coefficients. This approach is named as zero order analytical (ZOA) method. The ZOA method later extended to general helical end mills and non-uniform pitch cutters [17] [18]. Bravo et al. [19], and Thevenot et al. [20] applied ZOA method to extract stability information for thin walled plate machining. Consequent studies showed that in the case of highly intermittent machining, ZOA method cannot account the additional instabilities come from the multiple harmon- ics of the tooth passing frequency. Multi-frequency (MF) solution which includes multiple terms of Fourier series has been investigated by Budak and Altintas [21], and then later Merdol and Altintas [22] extended it by discussing on the additional instability regions. Temporal finite element method (TFEA) proposed by Bayly et al. [23] and semi discretization introduced by Insperger and Stepan [24] are time domain based chatter stability prediction approaches which are also accurate in intermittent cutting conditions. The cause of additional stability regions explicitly investigated by Bayly et al. [23], and then later by Insperger et al.[25]. The instability in point A is associated with Hopf (Quasi-periodic) bifurcations and B with Flip (period doubling) bifurcations in Figure 2-8. In Hopf bifurcations, two complex conjugate eigenvalues leave the unit circle and the frequency in Hz associated with it is expressed as: cQP ω 2π 60 kn f , (2.12) Chapter 2. Literature Review 14 Figure 2-8: Chatter stability lobe and movement of eigenvalues in discrete domain for points A (Hopf bifurcation) and B (Flip bifurcation) where cω is chatter frequency, n is the spindle speed in rpm and k is the harmonics, i.e. 0, 1, 2,...k . The Flip bifurcations are the reasons of additional lobes in intermittent cutting, where only one eigenvalue on the real axis moves out of the unit circle from -1. The frequency of the Flip bifurcations is defined as: P2 120 60 n kn f ’ (2.13) where frequency is related to spindle’s angular speed. Effect of helix angle on the chatter stability in intermittent cutting conditions is modeled in frequency domain using multi frequency solution by Zatarain et al. [26] and chatter islands caused by helix angle are discussed. Similar analysis has been presented by Patel et al. [27] using TFEA method. Non-uniform pitch and non-uniform helix cutters are considered by Sims et al. [28] using semi discretization method. Dombovari et al.[29], also solved stability of serrated cutters by semi discretization. Chapter 2. Literature Review 15 2.5 Process Damping in Milling The heat resistant materials such as titanium and steel alloys that are frequently used in thin walled components can only be cut at low cutting speeds. When machining these materials at low speeds, it is known that the chatter stability limit increases due to additional damping caused by the process. Finding a mathematical model to represent this behavior has been searched significantly, and several theories are created to explain it in physical sense. In this thesis, flank indentation approach is adopted. According to this theory, at low speeds vibration marks imprinted on the surface becomes narrower and hence workpiece flank rubs against it that causes the additional damping as illu- strated in Figure 2-9. Figure 2-9: Tool flank’s indentation into wavy surface in low cutting speeds Sisson and Kegg [30] explained the low speed cutting behavior by analyzing the tool flank rubbing to wavy inner workpiece surface. They qualitatively searched the effects of parameters such as clearance angle γ , feed rate and tool wear. Wu [31] claimed that at low cutting speeds plowing forces are dominant due to indentation of the tool edge into workpiece. Since the edge radius is always present, plowing forces are exerted: Chapter 2. Literature Review 16 x sp y x F K V F F , (2.14) where xF is the friction force in feed direction proportional to material dependent on indentation constant spK and indented volume V ; yF is the normal force which is proportional to dry Coulomb friction constant . Chiou and Liang [32], utilized Wu’s indentation approach and approximated the indented volume when only effect of tool wear is considered for small vibration amplitudes: 2 w c2 L a V x v , (2.15) where wL is the average flank wear length of the cutter, cv is the cutting velocity of the edge in cut, x is the vibration velocity of the surface in radial direction, and a is axial depth of cut. In their work, material dependent indentation coefficient is identified from static indentation tests of the tool edge. Eynian and Altintas [33] applied the same approach to milling problem, and found matching experimental results. In recent studies, indented volume is approximated with more comprehensive models that account clear- ance angle, edge radius, and amplitudes of the wavy surface. Budak and Tunc [34], calculated the indented volume geometrically over a period, and took the average of energy dissipated to determine average process damping coefficient. Their model included edge radius, separation angle due to indentation, and clearance angle contribu- tions to indented volume. Ahmadi and Ismael [35], similar to Budak’s approach, utilized average energy dissipation method and they also included the effect of vibration wave amplitudes imprinted on the surface. They calculated the stability lobes for maximum amplitude and minimum amplitude of vibrations to determine a region that chatter can happen in low cutting speeds. They analyzed this region by changing feed rate, hence the amplitude of vibrations, to validate their model. Chapter 3. Dynamics of Peripheral Milling 17 3. Dynamics of Peripheral Milling 3.1 Introduction Machining thin walled components is a challenge due to undesired static and dy- namic surface location errors and unstable chatter vibrations. In order to mathematically model the problem, the relative motion between flexible workpiece and slender end mill is required. This can be achieved by formulating the periodic milling forces, obtaining the dynamics of the structures, and associating the external loading with the dynamics and the laws of cutting mechanics. In this chapter, the discretized representation of peripheral milling is discussed by considering varying structural properties and cutting mechanics along the cutter axis. The differential milling forces and discrete structural dynamics are presented. The static and dynamic 3-axis cutting forces are formulated using linear edge force model [4][36]. The process damping force is included by considering the tool flank indentation-rubbing approach as presented in the literature [31][32]. Contrary to common methods on modeling the dynamics of the structures via finite element analysis [8][37][19][20], the structural dynamics are obtained through experimental modal analysis. The equations of motion of each discrete point are later collected and subsequently, modal space represen- tation is used to reduce the number of equations in order to shorten the computation time without sacrificing the overall solution accuracy. Finally, modal equations of motion in time domain are expressed for peripheral milling. Chapter 3. Dynamics of Peripheral Milling 18 3.2 Discretized Representation of Process Dynamics In the peripheral milling of thin walled structures, the cutter is engaged with the workpiece along its axis. The end mill geometry, structural dynamics, cutting conditions such as radial depth of cut may vary continuously as axial position changes. Thus, discretizing axial depth of cut into sufficient number of points in order to account changing parameters accurately, and constructing the equation of motion considering the effect of each individual point are necessary. In the following subsections, formulation of 3-axis milling forces and equation of motion in frequency domain with the addition of identified structural parameters are provided in discretized form. 3.2.1 Forces in Milling For ideal axially symmetric end mill, the forces are generated when the tool’s flutes remove material. Engagement conditions depend on cutter geometry, tool motion and selected axial and radial depth cuts. After engagement takes place, process mechan- ics affect the milling forces. In this section, geometrical relations regarding the general end mill, formulated by Engin and Altintas [36], are given followed by the comprehen- sive cutting mechanics model for three axis milling. General end mills may have varying shape, helix, and pitch angle along the cutter axis. In order to assess each cutting edge’s contribution, the cutter axis is discretized into l number of points as demonstrated in Figure 3-1a. The number of flutes is denoted by N , tooth number is indicated by j , and axial elevation from tip of the tool is represented by z . The relation between z and the discrete points is: z kdz , (3.1) where dz is the differential axial depth of cut and k is the axial index of the cutting point. Chapter 3. Dynamics of Peripheral Milling 19 a) b) Figure 3-1: a) Axially discretized general end mill geometry, b) milling forces and angle convention Instantaneous immersion angle, ( , )j t z , is defined to evaluate cutting edge’s position in rotating frame by measuring the angle from fixed y-axis in clockwise direc- Chapter 3. Dynamics of Peripheral Milling 20 tion as shown in Figure 3-1b. Instantaneous immersion angle for tooth j at axial eleva- tion z is : 1 p, 0 ( , ) ( ) ( ) ψ ( ) j j e j e t z t z z , (3.2) where tooth j ’s angular position depends on cutter rotation, ( )t , lag angle, ψ ( )j z , and the summation of the pitch angles, 1 p, 0 ( ) j e e z in which the indexing is one less from the tooth number and zero term is set p,0( ) 0.z As spindle rotates with speed n [rev/min], the angular position of the cutter varies as: 2 ( ) 60 n t t . (3.3) Pitch angle, p, ( )j z , is the angular spacing between adjacent flutes j and 1j at axial elevation z . For three cases stated below, the pitch angles are given: 0p, p, 0 p, p, 1 2 uniform helix and pitch ( ) uniform helix,non-uniform pitch ( ) ( ) ( j j j j j j N z z z z ) non-uniform helix and pitch , (3.4) where 0 p, j is the pitch angle between flutes j and 1j at the tip of the tool. The lag angle ( )j z is the angle between the cutting edge and the end point (at the tool tip) of flute j due to helix of the cutter. It is formulated along the cutting edge for end mills with constant lead or constant helix in [38]. For cylindrical end mills, lag angle is calculated as: 2 tan( ( )) ( ) j j z i z z D , (3.5) where given axial elevation, the cutter diameter is denoted by ( )D z , and the helix angle is represented with ( )ji z . The diameter of the general end mill is formulated in [38] with Chapter 3. Dynamics of Peripheral Milling 21 respect to axial position. For uniform helix cutters, the helix angle stays same for each flute and axial element. However, for non-uniform helix cutters, the helix may vary for each flute and/or axial level and this is illustrated in Appendix C. Radial, tangential and axial directions in rotating frame define the milling forces for general end mills as demonstrated in Figure 3-2. The instantaneous differential forces corresponding to axial position are expressed as: rta rta, 1 ( , ) ( ( , )) ( , ) N j j j d t z g t z d t z F F , (3.6) where T rta, r, t, a,( , ) { F ( , ) F ( , ) F ( , )}j j j jd t z d t z d t z d t zF are differential forces in radial, tangential and axial forces for each flute, respectively. The unit step function ( ( , ))jg t z is used to consider whether the cutting edge is in cut or not. If instantaneous immersion angle ( , )j t z is in between predetermined start st and exit ex angles, the cutting edge is removing material, hence generating forces: Total Milling Force, F A. Cutting Forces, Fc A1. Static Cutting Forces, Fcs A2. Dynamic Cutting Forces (Regeneration), Fcd B. Edge Forces, Fe B1.Static Edge Forces, Fes B2. Dynamic Edge Forces (Process Damping), Fed Figure 3-2: Milling Forces Scheme Chapter 3. Dynamics of Peripheral Milling 22 st ex1 if ( , ) ( ( , )) 0 otherwise j j t z g t z . (3.7) Milling forces consist of shearing of material at the rake face, and edge forces caused by the plowing between flank face of the tool and the finished surface as indi- cated in Figure 3-2. Intended material removal amount and effects coming from structural flexibilities constitute both forces as sum of static and dynamic components. Static part is constructed owing to the general cutting mechanics, treating the workpiece and tool as rigid bodies. In dynamic force formulation, the interaction between structural vibrations and milling forces is considered. The general form in Eq. (3.6) is separated into cutting and edge forces: c e rta rta, rta, 1 ( , ) ( ( , ))( ( , ) ( , )) N j j j j d t z g t z d t z d t z F F F , (3.8) where differential cutting forces c c c c T rta, r, t, a,( , ) { F ( , ) F ( , ) F ( , )}j j j jd t z d t z d t z d t zF and differential edge forces are e e e e T rta, r, t, a,( , ) { F ( , ) F ( , ) F ( , )}j j j jd t z d t z d t z d t zF . Cutting forces are proportional to the chip thickness and cutting force coeffi- cients. Static and dynamic forces at axial elevation z and flute number j are expressed for general end mill respectively, using the linear edge force theory [4] [39]: cs c rta, rta( , ) ( , ) ( )j jd t z h t z dS zF K , (3.9) and cd c rta, rta d,( , ) ( , ) ( )j jd t z h t z dS zF K , (3.10) where cutting force coefficient matrix is c T rta rc tc ac{ }K K KK . Cutting force coeffi- cients can be identified mechanistically from the measured cutting forces by changing feed rate and fitting a linear curve [40], or applying orthogonal to oblique transforma- Chapter 3. Dynamics of Peripheral Milling 23 tion as described in [4]. Differential contact length is denoted by ( )dS z and it is ciated with differential axial depth by: ( ) / sin ( )dS z dz z , (3.11) where ( )z is the axial immersion angle indicating the angle between cutter axis and the normal of the cutting edge at axial elevation z . Axial immersion angle can be calculated for tip, arc and tapered sections of the cutting edge for general end mills as presented in [38]. For cylindrical end mills, it is always ninety degrees which makes differential contact length and differential axial depth equal. Figure 3-3 : Demonstration of dynamic and static chip thicknesses The static chip thickness, ( , )jh t z , and undulations generating the dynamic chip thickness are demonstrated in Figure 3-3. When the feed in x direction is projected on to radial direction of the cutter, static chip thickness is derived [2]: ( , ) sin ( )sin ( , )j jh t z c z t z . (3.12) Chapter 3. Dynamics of Peripheral Milling 24 The feed per tooth [mm/rev/tooth] is indicated by ( )c z , which is constant for cutters with regular pitch and regular helix. For non-uniform pitch and non-uniform helix cutters, it is calculated via the pitch angle between adjacent flutes [17]: p, ( )Feedspeed[mm/min] ( ) 2 j z c z n . (3.13) Substituting Eqs. (3.11) and (3.12) into Eq. (3.9) gives the differential static cutting force for flute j : cs c rta, rta( , ) sin ( , )j jd t z c t z dz F K . (3.14) As milling forces are applied on the flexible end mill and the thin wall, relative displacements occur, predominantly in the principal vibration directions. These vibra- tions cause undulated surface which changes the overall chip thickness. Regarding the regeneration phenomenon introduced by Tobias [41] and Tlusty [14], the dynamic chip thickness is comprised of instantaneous vibrations and vibration marks imprinted on the surface from the previous tooth. The difference between current and previous vibrations at axial elevation z , assuming that the principal vibration directions are in Cartesian coordinates is: ( , ) (( ( )), ) ( , ) x( , ) y( , ) z( , )jt z t T z z t z t z t z t z q q q i j k , (3.15) where T( , ) {x( , ) y( , ) z( , )}t z t z t z t zq is current vibrations vector in displacement form, (( ( )), )jt T z zq is the vibration marks left by the previous tooth. Time delay, indicated by ( )jT z , provides the instant that previous tooth left vibration marks on the surface. Using pitch angle and time passed in one full revolution, the delay is derived: p, ( ) 60 ( ) 2 j j z T z n . (3.16) The dynamic chip thickness is obtained when vibrations are projected on the radial direction of the cutting edge [2]: Chapter 3. Dynamics of Peripheral Milling 25 d, ( , ) x( , )sin ( )sin ( , ) y( , )sin ( )cos ( , ) z( , )cos ( )j j jh t z t z z t z t z z t z t z z . (3.17) Substituting Eqs. (3.11), (3.17) to Eq. (3.10), the differential dynamic cutting force for flute j is expressed: cd c rta, rta( , ) {sin ( , ) cos ( , ) cot ( )} ( , )j j jd t z t z t z z t z dz F qK . (3.18) Wu [31] states that the contact pressures are exerted on the flank face as a result of tool indentation into workpiece, when finite edge radius exists. Since cutting edge of the tool is not perfectly sharp, the edge or plowing forces are always available. Moreo- ver, vibration marks on the surface affect the amount of material pressed under the tool. Figure 3-4 shows static displacement owing to edge radius and dynamic displacement due to wavy surface. Figure 3-4 : Indentation of flank due to edge radius and undulated surface Chapter 3. Dynamics of Peripheral Milling 26 The static plowing forces can be calculated through extensive FE models deter- mining the chip separation and elastically deformed volume under the tool flank [42]. Since it is out of this thesis’ scope, the static edge forces are approximated by linear edge force theory for general end mills [4][36]: es e rta, rta( , ) ( )jd t z dS z F K , (3.19) where e T rta re te ae{ }K K KK are the edge force coefficients [40]. The dynamic edge forces are exerted against the undulated surface, and they are assumed to be a function of the elastically deformed volume. The instantaneous differen- tial dynamic edge force in radial direction is expressed as [31]: ed d r, spF ( , ) ( , )j jd t z K dV t z , (3.20) and the differential dynamic edge force in tangential direction is approximated by dry friction model given by Coulomb: ed ed t, r,F ( , ) F ( , )j jd t z d t z , (3.21) where spK is the material dependent contact force coefficient. The contact force coeffi- cient can be calculated through tool-material indentation tests as reported in [32][43], or by using fast tool servo mechanism that eliminates regeneration forces but stimulates process damping forces as discussed in [44]. The dynamic Coulomb friction coefficient is , and d ( , )jdV t z is instantaneous differential dynamic volume indented into the surface, which is calculated from: d d( , ) ( , )j jdV t z dA t z dS . (3.22) The instantaneous indented area in radial and tangential plane is denoted by d ( , )jdA t z . Slope of the undulated surface, clearance angle , geometry of the cutting edge such as edge chamfer, edge radius and tool wear length are the reported factors affecting the area [31][32][34][45][46][47]. In this thesis, the approach of Chiou and Liang [32], who considered the tool wear effect, is used in evaluating the indented area: Chapter 3. Dynamics of Peripheral Milling 27 2 d w c ( , ) ( , ) ( , ) 2 ( , ) j j L dA t z p t z r t z v t z , (3.23) where wL is the average flank wear length of the cutter, c ( , )v t z is the cutting velocity of the edge in cut, ( , )jr t z is the vibration velocity of the surface in radial direction, and ( , )p t z is the unit step function accounting the instants that indentation due to wavy surface took place. Since the dynamic edge force incorporates the vibration velocity and tries to damp out the surface vibrations, it is called process damping force. The indirect proportionality to cutting speed, which comes from the slope of the surface, represents the increasing stability behavior at low cutting speeds. The cutting velocity is expressed as: c ( ) ( , ) ( ) 2 D z v t z t . (3.24) The surface vibration velocity is obtained when velocities in x,y,z are projected to the radial direction: ( , ) x( , )sin ( )sin ( , ) y( , )sin ( )cos ( , ) z( , )cos ( )j j jr t z t z z t z t z z t z t z z . (3.25) Several methods have been proposed for the unit step function calculation based on the instant of flank contact with workpiece [32][48]. Eynian and Altintas [49] sug- gested that since static indentation, in volume wise s ( , )jdV t z , is happening at all times, the dynamic edge forces are applied against the surface vibrations regardless of uphill or downhill motion as shown in Figure 3-5. In uphill motion, the indentation amount decreases and pushes tool less to bring it close to equilibrium; whereas in downhill motion, the indentation increases and more contact force is applied to push the tool to equilibrium position. The unit step function is set: ( , ) 1p t z . (3.26) Chapter 3. Dynamics of Peripheral Milling 28 Figure 3-5: Indentation of flank at several positions of the wavy surface The process damping force is derived when Eqs. (3.22)- (3.26) are substituted into Eqs. (3.20) and (3.21) consecutively: 2 sp wed rta, 1 ( , ) {sin ( , ) cos ( , ) cot ( )} ( , ) ( ) ( ) 0 j j j K L d t z t z t z z t z dz t D z F q (3.27) After defining each component of the Eq. (3.8), the forces in radial, tangential and axial directions are projected to Cartesian coordinates by the transformation matrix T [39]: xyz rta( , ) ( , ) ( , )d t z t z d t zF T F , (3.28) where the force vector is T xyz x y z( , ) { F ( , ) F ( , ) F ( , )}d t z d t z d t z d t zF . The transforma- tion matrix, transforming the rotating tool coordinates to Cartesian tool coordinates, is expressed as: sin ( , )sin ( ) cos ( , ) sin ( , )cos ( ) ( , ) cos ( , )sin ( ) sin ( , ) cos ( , )cos ( ) cos ( ) 0 sin ( ) j j j j j j t z z t z t z z t z t z z t z t z z z z T . (3.29) When Eq. (3.14), (3.18), (3.19), (3.27) are substituted into Eq. (3.28), the overall diffe- rential milling force in Cartesian tool coordinates for axial elevation z is formulated: Chapter 3. Dynamics of Peripheral Milling 29 xyz 1, 2, 1, 2, 1 ˆ ˆ ˆ ˆ( , ) ( ( , ) ( , ) ( , ) ( , ) ( , ) ( , )) N j j j j j d t z t z t z t z t z t z t z dz F A A q B B q , (3.30) The coefficient matrices for a single axial element are formulated when matrix products are obtained: T 1, 1,x 1,y 1,z 1,x rc ac tc 1,y rc ac tc 1,z ˆ ( , ) ( ( , )){ ( , ) ( , ) ( , )} ( , ) ( sin ( ) cos ( ))(cos 2 ( , ) 1) sin 2 ( , ) 2 ( , ) ( sin ( ) cos ( ))sin 2 ( , ) (1 cos 2 ( , )) 2 ( , ) j j j j j j t z g t z a t z a t z a t z c a t z K z K z t z K t z c a t z K z K z t z K t z a t z A rc acsin ( , )( cos ( ) sin ( ))jc t z K z K z , (3.31) Static Cutting Force Dynamic Cutting Force Static Edge Force Dynamic Edge Force (P. Damping) Chapter 3. Dynamics of Peripheral Milling 30 2,xx 2,xy 2,xz 2, 2,yx 2,yy 2,yz 2,zx 2,zy 2,zz 2,xx rc ac tc 2,xy ( , ) ( , ) ( , ) ˆ ( , ) ( ( , )) ( , ) ( , ) ( , ) ( , ) ( , ) ( , ) 1 ( , ) ( sin ( ) cos ( ))(cos 2 ( , ) 1) sin 2 ( , ) 2 ( , j j j j a t z a t z a t z t z g t z a t z a t z a t z a t z a t z a t z a t z K z K z t z K t z a t z A rc ac tc rc 2,xz tc ac 2,yx rc ac tc 1 ) ( sin ( ) cos ( ))sin 2 ( , ) (1 cos 2 ( , )) 2 cos ( )sin ( , ) ( , ) ( cot ( , ) cos ( )) cot ( )sin ( , ) 1 ( , ) ( sin ( ) cos ( ))sin 2 ( , ) ( 2 j j j j j j K z K z t z K t z K z t z a t z K t z K z z t z a t z K z K z t z K 2,yy rc ac tc rc 2,yz tc ac 2,zx rc cos 2 ( , ) 1) 1 ( , ) ( sin ( ) cos ( ))(1 cos 2 ( , )) sin 2 ( , ) 2 cos ( )cos ( , ) ( , ) ( tan ( , ) cos ( )) cot ( ) cos ( , ) ( , ) sin ( , )( j j j j j j j t z a t z K z K z t z K t z K z t z a t z K t z K z z t z a t z t z K ac 2,zy rc ac 2,zz ac rc cos ( ) sin ( )) ( , ) cos ( , )( cos ( ) sin ( )) ( , ) cos ( )( cot ( )) j z K z a t z t z K z K z a t z z K K z , (3.32) Chapter 3. Dynamics of Peripheral Milling 31 T 1, 1,x 1,y 1,z 1,x re ae te 1,y re ae te 1,z re ae ˆ ( , ) ( ( , )){ ( , ) ( , ) ( , )} ( , ) ( cot ( ))sin ( , ) cos ( , ) / sin ( ) ( , ) ( cot ( ))cos ( , ) sin ( , ) / sin ( ) ( , ) cot ( ) j j j j j j t z g t z b t z b t z b t z b t z K K z t z K t z z b t z K K z t z K t z z b t z K z K B , (3.33) 2,xx 2,xy 2,xz2 sp w 2, 2,yx 2,yy 2,yz 2,zx 2,zy 2,zz 2,xx 2,xy ( , ) ( , ) ( , ) ˆ ( , ) ( ( , )) ( , ) ( , ) ( , ) ( ) ( ) ( , ) ( , ) ( , ) 1 ( , ) sin ( )(cos 2 ( , ) 1) sin 2 ( , ) 2 1 ( , ) 2 j j j j b t z b t z b t z K L t z g t z b t z b t z b t z t D z b t z b t z b t z b t z z t z t z b t z B 2,xz 2,yx 2,yy 2, sin ( )sin 2 ( , ) (1 cos 2 ( , )) ( , ) cos ( )sin ( , ) cot ( ) cos ( , ) 1 ( , ) sin ( )sin 2 ( , ) (cos 2 ( , ) 1) 2 1 ( , ) sin ( )(1 cos 2 ( , )) sin 2 ( , ) 2 j j j j j j j j z t z t z b t z z t z z t z b t z z t z t z b t z z t z t z b yz ( , ) cos ( )cos ( , ) cot ( )sin ( , )j jt z z t z z t z (3.34) 2,zx 2,zy 2,zz ( , ) cos ( )sin ( , ) ( , ) cos ( )cos ( , ) 1 ( , ) cos 2 ( ) / sin ( ) 2 j j b t z z t z b t z z t z b t z z z . Chapter 3. Dynamics of Peripheral Milling 32 3.2.2 Equations of Motion in Laplace Domain In order to associate forces with structural dynamics, inertial, elastic and damp- ing characteristics of the workpiece and the tool must be identified separately. The structure’s dynamic properties correspond to amplitudes and phases of the calculated frequency response function which is comprised of measured force (input) and motion (output). Accurate representation of the dynamics can be obtained by scanning sufficient range of excitation frequency with modal hammer impact or shaker. Figure 3-6 : Discretized representation of the plate and the end mill In theory, the thin walled element is most flexible in its thickness direction ( y ) and in most of the cases; equations of motion in the other two orthogonal ( x,z ) direc- tions are neglected. The cutter is usually considered rigid compared to thin plate. In this Chapter 3. Dynamics of Peripheral Milling 33 thesis, for generality, three orthogonal coordinates which are assumed to lie at principal vibration directions are considered, and later simplifications are made for specific applications. In order to model the continuous dynamic nature of the structures, similar to milling force approach, spatial discretization is applied. The thin walled plate is divided into p w modal nodes and cutter into 1l where the differential axial depth, dz kept same for both structures as presented in Figure 3-6. The nodes that are in contact have the same number in workpiece and tool side. Assuming that the transfer functions are identified for each node, the steady state representation of equations of motion for arbitrary structure is expressed in Laplace domain: xx xy xz x xyz yx yy yz y zx zy zz z ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) s s s s s s s s s s s s s ss s s s G G G F x G F Q G G G F y zG G G F , (3.35) where receptance (displacement) transfer function matrix is denoted by ( )sG , milling force vector is xyz ( )sF and physical displacement vector is ( )sQ . The sub-matrices in transfer function matrix are defined as ( )ji sG , where dis- placement (response) direction subscript is x,y,zj and force direction (excitation) subscript is x,y,zi . When j i , the term is called direct transfer function for the specific direction; when j i , the term is called cross transfer function, accounting the geometric coupling effects coming from the other orthogonal directions. The transfer functions in arbitrary directions are identified through measurement for both tool and workpiece, respectively: Chapter 3. Dynamics of Peripheral Milling 34 1 1,t 1 2,t 1 ,t 1 ,t 2 1,t 2 2,t 2 ,t 2 ,t ,t 1,t ,t ,t 1,t ,t ,t ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) j i j i j iq j il j i j i j iq j il ji jqi jqiq jqil jli jliq jlil l l G s G s G s G s G s G s G s G s s G s G s G s G s G s G s G , (3.36) and ,wp 1 1,wp 1 2,wp 1 ,wp 1 ( ),wp 2 1,wp 2 2,wp 2 ,wp 2 ( ),wp 1,wp ,wp ( ),wp ( ) 1,wp ( ) ,wp ( ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ji j i j i j iq j i p w j i j i j iq j i p w jqi jqiq jqi p w j p w i j p w iq j p s G s G s G s G s G s G s G s G s G s G s G s G s G s G G ) ( ),wp ( ) ( )( )w i p w p w p ws (3.37) where looking at the form ( )jaibG s , the subscript a indicates the measurement node, and b represents the node where force is applied. For linear elastic bodies, Maxwell’s reciprocity theory holds [50]: ( ) ( )jaib ibjaG s G s , (3.38) therefore measuring just a column or a row of Eqs. (3.36) and (3.37) is enough to construct whole symmetrical matrix. A sub-vector in force vector, iF , and sub-vector in displacement vector, i , for x,y,zi are expressed, respectively: T ,1 ,2 ,( ) {F ( ) F ( ) F ( ) }i i i i qs s s sF , (3.39) T 1 2( ) { ( ) ( ) ( ) }qs i s i s i s i , (3.40) where the sizes are 1l and ( ) 1p w for the tool and workpiece, respectively. Chapter 3. Dynamics of Peripheral Milling 35 3.3 Modal Space Formulation In previous section the equations of motion are presented in Laplace domain. Due to discretized representation, the number of equations to solve can be plenty de- pending on the axial engagement length. It is known that, analyses in dynamics consume computational effort since the equations need to be evaluated repeatedly in specified time interval [51]. In order to reduce the amount of computation, the reduced degrees of freedom (DOF) are used instead of full DOF. Although, this approach brings errors by disregarding information coming from eliminated DOF, or in other words natural frequencies (modes), in practice most applications have dominant dynamics at certain frequency range. Consequently, selecting sufficient number of DOF for certain frequen- cy bandwidth decreases the computation time dramatically without affecting the solution accuracy. The frequency response function obtained from the measurement with real mode shape assumption has the form [50]: , , 2 2 2 1 , , 1 1 ( ) 2 m ja ib jaib R R jaib n n jaib u u G s s M s s K , (3.41) where , R R jaib jaibM K are residual mass and stiffness terms respectively, which are included to correct truncation errors coming from the elimination of low and high frequency modes. , ,,ja ibu u are the mass normalized mode shapes for nodes in assigned directions ja and ib ; , ,n are the natural frequency and the modal damping ratio for specific mode λ , respectively. The displacement is formulated when each term in Eq. (3.41) superposed: * rbody st 1 ( ) ( ) ( ) ( ) m s s s s Q Q Q Q , (3.42) Chapter 3. Dynamics of Peripheral Milling 36 where rbody ( )sQ implies the rigid body motion behavior at low frequencies [52], *( )sQ is m number of identified modes’ contribution to displacement, and st ( )sQ represents the displacement correction for high frequency modes, which can be associated with static correction in finite element analysis [51]. It is known that, for the machining processes, the instability condition and the detrimental forced vibrations happen around the dominant natural frequencies of the structure. Selecting sufficient number of domi- nant modes and neglecting the residual terms, if they have small amplitude in considered frequency range, can be applied without affecting the accuracy. The total displacement in Eq. (3.42) is approximated as the superposition of the modal displacements in physi- cal coordinates as: * 1 ( ) ( ) m s s Q Q . (3.43) In practice number of selected modes is smaller than the number of DOF such as t 3m l and wp 3( )m p w , where tm is the number of tool’s natural frequencies, and wpm is number of workpiece’s natural frequencies. Benefiting from the linear modal superposition shown in Eq. (3.43), a modal space can be defined to collect all the nodal responses to a modal coordinate, for specific mode λ . This approach reduces the number of equations from 3l and 3( )p w to tm and wpm for tool and workpiece, respectively. Applying real eigenvalue analysis, Eq. (3.43) is arranged by defining a modal space: ( ) ( )s sQ U , (3.44) where t T t t,1 t,2 t,( ) { ( ) ( ) ( )}ms s s s is modal displacement for the tool, and wp T wp wp,1 wp,2 wp,( ) { ( ) ( ) ( )}ms s s s is modal displacement for the work- piece. The mass normalized modal (mode shape) matrices, which are related with the eigenvectors of the eigenvalue problem, are expressed for tool and workpiece, respec- tively: Chapter 3. Dynamics of Peripheral Milling 37 t t t t t t t t x1,1,t x1,2,t x1, ,t x ,1,t x ,2,t x , ,t x ,1,t x ,2,t x , ,t y1,1,t y1,2,t y1, ,t y ,1,t y ,2,t y , ,tt y ,1,t y ,2,t y , ,t z1,1,t z1,2,t m q q q m q m l l l m m q q q m q m l l l m u u u u u u u u u u u u u u u u u u u u U t t t t t z1, ,t z ,1,t z ,2,t z , ,t z ,1,t z ,2,t z , ,t 3 m q q q m q m l l l m l m u u u u u u u , (3.45) and tÛ Chapter 3. Dynamics of Peripheral Milling 38 wp wp wp wp wp x1,1,wp x1,2,wp x1, ,wp x ,1,wp x ,2,wp x , ,wp x( ),1,wp x( ),2,wp x( ), ,wp y1,1,wp y1,2,wp y1, ,wp y wp m q q q m q m p w p w p w m m u u u u u u u u u u u u u U wp wp wp wp wp ,1,wp y ,2,wp y , ,wp y( ),1,wp y( ),2,wp y( ), ,wp z1,1,wp z1,2,wp z1, ,wp z ,1,wp z ,2,wp z , ,wp q q q m q m p w p w p w m m q q q m u u u u u u u u u u u wp wp wp z( ),1,wp z( ),2,wp z( ), ,wp 3( ) q m p w p w p w m p w m u u u , (3.46) where each column of the matrices gives mass normalized displacements in x,y,z directions for specific mode. The q number of nodes is highlighted to indicate the nodes that are in cut. The receptance transfer function matrix is expressed in terms of mass normalized mode shapes and modal parameters, with symmetric mass and stiffness matrices ap- proach [50]: 2 2 1 T( ) ( 2 )n ns s s G I ζω ωU U , (3.47) where I is the identity matrix, ζ is the diagonal damping ratio matrix, and nω is the diagonal natural frequency matrix. When Eqs. (3.44) and (3.47) are substituted in Eq. (3.35), the following form appears: 2 2 1 T xyz( 2 ) ( ) ( )n ns s s s I ζω ω FU U U . (3.48) wpÛ Chapter 3. Dynamics of Peripheral Milling 39 Eq. (3.48) is pre-multiplied by 1 U (if mass matrix is introduced, benefiting from the orthogonality properties, T U M can be also used for pre-multiplication [52]): 2 2 1 T xyz( 2 ) ( ) ( )n ns s s s I ζω ω FU , (3.49) which presents the modal domain model with tm number of equations to solve for tool, and wpm number of equations for workpiece. The time domain model is obtained when Inverse Laplace Transform is applied to Eq. (3.49): 2 T xyz( ) 2 ( ) ( ) ( )n nt t t t ζω ω F U . (3.50) In general, full size modal matrices are used for modal space representation. However, in this particular problem, only the nodes (1, , )q are in cut for both workpiece and the cutter as demonstrated in Figure 3-6. Although, other nodes vibrate due to off diagonal terms in transfer functions, they do not contribute to cutting forces via regeneration and process damping. The effective terms in force vector, in physical domain are identified using Eqs. (3.1), (3.28) and (3.39) for x,y,zi : T i ˆ ( ) { F ( , ) F ( , )}i it d t dz d t qdzF . (3.51) When the equations are expressed in modal domain, truncating the mode shape vectors after node number q gives the same result with the full modal matrix approach for the force side. The truncated modal matrices, which are shown in Eqs. (3.45) and (3.46), are denoted by t wp ˆ ˆ,U U for tool and workpiece, where the matrix sizes are t3q m and wp3q m accordingly. The spatial terms in physical coordinates in the force vector are transformed to modal space by ˆ( ) ( )t t Q U and ˆ( ) ( )t tQ U . However, special care is required Chapter 3. Dynamics of Peripheral Milling 40 for non-uniform helix cutters where delay changes with axial position. The modal coordinate representation of non-uniform helix cutter’s delayed term is: T T T 1 2x1,1 x1,2 x1, T x2,2 x2,2 x2, 1 2 x ,1 x ,2 x , 1 2 ˆ( ) ( ) [ ( ( )) ( ( )) ( ( ))][ ] [ ] [ ( (2 )) ( (2 )) ( (2 ))] [ ] [ ( ( )) ( ( )) ( ( j j m jm m j j m j q q q m j j m j t t t T dz t T dz t T dzu u u u u u t T dz t T dz t T dz u u u t T qdz t T qdz t T q Q U T T y1,1 y1,2 y1, 1 2 T y ,1 y ,2 y , 1 2 T z1,1 z1,2 z1, 1 2 ))] [ ] [ ( ( )) ( ( )) ( ( ))] [ ] [ ( ( )) ( ( )) ( ( ))] [ ] [ ( ( )) ( ( )) ( ( ))] m j j m j q q q m j j m j m j j m j dz u u u t T dz t T dz t T dz u u u t T qdz t T qdz t T qdz u u u t T dz t T dz t T dz T y ,1 y ,2 y ,1 1 2 3 1 [ ] [ ( ( )) ( ( )) ( ( ))]q q q j j m j q u u u t T qdz t T qdz t T qdz , (3.52) where {1, , }r q is the axial index of the cutting point. ,j rT is the delay value for the corresponding axial elevation and tooth number j , as calculated in Eq. (3.16). The sign indicates row by column multiplication as shown in Eq. (3.52). The detailed discussion on mapping of each term in semi discretization method is discussed in chapter 4. The regenerative displacement and the velocity in forcing function are calculated regarding the relative motion between tool and workpiece as discussed by Bravo et al. [19]. The relative displacement for axial position r is shown in Figure 3-7. In order to define relative regenerative displacement and relative vibration velocity in modal coordinates, the following variables are introduced: t t wp wp t t t t , wp wp wp wp , ˆ ˆ( ) ( ) uniform helix ( ) ˆ ˆ ˆ ˆ( ) ( ) ( ) ( ) non-uniform helixj r j r t t t t t T t t T Ψ U U U U U U , (3.53) and t t wp wp ˆ ˆ( ) ( ) ( )t t t Ψ U U . (3.54) Chapter 3. Dynamics of Peripheral Milling 41 Figure 3-7: Relative displacement for axial point r Eq. (3.30) is expanded for all nodes and substituted into Eq. (3.50). Consequent- ly, the time domain modal equations of motion for tool and workpiece are formulated: 2 t t ,t t ,t t T t 1, 2, 1, 2, 1 ( ) 2 ( ) ( ) ˆ ( ( ) ( ) ( ) ( ) ( ) ( )) n n N j j j j j t t t t t t t t t dz ζ ω ω A A Ψ B B Ψ U , (3.55) and 2 wp wp ,wp wp ,wp wp T wp 1, 2, 1, 2, 1 ( ) 2 ( ) ( ) ˆ ( ( ) ( ) ( ) ( ) ( ) ( )) n n N j j j j j t t t t t t t t t dz ζ ω ω A A Ψ B B Ψ U , (3.56) Chapter 3. Dynamics of Peripheral Milling 42 where 1, 2, 1, 2,( ), ( ), ( ), ( )j j j jt t t tA A B B are expanded and arranged forms of corresponding 1, 2, 1, 2, ˆ ˆ ˆ ˆ( , ), ( , ), ( , ), ( , )j j j jt z t z t z t zA A B B terms in Eq. (3.30). They are expressed as: 1, T 1,x 1,x 1,y 1,y 1,z 1,z 3 1 ( , ) ( ( )) { ( , ) ( , ) ( , ) ( , ) ( , ) ( , )} j j q t z t a t dz a t qdz a t dz a t qdz a t dz a t qdz A g , (3.57) 2, 2,xx 2,xy 2,xz 2,xx 2,xy 2,xz 2,yx 2,yy 2,yz 2,yx 2,yy 2,yz 2, ( ) ( ( )) ( , ) 0 ( , ) 0 ( , ) 0 0 0 0 0 ( , ) 0 ( , ) 0 ( , ) ( , ) 0 ( , ) 0 ( , ) 0 0 0 0 0 ( , ) 0 ( , ) 0 ( , ) j jt t a t dz a t dz a t dz a t qdz a t qdz a t qdz a t dz a t dz a t dz a t qdz a t qdz a t qdz a A g zx 2,zy 2,zz 2,zx 2,zy 2,zz 3 3 ( , ) 0 ( , ) 0 ( , ) 0 0 0 0 0 ( , ) 0 ( , ) 0 ( , ) q q t dz a t dz a t dz a t qdz a t qdz a t qdz , (3.58) 3 1 1, T 1,x 1,x 1,y 1,y 1,z 1,z ( ) ( ( )) { ( , ) ( , ) ( , ) ( , ) ( , ) ( , )} q j jt t b t dz b t qdz b t dz b t qdz b t dz b t qdz B g , (3.59) Chapter 3. Dynamics of Peripheral Milling 43 2, 2,xx 2,xy 2,xz 2,xx 2,xy 2,xz 2,yx 2,yy 2,yz 2,yx 2,yy 2,yz 2, ( ) ( ( )) ( , ) 0 ( , ) 0 ( , ) 0 0 0 0 0 ( , ) 0 ( , ) 0 ( , ) ( , ) 0 ( , ) 0 ( , ) 0 0 0 0 0 ( , ) 0 ( , ) 0 ( , ) j jt t b t dz b t dz b t dz b t qdz b t qdz b t qdz b t dz b t dz b t dz b t qdz b t qdz b t qdz b B g zx 2,zy 2,zz 2,zx 2,zy 2,zz 3 3 ( , ) 0 ( , ) 0 ( , ) 0 0 0 0 0 ( , ) 0 ( , ) 0 ( , ) q q t dz b t dz b t dz b t qdz b t qdz b t qdz . (3.60) The diagonal unit step function matrix is defined as: 3 3 ( ( )) ( ( , )) 0 0 0 0 0 0 ( ( , )) ( ( , )) 0 0 0 0 0 0 ( ( , )) ( ( , )) 0 0 0 0 0 0 ( ( , )) j j j j j j j q q t g t dz g t qdz g t dz g t qdz g t dz g t qdz g 0 0 0 0 0 0 . (3.61) The workpiece and tool equations are coupled to each other due to relative rege- nerative displacement and velocity terms in the force function. The Eqs. (3.55) and (3.56) are combined by expanding the matrices: 2 t ,t ,t tt t 2 wp ,wp wpwp wp ,wp T t te 1, 2, 1, 2, t wpT wpwp 2 ( )( ) ( ) 2 ( )( ) ( ) ˆ ( ) ˆ ˆ( ) ( ) ( ) ( ) ( ) ˆ ( ) n n n n j j j j tt t tt t t t t t t t t ζ ω 0 ω 0 0 ζ ω 0 ω A A Ψ B B U U U U 1 N j dz , (3.62) Chapter 3. Dynamics of Peripheral Milling 44 where expanded relative regeneration term is defined similar to Eq. (3.53): t t wp wp e t ,t t wp t wp wp wp , ( ) ˆ ˆ uniform helix ( ) ( ) . ( )( ) ˆ ˆ ˆ ˆ non-uniform helix ( ) ( ) j r j r t t t t Tt t t T Ψ U U U U U U (3.63) Eq. (3.62) gives the equations of motion of peripheral milling in modal domain. It is simplified as: 2 T t cs cd cd es ed 1 2 TT wp ( ) 2 ( ) ( ) ˆ ( ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )) ˆ n nt t t t t t t t t t t dz ζω ω U f f f f f U , (3.64) where the expanded displacement matrix is T t wp( ) { ( ) ( )}t t t . The static cutting force component is: cs 1, 1 ( ) ( ) N j j t t Af . (3.65) The dynamic cutting force component due to the current vibrations is: cd1 2, t wp 1 ˆ ˆ( ) ( ) N j j t t Af U U , (3.66) The dynamic cutting force component due to the previous vibrations is defined for three cases: cd 2 cd cd 2 T 2, 1 cd 2, , , 1 1 ˆ ( ) ( ) uniform helix and pitch ˆ( ) ( ) ( ) ( ) uniform helix, non-uniform pitch ˆ ( ) ( ) non-uniform helix and pitch N j j j qN j r j r j r t t T t t t t T t t T f f f f , (3.67) Chapter 3. Dynamics of Peripheral Milling 45 where cd 2, , 2, t wp ˆ ˆ ˆ( ) ( )j r j r t t Af U U . (3.68) The expanded modal matrix for axial element r is: t wp t wp t wp t wp xr,1,t xr, ,t x ,1,wp x , ,wp yr,1,t yr, ,t y ,1,wp y , ,wpt wp zr,1,t zr, ,t z ,1,wp z , ,wp 3 ( ) ˆ ˆ m r r m m r r m r m r r m q m m u u u u u u u u u u u u 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 U U . (3.69) The other dynamic cutting force components are calculated as cd cd2, 2, , 1 ˆ ˆ( ) ( ) q j j r r t t f f and cd cd 2 2, 1 ˆ ˆ( ) ( ) N j j t t f f . The static edge force component is: es 1, 1 ( ) ( ) N j j t t Bf . (3.70) The dynamic edge force or process damping force component is: ed 2, t wp 1 ˆ ˆ( ) ( ) N j j t t Bf U U . (3.71) Chapter 4. Time Domain Simulation 46 4. Time Domain Simulation 4.1 Introduction The methods to solve equations of motion for the peripheral milling process are dis- cussed in this chapter. Semi discretization method proposed by Insperger and Stepan [24] for delayed differential equations with periodic coefficients is extended by includ- ing the static and process damping forces for the time domain simulation. Moreover, for stable cutting conditions, surface location error method proposed by Schmitz [10] is used accounting the steady state response of the differential equation. 4.2 Semi Discrete Time Domain Solution The equations of motion derived in previous chapter are delayed differential equ- ations (DDE) with periodic coefficients, where the solution cannot be given in closed form [53]. In order to solve this type of equations, researchers have considered true kinematics model in time domain where they included inner and outer modulations by keeping the past and current surface data and considering the nonlinearities [7] [54]. This model is then solved employing explicit Runge-Kutta methods. However, storing all the geometric data is time consuming. Furthermore, stability of the time stepping scheme and order of accuracy enforce programmers to take small time steps, which adds on to the computation time. Insperger and Stepan [24] proposed semi discretization (SD) method, where they solved the stability of the linear milling equation semi analytically by just taking regenerative force term into account. The objective of SD technique is to approximate delayed and time dependent spatial terms at each discrete time step by piecewise constant values, while keeping actual time domain terms in the original form [53]. In this section, DDE representing the peripheral milling dynamics is solved using time domain SD method by considering all reported force components, while disregard- ing nonlinearities like runout, tool jump and change of radial immersion due to deflections. Chapter 4. Time Domain Simulation 47 The modal equation of motion in time domain was derived in section 3.3: 2 T t cs cd cd es ed 1 2 TT wp ( ) 2 ( ) ( ) ˆ ( ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )) ˆ n nt t t t t t t t t t t dz ζω ω U f f f f f U (4.1) The equation is transformed to state space form by dropping the order of DDE one and doubling the number of equations: T( ) ( ) ( ) ( ) ( ) ( )t t t t t t Θ L Θ R Θ S , (4.2) where the motion vector is: t wp2( ) 1 ( ) ( ) ( ) m m t t t Θ , (4.3) the system matrix is: t wp t wp T T t t2 cd ed 1T T wp wp 2( ) 2( ) ˆ ˆ( ) ( ) 2 ( ) ˆ ˆn n m m m m t t dz t dz 0 I L ω ζω U U f f U U , (4.4) the static force matrix: t wp T t cs es T wp 2( ) 1 ˆ( ) ( ( ) ( )) ˆ m m t t t dz 0 S U f f U . (4.5) Chapter 4. Time Domain Simulation 48 The product of retarded matrix and delayed motion vector has 3 different forms: T 1 , , 1 1 ˆ ( ) ( ) uniform helix and pitch ˆ( ) ( ) ( ) ( ) uniform helix, non-uniform pitch ˆ ( ) ( ) non-uniform helix and pitch N j j j qN j r j r j r t t T t t t t T t t T R Θ R Θ R Θ R Θ (4.6) T t cd 2T wp T t cd 2,1 T wp ,T t cd 2,1 1 T wp ( )ˆ ˆ ( ) ( ) ˆ ( ) ˆ ˆ ( ) ( ) ˆ ( ˆ ˆ ( ) ˆ N j j jj qN j jj r t T t dz t T t T t dz t T t T t dz 0 0 0 0 0 0 0 0 0 U f U U f U U f U , ) ( ) r j rt T , where cd cd cd 2 2, 2, , ˆ ˆ ˆ, ,j j rf f f are discussed in section 3.3. The most general case , , 1 1 ˆ ( ) ( ) q N j r j r r j t t T R Θ is considered for the following equations. Chapter 4. Time Domain Simulation 49 Figure 4-1: Approximation of the delayed term by discretization SD method starts the formulation by discretizing the delay term into finite number of points as illustrated in Figure 4-1. The time interval is expressed as: , 0,1,2, ,it i t i k . (4.7) The time step t is a significant factor that affects the accuracy of the numerical solu- tion. Two considerations need to be present when transforming the system to discrete domain. The time step should be small enough to capture all vibrations without having the aliasing problem: ,maxn t , (4.8) where ,maxn is the highest natural frequency of the structure. Moreover, the time varying modal forces or modal force components should be sufficiently represented in discrete domain. Figure 4-2 compares the effect of different number of discretization points on the time varying modal force component. Sudden changes, which mainly occur Chapter 4. Time Domain Simulation 50 in intermittent cutting conditions, cannot be captured by the method, if there is no point to represent it. This brings sizable numerical error to the solution as discussed in [55]. Figure 4-2: Discretization of an element in dynamic modal force component matrix with different number of points The lower limit of the time step is confined by simulation accuracy and the upper limit is confined by computational limitations. The number of arithmetic operations is roughly proportional to 4 1 t ; therefore compromise needs to be made between order of the error and the simulation time. The integer number of discretized points is calculated as discussed in [29]: maxint( / / 2)k T t p , (4.9) where maxT is the maximum delay in the system, p is the polynomial order of the approximation. It has been reported that the order of approximation of delayed term affects the convergence to exact solution [56][57]. Moreover, when non-uniform helix and pitch cutters are concerned, the delay term can be in between two discretized points. 60 65 70 75 80 85 90 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 10 6 Cutter Rotation Angle [deg] D ir e c ti o n a l M o d a l F o rc e C o m p o n e n t] k=24 k=60 k=119 Chapter 4. Time Domain Simulation 51 In order to increase solution’s accuracy first order interpolation of the delay is reported to be convenient [56], which has the form: , , , , , , , 1 ( ) ( 1 ) ( ) ( ) ( ) j r j r j r j r j r j r j r i k i k t T i k t t T i k t t T t t t t Θ Θ Θ , (4.10) where , ,int( / / 2)j r j rk T t p is mapping the delayed term coming from particular tooth and axial elevation for ,2 j rk k . SD method approximates time dependent matrices in each interval as constant values: 1 1 1 , , , , 1 1 1ˆ ˆ ˆ( ) ( ) , ( ) ( ) , ( ) ( ) i i i i i i t t t i i j r i j r j r i i i t t t t t dt t t dt t t dt t t t L L L R R R S S S . (4.11) The differential equation in Eq. (4.2) can be treated as constant ordinary differential equation within the interval, 1i it t t , when the general solution is attained: . .( ) , , , 1 10 ˆ( ) (0) ( )i i t qN t t j r i j r i j r t e e T d L L Θ Θ R Θ S . (4.12) Eq. (4.10) is substituted into Eq. (4.12) and it is assumed that within the interval, 1i it t t , 1 i L , exists and , , ˆ j r iR stays constant. The following equation is obtained when 1( )it Θ is written in terms of ( )itΘ , similar to Insperger and Stepan [24]: , ,1 ,0, , 1 ,1, , 1 1 1 1 ( ) ( ) ( ) ( ) j r j r q qN N i i i i j r i k i j r i k i j r j r t t t t Θ CΘ D Θ D Θ E (4.13) where .i t i e LC , (4.14) .1 2 1,0, , , , , , 1 ˆ( ( 1) ) ( )i ti j r i i j r j r i j r iT k t e t L D L L L I R , (4.15) Chapter 4. Time Domain Simulation 52 .1 2 1,1, , , , , , 1 ˆ( ) ( )i ti j r i i j r j r i j r iT k t e t L D L L L I R , (4.16) . 1( )i t i i ie LE I L S . (4.17) For non-uniform pitch and uniform helix, and uniform pitch and helix cutters, , ,i e jD and ,i eD for 0,1e are introduced, respectively. The terms are defined as , , , , , 1 q i e j i e j r r D D and , , , 1 N i e i e j j D D . The Eq. (4.14) is mapped using the form 1i i i i Θ ZΘ W : , , , 1 ,0, , ,1, , ( 1) ( 1) ( )( ) .( ) .. .. .. ( ). ( ) ( ) j r j r j r ii i i j r i j r i i k i k i k tt t t t t ΘΘ C 0 0 D D E Θ I 0 0 0 0 0 I 0 0 0 Θ0 0 0 I 0 0 Θ Θ0 0 0 0 I 0 . . . i 0 0 0 , (4.18) where given the process mechanics and information on past vibrations, the solution of next time step can be calculated. For uniform pitch and helix cutters, there is only one delay. The matrices asso- ciated with the delayed term are located at the first row, k th and ( 1)k th columns of the transition matrix. When multiple delays exist due to non-uniform helix and pitch angles, the delayed matrices are accumulated at first row and corresponding ,j rk th and ,( 1)j rk th columns. SD solution requires t wp2 ( )k m m number of initial conditions to start the time domain simulation. Similar to models considering the true kinematics, at the beginning, Chapter 4. Time Domain Simulation 53 all initial conditions are set to zero. For the first iterations, the only excitation component affecting the solution is forced vibration term iE . Then algorithm starts to assign pre- vious solutions obtained to corresponding delay elements of the motion vector iΘ , thus regenerative vibration term becomes effective. After solving for the modal displace- ments and velocities for determined time interval, the physical displacements for nodes that are in cut are calculated by using transformation from modal space to physical space as: t wp ˆ ( ) ˆ ˆ ( ) ˆ ( ) t t t Q Θ Q U U . (4.19) The size of the coefficient matrix iZ is t wp t wp2( 1)( ) 2( 1)( )k m m k m m .In order to calculate the next time step, multiplication of first row and the motion vector is sufficient. Moreover, the velocity terms in retarded matrix which have no effect on the solution, can be eliminated and the transition matrix size can be reduced to t wp t wp( 2)( ) ( 2)( )k m m k m m as discussed in [58]. Some other efficient calcula- tion methods for SD algorithm had been proposed by Henninger and Eberhard [55]. After finding the velocities and displacements for each node, differential dynam- ic force can be calculated by using Eq. (3.18), and differential process damping force can be obtained from Eq. (3.27). 4.2.1 SD Time Domain Solution Example The SD time domain solution is coded for different types of end mills and validi- ty of the method has been analyzed by comparing it to commercial software (CUTPRO). Standard cutter’s dynamics is tested for stable and unstable cutting cases. The multi point mechanics described in the previous sections is lumped into a single point to exactly mimic the commercial software. The methodology of single point approach is explained in the Appendix A. Chapter 4. Time Domain Simulation 54 Among the three different cutter types tested, standard uniform pitch and uniform helix cutter example is chosen from Turner et al. [59] due to good experimental match with the simulations. In Table 4-1, the dimensions of the cutter are specified. Table 4-1: Parameters for regular cylindrical end mill Tool Diameter [mm] Teeth Helix 1 [deg] Helix 2 [deg] Pitch 1 [deg] Pitch 2 [deg] ST3 16 4 30 30 90 90 The modal parameters given in Table 4-2 are experimentally identified for the tool tip in x and y directions, while assuming that the workpiece is considerably rigid compared to cutter. Three modes in each direction are extracted. The axisymmetric geometry of the cutter brings similar frequency modes with mode shape vectors in different directions. It is assumed that the modes in each direction are orthogonal to each other, hence no cross talking exists between the directions. When each x, y and z direc- tion mass normalized mode shape value for the tool tip is placed side by side, lumped mode shape matrix (modal matrix) is constructed: t tip 3 6 3.543 1.512 1.303 0 0 0 0 0 0 3.653 1.616 1.483 0 0 0 0 0 0 U . (4.20) Table 4-2: Modal parameters of the regular cylindrical end mill Tool and Direction Natural Frequency [Hz] Effective Stiffness [N/m] Damping Ratio ST3 (x) 2062, 2400, 2956 1.337e7, 9.944e7, 2.031e8 0.0206, 0.0183, 0.0151 ST3 (y) 2063, 2388, 2957 1.259e7, 8.624e7, 1.57e8 0.0199, 0.0175, 0.0150 Chapter 4. Time Domain Simulation 55 The material is Al 6061 where mechanistically determined cutting force coeffi- cient vector is c T 2 rta {160 400 0} N/mmK and edge force coefficient vector is e T rta {30 26 0} N/mmK . The process damping is neglected for the selected speeds. The radial immersion is set as 5 mm. The stability chart is for low immersion down milling in Figure 4-3. In order to analyze stable and unstable behavior, 7 mm axial depth of cut (ADOC) with 7500 rpm spindle speed and 7 mm ADOC with 6800 rpm are simulated. The displacement responses in Figure 4-3c,e show that at stable cutting condition CUTPRO and SD results match well. The difference between the results are attributed to circular motion, no runout and change of radial immersion assumptions made by SD method. As vibration magnitude - radial immersion ratio becomes higher; the discrepancy becomes apparent as demonstrated in Figure 4-3b, d. When the system chatters, SD result grows quicker compared to true kinematics model. In CUTPRO simulation, vibrations directly affect the radial immersion, hence the cutting forces, so that the system tries to damp out the grow rate of the oscillations as illustrated in Figure 4-4a. It can be concluded that SD time domain simulation works accurately at less oscillatory stable cutting conditions, and determining the chatter, but it fails at cases where large vibration displacements occur. Chapter 4. Time Domain Simulation 56 a) b) c) d) e) Figure 4-3: a) Stability chart for standard cutter [59] (with the permission of D. Merdol) b)-d) time domain simulation in x and y directions respectively for 7 mm ADOC and 6800 rpm c)-e) simula- tions in x and y directions respectively for 7 mm ADOC and 7500 rpm 0 0.005 0.01 0.015 0.02 0.025 -0.2 -0.1 0 0.1 0.2 0.3 Time [s] D is p la c e m e n t [m m ] X Direction Displacement Semi Discretization CUTRPO 0 0.005 0.01 0.015 0.02 0.025 -0.01 0 0.01 0.02 0.03 Time [s] D is p la c e m e n t [m m ] X Direction Displacement 0 0.005 0.01 0.015 0.02 0.025 -0.2 -0.1 0 0.1 0.2 0.3 Time [s] D is p la c e m e n t [m m ] Y Direction Displacement 0 0.005 0.01 0.015 0.02 0.025 0 0.02 0.04 0.06 Time [s] D is p la c e m e n t [m m ] Y Direction Displacement Chapter 4. Time Domain Simulation 57 a) b) Figure 4-4: Resultant cutting forces on x-y plane a) unstable cutting b) stable cutting Frequency spectrums of both cases give insight about the modes that cause chat- ter, and also checks semi discretization method’s validity at determining unstable poles. Figure 4-5a shows the unstable x direction force’s spectrum, and due to difference between force amplitudes, it is presented in logarithmic scale. The main chatter fre- quency at 2119 Hz and its sub-super harmonics due to tooth passing frequency (453 Hz), which are 1666 Hz, 2569 Hz and 3022 Hz, are identified similarly in both methods. In stable cutting condition as illustrated in Figure 4-5b, the tooth passing frequency at 500 Hz and its harmonics are the main excitation sources. 0 0.005 0.01 0.015 0.02 0 500 1000 1500 2000 2500 Time [s] F o rc e [ N ] Resultant Force on XY Plane Semi Discretization CUTPRO 0 0.02 0.04 0 100 200 300 400 500 600 Time [s] F o rc e [ N ] Resultant Force on XY Plane Chapter 4. Time Domain Simulation 58 a) b) Figure 4-5: FFT of force data a) unstable cutting b) stable cutting In summary, it is noted that the SD accurately solves the response in stable cutting conditions. The advantage of semi discrete time domain solution is that it calculates the whole response by just doing one matrix multiplication as given in Eq. (4.19). If the number of degrees of freedom is high, the simulation time does not increase excessively as it happens in current true kinematics models. On the other hand, if the user needs the accurate behavior where nonlinearities come in, current time domain models should be used. 4.3 Surface Location Error Forced vibrations and static deflections may occur during the stable cutting process. The cutter deviates from the intended surface finish and leaves undercut or overcut regions called the surface location error (SLE) as illustrated in Figure 4-6. This problem necessitates additional finishing passes, or may even cause to scrap the whole part due to tolerance limitations. 0 1000 2000 3000 4000 10 2 10 4 10 6 Frequency [Hz] F o rc e [ N ] FFT of X Direction Force 0 500 1000 1500 2000 2500 0 2 4 6 8 x 10 5 Frequency [Hz] F o rc e [ N ] FFT of Resultant Force Semi Discretization CUTPRO453 Hz 2119 Hz 2569 Hz 3022 Hz 1666 Hz 500 Hz 1000 Hz Chapter 4. Time Domain Simulation 59 Figure 4-6: Surface location error a) undercut example b) overcut example The surface form errors had been solved in time domain with exact kinematics model as presented in [8][54]. Schmitz [10] extracted the SLE from the steady state response by solving the differential equations in frequency domain. This eliminated lengthy convolution operation used in time domain models. In this section SLE is formulated according to the differential equation derived for the peripheral milling process in chapter 3. The flowchart is demonstrated in Figure 4-7. Chapter 4. Time Domain Simulation 60 Figure 4-7: Flowchart of SLE method In stable cutting, the dynamic forces are assumed to have no influence. The peripheral milling equation becomes constant coefficient ordinary differential equation: T t2 cs es T wp ˆ ( ) 2 ( ) ( ) ( ( ) ( )) ˆn n t t t t t dz ζω ω U f f U . (4.21) Input static cutting and edge forces (which includes cutting force coefficients and cutting conditions), cs es( ) ( )t tf , f , structural parameters of tool and workpiece, T Tt wp ˆ ˆ, ,nζ ω U ,U , set spindle speed (varying) Express forces cs esf , f using Fourier series to transform them to frequency domain. Solve for Fourier coefficients either numerically or analytically. Express modal dynamics in frequency domain via transfer function ( )iM . Calculate the displacement by taking the product of force and transfer function in magnitude-phase form. Write the solution for displacement using Fourier series. Express the displacement in physical domain by taking modal space transformation and find the displacements for each point when ( , )j t z equals to 0 (up milling) or π (down milling). Store the SLE value for corresponding axial elevation. Change spindle speed, n Chapter 4. Time Domain Simulation 61 The time dependent components of the force are expressed by exponential Fourier series in order to have compact form: p T it cs es s s T wp ˆ ( ( ) ( )) ( ) e ˆ h h t h h t t dz t U f f f f U , (4.22) where h is the index of harmonics and i is the imaginary unit. The Fourier coefficients are defined as: p T -its cs es T wp0 ˆ 1 ( ( ) ( ))e ˆ h t h dz t t dt U f f f U . (4.23) The least common period, , and frequency, pω , for different cutter geometries are defined in Table 4-3. Table 4-3: Least common period and frequency for different cutters Pitch and Helix Angles Least common period (sec) Least common fre- quency ( / )p rad s Uniform pitch and helix angle 60 Nn p 2 60 n N Alternating pitch and/or alternating helix angle i.e. p,1 p,2 p,3 p,4, , , ε : number of different angles 60 ( / )N n p 2 60 n N Random pitch and/or helix angle, axially changing helix angle 60 n p 2 60 n The rotation angle is dependent on the spindle speed and time, thus change of va- riable is applied to the integration. Since 2 2 ( ) and 60 60 n n t t d dt , Eq. (4.23) is redefined as: Chapter 4. Time Domain Simulation 62 p s s T -i ts cs es T s wp0 ˆ 1 ( ( ) ( ))e ˆ h h dz d U f f f U , (4.24) where is the rotation angle, 2 / 60 s n is the spindle frequency, and s is the angle corresponding to one period of the forcing function. For complex tools, such as variable pitch, variable helix and non cylindrical cutters, obtaining coefficients via Fast Fourier Transform (FFT) will be faster compared to analytical solution as number of mathemati- cal operations is small. In this section both numerical and analytical approaches are considered. Existence of switching function ( ( ))t j g makes it necessary to redefine the Fourier series for the analytical solution. Eq. (4.23) has the form similar to Schmitz [10]: s T t is cs es , ,T 1wp ˆ ( ) ( + ) e ˆ N h h t h j h j j h t dz U f F F U , (4.25) where is the sign for element wise Hadamard product ( , , ,( ) A Bi j i j i j A B ),and vector containing the exponential terms is expressed as: 1 p p, 0 1 p p, 0 1 p p, 0 1 p p, 0 p i ( ) ψ ( ) i ( ) ψ ( ) i ( ) ψ ( ) , i ( ) ψ ( ) i e e e ( ) e e j e j e j e j e j e j e j e j e h t dz dz h t dz qdz h t dz dz h j h t dz qdz h t t e 1 p, 0 1 p p, 0 ( ) ψ ( ) i ( ) ψ ( ) 3 1 e j e j e j e j e dz dz h t dz qdz q . (4.26) Chapter 4. Time Domain Simulation 63 The Fourier coefficient vectors have the form: cs cs cs cs cs cs cs T , , ,x , ,x , ,y , ,y , ,z , ,z 3 1{ ( ) ( ) ( ) ( ) ( ) ( )} ,h j h j h j h j h j h j h j qF dz F qdz F dz F qdz F dz F qdz F (4.27) es es es es es es es T , , ,x , ,x , ,y , ,y , ,z , ,z 3 1{ ( ) ( ) ( ) ( ) ( ) ( )} .h j h j h j h j h j h j h j qF dz F qdz F dz F qdz F dz F qdz F (4.28) Each element in Eqs. (4.27) and (4.28) are calculated for one full spindle revolution and integration limits are simplified since switching function becomes non zero between start and exit angles of the cut: ex st 2 cs -i -i , , 1, 1, 0 1 1 ( ) ( , )e ( , )e 2 2 h h h j i i iF z a z d a z d , (4.29) ex st 2 es -i -i , , 1, 1, 0 1 1 ( ) ( , )e ( , )e 2 2 h h h j i i iF z b z d b z d , (4.30) where x,y,zi . The analytical solutions of the integrals are given in Appendix B. Eq. (4.22) can be simply solved for modal displacement in frequency domain by taking the product of transfer function and the modal force: s(i ) (i ) (i ) M Γf , (4.31) where the modal transfer function is: 2 2 1(i ) ( 2 i )n n M I ζω ω . (4.32) In frequency domain the multiplication is expressed in magnitude-phase form as: s(i ) (i ) (i ) (i ) Γ M Mf , (4.33) where the magnitudes are denoted by and the phases are represented by . The magnitude of the transfer function is: 2 2 2 2 1(i ) (( ) (2 ) )n n M I ω ζω , (4.34) and the phase of the transfer function is: Chapter 4. Time Domain Simulation 64 1 2 2 2 (i ) tan n n ζω M I ω . (4.35) Using Fourier series, solution in Eq. (4.34) is expressed in time domain for numerical and analytical solutions, respectively: p p, i( (< ( ))s p,( ) ( ) e h h h t i h h h t i M Γ Mf , (4.36) p,s T i(< ( ))t ics es , , p,T 1 1wp ˆ ( ) ( + ) e ( ) e ˆ h N h ih t h j h j h h j h t dz i M Γ M U F F U . (4.37) The physical displacements are obtained using the modal space transformation: t wp ˆ ˆ ˆ( ) ( )t t Q ΓU U , (4.38) where ˆ ( )tQ is the relative displacement for all points in cut. SLE is found for all points when corresponding immersion angle at determined axial elevations equal to zero for up milling, or for down milling. For cases where only the tool tip deflection is important, using spatial coordi- nates for the SLE calculation will be faster compared to modal solution described above. Moreover, in spatial domain, user can use the measured transfer function directly and eliminate the approximations made in curve fitting. The single point spatial domain solution is given in the Appendix A. 4.4 SLE Example In this section a theoretical example of surface location error calculation is given. The Al 7075 plate with 10.5 mm thickness, which is used in roughing chatter tests in chapter 5, is taken as a case scenario in order to apply multipoint approach. Along its z direction, it has changing dynamics due to different effective masses and stiffnesses. Chapter 4. Time Domain Simulation 65 Figure 4-8: Front view of the Al 7075 plate, with four measurement points indicated along its right corner The plate’s dimensions are given in Figure 4-8. The four points along z direction are indicated with thick dots and the distances are also given. The corresponding mass normalized mode shape values for those points and modal parameters are given in Table 4-4 for the y direction. The other directions of the plate are assumed to be very rigid compared to y direction. The reader is advised to check chapter 5 for more detailed analysis of the plate. Table 4-4: Modal parameters of Al 7075 plate along its y direction Mode (Hz) Modal Damping Ratio (%) Mode shape values for points: 1,2,3,4 731.17 1.00 3.315, 2.030, 0.9030, 0.266 1718.54 0.22 4.0011, 3.3502, 1.9740, 0.767 4252.81 1.30 1.8855, -1.3957, -2.8167, -1.5987 A regular two fluted cylindrical solid cutter with sharp edge is selected for the finishing operation. The cutter parameters are given in Table 4-5. The flexibility of the cutter is Chapter 4. Time Domain Simulation 66 assumed to be rigid compared to workpiece, thus its dynamics are not included in this example. Table 4-5: Parameters of cylindrical cutter Diameter [mm] Teeth Helix 1 [deg] Helix 2 [deg] Pitch 1 [deg] Pitch 2 [deg] 20 2 30 30 180 180 The modal matrix for the multipoint case is expressed as: 1 1 1 wp 1 1 1 1 1 0.266 0.767 1.5987 0.9030 1.9740 2.8167ˆ 2.030 3.3502 1.3957 3.315 4.0011 1.8855 q q q q q q q q 0 0 0 0 0 0 U 1 3 3q q , (4.39) where q indicates the points in cut. The levels that are already measured are written in the modal matrix, the points in between can be calculated using linear interpolation regarding the differential depth of cut value taken by the user. The levels that need to be considered are related with the axial depth of cut, thus the modal matrix should be updated for different axial depth of cuts. The radial depth of cut is 0.5 mm; the cutting and the edge force coefficients are c T 2 rta {168 796 222} N/mmK and e T rta {30.8 27.7 1.4} N/mmK , respectively. In order to apply SLE algorithm for given conditions, the stability of the operation should be known. The multipoint semi discrete chatter stability simulation for down milling is conducted for one stability pocket as shown in Figure 4-9. The details of semi discrete time domain chatter simula- Chapter 4. Time Domain Simulation 67 tion are given in chapter 5. The stable cutting region at 10 mm ADOC, and 16400 rpm to 17200 rpm is taken for the SLE analysis. Figure 4-9: Chatter stability for low immersion down milling operation and the region for SLE calculation is indicated with dashed line. The cutting forces for 10 mm ADOC can be calculated by integrating the contri- bution of each axial level and flute. It should be noted that since the vibration is assumed to take place only in y direction, the cross interactions between y and two other direc- tions are zero, therefore only y force is calculated as shown in Figure 4-10a. The finishing cut is highly intermittent, thus a wide frequency range scanned by tooth passing frequency harmonics as illustrated in Figure 4-10b. The excitation frequencies are expressed in terms of spindle frequency harmonics. The force is assumed to be constant for the same axial depth of cut, neglecting the effect of regeneration and change of radial immersion for the stable cutting case. The frequency content of force should be adjusted according to the spindle frequency. 1.4 1.45 1.5 1.55 1.6 1.65 1.7 1.75 1.8 x 10 4 0 5 10 15 20 25 30 Spindle speed [rev/min] C ri ti c a l a x ia l d e p th o f c u t [m m ] Multipoint Chatter Stability The region for SLE Chapter 4. Time Domain Simulation 68 a) b) Figure 4-10: a) Force in y direction for one spindle rotation, b) Fast Fourier Transform (FFT) of the force signal The algorithm for SLE calculation is shown in Figure 4-11. After getting the dif- ferential forces for each axial level, the numerical way of calculating Fourier coefficients is adopted for speed concerns, thus FFT with number (powers of 2) of points is applied. Consequently, modal transformation is used to get modal forces. The relevant frequencies of force spectrum are employed to construct (i )M . The attention should be given to take the complex conjugate mirror of (i )M after 1 2 points. In frequency domain, modal forces and the transfer function are multiplied to get the modal displace- ments. Applying Inverse Fast Fourier Transform (IFFT), the modal displacements in time domain is attained, and then SLE calculation can be conducted. For the instants and axial depths that the tool leaves the cut, corresponding modal displacement is retrieved, and consequently physical displacement is obtained by multiplying with the modal matrix vector of axial level that is in contact. After calculating the errors for each axial level, spindle speed is varied and new calculation is done. For 81 spindle speed points, 201 differential axial elements, and 512 Fourier points ( υ ), the total elapsed time is 2.493 seconds with a 2.4 GHz Quad CPU speed and 3.5 GB RAM PC. 0 100 200 300 0 50 100 150 200 250 F o rc e i n y d ir e c ti o n [ N ] Rotation of the spindle [deg] 0 5 10 15 20 0 5 10 15 20 25 30 35 40 45 F o rc e i n y d ir e c ti o n [ N ] Frequency in terms of cutter rotation frequency Chapter 4. Time Domain Simulation 69 Figure 4-11: Fast SLE calculation algorithm The SLE graph for corresponding spindle speeds and axial depths is given in Figure 4-12. It is observed that when the harmonics of tooth passing frequency starts to excite the mode at 1718.54 Hz, the errors become dominant. Since change of radial immersion with vibrations is not included, the values that are even bigger than radial depth of cut do exist in the resonance region. The positive values of error means under- cut, and negative values are overcut for down milling operation. The detailed SLE along the axial depth for 16580 rpm is given in Figure 4-12. After 16580 rpm, the overcut happens in finishing pass which is not tolerable; hence the best operating speeds are in between 16400 rpm and 16580 rpm. Obtain y component of forces for each differential axial depth, take FFT of differential forces Apply modal transformation to differential forces Calculate M (iω) for each regarding point of force FFT Multiply force and transfer function, take the IFFT to obtain modal displacements Find the instants and axial level that cutter leaves the surface, calculate physical displacements and store SLE Change spindle speed Chapter 4. Time Domain Simulation 70 Figure 4-12: Three dimensional SLE graph for low immersion down milling of thin plate and SLE at 16580 rpm 0 50 100 150 0 1 2 3 4 5 6 7 8 9 10 SLE at 16580 rpm A x ia l D e p th f ro m t h e t ip [ m m ] Surface Location Error [micrometer] Chapter 5. Chatter Stability in Milling 71 5. Chatter Stability in Milling 5.1 Introduction In this chapter time domain semi discretization and frequency domain formulations are used to predict chatter stability diagrams. The algorithms are applied for real case scenarios for thin wall machining and low speed cutting. Two examples are given in thin wall machining to illustrate the variation of dynamics as material is removed from aluminum plate. Due to complexities in low speed cutting, stiff workpiece and cutter are used for the tests to identify low frequencies effect on chatter. 5.2 Semi Discrete Time Domain Chatter Stability The time periodic delayed differential equation was solved for each time step using semi discretization (SD) method in chapter 4. In this section, the stability of the process is investigated semi analytically by applying extended Floquet theory due to presence of time periodic coefficients in the differential equation [53]. In this way, the stability of the process is obtained quicker with quantitative measures, contrary to other time domain simulations which contain qualitative decision making algorithms. It is assumed that static forces have no effect on the stability since change of radial immersion is neglected throughout the formulation. The current state transition equation is modified and written in the form of 1i i i Θ ZΘ : , , , 1 ,0, , ,1, , ( 1) ( 1) ( )( ) .( ) .. .. .. ( ). ( ) ( ) j r j r j r ii i i j r i j r i i k i k i k tt t t t t ΘΘ C 0 0 D D Θ I 0 0 0 0 0 I 0 0 0 Θ0 0 0 I 0 0 Θ Θ0 0 0 0 I 0 . (5.1) If the system is periodic, ( ) ( )t t T Θ Θ , then Floquet theory can be applied to the problem. The Floquet transition matrix for the discrete system is calculated as: Chapter 5. Chatter Stability in Milling 72 1 1 0m mΦ Z Z Z Z , (5.2) where int( / / 2)m t p is the number of points taken to constitute least common period . The least common periods for different cutters are defined in the Table 5-1. Table 5-1: Least common period and frequency for different cutters Pitch and Helix Angles Least common period Least common frequency p Uniform pitch and helix angle 60 Nn p 2 60 n N Alternating pitch and/or alternating helix angle i.e. p,1 p,2 p,3 p,4, , , ε : number of different angles 60 ( / )N n p 2 60 n N Random pitch and/or helix angle, axially changing helix angle 60 n p 2 60 n The stability of the process is assessed by plotting the Floquet transition matrix’s eigenvalues, or as formally called the characteristic multipliers, in the complex plane. If the eigenvalues are inside the unit circle, then the process is stable. The critically stable point is obtained by changing the axial depth of cut inside the algorithm for one spindle speed and stability lobes are constructed when all the corresponding axial depths and spindle speeds are plotted. In here several computational considerations are worth to mention. Firstly, due to large size of the coefficient matrix iZ , calculation of Floquet transition matrix is time consuming. Use of sparse matrix calculation algorithms due to zeros available in the matrix can bring apparent reduction to simulation time. Moreover, the number of discre- tized points is proportional to least common period and period is dependent on the Chapter 5. Chatter Stability in Milling 73 spindle speed. Varying the number of points regards to spindle speed, increases the simulation time at high speed end without the loss of accuracy. 5.3 Frequency Domain Chatter Stability Frequency domain based methods do not require modal identification and they are computationally advantageous to use for chatter stability prediction. On the other hand, for intermittent cutting conditions, determining the number of excitation frequency harmonics to consider, and avoiding the truncation errors are complications the frequen- cy domain solution brings within. Similar to SD time domain stability method, only the dynamic forces have effect on the regenerative vibrations. The dynamic part of the modal milling force is defined as: T tT d cd cd ed xyz 1 2 TT wp ˆ ( ) ( ( ) ( ) ( ) ( ) ( ) ( )) ˆ t t t t t t t dz F U U f f f U . (5.3) The time periodic coefficients such as cd cd ed 1 2( ), ( ), ( )t t tf f f are expressed by Fourier series [2]: cd cd cd cd ed ed 1 1 2 2( ) , ( ) , ( ) p p pih t ih t ih t h h r h h h t e t e t e f f f f f f , (5.4) where the Fourier coefficients are 0 1 ( ) p ih t h h t e dt f f for cd cd ed 1 2{ , , }h h h hf f f f . The common frequency p and the common period are defined in Table 5-1. The harmon- ics are denoted by h . The Fourier coefficients cd cd ed 1 2{ , , }h h h hf f f f are given in [2]. The displacement is expressed in terms of transfer function and dynamic force [22]: T d xyz( ) ( ) ( )t i t M FU . (5.5) Chapter 5. Chatter Stability in Milling 74 Eqs. (5.4) and (5.5) indicate that the forces consist of the excitation frequency p and chatter frequency c due to vibration term ( )t at the marginally stable condition [2]. The time domain equation is transformed into frequency domain when Fourier series expansion is used as shown by Merdol and Altintas [22]. Due to computational difficul- ties, the infinite number of Fourier harmonics are truncated and finite sized Toeplitz matrix is constructed: (0,0) ( 1,1) (1, 1) o o o (1,0) (0,1) (2, 1) o o o ( 1,0) ( 2,1) (0, 1) c o o o (( ), ) o 0 1 -1 - ( , ) h h p p k k k m m p M M M M M M Η M M M M 0 1 1 - h h m m . (5.6) The oriented transfer function ( , )c p Η is t wp t wp( )(2 1)-by-( )(2 1)h hm m m m m m matrix where, cd cd ed 1 c 2 c c c{ ( , ), ( , ), ( , )} ( , )p p p p Η Η Η Η and hm is the maximum number of harmonics. The elements inside the matrix have the form: (( ), ) o ( ) p k k p k ki M Mf , (5.7) where multi frequency excitation is represented with chatter frequency and excitation frequency by ck pk for 0, 1, 2,k Due to harmonics of the excitation frequency, the modal transfer function ( )kiM is defined for positive and negative frequency domain. The accuracy of the solution depends on how well the forces are represented by the truncated Toeplitz matrix. With sufficient amount of harmonics, the error can be minimized. The intermittency of the forces determines the maximum number of harmon- ics. For continuous cutting forces, zero order solution described by Budak and Altintas [16] gives good results. The Toeplitz representation of zero order solution is obtained by just taking the average Fourier term: Chapter 5. Chatter Stability in Milling 75 (0,0) c o( ) Η M , (5.8) where the excitation frequency no longer affects the chatter equation, hence oriented transfer function becomes (0,0) o 0 c( )i M Mf . The equation of motion now can be expressed in frequency domain for critically stable case: c T d xyz c T t cd cd ed T d 1 c 2 c c c xyz cT wp ( ) ˆ ( ( , ) ( , ) ( , )) ( ) ˆ i T p p p i e i dz i F Η Η Η F U U U U ’ (5.9) which can be simplified as an eigenvalue problem: c T t cd cd ed T d 1 c 2 c c c xyz cT wp ˆ ( ( , ) ( , ) ( , )) ( ) 0 ˆ i T p p pe i dz i I Η Η Η F U U U . (5.10) The characteristic equation for peripheral milling is formulated for the non trivial solutions of the eigenvalue problem: c T t cd cd ed 1 c 2 c cT wp ˆ det ( ( , ) ( , ) ( , )) 0 ˆ i T p p pe i dz I Η Η Η U U . (5.11) The product of delay term and the corresponding oriented transfer function is defined for three cases: c cc c , cd 2 c cd cd 2 2 c 1 cd 2 c 1 ˆ ( , ) uniform helix and pitch ˆ( , ) ( , ) uniform helix, non-uniform pitch ˆ ( , ) non-uniform helix and j j r i T p N i Ti T c p p j q i T p r e e e e Η Η Η Η 1 pitch N j . (5.12) Chapter 5. Chatter Stability in Milling 76 The characteristic equation’s poles are the zeros of the system. In order to deter- mine the stability, the location of the zeros should be assessed in the complex plane. Since the presence of the delay term brings infinitely many zeros, direct calculation is not viable. Spindle speed dependency of multi-frequency approach and process damping term make analytical stability solution inapplicable. For generality purposes Nyquist stability criterion is used in this thesis. In the literature, Nyquist stability criterion had been used in the calculations by Minis and Yanushevsky [15] based on the findings of Wereley [60] for the linear time periodic systems. Moreover, the validity of the method for multiple delays had been proved by Freedman et al. [61]. Figure 5-1: Nyquist plot of the characteristic equation The Nyquist stability criterion sweeps certain frequency range and looks clock- wise encirclements of the complex plane origin as illustrated in Figure 5-1. If it does encircle the origin, the criterion says the system is unstable. Similar to semi- discretization, critically stable points are searched by changing the ADOC and spindle speed. -5 0 5 10 15 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 Real Part of the Characteristic Equation Im a g in a ry P a rt o f th e C h a ra c te ri s ti c E q u a ti o n Chapter 5. Chatter Stability in Milling 77 5.4 Chatter Tests Chatter stability tests are conducted to validate the developed theoretical models. In the first example low immersion peripheral milling is conducted for low axial depths of cut. In the second example, process damping model is tested by cutting rigid steel block, where low frequency column mode and high frequency spindle modes are the major vibration sources in the system. In both examples the stability predictions are done by semi discretization and traditional zero order analytical (ZOA) frequency domain solution. 5.4.1 Thin Wall Roughing Test for Chatter Analysis Aluminum 7075 plate is cut as a thin wall with a thick bottom to minimize damp- ing and stiffness nonlinearities coming from the clamping. The dimensions of the plate are given in Figure 5-2. The thick bottom of the plate is clamped on the vise with a certain torque for repeatable measurements. The vise is fixed on the Kistler 9255B dynamometer for force measurements in three axes. Figure 5-2: Side and front views of the plate with dimensions and measurement locations Chapter 5. Chatter Stability in Milling 78 The dynamics of the plate are obtained through experimental modal analysis. For the roughing measurements, small size Dytran 3225F1 (1.284 g mass) accelerometer and medium size Dytran 5800B4 impact hammer are used. The most compliant point in y direction (along the thin section) is approximately hundred times more flexible than the tip of the workpiece in x direction (along the feed), thus x direction is considered to be rigid. In order to get the first three dominant modes that found to be causing vibrations, the frequency bandwidth of the measurement is set to 5000 Hz, and steel tip is selected for that range. Due to changing dynamics of the plate along feed and axial directions of the machine coordinate system, 11 points are marked and tested by keeping the impact position same and moving the accelerometer. It should be noted that the number of measurement points, and their locations may vary with respect to the modes and mode shapes of the workpiece. Due to dominant first bending and torsional modes seen on this plate, two points are placed at the corners, and one point placed in the middle along the feed direction in order to see torsional behavior; axial level is discretized by 4 points on the edges and 3 points on the middle region to measure both torsional and bending behavior along the z axis. The high flexibility of the workpiece impedes impact testing from the tip of the structure due to multiple hits. Down left corner is selected for the hitting location as illustrated with the green dot in Figure 5-2, where all the modes in the set frequency range are properly excited. The errors coming from the accelerometer mass are eliminated using method proposed by Cakar and Sanliturk [62] as presented in Appendix D. After calculating the modified frequency response functions for each point, modal parameters of well separated modes are identified globally using linear least squares fit algorithm in CUTPRO. The modes and the mode shapes are given in Table 5-2. Chapter 5. Chatter Stability in Milling 79 Table 5-2: Modal Parameters of the 10.5 mm thick plate Mode (Hz) Modal Damping Ratio (%) Mode shape 731.17 (1 st Bending mode) 1.00 1718.54 (1 st Torsional Mode) 0.22 4252.81 (2 nd Bending mode) 1.30 Chapter 5. Chatter Stability in Milling 80 In Figure 5-3a-c, the change of dynamics in feed direction and axial directions are demonstrated with the calculated direct frequency response functions (FRF) at the points 1, 2, 3 and 4 as marked in Figure 5-2. Using the relation 2 2 ,wp , ,wp y, , ,wp/ij n j i jk u between modes and mass normalized mode shapes, effective stiffness ,wpijk is calculated, where i is the mode number and j is the point number. Similar to modal stiffness approach used in commercial software packages, it indicates the flexibility of the mode at picked locations; hence it is inversely proportional to chatter stability limit. The Figure 5-3a shows that the first bending and torsional modes become 42% and 14% stiffer respectively, when descended 11 mm from the tip of the workpiece. Along the feed direction, the change in torsional mode is obvious as stiffness becomes 364 times more when moved from point 1 to 2. In the tested region, each location has different effective stiffness values, thus many stability lobes should be calculated separately to find the global stability. In order the overcome this time consuming occasion, stability lobes are constructed for the worst case by selecting most compliant points for each separate mode as given in Table 5-8. A 4 fluted helical end mill with 19.05 mm diameter is used for the tests. The he- lix angle of the tool is 35 . The tool tip FRFs at x and y directions are given in Figure 5-11. The low frequency modes are from the spindle dynamics, and high frequency modes are bending modes of the tool in x and y direction. The tool is relative- ly stiff compared to workpiece, but its dynamics are added for complete analysis of the process. The modal parameters of the tool are given in Table 5-8. Chapter 5. Chatter Stability in Milling 81 a) b) c) d) Figure 5-3: Frequency response functions of workpiece at different locations and tool tip a) first bending mode of the workpiece, 6 611,wp 12,wp1.93 10 N/m, 2.18 10 N/m,k k 6 6 13,wp 14,wp2.61 10 N/m, 2.74 10 N/mk k b) first torsional mode of the workpiece, 6 9 6 21,wp 22,wp 23,wp7.52 10 N/m, 2.34 10 N/m, 8.58 10 N/m,k k k 9 24,wp 1.92 10 N/mk c) second bending mode of the workpiece 831,wp 1.97 10 N/m,k 8 31,wp 1.97 10 N/m,k 7 32,wp 6.90 10 N/m,k 8 8 33,wp 34,wp2.40 10 N/m, 1.67 10 N/mk k d) tool’s FRFs at x and y directions 680 700 720 740 760 780 -1 0 1 x 10 -5 R e a l p a rt [ m /N ] 680 700 720 740 760 780 -4 -2 0 x 10 -5 Frequency in HzI m a g in a ry p a rt [ m /N ] Point 1 Point 2 Point 3 Point 4 1680 1700 1720 1740 -2 0 2 x 10 -5 R e a l p a rt [ m /N ] 1680 1700 1720 1740 -4 -2 0 x 10 -5 Frequency in HzI m a g in a ry p a rt [ m /N ] 4200 4250 4300 -5 0 5 x 10 -7 R e a l p a rt [ m /N ] 4200 4250 4300 -5 -4 -3 -2 -1 x 10 -7 Frequency in HzI m a g in a ry p a rt [ m /N ] 500 1000 1500 2000 2500 -1 0 1 2 x 10 -7 R e a l p a rt [ m /N ] 500 1000 1500 2000 2500 3000 -3 -2 -1 0 x 10 -7 Frequency in Hz Im a g in a ry p a rt [ m /N ] Tool y Tool x Chapter 5. Chatter Stability in Milling 82 Table 5-3: Modal parameters for tool and workpiece Structure-Direction Mode (Hz) Effective Stiffness (N/m) Modal Damping Ratio (%) Tool x 610.93, 854.97, 2379.92 1.25e8, 1.51e8, 1.13e8 4.48, 2.25, 1.94 Tool y 2377.82 1.25e8 1.41 Workpiece y 731.17, 1718.54, 4252.81 1.93e6, 7.52e6, 6.90e7 1.00,0.22,1.30 Five axis Mori Seiki NMW5000 high spindle speed machine is used in chatter tests. The experimental setup is demonstrated in Figure 5-4. Cutting coefficients for aluminum 7075 are taken as c T 2 rta {168 796 222} N/mmK . The radial depth is set to 2 mm, where the radial immersion is 10.50%. The feed per tooth is 0.1 mm. The force and sound measurements are collected. The tested conditions are given in Table 5-4. Figure 5-4: Setup for stability test Chapter 5. Chatter Stability in Milling 83 Table 5-4: Tested conditions and results Spindle Speed (rpm) Axial Depth of Cut (mm) Distance from the workpiece tip (mm) Chatter (visual inspection)-frequency 9250 0.5 0 Marks on the edges (1727 Hz) 10300 0.5 0.5 No 10300 1 1 Slight marks on the edges (1764 Hz) 11000 2 2 No 11400 1 4 Marks on whole path (740 Hz) 7500 1 5 Marks on the edge (1785 Hz) 7000 1 6 Marks on the edges (1746 Hz) 7900 1 7 No 8350 1.5 8 No 8640 1.5 9.5 No Inputting the dynamic parameters and cutting conditions provided, stability lobes are plotted using both traditional zero order frequency domain analytical solution and semi discretization method. In order to be consistent with the current single point approach used in analytical solution, dynamics and forces are lumped into a single point in semi discretization approach as discussed in Appendix A. Since the axial depth of cut is very low, this approach does not bring excessive errors due to varying dynamics along the tool axis direction. Assuming all modes are orthogonal to each other, the lumped mode shape matrix for the most flexible scenario is expressed as: t wp 3 7 0.3432 0.4372 1.4085 0 0 0 0 0 0 0 1.3362 3.3061 3.9385 3.2168 0 0 0 0 0 0 0 U U . (5.13) Chapter 5. Chatter Stability in Milling 84 The predictions and the experimental results showed satisfactory match in tested conditions. The high helix angle and four fluted cutter provided good continuity in the cutting force, thus semi discretization and zero order frequency domain solutions gave similar results except the Flip bifurcation lobe at 8000 rpm. Figure 5-5: Predicted and measured chatter stability results for the most flexible condition The measured forces and simulated displacements by semi discrete time domain solution for y direction are compared for three points denoted by A, B, C as in Figure 5-5. The frequency spectrums are calculated in order to correlate between measurements and simulations. The force measurement at point A, demonstrated in Figure 5-6a, shows that the force grows at the start and the end of the cutting process, damps out at the middle. This behavior corresponds to torsional mode of the plate. The marks on the edges of the workpiece as reported in Table 5-4, validates the chatter caused by torsional mode. The frequency spectrum of the force applied on the edge shows a dominant peak at 1725 Hz, which is independent from tooth passing frequency 617 Hz and its harmon- ics. Using top corner location’s dynamics, semi discrete time domain simulation is conducted for vibration displacements. The edge force coefficients are taken as 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 x 10 4 0 0.5 1 1.5 2 2.5 3 3.5 4 Spindle speed [rev/min] C ri ti c a l a x ia l d e p th o f c u t [m m ] Semi Discretization Zero Order Frequency Domain Experiment: Chatter Experiment: No chatter A B C Chapter 5. Chatter Stability in Milling 85 e T rta {30.8 27.7 1.4} N/mmK . In Figure 5-6c, the time domain simulation results show agreement with the experiment predicting the growth of the vibrations in y direc- tion. The calculated chatter frequency, 1727 Hz, illustrated in Figure 5-6d, complies with the measured frequency. a) c) b) d) Figure 5-6: Y direction measurement and simulation for point A, a) Force measurement b) FFT of force measurement c) simulated vibration displacement d) FFT of displacement The measurement for the stable cutting condition at point B shows that the force is steady along the feed direction. The axial depth of cut is quadrupled respect to point A, but the forces are ten times more than the steady region forces of point A. This should be attributed to heavy forced vibrations due to excitation of first bending mode at 731.17 0 0.5 1 1.5 -600 -400 -200 0 200 400 Time [s] F o rc e [ N ] 0 0.02 0.04 0.06 0.08 0.1 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 Time [s] D is p la c e m e n t [m m ] 500 1000 1500 2000 2500 0 20 40 60 80 100 120 Frequency [Hz] F o rc e [ N ] 500 1000 1500 2000 2500 0.005 0.01 0.015 0.02 0.025 0.03 D is p la c e m e n t [m m ] Frequency [Hz] 1725 Hz 1727 Hz 617 Hz 617 Hz Chapter 5. Chatter Stability in Milling 86 Hz by tooth passing frequency at 733.33 Hz. The damped forced vibration simulation is presented in Figure 5-7c. At steady state, the displacement, disregarding the change of force due to immersion angle in the simulation, approaches to 1.4 mm which is 70% of the radial depth. It can be concluded that chatter free high axial depth of cuts can be achieved by running the system at resonant frequencies of the structures with high forces and surface location errors as consequences. a) c) b) d) Figure 5-7: Y direction measurement and simulation for point B, a) Force measurement b) FFT of force measurement c) simulated vibration displacement d) FFT of displacement The point C is an example of chatter triggered by the first bending mode. The force is constantly growing throughout the feed as shown in Figure 5-8a. The chatter marks are left on the whole tool path. The frequency spectrum of the selected data as presented in Figure 5-8b indicates that chatter occurs at 740 Hz, and the tooth passing frequency at 0 0.5 1 -2000 -1000 0 1000 2000 Time [s] F o rc e [ N ] 0 0.02 0.04 0.06 0.08 0.1 -2 -1 0 1 2 Time [s] D is p la c e m e n t [m m ] 0 500 1000 1500 2000 2500 0 500 1000 1500 Frequency [Hz] F o rc e [ N ] 0 500 1000 1500 2000 2500 0 0.5 1 1.5 D is p la c e m e n t [m m ] Frequency [Hz] 732 Hz 1464 Hz 733 Hz Chapter 5. Chatter Stability in Milling 87 760 Hz is also amplified by the forced vibrations. Availability of these two close fre- quencies creates undulated behavior in force due to “beat” phenomenon. The similar response is identified in semi discrete time domain simulation as shown in Figure 5-8c. The stability solution (Figure 5-5) recognizes this speed and axial depth combination as critically stable point with chatter frequency very close to the first bending mode. In order to identify the effect of chatter, the frequency spectrums of the time domain solutions with and without regeneration term are plotted together as shown in Figure 5-8d. The contribution of chatter vibration is clearly seen as the peak at 731.5 Hz is amplified with the regeneration. In here, the possible reason behind the small difference between simulated and measured chatter frequencies should be mentioned. The cut at point C is conducted after removing 4 mm material axially as reported in Table 5-4. Significant decrease in mass of the structure and comparatively slight change in stiffness due to stepped geometry tends to increase the frequency of first bending mode. Moreo- ver, mass normalized mode shapes which are indirectly proportional to chatter stability increases as mass become less. For tests with higher axial depth of cuts, the changes become dominant impairing the accuracy of the stability predictions done for the blank material. Chapter 5. Chatter Stability in Milling 88 a) c) b) d) Figure 5-8: Y direction measurement and simulation for point C, a) Force measurement b) FFT of force measurement c) simulated vibration displacement d) FFT of displacement 5.4.2 Thin Wall Finishing Test for Chatter Analysis The same Al 7075 part is cut to 5 mm thickness and 9 points are determined on the plate for modal testing. The dimensions of the plate and the measured locations are shown with red points, and fixed hitting location is given with green point. Since the thin wall’s mass is relatively low, PCB laser vibrometer is used for the tests. 0 0.5 1 -1000 0 1000 2000 Time [s] F o rc e [ N ] 0 0.1 0.2 0.3 0.4 -0.4 -0.2 0 0.2 0.4 0.6 Time [s] D is p la c e m e n t [m m ] 600 700 800 900 0 500 1000 1500 Frequency [Hz] F o rc e [ N ] 600 700 800 900 0 0.05 0.1 0.15 0.2 D is p la c e m e n t [m m ] Frequency [Hz] With regeneration Without regeneration 740 Hz 760 Hz 731.5 Hz 760 Hz Chapter 5. Chatter Stability in Milling 89 Figure 5-9: Side and front views of the plate with designated modal test points The frequency range is set to 5000Hz, and 4 modes are identified. Severe changes observed in the dynamics of the structure compared to the 10.5 mm thick plate. The natural frequencies, modal damping ratios are dropped, mass normalized mode shape values increased. Identified modal parameters and effective stiffness values correspond- ing to each mode’s most compliant point are presented in Table 5-5. Two fluted finishing tool with 10 mm diameter is used for the tests. The helix angle of the tool is 30 . The dynamic parameters of the tool tip at x and y directions are given in Table 5-6. Chapter 5. Chatter Stability in Milling 90 Table 5-5: Modal Parameters of the 5 mm thick plate Mode (Hz) Modal Damping Ratio (%) Most Flexible Effective Stiffness (N/m) Mode shape 440.30 (1 st Bending mode) 0.12 3.04e5 972.87 (1 st Torsional mode) 0.03 8.78e5 2549.18 (Bending- torsional mode) 0.78 8.64e6 3584.49 (2 nd Torsional mode) 0.196 1.40e7 Chapter 5. Chatter Stability in Milling 91 Table 5-6: Modal parameters of the tool Direction Mode (Hz) Effective Stiffness (N/m) Modal Damping Ratio (%) Tool x 2455.20, 3244.71, 3826.76 5.52e7, 2.63e7, 1.41e7 1.78, 1.75, 3.06 Tool y 2478.73, 3226.11, 3810.10 4.99e7, 1.42e7, 2.09e7 1.58, 4.17, 2.58 The radial depth is set to 0.5 mm which corresponds to 5% radial immersion. The feed per tooth is 0.1 mm. The tested conditions are given in Table 5-7. The chatter condition is inspected visually from the cut surface, and the identified dominant frequencies are noted. A wide frequency range is excited by the tooth passing frequency and its multiple harmonics. Due to high number of modes available in the structure and in the measure- ment system (dynamometer), detecting chatter through frequency spectrum is a challenge. The stability lobes are plotted with analytical zero order frequency domain and semi discretization methods with single point approach. The mass normalized mode shape matrix is expressed as: 3 10 2.08 3.97 6.39 0 0 0 0 0 0 0 0 0 0 2.20 5.38 5.23 5.02 6.52 5.45 6.02 0 0 0 0 0 0 0 0 0 0 . (5.14) The theory and the experiments show good match around the tested stability pocket as shown in Figure 5-10. The stable cut at 12600 rpm created high forced vibrations with undercut surface finish, thus decreasing the radial immersion and increasing the stability at that speed. This situation is observed at 3 mm axial depth of cut where slight chatter marks were seen on the tool path. Another important point to discuss is the availability of chatter frequency around 3050-3150 Hz which is not detected by the modal tests done at the static condition. However, tool’s bending modes at 3244 and 3226 can be de- Chapter 5. Chatter Stability in Milling 92 creased to those frequencies due to gyroscopic effects at high speeds and variations in spindle dynamics. Table 5-7: Tested conditions for finishing operation Spindle Speed (rpm) Axial Depth of Cut (mm) Distance from the workpiece tip (mm) Chatter (visual inspection)-frequency 12600 1 0 No 11800 0.5 1 Slight marks on the edges (3070 Hz) 9400 0.5 1.5 Chatter marks on the edges 9600 0.5 2 Slight marks on whole path (458 Hz) 12600 1.5 2.5 No 12600 2 4 No 12600 3 6 Slight marks 133380 1 9 Marks on edges (3121 Hz) 14000 0.5 10 Heavy chatter (466 Hz) Chapter 5. Chatter Stability in Milling 93 Figure 5-10: Predicted stability lobes for the most flexible case and measured points are presented for 5% radial immersion thin wall machining 0.9 1 1.1 1.2 1.3 1.4 1.5 x 10 4 0 0.5 1 1.5 2 2.5 3 3.5 Spindle speed [rev/min] C ri ti c a l a x ia l d e p th o f c u t [m m ] Semi Discretization Zeroth Order Frequency Domain Experiment: Chatter Experiment: No Chatter Experiment: Uncertain Chapter 5. Chatter Stability in Milling 94 5.4.3 Low Speed Cutting Test for Rigid Workpiece AISI1045 rigid steel block is cut by 7 fluted inserted face mill, with a 50 mm di- ameter, and 0.8 mm nose radius on the inserts. The cutting force coefficients are identified mechanistically, by changing feed, measuring forces and fitting a line to the measured points. The cutting force coefficient vector is c T 2 rta {970 2184 296} N/mmK . The dynamics of the cutter is measured with PCB 353B31 model big size accelerometer and Kistler mid-size hammer to capture low frequency dynamics properly. Approximately 1000 N force is applied to excite spindle and machine column modes. Two orthogonal directions, namely x and y, are assumed to be the principal directions and frequency response functions are measured directly at the tool tip. The modal parameters are given in Table 5-8. Table 5-8: Modal parameters at the tool tip Measurement Direction Natural Frequency [Hz] Effective Stiffness [N/m] Damping Ratio (%) x 104.8, 131.3, 195.0, 257.2, 324.4, 419.6, 584.4, 739.4, 878.7 3.21e8, 1.02e9, 6.90e8, 2.65e9, 6.87e8, 5.64e8, 3.30e8, 1.83e8, 3.63e8 5.6, 6.7, 3.5, 6.5, 6.2, 1.8, 3.1, 3.3, 3.3 y 77.4, 128.7, 250.3, 404.8, 544.1, 792.0, 869.0 1.41e8, 1.88e9, 4.09e8, 3.39e8, 4.51e8, 3.02e8, 1.49e8 2.7, 1.9, 6.9, 3.9, 4.5, 3.8, 3.1 The low frequency mode in y direction at 77.4 Hz is coming from the machine head, the modes between 500 and 900 Hz belong to spindle itself. As shown in Figure 5-11, low frequency mode and high frequency spindle mode have similar flexibilities; where one expect to see similar contribution to chatter stability from both modes. Chapter 5. Chatter Stability in Milling 95 Figure 5-11: Two orthogonal direction FRFs measured at the tool tip The process damping parameters for 1045 steel are taken from Eynian et al. [33], as 13 3 sp 4.0 10 N/mK , 0.3 and w 130 μmL . For half immersion down milling stability lobes are constructed with three different approaches and stability is tested for the first two lobes of the low frequency mode at 77.4 Hz. In Figure 5-12, the traditional frequency domain approach shows that the lobes coming from spindle and the mode of column have similar flexibilities, thus giving the similar critical axial depths. The first two lobes of low frequency mode are at 400 and 200 rpm. When process damping is added on the model, these two lobes remain but the spindle modes damp out. The exponential increase of the spindle lobes is demonstrated clearly between 650-400 rpm. The difference between semi discretization and zero order frequency domain methods with process damping included is due to time varying elements in regenerative and process damping force matrices. It is observed that for this case semi discretization shows faster damping behavior compared to zero order frequency domain. 0 200 400 600 800 1000 1200 -1 -0.5 0 0.5 1 x 10 -7 R e a l p a rt [ m /N ] X direction FRF Y direction FRF 0 200 400 600 800 1000 1200 -1.5 -1 -0.5 0 x 10 -7 Frequency in Hz Im a g in a ry p a rt [ m /N ] Chapter 5. Chatter Stability in Milling 96 Figure 5-12: Chatter stability prediction using different algorithms and test results 5.4.4 Summary of Experiments In summary three tests are conducted to check the accuracy of algorithms. The thin wall machining tests showed the dominance of workpiece dynamics, and low speed cutting test showed the significance of low frequency modes due to process damping. All models are simulated with single point approach since axial depth of cut limits are very low for thin wall machining tests, and in low speed the dynamics of the structure does not vary significantly along axial depth of cut. The multi point formulation which takes variation of structural parameters along contact axis can be verified if a setup that consists of stiffer workpiece along with change of dynamics along axis. Machining titanium thin walled components can be an example to this, however due to significant damping at low speeds, process damping coefficient should be identified considering all the geometric factors affect process damping. 250 300 350 400 450 500 550 600 650 2 4 6 8 10 12 14 16 18 Spindle speed [rev/min] C ri ti c a l a x ia l d e p th o f c u t [m m ] Semi Discretization Z.O. Frequency Domain Z.O. Frequency Domain (without process damping) Experiment: Chatter Experiment: No Chatter Experiment: Uncertain Chapter 6. Conclusions 97 6. Conclusions The dynamics and mechanics of thin wall machining are studied in this thesis. A general force model for milling is utilized to predict mechanics for different cutter geometries at high and low speed machining operations. Interactions between the cutting forces and structural dynamics are studied, and second order linear periodic delayed differential equation is derived. This delayed differential equation is then solved using time domain semi discretization method to predict workpiece vibrations, cutter vibra- tions, and dynamic milling forces. For stable cutting cases, the periodic delayed differential equation is simplified to constant coefficient ordinary differential equation, and fast frequency domain solver for steady state vibrations is developed to predict surface location errors. The same delayed differential equation is then solved for chatter stability using semi discrete time domain stability and frequency domain stability prediction algorithms. Sample chatter tests are conducted for intermittent milling of thin plates at roughing and finishing stages, and half immersion low speed cutting of rigid steel block. The contributions of this thesis are summarized as follows: A general dynamic force model for cylindrical end milling is developed which considers process damping, non-uniform pitch and helix angles. The stability of the system is solved both in frequency and semi-discrete time domains. Cutter and workpiece engagement is discretized into differential points along the tool axis to include effect of changing dynamics of structures and non-uniform tooth spacing and helix angles. As opposed to solving the system dynamics in computationally costly local coordinate system, the discretized points are transformed into modal domain to reduce the size of equations. Semi discrete time domain solution is developed in or- der to predict dynamic forces, and tool and workpiece vibrations in modal domain. Chapter 6. Conclusions 98 The existing surface location error prediction algorithm which considered only the tool tip dynamics is extended to include multipoint contacts while considering the varying dynamics of workpiece and cutter along engagement axis. Semi discrete chatter stability and frequency domain solver which uses Nyquist stability criterion are extended for end mills with non-uniform pitch and helix. Future research in this area is recommended in order to advance mathematical modeling of the process, and improve the computational efficiency. In this thesis, semi discretization approach is utilized for time domain simulation. Being an accurate and stable solver, it is known that semi discretization can consume great amount of time, if the size of the equations being solved are huge. Full discretization methods need to be investigated to obtain more efficient solver. However, attention should be given to stability of the time domain solver such as A stability that considers time domain solv- er’s influence on state matrix, P stability for influence on delayed state matrix, and L stability which accounts fictitious vibrations of the solver. Furthermore, a comprehen- sive process damping model is required to include effects of clearance angle, edge radius and magnitude of vibrations imprinted on the surface as reported in the literature review. Especially for titanium alloys, a comprehensive study on process damping effects on stability and plowing forces should be conducted to utilize the thin wall machining theory presented in this thesis. Moreover, in case of process dynamics, it is known that as material removed from the thin webs, the dominant modes vary. A computationally effective solver employing plate theory can be developed. This would eliminate the excessive modal testing after every pass the cutter makes. Finally, the models developed in this thesis are only verified for few examples. The effect of multipoint approach for high axial depth of cuts, measurement of surface location errors, validation of non- uniform pitch and helix cutters can be carried forward as future research topics. 99 Bibliography [1] M. E. Martelotti, "Analysis of the Milling Process," Transactions of the ASME, vol. 63, p. 677, 1941. [2] S.D. Merdol, "Virtual Three-Axis Milling Process Simulation and Optimization," University of British Columbia, Vancouver, Ph. D. Thesis 2008. [3] F. Koenigsberger and A. J. P. Sabberwal, "An Investigation into the Cutting Force Pulsations During Milling Operations," International Journal of Machine Tool Design and Research, vol. 1, pp. 15-33, 1961. [4] E. Budak, Y. Altintas, and E.J. A. Armarego, "Prediction of Milling Force Coefficients from Orthogonal Cutting Data," Journal of Manufacturing Science and Engineering, vol. 118, no. 2, pp. 216-224, 1996. [5] W.A. Kline, R.E. DeVor, and I.A. Shareef, "The prediction of surface accuracy in end milling," Transactions of ASME, vol. 104, pp. 272-278, August 1982. [6] J. W. Sutherland and R. E. Devor, "An Improved Method for Cutting Force and Surface Error Prediction in Flexible End Milling Systems," Transactions of the ASME, Journal of Engineering for Industry, vol. 108, pp. 269-279, 1986. [7] D. Montgomery and Y. Altintas, "Mechanism of Cutting Force and Surface Generation in Dynamic Milling," Transactions of ASME, vol. 113, pp. 160-168, May 1991. [8] E. Budak, "Mechanics and Dynamics of Milling Thin Walled Structures ," The University of British Columbia, Vancouver, Ph.D. Thesis 1994. [9] E. Budak and Y. Altintas, "Modeling and avoidance of static form errors in peripheral milling of plates," Int. J. Mach. Tools Manufact., vol. 35, no. 3, pp. 459- 476, 1995. [10] T. Schmitz and Mann B.P., "Closed-form solutions for surface location error in milling," International Journal of Machine Tools & Manufacture, vol. 46, pp. 1369- 1377, 2006. [11] T.L. Schmitz and K.S. Smith, Machining Dynamics, Frequency Response to Improved Productivity, 1st ed. US: Springer , 2009. 100 [12] D. Bachrathy, T. Insperger, and G. Stepan, "Surface properties of the machined workpiece for helical mills," Machining Science and Technology, vol. 13, pp. 227- 245, 2009. [13] S.A. Tobias and W. Fishwick, "Theory of regenerative machine tool chatter," The Engineer, vol. 205, pp. 16-23, 1958. [14] J., Polacek, M. () Tlusty, "The stability of machine tools against self excited vibrations in machining," International Research in Production Engineering,ASME, pp. 465-474, 1963. [15] I. Minis and R. Yanushevsky, "A new theoretical approach for the prediction of machine tool chatter in milling," Journal of Engineering Industry, vol. 115, no. 1, pp. 1-8, February 1993. [16] Y. Altintas and E. Budak, "Analytical prediction of stability lobes in milling," CIRP Annals-Manufacturing Technology, vol. 44, no. 1, pp. 357-362, 1995. [17] S. Engin, "Mechanics and Dynamics of Milling with Generalized Geometry.," Istanbul Technical University, Istanbul, Ph.D. Thesis 1999. [18] Y. Altintas, S. Engin, and E. Budak, "Analytical stability prediction and design of variable pitch cutters," J. Manuf. Sci. Eng., vol. 121, no. 2, pp. 173-178, May 1999. [19] U. Bravo, O. Altuzarra, L.N. Lopez de Lacalle, J.A. Sanchez, and F.J. Campa, "Stability limits of milling considering the flexibility of the workpiece and the machine," International Journal of Machine Tools & Manufacture, vol. 45, pp. 1669-1680, 2005. [20] V. Thevenot, L. Arnaud, G. Dessein, and G. Cazenave-Larroche, "Influence of Material Removal on the Dynamic Behavior of Thin Walled Structures in Peripheral Milling," Machining Science and Technology, vol. 10, pp. 275-287, 2006. [21] E. Budak and Y. Altintas, "Analytical prediction of chatter in milling-Part I:general formulation," Journal of Dynamic Systems, Measurement and Control , Transactions of the ASME, vol. 120, no. 1, pp. 22-30, 1998. [22] S.D. Merdol and Y. Altintas, "Multi frequency solution of chatter stability for low 101 immersion milling," Journal of Manufacturing Science and Engineering, vol. 126, no. 3, pp. 459-466, August 2004. [23] B.P., Mann,B.P. Bayly, T.L. Schmitz, D.A. Peters, G. Stepan, and T. Insperger, "Effects of radial immersion and cutting direction on chatter instability in end milling," American Society of Mechanical Engineers, vol. 13, pp. 351-363, 2002. [24] T. Insperger and G. Stepan, "Semi-discretization method for delayed systems," International Journal For Numerical Methods in Engineering, vol. 55, pp. 503-518, 2002. [25] T., Mann,B.,Surmann, T., Stepan,G. Insperger, "On the chatter frequencies of milling process with runout," International Journal of Machine Tools and Manufacture, vol. 48, pp. 1081-1089, February 2008. [26] M., Munoa,J., Peigne, G., Insperger, T. Zatarain, "Analysis of the influence of mill helix angle on chatter stability," CIRP Annals-Manufacturing Technology, vol. 55, no. 1, pp. 365-368, 2006. [27] B.R., Man,B.P.,Young,K.A. Patel, "Uncharted islands of chatter instability in milling ," International Journal of Machine Tools and Manufacture, vol. 48, pp. 124-134, 2008. [28] N.D., Mann, B., Huyanan, S. Sims, "Analytical prediction of chatter stability for variable pitch and variable helix milling tools," Journal of Sound and Vibration, vol. 317, no. 3-5, pp. 664-686, 2008. [29] Z. Dombovari, Y. Altintas, and G. Stepan, "The effect of serration on mechanics and stability of milling cutters," International Journal of Machine Tools & Manufacture, vol. 50, pp. 511-520, 2010. [30] T.R. Sisson and R.L. Kegg, "An explanation of low speed chatter effects," ASME Journal of Engineering for Industry, vol. 91, pp. 951-958, November 1969. [31] D.W. Wu, "A New Approach of Formulating the Transfer Function for Dynamic Cutting Process," ASME J. Eng. Ind., vol. 111, pp. 37-47, 1989. [32] R.Y. Chiou and S.Y. Liang, "Chatter Stability of a Slender Cutting Tool in Turning with Wear Effect," International Journal of Machine Tools and Manufacture, vol. 102 38, pp. 315-327, 1998. [33] M. Eynian and Y. Altintas, "Analytical chatter stability of milling with rotating cutter dynamics at process damping speeds," ASME Journal of Manufacturing Science and Engineering, vol. 132, no. 2, April 2010. [34] E. Budak and L.T. Tunc, "A New Method for Identification and Modeling of Process Damping in Machining," Journal of Manufacturing Science and Engineering, vol. 131, 2009. [35] K. Ahmadi and F. Ismail, "Analytical stability lobes including nonlinear process damping effect on machining chatter," International Journal of Machine Tools and Manufacture, vol. 51, no. 4, pp. 296-308, January 2011. [36] S. Engin and Y. Altintas, "Mechanics and dynamics of general milling cutters.Part I: helical end mills," International Journal of Machine Tools & Manufacture, vol. 41, pp. 2195–2212, 2001. [37] S. Ratchev, S. Liu, A. Huang, and A.A. Becker, "Milling error prediction and compensation in machining of low-rigidity parts," International Journal of Machine Tools & Manufacture, vol. 44, pp. 1629-1641, 2004. [38] D. Merdol and Y. Altintas, "Virtual Simulation and Optimization of Milling Operations-Part I: Process Simulation," Journal of Manufacturing Science and Engineering, vol. 130, October 2008. [39] Y. Altintas and P. Lee, "Mechanics and dynamics of ball end milling," Transactions of ASME, Journal of Manufacturing and Science Engineering, vol. 120, pp. 684- 692, 1998. [40] H. J. Fu., R. E. DeVor, and S. G. Kapoor, "A mechanistic model for the prediction of the force system in face milling operations," Transactions of ASME, Journal of Engineering for Industry, no. 106, pp. 81-88, 1984. [41] S.A. Tobias, Machine Tool Vibration.: Blackie and Sons Ltd., 1965. [42] Y.C. Yen, A. Jain, and T. Altan, "A finite element analysis of orthogonal machining using different tool edge geometries," Journal of Material Processing Technology, vol. 146, no. 1, pp. 72-81, 2004. 103 [43] M. Eynian and Y. Altintas, "Chatter Stability of General Turning Operations with Process Damping," Transactions of the ASME, Journal of Manufacturing Science and Engineering, vol. 131, 2009. [44] Y. Altintas, M. Eynian, and H. Onozuka, "Identification of dynamic force coefficients and chatter stability with process damping," CIRP Annals- Manufacturing Technology, vol. 57, pp. 371-374, 2008. [45] E. Budak and L.T. Tunc, "Identification and Modeling of Process Damping in Turning and Milling Using a New Approach," Annals of CIRP, vol. 1, 2010. [46] Ramin Rahnama, Mozhdeh Sajjadi, and Simon S. Park, "Chatter suppression in micro end milling with process damping," Journal of Materials Processing Technology, vol. 209, pp. 5766-5776, 2009. [47] K. Ahmadi and F. Ismail, "Machining chatter in flank milling," International Journal of Machine Tools & Manufacture, vol. 50, pp. 75-85, 2010. [48] V. Sellmeier and F., Denkena B. Hackeloeer, "Process Damping in Milling- Measurement of Process Damping Forces for Chamfered Tools by Means of Electromagnetically Guided Spindle". [49] M. Eynian, "Chatter Stability of Turning and Milling with Process Damping," University of British Columbia, Vancouver, Ph. D. Thesis 2010. [50] D.J. Ewins, Modal Testing: theory, practice and application, 2nd ed. Baldock, Hertfordshire, England: Research Studies Press Ltd., 2000. [51] R.D. Cook, D.S. Malkus, M.E. Plesha, and R.J. Witt, Concepts and Applications of Finite Element Analysis, 4th ed. United States: John Wiley & Sons, Inc., 2002. [52] M. Geradin and D. Rixen, Mechanical Vibrations: Theory and Application to Structural Dynamics, 2nd ed. West Sussex, England: John Wiley & Sons Ltd., 1997. [53] T. Insperger, "Stability Analysis of Periodic Delay-Differential Equations Modeling Machine Tool Chatter," Budapest University of Technology and Economics, Budapest, Ph.D. Thesis 2002. [54] M. Campomanes, "Dynamics of Milling Flexible Structures," University of British 104 Columbia, Vancouver, M.A.Sc. Thesis 1998. [55] C. Henninger and P. Eberhard, "Improving the computational efficiency and accuracy of the semi-discretization method for periodic delay-differential equations," European Journal of Mechanics A/Solids, vol. 27, pp. 975-985, 2008. [56] T. Insperger, G. Stepan, and J. Turi, "On the higher-order semi discretizations for periodic delayed systems," Journal of Sound and Vibration, vol. 313, pp. 334-341, 2008. [57] T. Insperger, "Full-discretization and semi-discretization for milling stability prediction: Some comments," International Journal of Machine Tools and Manufacture (in press), 2010. [58] T. Insperger and G. Stepan, "Updated semi-discretization method for periodic delay-differential equations with discrete delay," International Journal for Numerical Methods in Engineering, vol. 61, pp. 117-141, 2004. [59] S. Turner, D. Merdol, Y. Altintas, and K. Ridgway, "Modelling of the stability of variable helix end mills," International Journal of Machine Tools & Manufacture, vol. 47, pp. 1410–1416, 2007. [60] Wereley, "Analysis and Control of Linear Periodically Time Varying Systems," Massachusetts Institute of Technology, Boston, Ph.D. Thesis 1991. [61] "Stability Criteria for a System Involving Two Time Delays," SIAM J. APPL. MATH., vol. 46, no. 4, pp. 552-560, August 1986. [62] O. Cakar and K.Y. Sanliturk, "Elimination of transducer mass loading effects from frequency response functions," Mechanical Systems and Signal Processing, vol. 19, no. 1, pp. 87-104, January 2005. Appendix A. Single Point Formulation of Milling Equation 105 Appendix A: Single Point Formulation of Milling Equation In this section multipoint approach which contains the varying dynamics of the structures along axial direction, is lumped to single point as used in traditional algo- rithms to calculate stability. In general, one tap test for each flexible direction is sufficient at the tool tip and aimed workpiece section, assuming that the dynamics do not change distinctively along the axial depth. Then forces are accumulated and exerted on this point to find dynamic response of both structures as presented in Figure A-1. This method fastens the simulation without losing much from the accuracy, if assumptions are valid throughout the conditions given. Figure A-1: Resultant force lumped on point k The single point approach starts with addition of the milling forces, which are formulated in chapter 3, to calculate resultant lumped force on a single point as ex- pressed: Appendix A. Single Point Formulation of Milling Equation 106 2 T t cs cd cd es ed 1 2 TT wp k ( ) 2 ( ) ( ) ( ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )) n nt t t t t t t t t t t dz ζω ω U f f f f f U . (A.1) The resultant forces contain loops that add the each flute’s and axial slice’s contribution to the force implied by 1,...,j N and 1,...,r q , respectively. The lumped static cutting force component is: 1,x 1 cs 1,y 1 1 1,z 1 ( ( , )) ( , ) ( ) ( ( , )) ( , ) ( ( , )) ( , ) q j r qN j j r q j r g t rdz a t rdz t g t rdz a t rdz g t rdz a t rdz f . (A.2) Lumped dynamic cutting force component due to the current vibrations is: cd 1 2,xx 2,xy 2,xz 1 1 1 2,yx 2,yy 2,yz 1 1 1 2,zx ( ) ( ( , )) ( , ) ( ( , )) ( , ) ( ( , )) ( , ) ( ( , )) ( , ) ( ( , )) ( , ) ( ( , )) ( , ) ( ( , )) q q q j j j r r r q q q j j j r r r j t g t rdz a t rdz g t rdz a t rdz g t rdz a t rdz g t rdz a t rdz g t rdz a t rdz g t rdz a t rdz g t rdz a f t wp k 1 2,zy 2,zz 1 1 1 . ( , ) ( ( , )) ( , ) ( ( , )) ( , ) N j q q q j j r r r t rdz g t rdz a t rdz g t rdz a t rdz U U Lumped dynamic cutting force component due to the previous vibrations is: cd 2 cd cd 2 T 2, 1 cd 2, , 1 ( ) ( ) uniform helix and pitch ( ) ( ) ( ) ( ) uniform helix, non-uniform pitch ( ) ( ) non-uniform helix and pitch N j j j N j j avg j t t T t t t t T t t T f f f f , (A.3) Appendix A. Single Point Formulation of Milling Equation 107 where the average delay due to non-uniform helix for the flute j is denoted by ,j avgT . The dynamic cutting force component for flute j is: cd 2, 2,xx 2,xy 2,xz 1 1 1 2,yx 2,yy 2,yz 1 1 1 2, ( ) ( ( , )) ( , ) ( ( , )) ( , ) ( ( , )) ( , ) ( ( , )) ( , ) ( ( , )) ( , ) ( ( , )) ( , ) ( ( , )) j q q q j j j r r r q q q j j j r r r j t g t rdz a t rdz g t rdz a t rdz g t rdz a t rdz g t rdz a t rdz g t rdz a t rdz g t rdz a t rdz g t rdz a f t wp zx 2,zy 2,zz 1 1 1 ˆ ˆ ( , ) ( ( , )) ( , ) ( ( , )) ( , ) q q q j j r r r t rdz g t rdz a t rdz g t rdz a t rdz U U , (A.4) and cd cd 2 2, 1 ( ) ( ) N j j t t f f . The lumped static edge force component is: 1,x 1 es 1,y 1 1 1,z 1 ( ( , )) ( , ) ( ) ( ( , )) ( , ) ( ( , )) ( , ) q j r qN j j r q j r g t rdz b t rdz t g t rdz b t rdz g t rdz b t rdz f (A.5) The lumped process damping force component is: ed 2,xx 2,xy 2,xz 1 1 1 2,yx 2,yy 2,yz 1 1 1 2,zx ( ) ( ( , )) ( , ) ( ( , )) ( , ) ( ( , )) ( , ) ( ( , )) ( , ) ( ( , )) ( , ) ( ( , )) ( , ) ( ( , )) ( q q q j j j r r r q q q j j j r r r j t g t rdz b t rdz g t rdz b t rdz g t rdz b t rdz g t rdz b t rdz g t rdz b t rdz g t rdz b t rdz g t rdz b f 1 2,zy 2,zz 1 1 1 , ) ( ( , )) ( , ) ( ( , )) ( , ) N j q q q j j r r r t rdz g t rdz b t rdz g t rdz b t rdz . (A.6) Appendix A. Single Point Formulation of Milling Equation 108 Only one point is measured via tap tests therefore mode shape matrix’s size needs to be reduced. The modified mode shape matrix for measured point k is: t wp t wp t wp t wp xk,1,t xk, ,t xk,1,wp xk, ,wp t wp yk,1,t yk, ,t yk,1,wp yk, ,wpk zk,1,t zk, ,t zk,1,wp zk, ,wp 3 ( ) m m m m m m m m u u u u u u u u u u u u U U . (A.7) The single point can be directly applied to semi discrete time domain and stabili- ty solutions. For surface location error, formulating the milling equation in spatial domain is faster than the modal domain method. The Fourier series expansion of the force is explained in chapter 3, Fourier coef- ficients are given in appendix B for analytical solution. Due to helix angle, the integral for each axial slice should be calculated for Fourier coefficients, which is computational- ly difficult as discussed in chapter 4. For fast algorithm, the following integral can be calculated numerically for Fourier coefficients: p s s -i s cs es xyz, s 0 1 ( ( ) ( ))e h h d F f f . (A.8) Frequency domain solution for displacement in spatial coordinates is obtained using single point approach as expressed: s xyz(i ) (i ) (i ) G F Q , (A.9) where the single point transfer function is defined: xx xy xz yx yy yz zx zy zz (i ) (i ) (i ) (i ) (i ) (i ) (i ) (i ) (i ) (i ) G G G G G G G G G G . (A.10) Correlating between frequency domain solution and Fourier series expansion as dis- cussed in chapter 4, the steady state response of the system is obtained: Appendix A. Single Point Formulation of Milling Equation 109 p p, i( (< ( ))s xyz, p,( ) ( ) e h h h t i h h h t i G Q F G . (A.11) Appendix B. Fourier Coefficients for Time Dependent Force Terms 110 Appendix B: Fourier Coefficients for Time Dependent Force Terms In chapter 4, elements in Fourier coefficient vector were defined as: ex st cs -i , , 1, 1 ( ) ( , )e 2 h h j i iF z a z d , (B.1) ex st es -i , , 1, 1 ( ) ( , )e 2 h h j i iF z b z d . (B.2) The Fourier coefficients for the static cutting force are: ex st ex st s cs 0, ,x rc tc ac cs 0, ,y rc tc ac cs 0, ,z rc ac ( ) sin ( ) sin(2 ) 2 cos(2 ) cos ( ) sin(2 ) 2 4 ( ) sin ( ) cos(2 ) sin(2 ) 2 cos ( ) cos(2 ) 4 ( ) ( cos( )cos ( ) cos( )sin ( )) j j j c F z K z K K z c F z K z K K z F z c K z K z ex t (B.3) If 1Nh ex st -i2 cs rc ac 1, ,z (e 2 i)( cos ( ) sin ( )) ( ) 4 j c K z K z F z (B.4) If 2Nh ex st ex st -i2 -i4 -i4 rc tccs , ,x -i2 -i4 ac -i4 -i2 -i4 rc tccs , ,y -i4 ac sin ( )(4 4e i+e i) (e 4 i) ( ) 16 cos ( )(4 4e i+e i) sin ( )(e 4 i) (4 4e i+e i) ( ) 16 cos ( )(e 4 i) h j h j K z Kc F z K z K z Kc F z K z (B.5) If 2,3,4,....Nh Appendix B. Fourier Coefficients for Time Dependent Force Terms 111 ex st 2 2 2 2 rc-i cs , ,x tc3 3 2 2 2 2 ac -i cs , ,y sin ( )(2 sin(2 ) cos(2 )i 4i+ i) ce ( ) cos ( )(2cos(2 ) sin(2 )i) 2 8 cos ( )(2 sin(2 ) cos(2 )i 4i+ i) ce ( ) N h h j N h j K z Nh N h N h F z K Nh z Nh N h Nh K z Nh N h N h F z ex st rc 2 2 2 2 tc3 3 ac -i cs , ,z rc2 2 sin ( )(2cos(2 ) sin(2 )i) cos ( )( 2 sin(2 ) cos(2 )i 4i i) 2 8 cos ( )(2cos(2 ) sin(2 )i) ce ( ) [ cos ( )(cos( ) 1 h N h h j K Nh z Nh K z Nh N h N h N h Nh K Nh z Nh F z K z Nh N h ex st acsin( )i) sin ( )(cos( ) sin( )i)]K z Nh (B.6) The Fourier coefficients for the static edge force are: ex st ex st ex st es 0, ,x re te ae es 0, ,y re te ae es 0, ,z re ae ( ) [ cos( ) (sin( ) / sin ( )) (cos( )cot ( ))] ( ) [ sin( ) (cos( ) / sin ( )) (sin( )cot ( ))] ( ) ( cot ( ) ) j j j F z K K z K z F z K K z K z F z K z K (B.7) If 1Nh ex st -i2 -i2 re te es 1, ,x -i2 ae -i2 -i2 re te es 1, ,y a e e i i 4 2 2sin ( ) 4sin ( ) ( ) e cot ( ) cot ( ) i 4 2 e e i i 2 4 4sin ( ) 2sin ( ) ( ) j j K K z z F z z z K K K z z F z K ex st ex st -i2 e es -i -i 1, ,z re ae cot ( ) e cot ( ) i 2 4 ( ) ( e cot ( )i e i)j z z F z K z K (B.8) 2,3,4,....Nh Appendix B. Fourier Coefficients for Time Dependent Force Terms 112 ex st -i re tees , ,x 2 2 ae -i re tees , ,y 2 2 sin( ) cos( )i cos( ) sin( )ie sin ( )( ) 1 (cos( ) sin( )i) cot ( ) cos( ) sin( )i sin( ) cos( )ie ( ) sin ( 1 N h h j N h h j Nh K Nh K zF z N h K Nh z Nh K Nh K F z N h ex st ex st ae -i es , ,z re ae ) ((sin( ) cos( )i) cot ( )) e ( ) ( cot ( )i i) N h h j z K Nh z F z K z K Nh (B.9) Appendix C. Variable Helix Cutter Representation 113 Appendix C: Variable Helix Cutter Representation Changing the helix angle for each flute, and/or changing the helix angle axially are types of non-uniform helix cutters. Non-uniform helix introduces different pitch angle at each axial segment of the cutter. In order to get the actual angle between the flutes, pitch angle for each axial segment is formulated: p,eff , p, 1( ) ( ) ( )j j j jz z z (C.1) The lag angle for corresponding to axial elevation is obtained when the equation for the helix curve is known. For instance, if the helix change is approximated linearly along the axial position, then following formula is derived for the cylindrical end mill: ,st ,e c ,st ,e ctan( ( 1) / ) tan( / ) ( ) ( ) j j j j j i k dz i L i kdz i L kz dz C D z (C.2) where ,e ,st,j ji i are exit and start helix angles respectively, cL is the length of cut and C is the sum of the lag angles coming from the previous axial discretizations. For different types of variable helix cutters, geometry of the helix needs to be known for the calcula- tion. Appendix C. Variable Helix Cutter Representation 114 Figure C-1 : Crest cut end mill with varying helix angle along the flute Appendix D. Accelerometer Mass Elimination 115 Appendix D: Accelerometer Mass Elimination When multiple FRFs are available, and structure’s behaviour throughout the complete geometry is significant, then global identification methods take place. Usually, either motion sensor or impact sensor get fixed to a point and the other one is moved on many points. Due to nonlinearities in the structure, and other effects like accelerometer mass, the curve fitting for all of the measured FRFs is a challenge. For instance, in Figure D-1, for torsional mode the variation of frequncy is seen throughout different measurement locations. This variation is found to be coming from accelerometer mass, where it affects structure excessively at the most flexible points of that mode. Figure D-1: Overlayed measurements indicating the torsional mode of the structure When curve fitting is done from the direct FRF location ( at the bottom left corner of the roughing plate), big errors arised at the left and right top corner’s curve fits as presented in Figure D-2, there is a huge error in the algorithm for torsional mode as frequency is shifted for that point. Appendix D. Accelerometer Mass Elimination 116 Figure D-2: Curve fit for the flexible top right corner, the pink line is the fit curve and the blue line is the measurement. The discrepancy at torsional mode at 1720 Hz is apparent. In order to overcome this problem, accelerometer mass elimination algorithm proposed by Cakar and Sanliturk [62] is applied on each curve fit. This formulation is more general, containing both direct and cross FRF mass elimination. The corrected direct FRF is expressed as: * 2 acc G G 1 m ω G ii ii ii ’ (D.1) where *Gii is the corrected FRF, Gii is the measured FRF, macc is the mass of the accelerometer and ω is the angular frequency. The corrected cross FRF is defined as: * 2 acc G G 1 m ω G ji ji jj , (D.2) where i is the accelerometer position and j is the hitting location. G jj in the denominator makes direct FRF measurement necessary for the mass cancellation. Since hitting from plate tip is not possible practically due to multiple hits, the authors Appendix D. Accelerometer Mass Elimination 117 recommend use of dummy mass during the measurement. However, this may even contaminate the measurement. Assuming that the system is linear, accelerometer mass only affects the diagonal terms in the global mass matrix and the conducted measurements are precise, following procedure can be applied. It is known that the natural frequency is ratio of effective mass and stiffness of the corresponding mode (similar to modal mass and stiffness) as: eff,2 n eff, k ω m j j , (D.3) where effective mass is the inverse of the square of the direct point FRF’s mass normalized mode shape: eff, 2 1 m =j jju . (D.4) When we add accelerometer to the structure at point j , it is assumed that only the mass is affected hence changed natural frequency is: eff,2 eff, acc k m +m j n j . (D.5) For the hitting location, since corrected direct FRF can be calculated from Eq. (D.1), the effective stiffnesses eff,k i for all modes are known for that point. Using the inverse calculation the modified (measured system’s) mode shape matrix can be calculated using: eff, 2 n k ω i iu . (D.6) The modified mode shape of the point where accelerometer is mounted, namely ju , can be obtained when the curve is fitted and residue terms are found (real part of the residue terms is equal to j iu u ). Using the cross measurement, direct measurement at that point is defined as: Appendix D. Accelerometer Mass Elimination 118 2 2 1 n, n, G G 2ζ ω ω m j j j i jj ji p p p p u u u u s s . (D.7) Then all the FRFs can be corrected using the measured accelerometer mass accm 1.2784g . The frequency variation and the new fit for the torsional mode are illustrated in Figure D-3 and Figure D-4, respectively. Figure D-3: The most of the frequencies are 1719 Hz and variation is around 2 Hz. Appendix D. Accelerometer Mass Elimination 119 Figure D-4: The curve is fitted from the direct FRF and since there is no shift, the errors are less at the top corners of the plate.