VIRTUAL THREE-AXIS MILLING PROCESS SIMULATION AND OPTIMIZATION by SÜKRÜ DORUK MERDOL BSc, Middle East Technical University, Ankara, Turkey, 2000 MASc, University of British Columbia, Vancouver, Canada, 2003 A THESIS SUBMITTED IN PARTIAL FULLFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April 2008 © Sükrü Doruk Merdol, 2008 Abstract The ultimate goal in the manufacturing of a part is to achieve an economic production plan with precision and accuracy in the first attempt at machining. A physics-based comprehensive modeling of the machining processes is a fundamental requirement in identifying optimal cutting conditions which result in high productivity rates without violating accuracy throughout the part production process. This thesis presents generalized virtual simulation and optimization strategies to predict and optimize performance of milling processes up to 3-axis. Computationally efficient mathematical models are introduced to predict milling process state variables such as chip load, force, torque, and cutting edge engagement at discrete cutter locations. Process states are expressed explicitly as a function of helical cutting edge - part engagement, cutting coefficient and feedrate. Cutters with arbitrary geometries are modeled parametrically, and the intersection of helical cutting edges with workpiece features are evaluated either analytically or numerically depending on geometric complexity. The dynamics of generalized milling operations are modeled and the stability of the process is predicted using both time and frequency domain based models. These algorithms enable rapid simulation of milling operations in a virtual environment as the part features vary. In an effort to reduce machining time, a constraint-based optimization scheme is proposed to maximize the material removal rate by optimally selecting the depth of cut, width of cut, spindle speed and feedrate. A variety of user defined constraints such as maximum tool deflection, torque/power demand, and chatter stability are taken into consideration. Two alternative optimization strategies are presented: pre-process optimization provides allowable depth and width of cut ii during part programming at the computer aided manufacturing stage using chatter constraint, whereas the post-process optimization tunes only feedrate and spindle speed of an existing part program to maximize productivity without violating physical constraints of the process. Optimized feedrates are filtered by considering machine tool axes limitations and the algorithms are tested in machining various industrial parts. The thesis contributed to the development of a novel 3-axis Virtual Milling System that has been deployed to the manufacturing industry. iii Table of Contents Abstract.......................................................................................................................................... ii Table of Contents ......................................................................................................................... iv List of Tables................................................................................................................................ vii List of Figures............................................................................................................................. viii Nomenclature ............................................................................................................................. xiii Acknowledgements ................................................................................................................... xvii 1. Introduction................................................................................................................................1 2. Literature Review ......................................................................................................................5 2.1. Overview..............................................................................................................................5 2.2. Kinematics of Milling ..........................................................................................................5 2.3. Modelling of Milling Forces ................................................................................................8 2.4. Process Simulation.............................................................................................................12 2.5. Process Optimization .........................................................................................................14 2.6. Dynamics of Milling and Chatter Stability ........................................................................17 3. Process Simulation ...................................................................................................................25 3.1. Introduction........................................................................................................................25 3.2. Kinematics of Milling ........................................................................................................25 3.3. Geometric Relations for Milling Cutters ...........................................................................27 3.3.1. Lag Angle Representation .......................................................................................29 3.3.2. Radius and Axial Immersion Angle Representations..............................................31 3.4. Mechanics of Milling - Generalized Kinematics ...............................................................33 3.5. Cutting Force Coefficient Identification ............................................................................38 iv 3.5.1. Orthogonal to Oblique Cutting Transformation ......................................................39 3.5.2. Mechanistic Model ..................................................................................................41 3.6. Cutter Engagement Feature (CEF) - Cutter Engagement Plane (CEP) .............................44 3.6.1. Determination of Integration Boundaries................................................................46 3.6.2. Performing Integration ............................................................................................52 3.6.3. Special Case: Flute Re-cut ......................................................................................56 3.7. Verification.........................................................................................................................58 4. Dynamics of Milling Operations.............................................................................................78 4.1. Introduction........................................................................................................................78 4.2. Dynamics of Milling ..........................................................................................................79 4.3. Stability of Milling.............................................................................................................86 4.3.1. Time Domain Based Stability Solution ...................................................................92 4.3.2. Frequency Domain Based Stability Solutions.........................................................97 4.3.2.1. Zero Order (ZO) Solution ........................................................................108 4.3.2.2. Generalized Zero Order (GZO) Solution .................................................109 4.3.2.3. Multi Frequency (MF) Solution ...............................................................113 4.3.2.4. Iterative Stability Solution and Eigenvalue re-analysis ...........................114 4.4. Comparison of Stability Methods ....................................................................................120 4.5. Summary ..........................................................................................................................127 5. Process Optimization .............................................................................................................129 5.1. Introduction......................................................................................................................129 5.2. Process Optimization .......................................................................................................130 5.3. Pre-Process Optimization.................................................................................................133 5.3.1. Stability Constraint................................................................................................133 5.3.2. Machine Tool Torque/Power Constraints ..............................................................137 5.3.3. Chip Thinning Constraint......................................................................................144 v 5.4. Verification for Pre-Process Optimization .......................................................................146 5.4.1. Case Study I...........................................................................................................146 5.4.2. Case Study II .........................................................................................................151 5.5. Post-Process Optimization ...............................................................................................158 5.5.1. Feedrate Optimization ...........................................................................................163 5.5.2. Feedrate and Spindle Speed Optimization ............................................................167 5.6. Filtering "Unrealizable" Feed Commands .......................................................................169 5.6.1. Classification of Feed Blocks................................................................................170 5.6.2. Filtering Rules .......................................................................................................173 5.7. Updating the Original CL File .........................................................................................181 5.8. Verification for Post-Process Optimization......................................................................183 6. Conclusions ............................................................................................................................190 6.1. Conclusions ......................................................................................................................190 6.2. Future Research Directions ..............................................................................................193 Bibliography ...............................................................................................................................195 vi List of Tables 3. Process Simulation Table 3.1 : Cutting force coefficients of Al 6061 ..........................................................................61 Table 3.2 : List of milling cutters used in experiments..................................................................67 Table 3.3 : List of Operations ........................................................................................................67 Table 3.4 : Average cutting force coefficients for Uddeholm Carmo - Carmo Cast......................68 4. Dynamics of Milling Operations Table 4.1 : Summary of different stability solution methods.......................................................115 Table 4.2 : Dynamic parameters [101].........................................................................................121 5. Process Optimization Table 5.1 : Integration boundaries for different regions ..............................................................140 Table 5.2 : Generalized coefficients for tangential force .............................................................141 Table 5.3 : Cutting force coefficients of Al 6061 ........................................................................153 Table 5.4 : All possible update configurations with proposed update rules ................................176 Table 5.5 : Assumed Characteristics of the Heidenhain Controller.............................................184 Table 5.6 : Other parameters of the machine tool used for verification ......................................184 Table 5.7 : Optimization Criteria for Different Operations and Cutting Tools............................185 Table 5.8 : Cycle time results of operations.................................................................................186 vii List of Figures 2. Literature Review Figure 2.1 : Methods of milling operations: (a) up milling; (b) down milling ................................6 Figure 2.2 : 3-axis milling examples (source: EdgeCAM [11]) ......................................................7 Figure 2.3 : Exact kinematics of a milling tool for calculation of the chip thickness ...................12 Figure 2.4 : A comprehensive process simulation and optimization flow chart by Altintas [26] .17 Figure 2.5 : Wave generation in turning ........................................................................................19 3. Process Simulation Figure 3.1 : Typical horizontal milling machine, different milling operations..............................26 Figure 3.2 : Coordinate systems in a milling operation .................................................................28 Figure 3.3 : (a) Generalized tool geometry [20], (b) angle convention .........................................32 Figure 3.4 : Mechanics and kinematics of 3-axis milling..............................................................34 Figure 3.5 : Orthogonal and oblique cutting geometries [64]........................................................39 Figure 3.6 : Cutter - workpiece intersection, CEF mapping .........................................................45 Figure 3.7 : Generalized CEF ........................................................................................................45 Figure 3.8 : Rectangular CEF; (a) possible intersection cases, (b) parameter definition ..............46 Figure 3.9 : Calculation of integration boundaries for general end mills ......................................47 Figure 3.10 : (a) Planar machining with a ball end mill, (b) example CEF for ball end milling...48 Figure 3.11 : Change of integration boundaries over time for the most general case ...................50 Figure 3.12 : Root bracketing ........................................................................................................51 Figure 3.13 : Variation of edge integrand in y direction for a 12 mm ball end mill ......................56 Figure 3.14 : (a) Low helix cutter (4-fluted) with small depth of cut, (b) high helix cutter (2fluted) with large depth of cut, (Inset) chip geometry..............................................56 Figure 3.15 : Cutting region at the front of the cutter....................................................................57 Figure 3.16 : Flute re-cut example, same depth, different helix angles.........................................59 viii Figure 3.17 : Multiple immersion milling, engagement map, simulated vs measured forces .......60 Figure 3.18 : (a) Face milling, (b) right side profiling, (c) top hole helical milling, (d) middle hole helical milling...........................................................................................................62 Figure 3.19 : Face milling a) measured, (b) simulated cutting forces ...........................................63 Figure 3.20 : Profile milling a) measured, (b) simulated cutting forces ........................................64 Figure 3.21 : Top hole helical milling a) measured, (b) simulated cutting forces .........................65 Figure 3.22 : Test part stock and final shape (a) isometric; (b) left side; (c) front; (d) top views .66 Figure 3.23 : Cutting force coefficient variation ...........................................................................70 Figure 3.24 : Actual versus commanded feedrate..........................................................................70 Figure 3.25 : T63 - Operation 3, finishing of the walls at z = 105.26 [mm] .................................72 Figure 3.26 : T63 - Operation 3, (a) engagement map at one point, (b) measured forces for two spindle revolutions at and around the engagement map in (a) at z = 105.26 [mm] .73 Figure 3.27 : T63-Operation 1, roughing of the top part at z=163.5 [mm] ...................................74 Figure 3.28 : T63 - Operation 1, roughing at z = 88.51 [mm] .......................................................75 Figure 3.29 : T20 - Operation 2, semi-finishing of the top............................................................76 Figure 3.30 : T20 - Operation 2, semi - finishing of door handle. area .........................................77 4. Dynamics of Milling Operations Figure 4.1 : Main sources of flexibilities in milling ......................................................................79 Figure 4.2 : Undulations on the chip, and dynamic chip thickness ...............................................81 Figure 4.3 : Dynamic chip thickness definition for general end mills...........................................82 Figure 4.4 : Multi degree-of-freedom milling dynamics ...............................................................86 Figure 4.5 : Approximation of the delayed term using the semi-discretization (SD) by ref. [73].94 Figure 4.6 : a) General end mill with various sections; b) Solution of stability for the bottom most part; c) Switching to second relative depth of cut for stability solution ................ 111 Figure 4.7 : Graphical representation of iterative root search .....................................................113 Figure 4.8 : Shift of eigenvalues when the frequency is changed ...............................................120 Figure 4.9 : Stability lobes obtained using various stability methods .........................................124 Figure 4.10 : Chatter frequencies obtained from frequency domain based solutions..................124 ix Figure 4.11 : Normalized X-direction eigenvector magnitudes of MF-5 ....................................125 Figure 4.12 : Eigenvalues of the transition matrix of the SD method .........................................126 Figure 4.13 : Comparison of exact and approximate root functions of one harmonic case ........127 5. Process Optimization Figure 5.1 : Manufacturing flow chart, current and proposed methodologies.............................131 Figure 5.2 : Design variables .......................................................................................................132 Figure 5.3 : Optimized NC code generation flow chart with the Pre-Process Optimization.......134 Figure 5.4 : Design variables that are important to suppress chatter vibrations ..........................135 Figure 5.5 : Chatter stability lobes for constant (a) radial width of cut, (b) axial depth of cut ...136 Figure 5.6 : Design space bounded by chatter stability ...............................................................136 Figure 5.7 : Important parameters for analytical solution ...........................................................138 Figure 5.8 : Cutter workpiece engagement, integration boundaries ............................................139 Figure 5.9 : (a) Case I, (b) Case II ...............................................................................................140 Figure 5.10 : Possible intersection scenarios ...............................................................................142 Figure 5.11 : Typical torque-power characteristics of a machine tool.........................................143 Figure 5.12 : Chip thickness for a face mill with a round insert..................................................144 Figure 5.13 : (a) Integration boundaries with various depths of cut, simulated tangential forces, (b) analytically obtained maximum tangential forces..................................................147 Figure 5.14 : Torque - Power curves of the machine tool............................................................148 Figure 5.15 : Radial and axial depths of cut curves for up-milling, feasible region....................149 Figure 5.16 : Stability lobes for full immersion milling case ......................................................150 Figure 5.17 : MRR curves of up-milling based on stability, torque and power limits.................150 Figure 5.18 : Radial and axial depths of cut curves for up-milling, feasible region....................150 Figure 5.19 : Gear box cover: (a) complete solid model (double sided machining), (b) solid model for machining front features, (c) real cut part with front features..........................151 Figure 5.20 : Transfer functions of 10 [mm] endmill in x and y directions (Shrink Fit).............152 Figure 5.21 : Transfer functions of 25 [mm] indexable cutter in x and y directions (Corogrip) .153 x Figure 5.22 : Stability lobes for full immersion milling (b=1) of (a) 10 [mm] end mill, (b) 25 [mm] indexable cutter. Cutting coefficients from Table 5.3 are used ..............................154 Figure 5.23 : Stability lobes for various immersion cases, 10 [mm] end mill. Cutting coefficients from Table 5.3 are used ..........................................................................................155 Figure 5.24 : MRR characteristics of 10 [mm] end mill based only on stability in Figure 5.23 with feed variation in Figure 5.26.d where =0.1 [mm/tooth.rev], and N=2 ...................155 Figure 5.25 : Mori Seiki SH403 - Torque/Power characteristics.................................................156 Figure 5.26 : Process planning charts for 10 [mm] end mill, n=20000 [rpm], N=2....................157 Figure 5.27 : Process planning charts for 25 [mm] indexable cutter, n=17000 [rpm], N=2........157 Figure 5.28 : Optimized NC code generation flow chart with Post-Process Optimization .........158 Figure 5.29 : Process optimization flow chart .............................................................................159 Figure 5.30 : Sample force variation for a uniform pitch cutter ..................................................160 Figure 5.31 : Two possible quadratic functions ...........................................................................161 Figure 5.32 : Inverse parabolic interpolation...............................................................................162 Figure 5.33 : CEFs and static tool deflection at surface generation point ...................................166 Figure 5.34 : Sample stability lobes.............................................................................................168 Figure 5.35 : Graphical representation of design space...............................................................170 Figure 5.36 : An example of optimized feedrate with sharp variations.......................................171 Figure 5.37 : Types of feed command used for filtering..............................................................172 Figure 5.38 : Allowable combination of different feed block types ............................................172 Figure 5.39 : Parameter definitions for a single feed block i.......................................................173 Figure 5.40 : Constant acceleration kinematic profiles ...............................................................174 Figure 5.41 : Update algorithm flow chart ..................................................................................177 Figure 5.42 : Feed filtering example with recursive updating .....................................................178 Figure 5.43 : Feed filtering example with recursive updating & MFC .......................................180 Figure 5.44 : Filtered optimum-feed-commands for 1 [g] acc/dec..............................................180 Figure 5.45 : CL file optimization routine ...................................................................................182 Figure 5.46 : Feedrate update strategy for linear versus circular segments.................................183 Figure 5.47 : Estimated Torque and Power Curves .....................................................................183 xi Figure 5.48 : Chip formation with (a) original, (b) optimized NC programs ..............................187 Figure 5.49 : Pictures of T63 from front/bottom after execution of original/optimized CL files188 Figure 5.50 : T20F - Part of the tool path from Operation 5, Semi-Finishing.............................188 Figure 5.51 : Pictures of T20F inserts after the (a) original, (b) optimized machining ...............189 xii Nomenclature a axial depth of cut a lim limiting axial depth of cut for chatter stability A closed form coefficients A(t) directional coefficient matrix b normalized width of cut B radial width of cut c feed per tooth det( . ) determinant operator dS length of the differential element defined along the cutting cutting edge dz height of the differential element defined along the cutter axis D tool diameter E delay term; youngs modulus f feed vector in the direction of resultant tool motion fc commanded feedrate fe end feedrate fs start feedrate fu unit feed vector in the direction of resultant tool motion F general process output (cutting force, torque, power,...) G overall transfer function matrix h uncut chip thickness xiii hd uncut dynamic chip thickness hr maximum number of harmonics i, j, k unit vectors in the principal directions x, y and z i0 helix angle I area moment of inertia I identity matrix j flute number K tc, K rc, K ac tangential, radial, and axial milling-cutting force coefficients K te, K re, K ae tangential, radial, and axial milling-edge force coefficients m modal mass n spindle speed nj unit vector perpendicular to the cutter body N number of teeth on the cutter p eigenvector P mode shape matrix R local radius of the cutter; nose radius x, y, z coordinate axes defining the cutter coordinate system X, Y, Z coordinate axes defining the machine tool coordinate system z lim, 0, z lim, 1 lower and upper axial integration limits of process outputs βa friction angle αn normal rake angle Δ vibration vector xiv δ tool deflection reflected on the surface ε phase shift between current and previous vibration marks ζ damping ratio η chip flow angle κ axial immersion angle λ eigenvalue λI imaginary part of the eigenvalue λR real part of the eigenvalue τ tooth period, time delay τs shear stress Φ oriented transfer function matrix φ rotation angle of the cutter φc shear angle φp pitch angle of the cutting tool φ st, φ ex radial start and exit angles in milling ϕ phase shift of the complex eigenvalue ψ lag angle ωc chatter frequency ωn natural frequency ωT tooth passing frequency T . ACIS transpose operator Andy Charles Ian’s System (solid modeling kernel) xv APT Automatically Programmed Tools B-rep Boundary representation CAM Computer Aided Manufacturing CEF Cutter Engagement Feature CL Cutter Location DDE Delayed Differential Equation FRF Frequency Response Function GZO Generalized Zero-Order HSK63A Hollow, Short, kegel - (the German word for Taper) and 63A is the size ICP Integrated Circuit Piezoelectric MATLAB MATrix LABoratory (software program) MFC Maximum Feed Correction MRR Material Removal Rate NC Numerical Control SD Semi-Discretization TD Time-Domain ZO Zero-Order xvi Acknowledgements It has been a long journey with the Manufacturing Automation Laboratory since I started as a graduate student. Every day, every hour, every moment of the past seven and a half years, I have truly enjoyed the ever-growing responsibility of working for a world-class academic person, Dr. Yusuf Altintas. Not only his academic excellence, but also his admirable stance on integrity have deeply influenced and guided me to successfully craft the upcoming stage of my professional life. It has been a great pleasure working, travelling and having discussions with him, and I would like to express my heartfelt gratitude to him for being generous with his time and knowledge in each step of my academic life. As much as I am thankful to the intellectual help, I am also grateful to have had the spiritual support of my family, friends and colleagues. I have better understood the importance of my family and their roles in gaining self-enlightenment, which tremendously helped me to overcome difficulties in my life and to make it an enjoyable journey. My dear mother, Mucella, my father, Erdal, and my sister, Ceylan, I am grateful to have you all. I also wish to express thanks to my close friends Kim Duffell, Mike Stewart, and Burak Sencer for sharing their hearts with me at different stages of my studies. We have spent a great deal of time together and connected with each other at various levels of our personal and social lives. Finally, I would like to dedicate this work to two beautiful members of the future generation, my beloved nieces: Lara and Mira. xvii Chapter 1 Introduction Milling operations are most widely used to more accurately produce net shapes in the manufacturing industry. Applications can be found in the milling of dies and molds, jet engine parts made of heat resistant alloys, aircraft fuselage and wing panels, and biomedical parts. The machining of small parts may technically only take a few minutes in overall mass production applications, while dies, molds and aerospace parts in fact may take several days of milling when accuracy and productivity rates become key determinants in the economic survival of the manufacturing industry. Both the accuracy and surface quality of machined parts are influenced by the relative vibrations between the machine tool and part structures, the selection of proper tool geometry and grade, the positioning accuracy, and the thermal stability of the computer controlled machine tools. Productivity within machining operations is measured by the volumetric material removal per unit time, which is desired to be as high as possible within the physical limitations of the machine and the processes. There is no mathematical model of integrated machining physics used in the industry. Process planners select cutting conditions, i.e. depths of cut, spindle speeds, and feedrates based on their accumulated experience, which is gained through trial and error over a period of years. When a part becomes considerably more costly, such as with an aircraft wing panel or a large stamping die, planners select the most conservative material removal rates in order to avoid scrapping of the 1 Chapter 1. Introduction parts caused by tool failure, vibrations, and overloading of the machine tool spindle. Nevertheless, it is a preferred approach to test the machining operations in a virtual environment and achieve an accurate and economically sound part that has been machined during the first trial on the factory floor. This scenario is ultimately possible if the physics behind the machining process is modeled and integrated to computer aided manufacturing systems where parts are graphically and virtually machined ahead of real-time production. While the geometric simulation of the part machining is already available, the physics behind machining process has not been developed with sufficient accuracy and computational efficiency to achieve a true Virtual Milling System. This thesis presents a novel 3-axis Virtual Milling System with integrated mechanics, dynamics and optimized process planning algorithms. The mechanics and dynamics of milling for arbitrary cutters and workpiece engagement conditions are modeled. First, the chatter vibration free spindle speeds, depth and width of cut are predicted with torque and power constraints of the machine. Closed form expressions are developed to predict cutting forces, torque, and power when the cutter-part engagement conditions can be mathematically defined. Cutters and engagements with arbitrary shapes are handled by the new, computationally efficient numerical methods. As the cutter travels along varying part features, the process states such as maximum forces, torque, power, chip load, tool deflection and vibrations are predicted. The spindle speed and feedrates are automatically adjusted to achieve highest material removal rate while respecting the physical limits of the machine and cutting tool. The developed algorithms are integrated to a prototype Virtual Machining System, and its effectiveness is demonstrated in machining a sample gear box cover and a stamping die for a passenger car door panel. 2 Chapter 1. Introduction The thesis is organized as follows: A review of related literature is presented in Chapter 2. The generalized mathematical model of parametric tool geometry and process mechanics are presented in Chapter 3. The generalized model uses parametric cutting tool geometry to predict resulting chip load force distribution along the cutting edge using closed form formulations as a function of arbitrary cutter-part intersection geometry. Process outputs such as cutting forces, bending moment, tool deflection, spindle torque and power consumptions are predicted using mechanics of 3-axis milling operations. Computationally efficient, closed-form formulations are derived for cutters with analytically defined cutting edges. When the cutter and part features have arbitrary geometries, computationally efficient numerical algorithms are employed. Chapter 4 is dedicated to generalized dynamics of 3-axis milling operations. The chatter stability of the milling operations are solved both in frequency and time domains. A detailed comparison of various stability models is also presented. Moreover, the analytical zero-order solution available in literature is generalized to model stability of end mills with complex geometries. In Chapter 5, virtual optimization of milling operations based on physics of metal cutting is discussed. First, the optimal spindle speed, depth and width of cut are identified by considering the chatter stability, torque and power limits of the machine to assist process planners during NC programming. Another optimization scheme that is proposed is based on the mechanics of cutting that varies along the tool path as the part geometry changes dynamically. This process is modeled at each tool location along the NC tool path; maximum values of chip, force, deflection, torque, and power are evaluated; chatter stability is assessed, and cutting parameters (feedrate and spindle speed) are automatically adjusted to avoid violating any of the constraints. Lastly, a filtering algo- 3 Chapter 1. Introduction rithm is developed to eliminate sudden changes in optimized feedrates that cannot be achieved due to both limited acceleration and deceleration of the motors. The thesis is concluded with a short summary of the main results and contributions. 4 Chapter 2 Literature Review 2.1. Overview Contrary to advancement in machine tool technology and milling tool development, cutting parameters such as spindle speed, feedrate, depth and width of cut are selected conservatively to avoid the risk of damaging costly workpieces and machine tools during manufacturing. There is an increasing demand for virtual process simulators that are capable of predicting performance measures such as cutting forces, tool deflection, workpiece/tool vibrations, spindle torque and power demands. Complex parts with many geometric features lead to dynamically changing engagement conditions during machining that require automatic selection of cutting parameters along the tool path to fully exploit the capacities of machine tool and milling cutter. Hence, virtual process simulation and optimization are essential components to maximize material removal rates while maintaining accuracy of the part in milling operations. The areas of kinematics of milling, dynamics of milling and chatter stability, milling process simulation and cycle time optimization are the foundations to achieve virtual machining studied in the thesis, and the related past research reported in the literature is reviewed in this chapter. 2.2. Kinematics of Milling Milling is an intermittent multi-point cutting operation, in which each tooth is in contact with the workpiece for a fraction of a spindle period. Tooth contact time varies as a function of width of cut, spindle speed and number of flutes. Unlike turning where chip thickness is constant, milling 5 Chapter 2. Literature Review tool follows a trochoidal path due to simultaneous feed and rotation motions leading to a variable chip thickness. Two different methods shown in Figure 2.1. are observed in planar milling: up (conventional) and down (climb) milling. In up milling, chip thickness increases as the cutter rotates and the direction of the cutter rotation opposes the feed motion. In down milling, on the other hand, the initial chip thickness, where tool enters into the workpiece, is at its maximum and decreases until it leaves the workpiece at zero chip thickness. The direction of the cutter rotation is same as the feed motion in down milling. The early research [1]-[7] focused on chip formation in milling and estimated power consumption of the spindle during machining. Martelotti modeled the kinematics of milling where a tooth followed a trochoidal path resulting in a complex chip model [8]. He showed that the chip thickness in planar milling can be approximated by a simpler expression when radius of the cutter is large compared to the feed: h ( φ ) = c ⋅ sin φ , (2.1) where h is the instantaneous uncut chip thickness, c is the feed per tooth and φ is the tooth immersion angle. Y φst= 0 Y Feed φ φex= π X X Feed (a) (b) Figure 2.1 : Methods of milling operations: (a) up milling; (b) down milling. 6 Chapter 2. Literature Review Sinusoidal chip thickness approximation in Eq (2.1) does not hold in 3-axis profiling and circular helical milling where the cutting tool is fed parallel to the cutter axis as well as perpendicular to it (see Figure 2.2). In order to calculate uncut chip thickness for general 3-axis milling applications, Fussell et al. [9] used the vector in the direction of travel to estimate how chip was generated in the case of a ball end mill. Projection of the travel vector onto the surface normal was used to calculate the uncut chip thickness. Lazoglu et al. [10]; on the other hand, separated conventional horizontal feed motion from downward/upward motion, and calculated uncut chip thickness contributions of each motion individually. Resultant uncut chip thickness was then obtained by combining chip thickness values from both motions. Figure 2.2 : 3-axis milling examples (source: Esprit [11]). A great selection of cutting tools are used in the industry to machine different parts. In the aerospace industry, cylindrical end mills are heavily involved in roughing and semi-finishing of aluminum ribs; whereas, helical ball and bull nose end mills are mostly preferred by the die and mold industry in machining sculptured surfaces. In the literature, milling mechanics had been developed for specific types of cutters such as face [12], cylindrical [13]-[15], ball end [16]-[18], and tapered ball end mills [19]. The most general approach was proposed by Altintas and Engin [20]. 7 Chapter 2. Literature Review In their model, standard Automatically Programmed Tools (APT) definition of a solid end mill was adopted and milling tools were successfully generalized by wrapping a helical flute around a parametrically defined cutter body. Using the exact kinematics of milling, the undeformed chip thickness distributions along flutes were accurately calculated in the presence of structural vibrations. 2.3. Modelling of Milling Forces Analytical and semi-analytical prediction of milling forces have been intensively researched in the past to avoid empirical techniques, which require an abundant amount of experimental data. In the late 1920s, Salomon [21] expressed the specific cutting pressure as an exponential function of the chip thickness based on the work done with a straight tooth cutter. Sabberwal and Koenigsberger [22][23] used similar exponential specific cutting pressures to model milling forces analytically. This model is known as mechanistic model and the instantaneous force acting at right angles to the chip area is calculated by the product of the chip area and the specific cutting pressure (K): F Cutting Force = K ⋅ a ⋅ h , (2.2) m where K = C ⋅ h is the specific cutting pressure, C and m are constants depending on the workpiece material and the milling cutter, a ⋅ h is the chip area, a is the depth of cut and h is the instantaneous uncut chip thickness. This approach has been adopted by many researchers [13][24]-[26] in the analysis of milling forces. Armarego [27][28] defined an alternative form of mechanistic modeling by separating shear deformation (cutting) from the edge effect (rubbing). Cutting coefficients were defined based on the classical oblique model [28][29][30], whereas edge forces were related to rubbing of the work 8 Chapter 2. Literature Review material on the flank face of the cutting edge. In this thesis, the linear edge force model by Armarego and Epp [27] is adopted to formulate the milling forces. The linear edge force model expresses the cutting force as a function of cutting coefficients, axial depth of cut and chip thickness as: F Cutting Force = K e ⋅ a + K c ⋅ a ⋅ h , (2.3) where K e , and K c are the edge and cutting force coefficients, respectively. In order to calculate cutting forces, cutting force coefficients need to be identified for the work material. One of the most widely studied models to determine cutting force coefficients is the mechanistic approach [13][24][31]. In this model, cutting forces are measured while keeping immersion constant and varying feedrates. Then, linear regression is applied to fit an approximate line to the average of measured cutting forces, which leads to identification of the unknown cutting force coefficients [30]. Since geometric properties of a milling cutter such as helix angle, rake angle, relief angle; workpiece material and other variables are all embedded in cutting coefficients, identified cutting coefficients become valid only for a specific cutter and a workpiece material. Another approach to determine cutting coefficients is based on an orthogonal to oblique cutting transformation. An orthogonal cutting database composed of shear stress, shear angle and friction coefficient identified from orthogonal cutting tests is transformed into three dimensional oblique milling mechanics [30]. The advantage of this method is that the coefficients obtained represent the true behavior of the material and local cutting edge geometry, and they are not limited to a specific tool geometry. Except for inserts having irregular rake face geometry with chip breaking grooves, this technique is applicable to end mills with arbitrary geometry, twist drills, boring and turning tools. 9 Chapter 2. Literature Review In an overview study, Smith and Tlusty [32] classified various models of the milling process into five groups: (1) the Average Rigid Force Static Deflection Model: This model takes material removal rate as a basis for calculation of process outputs. Cutting force, static tool deflection, torque and power are assumed to be linearly dependant on material removal rate. This model is the simplest and the least accurate as there is no simple relationship between material removal rate and process outputs. (2) the Instantaneous Rigid Force Model: Force is computed on incremental sections of the helical cutting edge and the resultant force is calculated using vectorial sum of all incremental forces. In this model, cutter deflection in response to force does not cause any change in force. Sabberwal and Koenigsberger [23] presented the first model in this group. (3) the Instantaneous Rigid Force, Static Deflection Model: This is similar to the previous model but static deflection of the tool is also considered under static loading. In this model, deflection of the cutter does not influence the uncut chip thickness, hence it does not affect the cutting forces. Form errors can be predicted at points where the helical flutes generate the finished surface. (4) the Instantaneous Force with Static Deflection Feedback Model: The deflection of the cutter is computed and its influence on the chip thickness and the force is considered by Sutherland et al. [13] and Armarego et al. [33]. (5) the Regenerative Force, Dynamic Deflection Model: This is the most accurate and complex model of the milling process. The equations of motion are solved in the time domain in order to calculate the vibration of the tool. These vibration terms are then used to calculate instantaneous dynamic chip thickness, which results in dynamic cutting forces. Term "regeneration" is used to 10 Chapter 2. Literature Review emphasize that the effect of previous time steps are taken into account. Tlusty and Ismail [34] presented the time domain simulation of helical end mills by including the structural dynamics of the system. Sutherland and DeVor [35] utilized the regeneration model and improved their static model by considering dynamic cutting force as a feedback. Montgomery and Altintas [36] and Altintas and Lee [37] contributed significantly to the prediction of chip formation using the exact kinematics of milling. In their model, surface and cutter locations were divided into small elements so that exact chip thickness was calculated by intersecting the tool and the current surface at each time step (see Figure 2.3). Surface and cutter locations were calculated using dynamics of workpiece and cutter; therefore, nonlinearities such as the tool jumping out of cut and the influence of vibrations could be easily incorporated into their general dynamic model. By integrating the process along each cutting edge in contact with the workpiece, cutting forces, vibrations, and dimensional surface finish were predicted. Smith and Tlusty [38] used time domain simulation to obtain peak-to-peak force ratios, which were later used as a criterion in identification of stability limits. Tlusty and McNeil [24] and Sutherland and DeVor [13] have also contributed to dynamic modeling of milling processes using the regeneration technique. In this thesis, two of the above models are adopted for milling simulations. The instantaneous rigid force - static deflection model (the model #3) type of approach is developed for fast and efficient process simulation of a complete part; whereas a detailed time domain simulation based on the regenerative force - dynamic deflection model (the model #5) is used to analyze a milling process with fixed cutting conditions (spindle speed, feedrate, depth and width of cut). 11 Chapter 2. Literature Review Y Exact Chip Thickness Workpiece X End Mill Current Surface (being generated) Previous Surface Combined Dynamics of Cutter, Workpiece, and Machine Tool Figure 2.3 : Exact kinematics of a milling tool for calculation of the chip thickness. 2.4. Process Simulation Mechanics of milling have been extensively studied in the past. It would be fair to state that the mechanics of milling to predict cutting forces and chip generation are well established when the cutter-part intersection boundaries are steady state and simple. Although metal cutting mechanics models still have room for improvements in modeling shear deformation and friction at the tool chip interface [39], the current challenge is to predict cutting forces, torque, power, form errors, temperature and vibrations along a tool path as part geometry, i.e. boundary condition, varies. There are two fundamental challenges in achieving a functional, machining process simulation system: identification of cutter-part intersection along a tool path at discrete feedrate intervals; and development of computationally efficient process models. There has been research effort in simulating the machining of parts in a CAM environment. Altintas and Spence [26] used constructive solid geometry to identify cutter-part intersection, and re- 12 Chapter 2. Literature Review formulated mechanistic milling process model to predict cutting forces; however, this approach was applicable only to planar operations along a straight line with simple cylindrical tools. Jerard and Fussel et al. [40][41][42] have contributed significantly in evaluating cutter-part engagement conditions with computationally efficient algorithms based on the Z-buffer method. In this method, the workpiece was broken into a set of evenly distributed discrete z direction vectors, with the length of each vector representing the depth of the workpiece. The intersection of the z vectors with the swept surfaces of the tool path envelope were used to update the workpiece. Cutting forces were then calculated using the removed material volume and other geometric information. Spence et al. [43][44] developed cutter-workpiece intersection for 2.5 D milling operations using constructive solid geometry, which was later transferred to the ACIS [45] Boundary representation (B-Rep) modeller [46]. Yip-Hoi et al. [47] presented feature based cutter-part intersection model for cylindrical end mills using ACIS solid modeler. Although more accurate cutterpart engagement surfaces can be obtained from solid modeling, computation time drastically increased due to surface-surface intersections. Weinert et al. [48][49] developed B-Rep based cutter-part intersection and mechanistic models to predict process forces in milling. Abrari and Elbestawi [50] developed closed-form formulations for ball and cylindrical end mills by projecting the chip area onto a set of reference coordinate planes; however, this method cannot be generalized if the cutter geometry gets more complex or if the kinematics of the process is not planar. In the feature-based machining area, other researchers mainly focused on the identification of geometric features on a part to identify number, type, and sequence of machining operations, and this area is not within the scope of this work. 13 Chapter 2. Literature Review In this thesis, cutter envelope - part intersection boundaries are identified by a commercially available algorithm [51]. Computationally efficient mathematical models are proposed to predict milling process state variables such as chip load, force, torque, and cutting edge engagement at discrete cutter locations. Process states are expressed explicitly as a function of helical cutting edge - part engagement, cutting coefficient, and feedrate. Cutters with arbitrary geometries are modeled parametrically, and the intersection of their helical cutting edges with workpiece features are evaluated either analytically or numerically depending on the geometric complexity. Process variables are computed efficiently for each cutting edge - part engagement feature through aggregation of oblique cuts distributed along the helical flute to predict total force, torque and power generated at each feedrate interval. In order to decrease computation time, closed form solutions are proposed to capture only the maximum and minimum values of process outputs eliminating the need for complete time history simulation. 2.5. Process Optimization It is difficult to achieve high productivity in milling as parts become more complicated with increased number of geometric features. Variation of a part geometry along the tool path results in constantly changing depth and width of cut; therefore, conservative machining parameters are preferred to prevent unforeseeable failures such as tool breakage and spindle stall. In one study [52], it was shown that the resultant force acting on a tool could increase by more than a factor of ten in a simple cornering cut. In order to increase productivity, process parameters should be assigned as the cutter-part intersection varies along the tool path. Bearing this objective in mind, many researchers have carried out volumetric analyses of the material removal rate (MRR) to schedule feedrate and to optimize the process. One of the early feedrate optimization schemes was 14 Chapter 2. Literature Review proposed by Wang [53]. He used solid modeling to calculate the swept volume, and adopted the specific energy method to calculate the spindle horse power requirement and the resultant force. Feedrate was then adjusted based on machinability constraints provided in handbooks. Since this method was based on pure volumetric analysis, it was not accurate in modeling the process. Ip et al. [54] proposed a fuzzy based MRR optimization approach to increase machining efficiency. In their model, an optimum feedrate for each cutting point in the machining of sculptured surfaces was calculated, and the cutting force was kept constant. The surface gradient and tool life were treated as fuzzy input variables. Commercial CAM packages are also based on volumetric analysis ignoring the physics of the metal cutting process. VeriCut’s Optipath [51] and MasterCam’s HiFeed [55] perform feedrate scheduling based on the amount of material removed per unit time. Assuming that cutting forces are proportional to the MRR, and changing the feedrate accordingly may result in not only tool breakage and form errors above tolerance, but also spindle stall due to excessive loading. Apart from volumetric approaches, some researchers concentrated on feedrate scheduling using physical constraints. In the early 1990s, Altintas [26] outlined a detailed flow chart of a comprehensive machining simulation and optimization scheme that considers not only the geometric factors but also mechanics and dynamics of milling, controller performance and feed drive dynamics as well as volumetric errors of a machine tool to compensate for machining errors and reduce cycle time by adjusting machining parameters (see Figure 2.4). Over the years, a large number of publications have been produced tackling problems related to individual areas; however, a systematic study of milling process simulation and optimization has not been conducted yet. Yazar et al. [56] used cutting force predictions for feedrate optimization in 3-axis milling and proposed 15 Chapter 2. Literature Review various rough milling strategies to reduce machining time and cost. Lim et al. [57] proposed cutting path adaptive feedrate strategy to increase the productivity of sculptured surface machining subjected to force and dimensional constraints. The finished surface was divided into grids, each of which represented a control point where chip thickness, force, and deflections were calculated. Finding acceptable feedrate limits at each control point, a feedrate map was constructed for NC programmer so that proper feed value and direction could be selected. Bailey et al. [58] reduced machining time by adjusting the feedrate based on maximum chip thickness and force constraints depending on the type of the operation, i.e., roughing, semi-finishing, or finishing. Lazoglu et al. [59] optimized feedrates by keeping the predicted cutting forces at a desired constant value during sculpture surface machining. They compared two different feedrate scheduling strategies, one based on the material removal rate (MRR) and the other based on constrained force. It was shown that the MRR based feedrate scheduling led to higher feedrates violating physical constraints such as cutting forces, and chip load. The geometric intersection between tool and workpiece was used to simulate cutting forces by Fussell et al. [60]; however, their model was still discrete in a sense that the cutter was divided into thin slices and the intersection of each slice was used to calculate cutting forces. Budak et al. [61] approached the optimization problem from a chatter stability point of view. They showed that the MRR variation with respect to the axial depth of cut became a bell like curve depending on the natural frequencies in orthogonal directions of the cutting, i.e., feed and normal directions. All of the proposed optimization methods presented thus far have focused on very specific milling operations with a dedicated constraint such as cutting force, chip load, stability or geometric volume removal. Hence, there is still a demand for a generalized optimization strategy that combines 16 Chapter 2. Literature Review different machining constraints including cutting forces, chip thickness, spindle torque-power, form errors on the workpiece and even stability of the system to determine the most efficient machining parameters. CAD MODEL NC Tool Path Cutter Geometry MACHINING DATABASE Spindle dynamics (FRF) Feed drive FRF and control Cutter-material cutting constants Kinematics of the machine Volumetric error model of the machine Cutter-part intersection geometry calculations FINAL PROCESS PLAN Optimized speed, feed, depth, width, error compensation Machining process simulation CAD BASED VIRTUAL MACHINING PROCESS SIMULATION SYSTEM MONITORING AND CONTROL DATA PATH PLANNER CL File improved path planning simulation results Strategy analysis M A C H I N E T O O L Peak force, torque, power, tracking error, modal frequencies Figure 2.4 : A comprehensive process simulation and optimization flow chart by Altintas [26]. 2.6. Dynamics of Milling and Chatter Stability When structural dynamics is taken into account, time-dependent and periodic cutting forces also become dependant on vibrations resulting in dynamic cutting forces. Depending on cutting conditions, dynamic cutting forces lead to chatter vibrations, which occur in milling due to dynamic characteristics of the machine tool, workpiece and cutting process; and can be best explained by a phenomenon called regeneration of waviness. The regeneration phenomenon was first explained by Tlusty [62] and Tobias [63]. The mechanism of chip regeneration is as follows: when one of the flutes is in cut, vibration is initiated due to an excitation force. The vibrating flute starts to 17 Chapter 2. Literature Review leave a wavy surface behind, which is then cut by the next flute, which also vibrates and leaves another wavy surface behind (see Figure 2.5.). The final chip thickness is determined by the difference between the two wavy surfaces. Depending on the phase shift between these waves, the maximum chip thickness may grow exponentially while oscillating at a chatter frequency that is close but not equal to a dominant structural mode in the system [64]. Chatter vibrations leave a wavy and poor surface behind leading to increased cutting forces which may destroy the workpiece, tool, tool holder and even the spindle. Like any other dynamical system, vibrations during machining can be stable, critically stable or unstable (chatter vibration) depending on the cutting conditions. There have been numerous efforts in the literature to map out stability charts for milling operations. There are some difficulties in the resolution of stability for milling. In milling, cutting forces rotate with the tool rotation and excite the structure in different orthogonal directions with varying magnitude. Time dependant - trigonometric coefficients, which orient cutting forces and change the direction of excitation, are called the directional coefficients. Tlusty [62] approached this problem by applying his orthogonal cutting stability formulation to milling with an average directional coefficient as well as an oriented transfer function. He proposed the chatter-free axial depth-of-cut as: –1 a lim = -------------------------------------------------- , b N 2K t ⋅ --- ⋅ ---- ⋅ G min ⋅ α d 2 (2.4) where K t is the cutting coefficient, b is the radial depth of cut, d is the diameter of the tool, N is the number of teeth, and the term G min ⋅ α is the real part of the transfer function oriented in the direction of chip thickness. Tobias [63] proposed a similar solution except that the cutting coeffi18 Chapter 2. Literature Review cient was expressed as a complex number with the dependency on effective cutting speed and chip thickness. Moreover, he developed stability lobes by relating the phase of vibrations to the spindle speed. Later, Merritt [65] used feedback control theory to verify the theory developed by Tlusty and Tobias and used the Nyquist stability criteria to obtain stability lobes. Although these theories were widely accepted and used by many researchers, they lacked accuracy for most of the milling operations because the two coupled directions were reduced to one directional flexibility (like turning operation) with the necessary adjustments on the transfer function (G). Workpiece a h0 n ε Workpiece h(t) y(t-T) y(t) Chuck Tool Tool Feed Direction h0 Figure 2.5 : Wave generation in turning. Minis and Yanushevsky [66] applied the theory of periodic differential equations and obtained the resulting characteristic equation of infinite order: 1 – λT det δ rk ⋅ I – --- ⋅ a ⋅ K t ⋅ ( 1 – e ) ⋅ W r – k ( λ + ikω T ) = 0 , 2 (2.5) where δ rk is the Kronecker delta, I is a (2m+1)-by-(2m+1) identity matrix, a is the depth of cut, K t is the cutting force coefficient in tangential direction, (r,k) are harmonic numbers, i.e r, k = 0, ± 1, ± 2, … ± m , ω T is the tooth passing frequency and W r – k is the oriented transfer function, which consists of directional milling coefficients and two dimensional transfer func- 19 Chapter 2. Literature Review tions. The equation above was truncated and numerically solved to determine the eigenvalues. Later, Budak and Altintas [67] solved Eq. (2.5) analytically, by considering the lowest Fourier mode (k=0) and λ = ± i ω c at the limit of stability (marginal stability) and obtained chatter freeaxial depth of cut as: ΛI 2 – 2 πΛ R a lim = ----------------- 1 + ⎛ -------⎞ , ⎝ Λ R⎠ NK t (2.6) where N is the number of teeth. Λ R and Λ I are the real and imaginary values of Λ , which is obtained as an eigenvalue of the following equation: det ( I + Λ ⋅ A ⋅ G ) = 0 , (2.7) where A ⋅ G is the oriented transfer function of the machine tool and workpiece in the feed and normal directions. Unlike Eq (2.4), the eigenvalue formulation in Eq (2.7) not only takes rotating cutting forces into account but it also allows dynamic excitation of orthogonal flexibilities individually. In recent years, stability of highly intermittent milling has been studied. When radial immersion is extremely small (e.g. less than 10% of the diameter), excitation forces (i.e. cutting forces) turn into impulses containing a large number of frequency components in their spectra. Stepan and Bayly et al. [68] studied intermittent machining by modeling dynamics of milling only in one direction. Two different sets of chatter frequencies were obtained from the criterion for stability obtained using the Floquet theory. The first set of frequencies referred to the Hopf bifurcation and was defined as: { ( ω H = ω c ± k ⋅ ω T ), k = 0, 1, 2… } , 20 (2.8) Chapter 2. Literature Review where ω c is the dominant chatter frequency and ω T = 2π ⁄ T is the tooth passing frequency. The second set referred to the Period Doubling Bifurcation leading to { ( ω PD = ( ω T ⁄ 2 ) ± k ⋅ ω T ), k = 0, 1, 2… } . (2.9) In addition to the Hopf Bifurcation, which is the type of bifurcation seen in traditional chatter, the Period Doubling (Flip Bifurcation) was shown to be a typical mechanism for stability loss in intermittent milling process. Davies et al. [69] also modelled intermittent milling by restricting structural flexibilities in the direction perpendicular to the feed (referred as Y axis) due to very low radial width of cut. In his model, the cutter was considered to be under free vibrations when it was out of cut, and delayed force vibrations within cut - since the ratio of time spent cutting to the spindle period was very small for low immersion milling. Combining these two vibration conditions - one becoming the boundary condition of the other - the forced periodic motion of the tool for intermittent milling was obtained. Stability solution revealed two types of vibrations in the system affecting the stability: the Hopf Bifurcation with the emergence of a traditional chatter frequency only and the Period Doubling (Flip Bifurcation) seen at odd multiples of the half-tooth passing frequency. Bayly et al. [70] proposed finite element analysis in time to solve stability in time domain for low immersion milling. The system was assumed to have a single-degree-of-freedom. Interrupted cutting was divided into two parts: in-cut and out-of-cut, similar to what was done in reference [69]. Time spent in-cut was divided into finite elements to model the interrupted cutting operations. A discrete system and the approximate solution was then matched with the analytical free response during out-of-cut. The complete, combined solution was cast in the form of a discrete map relating position and velocity at the beginning and end of each element to the corresponding values one period earlier. The eigenvalues of the linearized map were used to 21 Chapter 2. Literature Review determine stability. Bayly et al. [71], later, extended their stability analysis to cover two-degreeof-freedom systems by expressing the equations of motion in a matrix-vector form. Stepan et al. [72][73] proposed an alternative analytical solution for stability prediction of general milling operations. In their model, the delay differential equation representing dynamics of milling was approximated by a series of piece-wise ordinary differential equations and stability charts were constructed by point-by-point investigation of the parameter domain, i.e., spindle speed and depth of cut. This model was called the Semi-Discretization approach and proved to be an efficient time domain based analytical solution for these delay differential equations. Olgac et al. [74] presented an analytical solution for the stability of linear time invariant time delay systems. Time delay in the equation of motion, which makes the system nonlinear, was replaced by a bilinear expression and the stability of the new characteristic equation of the system was checked by using the Routh Hurwitz criterion while a range of depth of cut and spindle speed was swept. Merdol and Altintas [75] revisited and solved the multi frequency form of the milling stability model presented by Budak and Altintas [67] for low immersion case, and proved that added stability lobes could equally well be predicted in the frequency domain. They focused on a single-degree-of-freedom system and compared the stability charts with the benchmark experimental tests presented earlier by Stepan and Bayly et al. [68]. When the dynamics of metal cutting are correctly modeled, the analysis of stability of a milling operation reduces to solving an eigenvalue problem as shown in Eq (2.5). Except for a special case (when only the lowest Fourier component of the directional coefficients is used), stability lobes are generated iteratively requiring numerous calculation of the eigenvalues. Unfortunately, the size of this eigenvalue problem increases considerably depending on the number of harmonics 22 Chapter 2. Literature Review taken into account in characterizing the directional coefficients; therefore, the standard calculation of eigenvalues between iteration steps becomes computationally expensive. Since the system defined by the characteristic equation varies slightly at each iteration step, eigenvalues and eigenvectors can be estimated providing an initial guess for the next eigenvalue solution. Eigenvalue re-analysis, the study of variation of eigenvalues and eigenvectors due to variations in system parameters, has long been studied by various researchers. Approximate models for eigenvalue calculations are constructed based on the following factors: accuracy of the approximation; efficiency of the method (computational load); and ease of implementation [76]. One widely accepted re-analysis technique is the perturbation method used to calculate modal changes. Classical perturbation analysis [77]-[79] focused on calculation of distinct eigenvalues. Although the accuracy of the approximation can be improved by addition of higher order terms, classic analysis works only for small perturbations. Moreover, the required computational expense increases sharply due to the more complex formulas for higher order perturbations, and the accuracy of the solution inevitably suffers as computation errors accumulate. In order to improve accuracy for larger perturbations, a number of algorithms were proposed based on the incorporation of the Rayleigh quotient into the perturbation analysis [80]-[83]. Larger perturbations can be divided into a number of small perturbations, and the re-analysis computations are carried out step-bystep as suggested by Lu et al. [83]. The step-by-step perturbation method starts to lose its accuracy due to accumulation of computation errors when there are too many perturbation steps; therefore, an iterative method must also work in tandem with the step-by-step method [83][84]. Algorithms presented by Lu et al. [83] are found to be sufficient for eigenvalue re-analysis and 23 Chapter 2. Literature Review applied directly to increase computational efficiency of the solution of the machine tool stability problem in this thesis. 24 Chapter 3 Process Simulation 3.1. Introduction The simulation of the manufacturing processes in a computer environment is the ultimate goal of virtual machining. Although significant work has been reported in modeling individual but simplified machining processes, generalized mathematical modeling of arbitrary operations - an essential step to predict part specific operations in virtual machining - has not been addressed yet. This chapter presents a generalized mathematical model that uses parametric cutting tool geometry covering all possible helical end mills to model mechanics of 3-axis milling operations. Generalizations are carried out analytically or numerically depending on the complexity of the process and cutting tool geometry. Machining outputs such as cutting forces, form errors, power and torque demands are simulated and compared against experimental data. 3.2. Kinematics of Milling Process simulation of milling depends on both the kinematics of the cutting tool and cutter workpiece intersection. The kinematics of the tool is directly correlated to the motion of the machine tool axes. In this sense, milling machines are classified based on motions delivered by different axes, which are attached either to the table or the spindle. In 3-dimensional work space, a maximum of three translations and three rotations can be achieved. The majority of the milling machines used by industry comprise 3, 4, and 5-axis machine tools. The convention used by machine tool builders to assign axes is as follows: 25 Chapter 3. Process Simulation • x and y axes are orthogonal in direction and lie in the plane of the workpiece surface, • z axis is the remaining Cartesian axis assigned perpendicular to the workpiece surface, • rotations around x, y, and z are named as A, B, and C axes, respectively. A typical 3-axis horizontal milling machine is shown in Figure 3.1. For this machine, x and y axes are attached to the spindle where the cutter is mounted and the z axis motion is delivered along the guideways of the table to which the workpiece is clamped. Since each axis represents an independent motion, velocities at the tool centre and the workpiece can be represented by velocities of the corresponding axes, or equivalently, all relative motions can be referenced to the cutter, assuming the workpiece is stationary. Various combinations of velocity vectors define the final surface profile generated on the workpiece. For example, face milling requires only planar velocities whereas drilling, like plunge milling can be accomplished only by assigning a motion into the workpiece in the z direction. Different milling operations and corresponding velocity vector configurations are shown in Figure 3.1. Figure 3.1 : Typical horizontal milling machine, different milling operations. 26 Chapter 3. Process Simulation In milling, the tool follows a trochoidal path [8] due to simultaneous rotation and feed motions. A continuously changing chip geometry leads to varying cutting forces that are periodic at either tooth or spindle period - depending on the arrangement of tooth spacing, i.e. regular versus variable. In a general milling operation, as illustrated in Figure 3.2, there are two different coordinate systems defined in the Cartesian space. The global system represents the Machine Tool Coordinate System (MTCS) denoted by (X-Y-Z). A machine tool uses this coordinate system to drive the cutting tool to a desired location during machining. On the other hand, there is another coordinate system attached on the cutter called the Cutter Coordinate System (CCS) denoted by (x-y-z). The x-axis of the CCS is always aligned with the instantaneous feed direction on the XY plane, (represented by f xy ) which dynamically changes as the cutter generates features on a part. For 3-axis milling, the MTCS and the CCS have identical third orthogonal directions, z and Z axes. The rotation matrix between two coordinate frames can be defined as: ⎧ X ⎫ ⎪ ⎪ ⎨ Y ⎬ = ⎪ ⎪ ⎩ Z ⎭ cos α – sin α 0 ⎧⎪ x ⎫⎪ sin α cos α 0 ⎨ y ⎬ . ⎪ ⎪ 0 0 1 ⎩ z ⎭ (3.1) Throughout this thesis, the CCS is used as a default coordinate frame to model kinematics and mechanics of milling, however, the above transformation can be used to switch coordinate frames without loss of generality. 3.3. Geometric Relations for Milling Cutters Since cutting mechanics require geometric information of an end mill as an input, geometric variables of an end mill have to be analyzed in detail before proceeding with the modeling of cutting 27 Chapter 3. Process Simulation forces. In the literature, generalized geometric models for milling cutters were previously studied by Altintas and Engin [20]. They used standard CAD/CAM convention to define any milling cutter shape via seven parameters, which are highlighted in boxes in Figure 3.3.a. Although these seven parameters are sufficient to define the envelope of a milling cutter, the cutting flutes need to be wrapped around the cutter body to complete the geometric modelling. z Workpiece y x Y Z Y MTCS S f xy CC α X X Figure 3.2 : Coordinate systems in a milling operation. When a cutter is rotated at a certain spindle speed, n [rpm], the angular position of the cutter varies over time as: 2πn φ ( t ) = ---------- ⋅ t . 60 (3.2) At any time step t [sec], the angular position of an element located at an axial height z on flute j can be calculated as: φ j ( t, z ) = φ j = φ ( t ) + φ pj – ψ ( z ) , j = 1, …, N , 28 (3.3) Chapter 3. Process Simulation where N is the number of teeth, ψ ( z ) , which is dependant on the cutter envelope geometry and helix type, is the lag angle at elevation z due to the local helix angle. For an end mill, φ pj is the pitch angle of flute j and is expressed as: φ pj = ( j – 1 )φ p , (3.4) φ p = 2π ⁄ N , (3.5) and when teeth are evenly spaced (even pitch cutter). Angular positions are measured from the y-axis in a clock-wise direction as shown in Figure 3.3.b. In the following sections, geometric relations of general end mills will be presented. Based on the conventions used in [20], the envelope of the cutter is divided into three geometric portions namely TIP, ARC, and TAP, and the local geometry is defined for each portion. Once all the geometric relations are known, the unit vector perpendicular to cutter body at axial elevation z and angular position φ j can be defined as: n j ( z ) = sin κ ( z ) ⋅ sin φ j ( z ) ⋅ i + sin κ ( z ) ⋅ cos φ j ( z ) ⋅ j – cos κ ( z ) ⋅ k , (3.6) where i, j, and k are unit vectors in the principal directions x, y and z, respectively. This unit geometry vector will be important to determine the chip thickness later in this chapter. 3.3.1. Lag Angle ( ψ ( z ) ) Representation As mentioned before, complete geometric modelling of an end mill can be obtained only when cutting flutes are considered. The lag angle accounts for the local helix angle by describing the relative angular position of any point along the cutting flute, with respect to the flute tip’s angular position. For a general end mill, the lag angle along each section is represented as follows: • Tip Section ( 0 < z ≤ M z ): 29 Chapter 3. Process Simulation ln ( z ⋅ cot α ) ⋅ tan i 0 ψ TIP ( z ) = ---------------------------------------------- . cos α (3.7) To prevent singularity at the tool tip, z value must start from a finite but very small number, for example from z = M r ⁄ 20 as stated in reference [20]. • Arc Section ( M z ≤ z ≤ N z ): ( ( z – M z ) + R – R z ) ⋅ tan i 0 ψ ARC ( z ) = ---------------------------------------------------------------- + ψ TIP ( M z ) . R (3.8) • Tapered Section ( N z ≤ z ≤ h ), ⎧ ln ( 1 + ( ( z – N ) ⁄ N ) ⋅ tan β ) ⋅ tan i z r 0 ⎪ ------------------------------------------------------------------------------------- + ψ ARC ( N z ) if ( β ≠ 0 ) & Constant Helix ⎪ sin β ψ TAP ( z ) = ⎨ (3.9) ( ( z – N z ) ⁄ N r ) ⋅ tan i 0 + ψ ARC ( N z ) if ( β = 0 ) & Constant Helix ⎪ ⎪ ( ( z – N z ) ⁄ N r ) ⋅ tan i s + ψ ARC ( N z ) Constant Lead ⎩ where i 0 is the helix angle, and i s , the nominal helix for a constant lead cutter, is defined as: i s = atan [ ( 2πN r ) ⁄ ( L lead . cos β ) ] , (3.10) where L lead is the lead of the helical flute. Boundaries of each zone, TIP, ARC, and TAP are given as follows: 2 2 2 2 R z tan α + R r + ( R – R r ) tan2 α + 2R z R r tan α + ( R – R z ) -, M r = --------------------------------------------------------------------------------------------------------------------------------------------tan2 α + 1 (3.11) M z = M r tan α , (3.12) 2 2 2 2 ( R r – u ) tan β + R z – ( R – R z ) tan2 β + 2R z ( R r – u ) tan β + R – ( R r – u ) -, N r = -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------tan2 β + 1 (3.13) N z = u + N r tan β , (3.14) 30 Chapter 3. Process Simulation where u = ( 1 – tan α tan β )D ⁄ 2 , 0 ≤ α < π ⁄ 2 , and β < π ⁄ 2 . Some examples on how seven parameters are used to obtain lag angles for some commonly used end mills are as follows: • Cylindrical End Mill (R= R z = α = β =0, R r =D/2) with helix angle i 0 and diameter D, ψ ( z ) = k i0 z , (3.15) where k i 0 = ( 2 tan i 0 ) ⁄ D . • Tapered Helical Ball End Mill ( R r = α =0, R z =R= R b0 ) with taper angle β , Ball End Section with ball radius R b0 : ψ b ( z ) = ( tan i 0 ⁄ R b0 ) ⋅ z for z ≤ z 1 , (3.16) Tapered Section with constant helix ( i 0 ): ψ ( z ) = c 1 ⋅ ln [ 1 + c 2 ⋅ ( z – z 1 ) ] + ψ b ( z 1 ) for z ≥ z 1 , (3.17) Tapered Section with constant lead ( L lead ): ψ ( z ) = c 3 ⋅ ( z – z 1 ) + ψ b ( z 1 ) for z ≥ z 1 , where c 1 = tan i 0 ⁄ sin β , c 2 = tan β ⁄ R 0 (3.18) c 3 = tan i s ⁄ R 0 , i s = atan ( 2πR 0 ⁄ ( L lead . cos β ) ) , R0 = Nr , z1 = Mz . 3.3.2. Radius ( R ( z ) ) and Axial Immersion Angle ( κ ( z ) ) Representations The cutter envelope at each section (TIP, ARC, or TAP) is determined by the local radius. For a general end mill, the radius along the axis of the cutter is denoted by R(z) and axial immersion angle, κ ( z ) , is defined between axis of the cutter, and a line perpendicular to the tangent of the cutter envelope at axial height z. For each section, two variables can be represented as follows: • Tip Section: 31 Chapter 3. Process Simulation R ( z ) = z ⁄ tan α , (3.19) κ(z) = α . (3.20) • Arc Section: R(z) = 2 2 R – ( Rz – z ) + Rr , (3.21) R ( z ) – Rr κ ( z ) = asin ⎛ ----------------------⎞ . ⎝ ⎠ R (3.22) R ( z ) = u + z tan β , (3.23) κ(z) = π ⁄ 2 – β . (3.24) • Tapered Section z y φ+ ! S T κ(z) φ(t) φj(z) β x D TAP P Rr C N Rz M ARC TIP h R α O Mz xy Mr Nr # φpj ψ(z) Nz L @ !: Reference Cutter Position @: Tip Position for Flute j #: Position of Element Located at Axial Height z on Flute j (a) (b) Figure 3.3 : (a) Generalized tool geometry, (b) angle convention. 32 Chapter 3. Process Simulation 3.4. Mechanics of Milling - Generalized Kinematics Process simulation of a milling operation requires extensive and accurate modeling of the material removed by the cutter at each time step. The most fundamental parameter in cutting mechanics is the chip thickness. A typical milling operation with a general end mill is shown in Figure 3.4. In 3-axis milling, the feed motion of the cutter can be represented in a vector form as: f = c ⋅ fu , (3.25) where c [mm/rev.tooth] is the feed per tooth and f u is the unit feed vector in the direction of resultant tool motion expressed as: f u = f xy ⋅ i + f z ⋅ k , f u = 1 . (3.26) As a convention, the feed vector in XY plane is aligned with the x-axis of the cutter; therefore, the resultant planar feed is represented as a combination of X and Y feed components f xy = 2 2 f x + f y . Consider that the cutter is divided into differential axial elements with height dz. When the tool moves in the feed direction, flute j generates an uncut chip thickness that can be calculated by projecting the feed vector f onto the geometry vector n j ( z ) via a dot product ( • ) of two vectors: h j ( φ, z ) = n j ( z ) • f = c ⋅ ( f xy ⋅ sin κ ( z ) ⋅ sin φ j ( z ) – f z ⋅ cos κ ( z ) ) . (3.27) As material is removed from the workpiece, forces are exerted on the tool. Three orthogonal differential cutting forces generated by each axial element are represented in a rotating coordinate frame, and depicted in Figure 3.4. In order to relate differential rotating cutting forces to the instantaneous chip thickness, the Linear Edge Force Model [27] is implemented: 33 Chapter 3. Process Simulation ⎧ ⎪ dF rj ( φ, z ) ⎪ ⎨ dF tj ( φ, z ) ⎪ ⎪ dF aj ( φ, z ) ⎩ ⎫ ⎧ ⎪ ⎪ K re ⎪ ⎪ ⎬ = ⎨ K te ⎪ ⎪ ⎪ ⎪ K ae ⎩ ⎭ ⎫ ⎧ ⎪ ⎪ K rc ⎪ ⎪ ⎬ ⋅ dS ( z ) + ⎨ K tc ⎪ ⎪ ⎪ ⎪ K ac ⎭ ⎩ ⎫ ⎪ ⎪ ⎬ ⋅ h j ( φ, z ) ⋅ dS ( z ) , ⎪ ⎪ ⎭ (3.28) where dS(z) is the differential contact length expressed as: dS ( z ) = dz ⁄ sin κ ( z ) , (3.29) and dz is the axial height of the differential element. K rc, K tc, K ac are cutting force coefficients due to the shear and K re, K te, K ae are edge cutting force coefficients due to the rubbing of the tool flank with the workpiece, in radial, tangential and axial directions, respectively. z z Differential Chip Chip(ψ,φ,κ) n κ(z) zj,1 . Cutting Edge( j) zj,0 . dz Cutting Zone nj(z) dz dFtj f y dFrj x dFaj y hj dS x φj x Figure 3.4 : Mechanics and kinematics of 3-axis milling. 2 The oblique cutting coefficients, K rc, K tc, K ac [N/ mm ], are either obtained from mechanistic calibration tests or transformed from an orthogonal cutting database that contains shear stress ( τ ), shear angle ( φ c ), friction angle ( β a ) and the edge cutting force coefficients [86]. Although, in 34 Chapter 3. Process Simulation general, the edge cutting forces are dependent on the chip thickness, the linear model - Eq (3.28) is still acceptable as a first order approximation [29]. The influence of size effect, chip thickness, local rake and oblique angle on shear stress, shear angle and friction angles are considered in detail later. Differential cutting forces in Cartesian tool coordinates (CCS) can be obtained by using coordinate transformation: ⎧ ⎪ dF xj ( φ, z ) ⎪ ⎨ dF yj ( φ, z ) ⎪ ⎪ dF zj ( φ, z ) ⎩ ⎫ ⎪ – sin κ ( z ) sin φ j ( z ) – cos φ j ( z ) – cos κ ( z ) sin φ j ( z ) ⎪ ⎬ = – sin κ ( z ) cos φ j ( z ) sin φ j ( z ) – cos κ ( z ) cos φ j ( z ) ⎪ cos κ ( z ) 0 – sin κ ( z ) ⎪ ⎭ ⎧ ⎪ dF rj ( φ, z ) ⎪ ⎨ dF tj ( φ, z ) ⎪ ⎪ dF aj ( φ, z ) ⎩ ⎫ ⎪ ⎪ ⎬ . (3.30) ⎪ ⎪ ⎭ By substituting Eq (3.28) into Eq (3.30), the differential cutting forces can be represented as xyz dF j xyz, j ( φ, z ) = A e xyz, j ⋅ dz + A c ⋅ c ⋅ dz , (3.31) T where xyz dF j ( φ, ⎧ ⎫ z ) = ⎨ dF xj ( φ, z ) dF yj ( φ, z ) dF zj ( φ, z ) ⎬ , ⎩ ⎭ xyz, j Ae xyz, j Ac ⎧ ⎪ cos φ j sin φ j ⎪ – sin φ j K re – -------------- K te – ------------- K ae sin κ tan κ ⎪ ⎪ sin φ j cos φ = ⎨ – cos φ K + ------------ K te – --------------j K ae j re ⎪ sin κ tan κ ⎪ 1⎪ ---------K –K ⎪ tan κ re ae ⎩ ⎧ ⎪ – sin φ j sin κK rc – cos φ j K tc – sin φ j cos κK ac ⎪ = ⎨ – cos φ j sin κ K rc + sin φ j K tc – cos φ j cos κ K ac ⎪ cos κ K rc – sin κ K ac ⎪ ⎩ 35 ⎫ ⎪ ⎪ ⎪ ⎪ ⎬, ⎪ ⎪ ⎪ ⎪ ⎭ ⎫ ⎪ fz ⎞ ⎪ ⎛ - . ⎬ ⋅ ⎝ f xy sin φ j – ---------tan κ⎠ ⎪ ⎪ ⎭ (3.32) (3.33) Chapter 3. Process Simulation Differential cutting forces have to be integrated within the axial limits defined by the portion of the respective flute in-cut. If these limits are denoted by z j, 0 and z j, 1 at which the flute j enters into and exits the engagement zone, respectively, the total cutting forces acting on flute j can be expressed in an integral form as: z j, 1 xyz Fj ( φ ) = xyz, j ∫ ( Ae xyz, j c ) dz , + Ac (3.34) z j, 0 T ⎧ ⎫ where = ⎨ F xj ( φ ) F yj ( φ ) F zj ( φ ) ⎬ . The integrations are carried out to reduce cut⎩ ⎭ ting force experienced by each flute into the following form: xyz Fj ( φ ) xyz Fj xyz, j ( φ ) = A0 xyz, j ( φ ) + A1 (φ) ⋅ c , (3.35) where xyz, j A0 xyz, j A1 ⎫ ⎪ ⎪ ∫z j , 0 ⎪ ⎬. T ⎪ z j, 1 ⎫ xyz, j ⎪ = A d z ⎬ ∫z j , 0 c ⎪ ⎭ ⎭ T ⎧ ⎫ = ⎨ A x, j A y, j A z, j ⎬ = 0 0 ⎩ 0 ⎭ ⎧ = ⎨ A x, j A y, j A z, j 1 1 ⎩ 1 z j, 1 xyz, j A e dz (3.36) And finally, total cutting forces at time t ( φ = φ ( t ) ) are calculated by simply adding contribution of each flute: N F xyz (φ) = ∑ j=1 N xyz Fj ( φ ) = ∑ j=1 N xyz, j A0 + xyz, j ∑ A1 ⋅ c. (3.37) j=1 Note that the representation of cutting forces in Eq (3.37) is convenient because it separates the geometry of the cutter and the mechanics of milling from the feedrate in a closed form. Since xyz, j coefficients A 0 xyz, j , and A 1 need to be calculated once for a constant cutter-workpiece 36 Chapter 3. Process Simulation engagement at time t, the effect of the feedrate on process outputs can be easily taken into account. There are additional important process outputs beside cutting forces in x-y-z directions, some of which are spindle torque and power, and the resultant cutting force on the tool. The differential torque on flute j is obtained by multiplying differential tangential force with local radius of the cutter at axial height z: dTorque j ( φ, z ) = R ( z ) ⋅ dF tj ( φ, z ) , t, j t, j (3.38) t, j t, j where dF tj ( φ, z ) = ( A e + A c ⋅ c ) .dz and A e = K te ⁄ sin κ , A c = K tc ⋅ ( f xy sin φ j – f z ⁄ tan κ ) from Eq (3.28). Similar to cutting forces, the overall torque is obtained by carrying out the integration between engagement boundaries and adding the contribution of all flutes: N Torque ( φ ) = N ∑ Torque j ( φ, z ) = j=1 trq z j, 1 ∑ ∫z j=1 t, j t, j R ( z ) ( A e + A c ⋅ c ) ⋅ dz (3.39) j, 0 trq = A0 ( φ ) + A1 ( φ ) ⋅ c Since power is equal to torque times spindle speed (n), it can be trivially calculated as: pwr pwr Power ( φ ) = n ⋅ T orque ( φ ) = A 0 ( φ ) + A 1 ( φ ) ⋅ c , pwr where A i (3.40) trq ( φ ) = n ⋅ A i ( φ ) and i=0,1. Unlike torque and power, the resultant cutting force has a nonlinear relationship with the feedrate. Resultant force is defined as vector summation of x-y forces which has a magnitude of F res ( φ ) = 2 2 ( Fx ( φ ) ) + ( Fy ( φ ) ) . 37 (3.41) Chapter 3. Process Simulation Taking the square of both sides and substituting Eq (3.37) for the x and y forces, the following form is obtained: res F res ( φ ) = res res 2 A0 ( φ ) + A1 ( φ ) ⋅ c + A2 ( φ ) ⋅ c , (3.42) where coefficients are derived as: ⎫ ⎪ ⎪ res x x y y , A1 = 2 ( A0 ⋅ A1 + A0 ⋅ A1 ) ⎬ ⎪ ⎪ res x 2 y 2 A2 = ( A1 ) + ( A1 ) ⎭ res x 2 y 2 A0 = ( A0 ) + ( A0 ) q and A p = q, j ∑ Ap (3.43) (p = 0,1; q = x,y). All of these additional process outputs will be analyzed j more closely when they are used as constraints for feedrate optimization in the optimization chapter. 3.5. Cutting Force Coefficient Identification Two types of cutting force identification schemes are presented in this section. The orthogonal to oblique cutting transformation is based on geometric parameters of the cutter and cutting edge. This method transforms experimentally identified material properties into cutter-geometry-dependent cutting force coefficients for a variety of milling cutters. On the other hand, evaluation of cutting force coefficients for tools possessing complex cutting edges (e.g., chip breakers) are conducted by the mechanistic method. This model tunes the cutting force coefficients for a specific cutter and cutting conditions, hence cannot be used in a generalized sense. Details of both methods are given in the following sub-sections. 38 Chapter 3. Process Simulation 3.5.1. Orthogonal to Oblique Cutting Transformation Helical end milling is an oblique cutting technique where the cutting edge is no longer perpendicular to the cutting speed but inclined with an inclination angle ( i 0 ), which is essentially equal to the helix angle for an end mill. Orthogonal and oblique cutting geometries are shown in Figure 3.5. There are three important parameters for the determination of the cutting coefficients, the shear stress ( τ s ), the shear angle ( φ c ) and the friction angle ( β a ). These parameters can be determined from orthogonal cutting tests for different rake angles ( α n ) and cutting velocities [87]. The shear angle is an angle between the cutting speed (V) and the shear plane. The friction angle is defined on the rake face and the cutting edge is assumed to be sharp without edge radius. Oblique Cutting Orthogonal Cutting η Chip-flow angle Vc Rake face Chip Chip Vc Tool Tool i0 Fa Workpiece Workpiece Flank face V,Ft V,Ft Fr Fr h i0 Cutting edge inclination angle b Figure 3.5 : Orthogonal and oblique cutting geometries [64]. Using orthogonal cutting forces, namely feed/radial ( F rc ) and tangential ( F tc ), material cutting parameters, the shear stress, the shear angle and the friction angle can be obtained as: F tc cos φ c – F rc sin φ c τ s = -------------------------------------------------- , h b ------------sin φ c 39 (3.44) Chapter 3. Process Simulation r c cos α r φ c = atan ⎛ ---------------------------⎞ , ⎝ 1 – r c sin α r⎠ (3.45) F fc β a = α r + atan ⎛⎝ -------⎞⎠ , F tc (3.46) where h is the uncut chip thickness, b is the width of cut, α r is the rake angle and r c is the chip compression ratio, which is defined as a function of uncut (h) and cut chip thickness ( h c ) as, r c = h ⁄ h c . The oblique cutting coefficients for each chip segment are obtained by applying the orthogonal to oblique transformation proposed by Budak, Altintas and Armarego [30], where the chip flow angle ( η ) is assumed to be equal to the local oblique angle ( i 0 ), i.e. η = i 0 , as proposed by Stabler [88]. While Stabler’s rule may not yield an accurate prediction of chip flow angle, it does not influence the prediction of cutting forces significantly [27]. The cutting force coefficients for oblique cutting can be expressed: τs cos ( β n – α n ) + tan i 0 tan η sin β n K tc = -------------- --------------------------------------------------------------------------------------sin φ n 2 2 2 cos ( φ n + β n – α n ) + tan η sin β n τs sin ( β n – α n ) K rc = --------------------------- --------------------------------------------------------------------------------------- , sin φ n cos i 0 2 2 2 cos ( φ n + β n – α n ) + tan η sin β n (3.47) τs cos ( β n – α n ) tan i 0 – tan η sin β n K ac = -------------- --------------------------------------------------------------------------------------sin φ n 2 2 2 cos ( φ n + β n – α n ) + tan η sin β n where tan β n = tan β a cos η . In the oblique cutting transformation, the orthogonal shear angle is assumed to be equal to the normal shear angle in oblique cutting ( φ n ≡ φ c ); the normal rake angle is equal to the rake angle in orthogonal cutting ( α n ≡ α r ) and the chip flow angle is equal to the oblique angle ( η = i 0 ). The edge cutting force coefficients in tangential and radial directions, i.e. K te, K re , are calculated from the orthogonal cutting data by simply finding the force axis intercept 40 Chapter 3. Process Simulation of the force-feedrate graph where the feedrate is zero. In the verification part, these edge cutting force coefficients are expressed as a function of rake angle and cutting speed; however, they can just as well be assumed to be constant without losing too much accuracy. The edge cutting force coefficient in the axial direction K ae ; on the other hand, is taken as zero, which is known to be very small in oblique cutting [28]. 3.5.2. Mechanistic Model Although the orthogonal to oblique transformation provides cutting force coefficients for a variety of milling cutter geometries using three orthogonal parameters (shear stress, shear angle, friction coefficient), it is not always possible to identify these parameters when cutting tools with complex cutting edges are used. Moreover, constructing an orthogonal cutting database is very time consuming and costly. An alternative method is called the Mechanistic model [30]. In this model, average cutting forces per tooth period are measured while the feedrate is varied during each experiment with constant immersion and axial depth of cut (a). Any effect of run-out on measurements can be eliminated by dividing the total force per spindle revolution by the total number of teeth. Average cutting forces measured for m different feedrates is represented as a vector of length 3m, FM xyz . Since the total material removed per tooth is constant, the helix angle can be ignored for calculating average milling forces [64]. Due to the periodic nature of cutting forces in milling, averaging can be carried out only for one period, i.e., pitch angle φ p . Finally, analytical average milling forces can be described as: F xyz N ⎞ φp ⎛ 1 1- φex xyz xyz ⎜ ⎟ dφ = ----= ------ ∫ g ( φ )F ( φ ) F ( φ ) ⋅ dφ , j j ⎟ φp 0 ⎜ ∑ φ p ∫φst j=1 ⎝j = 1 ⎠ 41 (3.48) Chapter 3. Process Simulation where φ st and φ ex are the cutter entry and exit angles, respectively and g ( φ j ) is the switching function defined as: ⎧ 1 if ( φ st ≤ φ j ≤ φ ex ) ⎫ g ( φj ) = ⎨ ⎬. ⎩ 0 otherwise ⎭ (3.49) Considering a cylindrical end mill with κ ( z ) = π ⁄ 2 and zero helix angle, total cutting forces acting on flute j can be obtained by simply multiplying differential cutting forces in Eq (3.31) by the depth of cut (a): xyz Fj xyz, j ( φ ) = ( Ae xyz, j xyz, j Moreover, the immersion angle, φ j ( z ) , in A e ⋅ c) ⋅ a . + Ac xyz, j and A c (3.50) reduces to the following form when the helix angle is not in effect; φ j ( z ) = φ j = φ + ( j – 1 )φ p . Using Eqs (3.32) and (3.33) with φ j = 1 ( z ) = φ and planar milling feed conditions ( f xy =1, f z =0), the average cutting forces in Eq (3.48) are calculated as: F xyz 0.25 ( K tc cos 2φ – K rc ( 2φ – sin 2φ ) ) ⋅ c + ( – K te sin φ + K re cos φ ) Na = ------- 0.25 ( K rc cos 2φ + K tc ( 2φ – sin 2φ ) ) ⋅ c – ( K re sin φ + K te cos φ ) 2π K ac ⋅ c ⋅ cos φ – K ae φ φ ex , (3.51) φ st T ⎧ ⎫ = ⎨ F F F ⎬ . The above equation can be expressed in an alternative form, ⎩ x y z ⎭ which isolates unknown cutting force coefficients: where F xyz F xyz (c) = T(c) ⋅ K , T (3.52) ⎧ ⎫ where K = ⎨ K re K te K ae K rc K tc K ac ⎬ , and the transformation matrix T(c) represents the ⎩ ⎭ effect of both cutting geometry and feedrate (c) and is given by: 42 Chapter 3. Process Simulation 0 – c ( 2Δp 1 – Δs 2 ) cΔc 2 0 4Δc 1 – 4Δs 1 Na T ( c ) = ------- – 4Δs 1 – 4 Δc 1 , 0 cΔc 2 c ( 2Δp 1 – Δs 2 ) 0 8π 0 0 – 4 Δp 1 0 0 4cΔc 1 where Δc 1 = cos φ ex – cos φ st , Δs 1 = sin φ ex – sin φ st , (3.53) Δc 2 = cos 2φ ex – cos 2φ st , Δs 2 = sin 2φ ex – sin 2φ st and Δp 1 = φ ex – φ st . Averages of the measured cutting forces at m different feedrates are equated to analytically obtained average cutting forces: FM xyz = T ⋅ K, (3.54) T ⎧ ⎫ where T = ⎨ T ( c 1 ) T ( c 2 ) . . T ( c m ) ⎬ is a 3m-by-6 matrix, and c i ’s (i=1,...,m) are various ⎩ ⎭ feedrates used during experiments. Applying the Least Squares method to Eq (3.54), unknown cutting coefficients are successfully identified using the pseudo inverse as: T –1 T K = ( T T ) T FM xyz . (3.55) Note that these cutting coefficients are only valid for a specific cutter, depth of cut, and radial immersion. For cutters other than cylindrical end mills, integrations involved in Eq (3.36) must be carried out either analytically or numerically due to the effect of the axial immersion angle, κ ( z ) . For example, for a ball end mill: ⎛ z – R b0 ⎞ 1 ----------------= R d z atan ⎜ --------------------------------⎟ = R b0 atan ( – cot ( κ ( z ) ) ) , b0 ∫ sin κ ( z ) ⎝ z ( 2R – z )⎠ (3.56) b0 and 1 - dz ∫ -----------------tan κ ( z ) = z ⋅ ( 2R b0 – z ) . 43 (3.57) Chapter 3. Process Simulation Note that for general end mills, full or half immersion cases have to be used for cutting force coefficient identification. This way integration limits (z’s) do not depend on angular rotation of the cutter and Eq (3.48) will still be valid. 3.6. Cutter Engagement Feature (CEF) - Cutter Engagement Plane (CEP) The intersection between the cutter and the workpiece determines the region where the tool actually removes material. The shape of the intersection zone varies as a function of tool and workpiece geometries as well as width and depth of cut. The simplest intersection shape occurs between a cylindrical cutter and a prismatic workpiece as shown in Figure 3.6. For better visualization purposes, the intersection area, which is originally wrapped around the cutter, is unfolded and mapped onto a plane that has a vertical axis defined as the axial depth of cut and the horizontal axis as the angular position defined from the y-axis in clockwise direction. This plane is called the Cutter Engagement Plane (CEP). Each rectangle that is part of the intersection area will be referred to as the Cutter Engagement Feature (CEF) throughout this thesis. Using a similar mapping technique, the flute can be mapped onto the CEP, which then will appear as a slanted line when the helix is constant. The CEFs can be generalized when an arbitrary cutting tool intersects with a workpiece shaped other than prismatic. In this case, the CEFs will no longer be rectangular but they will be bounded by free form curves. An example of a generalized CEF is shown in Figure 3.7. In order to simulate cutting forces through Eq (3.37), first, boundaries of the engagement domain at which each flute enters and exits have to be determined, then, the integration is carried out over this range. Both calculations are performed either analytically or numerically. Although former is 44 Chapter 3. Process Simulation preferred due to the smaller computational burden, deriving a generalized analytical relation is not always possible. z y φ z x R CEP MAP 0 Flu CEF te CEF CEF φR πR 0 D L FO φ [rad] π UN Feed Direction Figure 3.6 : Cutter - workpiece intersection, CEF mapping. te Axial Depth of Cut Flu CEF 0 Immersion Angle π Figure 3.7 : Generalized CEF. It is important to state that the calculation of the engagement zone, i.e. intersection between cutter and workpiece, is not among one of the objectives of this research, it is rather taken as an input from a CAD system. In addition to available commercial packages that provide cutter/workpiece intersection maps, there is also considerable amount of research conducted in this field. Takata [89] used the Z-buffer method to obtain the cutter/workpiece intersection for machining force pre45 Chapter 3. Process Simulation diction. In many other studies, the Z-buffer method has been extensively used for machining process simulation [40][90][91]. Spence et al. [31] developed cutter-workpiece intersection for 2.5 D milling operations using constructive solid geometry, which was later transferred to the ACIS Boundary representation modeller [45][46]. 3.6.1. Determination of Integration Boundaries Integration boundaries to calculate machining state variables are obtained by calculating intersection points between each flute and engagement boundary. The simplest CEF, a rectangle, was studied earlier by Altintas et al. [26]. In their research, axial integration boundaries were obtained depending on the position of the flute, however, this method was limited to cylindrical end mills. In this research, the method in reference [26] is generalized so that integration boundaries for any cutter geometry can be determined. z z 0 2 4 zlim,1 3 zj,1 1 π 0 φ zlim,0 (zj,0) 0 (a) φj,1 φst φj,0 φex π φ (b) Figure 3.8 : Rectangular CEF; (a) possible intersection cases, (b) parameter definition for CEFs. Figure 3.8.a shows five distinct cases for flute-CEF intersection, one of which is detailed in Figure 3.8.b. Parameters bounding the CEF, i.e z lim, 0 , z lim, 1 , φ st , φ ex are provided by a CAD system hence known. For a given tool position φ = φ ( t ) , angular limits of the flute are determined 46 Chapter 3. Process Simulation both at the bottom and the top of the CEF using Eq (3.3), i.e., φ j, 0 = φ j ( z lim, 0 ) and φ j, 1 = φ j ( z lim, 1 ) . Based on these limits, a decision mechanism shown in Figure 3.9 is used to obtain axial positions between which the flute is in-cut, i.e. integration boundaries: φj,0>φst & φj,0<=φex NO ( ) YES ( ) φj,0>φex & φj,1<φex zj,0 = z(φj,0) zj,0 = z(φex) φj,1>φst & φj,1<=φex φj,1>φst 4 zj,1= z(φj,1) zj,1= z(φst) zj,0= z(φj,1) zj,1= z(φst) 1 0 zj,0 = 0.0 zj,1 = 0.0 NO ( ) YES ( ) NO ( ) YES ( ) NO ( ) YES ( ) 2 3 Figure 3.9 : Calculation of integration boundaries for general end mills. In Figure 3.9, z ( φ j ) represents a function that calculates the axial height z for a given angular position of flute j. Mathematically, it is equivalent to solving for z for a given angular position, φ j , in Eq (3.3). For a cylindrical end mill, the lag angle, ψ ( z ) = z tan i 0 ⁄ ( D ⁄ 2 ) , is substituted and z is simply solved as: ( φ ( t ) + φ pj – φ j ) ⋅ D z ( φ j ) = ------------------------------------------------- (Cylindrical End Mill) , 2 tan i 0 (3.58) where i 0 is the helix angle, and D is the diameter of the tool. The above solution is analytical because rectangular CEF imposes the condition that engagement boundaries z j, 0, z j, 1 ∈ [ z lim, 0, z lim, 1 ] can only vary along either φ j = φ st or φ j = φ ex line; however, such a useful property does not hold when more complex cutters like ball end mills are considered. 47 Chapter 3. Process Simulation 9 y D = 19.05 mm yos = 2 mm w = 8 mm 8 7 φst(z) + R(z) z yos w x z [mm] 6 workpiece 5 φbnd(z) = φst(z) φbnd(z) = φex(z) 4 3 φex(z) 2 CEF 1 y x 0 60 120 180 φ [deg] (a) (b) Figure 3.10 : (a) Planar machining with a ball end mill, (b) example CEF for ball end milling. D Consider a prismatic workpiece being machined by a ball end mill as shown in Figure 3.10.a, where w is the width of the workpiece and y os is the offset from its centre line. Since the cutter radius varies at each elevation along the z-axis, the entry and exit angles do not remain constant, i.e. φ st = φ st ( z ) and φ ex = φ ex ( z ) ; therefore, the CEF - shown in Figure 3.10.b - will not be a rectangle unless it is a full immersion case. Variation of immersion along cutter axis is not limited to ball end mills, but it is also true for any kind of general end mill. For a known workpiece width and centre line offset, start and exit angles of immersion at each axial point z or better known as the boundaries of the CEF, φ bnd ( z ) = { φ st ( z ), φ ex ( z ) } , can be calculated using local cutter radius as follows: 48 Chapter 3. Process Simulation Start Angle: yval = yos + w/2, bnd = st Exit Angle: yval = yos - w/2, bnd = ex yval >= 0 YES ( ) NO ( ) R(z) <= yval YES ( ) φbnd (z) = 0 NO ( ) φbnd (z) = cos (yval / R(z)) -1 R(z) <= |yval| YES ( ) NO ( ) φbnd (z) = π -1 φbnd (z) = π/2 + sin (|yval| / R(z)) If the workpiece is not prismatic but has a profile on its surface (Figure 3.7), then the intersection between the cutter and workpiece might not be obtained analytically as described above. In the most general case, the CEF is bounded by a curve referred as φ bnd ( z ) , and the variation of the integration boundaries as tool rotates is depicted in Figure 3.11. 49 Chapter 3. Process Simulation zmax Advancement of the flute as cutter rotates Integration Boundaries φbnd (z) z zmin π 0 Immersion Angle Immersion Angle Immersion Angle Figure 3.11 : Change of integration boundaries over time for the most general case. Assume that at a certain time step, flute j intersects with the CEF, φ bnd ( z ) curve, at axial height z root , which leads to the following relation at the intersection point: φ bnd ( z root ) = φ j ( z root ) = φ ( t ) + ( j – 1 )φ p – ψ ( z root ) . (3.59) This equation might not be solved for z root analytically due to the presence of non-linear terms both in φ bnd ( z ) and ψ ( z ) . Redefining Eq (3.59) in a residue form: g ( z ) = φ bnd ( z ) – φ ( t ) – ( j – 1 )φ p + ψ ( z ) , (3.60) roots satisfying g ( z root ) = 0 can be sought numerically. Brent’s numerical root finding algorithm [92] is selected due to its high convergence rate. This algorithm first brackets the root(s) within the domain of minimum and maximum axial points. A root is called "bracketed" between z 0 and z 2 when g ( z 0 ) and g ( z 2 ) have opposite signs (Figure 3.12). To account for all the roots, the domain is divided into a number of segments and all possible roots are bracketed. Next, algorithm fine tunes the roots which lie between brackets. 50 Chapter 3. Process Simulation g(z) zroot,1 z0 z1 z 2 zroot,2 z Figure 3.12 : Root bracketing. Brent’s method utilizes two methods for its fine tuning: the bisection method and inverse quadratic interpolation. The bisection method evaluates the function value at the midpoint of the bracket, checks its sign against g ( z 0 ) and g ( z 2 ) , and replaces the new point on the side that has the same sign as g ( z 0 ) or g ( z 2 ) . The method stops when two limiting brackets are within an error bound specified by the user. Since signs of the limits are always checked at each iteration step, this method is guaranteed to converge; however, at the cost of speed, i.e. linear convergence. On the other hand, inverse quadratic interpolation uses three points to fit an inverse quadratic func2 tion (z = g( y )) and determines the next root estimate at y=0 [92] as: g0 g1 z2 g1 g2 z0 g2 g0 z1 z = ------------------------------------------- + ------------------------------------------- + ------------------------------------------- , ( g2 – g0 ) ( g2 – g1 ) ( g0 – g1 ) ( g0 – g2 ) ( g1 – g2 ) ( g1 – g0 ) (3.61) where g i = g ( z i ) , i = 0, 1, 2, and z 1 is selected such that z 0 < z 1 < z 2 . Although a quadratic method has a higher convergence rate, its success depends solely on the smoothness of the function. In conclusion, Brent’s method successfully combines the fast convergence of inverse quadratic interpolation with the reliability of the bisection method, by switching between these two methods depending on the local behavior of the function g. It should also be mentioned that the 51 Chapter 3. Process Simulation Brent’s method does not need the derivative information of function g, hence complexity of g is not an obstacle when the solution is sought for complex CEFs. 3.6.2. Performing Integration Having calculated the boundaries of the integral, machining state variables - for example cutting forces in Eq (3.35) can be acquired by taking the integral along the cutting edge in contact with workpiece. Differential cutting forces in Eq (3.31) contain two geometric terms varying as a function of axial height z: axial immersion angle κ ( z ) , and lag angle ψ ( z ) . Except for cylindrical end mills, one or both of these terms bring non-linearity into the integrand. Moreover, cutting force coefficients, K rc, K tc, K ac , generally vary along the cutter axis for many complex cutters such as ball end and tapered ball end mills acting as an alternative source of nonlinearity. A numerical integration algorithm is indispensable for calculating machining state variables to handle cases with nonlinearities. The integration problem is defined as: z root, 2 I = ∫z g ( z ) dz , (3.62) root, 1 where g(z) is expressed by Eqs (3.32) and (3.33) depending on type of the machining state variable. Exact solution to Eq (3.62) can be approximated numerically using Trapezium rule as: z root, 2 ∫z 0 g ( z ) dz ≈ I k , (3.63) root, 1 k 0 where I k represent the Trapezoidal approximation with 2 panels, i.e., step size of k δ k = ( z root, 2 – z root, 1 ) ⁄ 2 , 2 and error estimate O ( δ ) : 52 (3.64) Chapter 3. Process Simulation 0 2 4 Ik = I ( 0 ) + c2 δ + c4 δ + … , (3.65) where I(0) is the approximate integration value, c 2 and c 4 are constants. A better approximation can be obtained by exploiting the error expansion. For example, consider two trapezoidal approximations with one has a half step size with respect to the other: 2 0 4 I1 = I ( 0 ) + c2 δ + c4 δ + … , (3.66) 2 4 0 I 2 = I ( 0 ) + 4c 2 δ + 16c 4 δ + … . (3.67) The first error terms are eliminated by multiplying Eq (3.66) by 4 and subtracting it from Eq 4 (3.67), an improved estimate is obtained with error estimate O( δ ): 0 1 I1 0 4I 1 – I 2 4 = ------------------ = I ( 0 ) – 4 c 4 δ + … . 3 (3.68) 1 The method of error elimination is called Richardson extrapolation and I 1 is called the first Richm ardson Extrapolant. In general, I k is the mth Richardson Extrapolant obtained by halving the step size, δ , at each stage, and the successive Richardson Extrapolants can be calculated by the following recursion: m m+1 Ik = m Ik + 1 m Ik + 1 – Ik -, + --------------------m 4 –1 (3.69) where m = 0, …, M , k = 1, …, m + 1 ,and M is the maximum number of refinements. The error term for this approximation becomes O( δ 2 m+1 ). The recursive formula can also be defined in a matrix form as follows: 53 Chapter 3. Process Simulation 0 I1 0 I2 0 I3 0 1 I1 1 I2 1 2 I1 2 . 3 I4 I3 I2 I1 . . . . . . . . . . . 0 Ik 1 Ik – 1 . . . . m I1 The numerical integration algorithm that employs both Richardson extrapolation and adaptive integration is called Romberg integration [92], which is used in this thesis to integrate differential cutting forces and other machining state variables when numerical integration due to non-linearities is required. In simulations, the maximum number of refinements, M, was limited to 20. Adaptive integration is achieved by simply refining step size, δ , until the difference between successive approximate integrals is less than some specified tolerance. In the matrix above, each column represents refinement of the sampling of the integrand as k is increased while m is kept constant. Romberg integration is quite powerful for sufficiently smooth integrands that do not posses any singular values over the interval including end points [92]. In this case, the integrands are the elements of force coefficient matrices (A) in Eqs (3.32) and (3.33), and among all, critical terms are the ones that have denominators comprised of either 1 ⁄ sin κ ( z ) or 1 ⁄ tan κ ( z ) . The singularity can appear only when κ ( z ) = 0 since κ ( z ) = π is not physically possible as seen in Figure (3.3). From Eq (3.22), κ angle of a ball end tool reduces to sin κ ( z ) = R ( z ) ⁄ R b0 , where R ( z ) = (3.70) z ( 2R b0 – z ) , R b0 is the ball radius, and from the geometry of ball end: cos κ ( z ) = ( R b0 – z ) ⁄ R b0 . 54 (3.71) Chapter 3. Process Simulation If the first tangent line (TIP) does not exist so that tool geometry starts with an arc then the sine of axial immersion angle varies considerably along the arc segment close to the tool tip, and becomes zero at the tip point (z = 0), which introduces an infinity in some integrands of the forces in Eq (3.32). In summary, any tool carrying a ball end part has to be analyzed carefully while numerically calculating cutting forces because some of the differential forces (integrands) are neither smooth nor non-singular at or near the tool tip. The following change of variable is proposed to relax singularity condition: z = R b0 ( 1 + sin ξ ) , (3.72) and it is substituted into radius: R(ξ) = 2 R b0 ( 1 + sin ξ ) ( 2R b0 – R b0 ( 1 + sin ξ ) ) = R b0 ( 1 – sin ξ ) = R b0 cos ξ . (3.73) Using the above conversion, new integrands are now obtained as: 1 1 ----------------dz → ----------- R b0 cos ξdξ = R b0 dξ , sin κ ( z ) cos ξ (3.74) 1 – sin ξ -----------------dz → -------------- R b0 cos ξdξ = – R b0 sin ξd ξ . tan κ ( z ) cos ξ (3.75) Note that at the tool tip, the relation ( z = 0 ) → ( ξ = 3π ⁄ 2 ) holds; therefore, singularity is avoided and computation time is immensely reduced due to faster convergence. The variation of edge integrand in y direction is shown in Figure 3.13 before and after the change of variable is applied to smooth out the integrand. 55 y,j Ae [N] ( x1e3 ) Chapter 3. Process Simulation 0.5 0.2 0 0.1 -0.5 0 -1.0 -0.1 -1.5 0 2 4 6 -0.2 -2.0 z [mm] -1.5 -1 ξ [rad] -0.5 0 Figure 3.13 : Variation of edge integrand in y direction for a 12 mm ball end mill. 3.6.3. Special Case: Flute Re-cut Under special circumstances, for example when depth of cut and/or helix angle on the cutter are very high, the cutting flute might enter into the same cutting region at several different axial locations at a fixed tool position. This special case, shown in Figure 3.14.b, has to be taken into consideration during process simulation for accurate calculation of cutting forces. Chip Geometry zmax (a) Flute Numbers Depth of Cut (a) Immersion Angle,( φ) #2 #1 zmax #2 #1 (b) Figure 3.14 : (a) Low helix cutter (4-fluted) with small depth of cut, (b) high helix cutter (2-fluted) with large depth of cut, (Inset) chip geometry. 56 Chapter 3. Process Simulation Consider a case when cutter-workpiece engagement is composed of multi CEFs as shown in Figure 3.15. The Cutter Engagement Plane (CEP) bounded between z ∈ [ z min = 0, z max ] and φ ∈ [ 0, π ] defines the possible Cutting Region at the front of the cutter, where z max is the highest of all CEF elevations for a given cutter and workpiece combination. The condition for a flute to re-enter into the cutting region can be derived using the lag angle ( ψ ), which is previously defined in Section 3.3.1 for each cutter geometry. In general, if the lag angle of a flute is bigger than 2 π for maximum height, z max , then it is concluded that very same flute exists in cutting region more than once at different axial locations. This is named as the re-cut condition. z z zmax CEF-3 CEF-3 Cutting Region CEF-1 0 CEF-1 CEF-2 φ [rad] π 0 CEF-2 φ [rad] π Figure 3.15 : Cutting region at the front of the cutter. During process simulation, flute re-cut condition can be detected by checking the number of occurrence: n occ ⎧ n pi ⎪ -----if n pi is even ⎪ 2 +1 = ⎨ , ⎪ ( n pi – 1 ) ⎪ -------------------- + 1 if n pi is odd 2 ⎩ (3.76) where ψ ( z max ) n pi = ceil ⎛ -----------------------⎞ , ⎝ ⎠ π 57 (3.77) Chapter 3. Process Simulation and the ceil(x) function returns a value representing the smallest integer that is greater than or equal to x. Flute re-cut will appear when n occ >1. Figure 3.16 presents different cutting configurations for same depth of cut but various helix angles. Four shaded zones (1,...,4) represent same cutting region that will overlap when z- φ plane is wrapped around body of the cutter. In case (a), flute helix is not high, therefore, it does not the re-cut the workpiece, i.e., the flute only resides in Zone 1. On the other hand, for the remaining cases, the helix is selected high enough to ensure that the flute enters into the cutting region more than once (Zones 1,2,...). In such cases, if any process output is to be calculated, then "Equivalent Flute-Positions" column shows how many more flute angles have to be used to consider re-cut condition. For example, when case (d) occurs, the flute re-cuts the workpiece three times, therefore, for a fixed tool position, forces acting on that flute are calculated by adding the contribution of forces for each occurrence of the same flute in the cutting zone. Similar analysis can be done when the CEP is defined at back of the cutter, i.e., φ ∈ [ π, 2π ] . 3.7. Verification Cutting experiments were conducted to verify the proposed model. 20 [mm] diameter cylindrical end mill with 30 [degree] helix was used to cut Aluminum Al 7050 with the following cutting 2 2 2 coefficients: K tc = 796 [N/ mm ], K rc = 169 [N/ mm ], K ac = 222 [N/ mm ], K te = 28 [N/mm], K re = 31 [N/mm], and K ae = 1.4 [N/mm]. The engagement domain is comprised of three different axial depths of cut (3, 7, 11 [mm]) as seen in Figure 3.17. The spindle speed and the feedrate used in tests were 500 [rpm] and 200 [mm/min], respectively. Simulation results are compared against measurements in Figure 3.17. 58 Chapter 3. Process Simulation Real Flute-Position Equivalent Flute-Positions z 3 2 Fl u 4 z 1 te -6π -5π -4π -3π -2π -π 0 π 2π φ = π 2π Section npi nocc a 1 1 b 2 2 c 3 2 d 4 3 φ (a) z 4 3 2 z 1 π 2π -6π -5π -4π -3π -2π -π φ z = π 2π φ + π 2π φ (b) z z 4 3 2 1 π 2π -6π -5π -4π -3π -2π -π φ = 3 2 -6π -5π -4π -3π -2π -π π 2π (c) z 4 z φ + π 2π z 1 π 2π φ Fflute z = z + = (d) φ π 2π ! idFflute φ π 2π + @ idFflute φ + π 2π + # idFflute + ... Figure 3.16 : Flute re-cut example, same depth, different helix angles. 59 4π φ Chapter 3. Process Simulation Measured Simulated Fx [N] 1000 500 0 Fy [N] -500 600 400 200 0 -200 Fz [N] 100 0 -100 -200 -300 0 180 360 Tool Rotation [deg] Axial Depth of Cut [mm] z MAP FORCE Cutter 11 7 3 Immersion Angle φ feed direction Figure 3.17 : Multiple immersion milling, engagement map, simulated versus measured cutting forces. A second set of experiments was conducted on an industrial part, a gear box cover. The part was made out of Al 6061 using two different sizes of milling cutters. The first cutter was a 10 [mm] end mill with 2 flutes, whereas the second one was a 25 [mm] indexable cutter with 2 inserts. Cutting coefficients, presented in Table 3.1, were identified using the mechanistic model, details of which were given in the "Mechanistic Model" section. Spindle speed values of 17000 [rpm] and 20000 [rpm] were used for 25 [mm] and 10 [mm] cutters, respectively. The maximum chip load on the cutter was maintained for varying widths of cut. According to the tool manufacturer’s cata- 60 Chapter 3. Process Simulation log, recommended chip load was 0.1 [mm/tooth.rev] for the 10 [mm] end mill, and 0.2 [mm/ tooth.rev] for the 25 [mm] indexable cutter. Tool Type 10 [mm] 25 [mm] Table 3.1 : Cutting force coefficients of Al 6061 K te K re K tc 2 [N/mm] [N/mm] [N/ mm ] 38.38 -9.792 403.396 21.34 23.141 662.296 K rc 2 [N/ mm ] 98.911 72.86 Various milling operations such as profile, face, and helical milling are shown on a piece part in Figure 3.18. For each operation, the measured experimental cutting forces are compared against simulated ones in Figures 3.19-3.21. Although measurements captured complete variation of cutting forces over time, cutting forces were simulated only to obtain minimum and maximum values for varying immersion and tool positions. The difference in process times and cutting force magnitudes between simulation and actual cut is due to the trajectory generation, and the limited maximum feed speed of the machine, which was set to 5000 [mm/min]. 61 Chapter 3. Process Simulation Figure 3.18 : (a) Face milling, (b) right side profiling, (c) top hole helical milling, (d) middle hole helical milling. 62 Chapter 3. Process Simulation Fx [N] 800 400 0 -400 -800 Fy [N] 800 400 0 -400 -800 0 5 10 15 20 25 time [s] (a) 30 35 40 45 Fx [N] 800 400 0 -400 -800 Fy [N] 800 400 0 -400 -800 0 5 10 15 time [s] (b) 20 25 30 Figure 3.19 : Face milling a) measured, (b) simulated cutting forces. 63 33 Fx [N] Fx [N] 1500 1000 500 0 -500 -1000 -1500 1500 1000 500 0 -500 -1000 -1500 Fy [N] 1500 1000 500 0 -500 -1000 -1500 Fy [N] Chapter 3. Process Simulation 1500 1000 500 0 -500 -1000 -1500 0 2 0 1.5 4 6 3.5 8 10 time [s] (a) 5.5 7.5 12 14 9.5 16 18 11.5 12 time [s] (b) Figure 3.20 : Profile milling a) measured, (b) simulated cutting forces. 64 Chapter 3. Process Simulation Fx [N] 2000 1000 0 -1000 -2000 Fy [N] 2000 1000 0 -1000 -2000 0 0.5 1 1.5 2 2.5 3 3.5 time [s] (a) Fx [N] 2000 1000 0 -1000 -2000 Fy [N] 2000 1000 0 -1000 -2000 0.1 0.3 0.5 0.7 0.9 time [s] (b) 1.1 1.3 1.5 1.7 Figure 3.21 : Top hole helical milling a) measured, (b) simulated cutting forces. In co-operation with Sandvik Coromant, Sweden, a die was machined to verify the proposed algorithms. The original part, a die for a car door panel, was first scaled down so that the part could be machined using a conventional machine tool. Steel casting of size 425x325x220 [mm] was provided by Uddeholm. Different views of both stock (casting) and final part are given in Figure 3.22.a-d. 65 Chapter 3. Process Simulation (a) (b) (c) (d) Figure 3.22 : Test part stock and final shape (a) isometric; (b) left side; (c) front; (d) top views. The operation plan for this part contained 9 steps with 8 different milling cutters. Cutters ranging from 63 [mm] to 4 [mm] were used to complete all operations. Important geometric parameters and serial numbers of the tools from Sandvik Coromant catalogue are presented in Table 3.2, and the operation list is provided in Table 3.3. The workpiece material was Uddeholm Carmo - Carmo Cast with a hardness of 270 [HB]. Cutting conditions for material modeling tests consisted of various depths of cut depending on ball or edge radius, and feedrates mostly selected from Sandvik Coromant’s Main Catalogue. Cutting forces were collected using CUTPRO’s data acquisition module, MALDAQ. Sampling frequency 66 Chapter 3. Process Simulation was set in such a way that forces were measured at least at each 3 [degree] angular rotation of the tool during cutting. Cutting tests were performed under dry conditions, and build-up edge was observed at high feedrates and depths of cut for some cutters. Such cutting tests were rejected from the measurement pool to eliminate the bias due to improper chip formation. Measured cutting forces were averaged out for approximately 20 [revolutions] of the cutter and the unknown cutting force coefficients were identified for each axial depth of cut. To reduce the number of coefficient sets to one per tool, coefficients were averaged over depth of cut and the final results are tabulated in Table 3.4 for each tool used in the tests. . Table 3.2 : List of milling cutters used in experiments. Tool No N H [deg] D [mm] SD [mm] Cutter Body Insert Type Shank Type Cutter Type Overhang T63 4 7 63 Arbor R300-06Q22-12L R300-1240M-KH/3040 Arbor Bull Nose 154 mm T20 2 -10 20 25 R216-20B25-050 R216-20 T3 M-M 1025 Weldon Ball End 50/51 mm T20F 2 0 20 25 R216F-20A25C-115 R216F-20 50 E-L P10A Cylindrical Ball End at the stepover T16 2 -10 16 20 R216-16B20-040 R216-16 03 M-M 1025 Weldon Ball End 50 mm T12 4 27 12 12 R216.44-12030-AI12G Cylindrical Ball End 37 mm T8 2 27 8 8 R216.42-08030-AI08G Cylindrical Ball End 27 mm T6 2 25 6 6 R216.42-06030-AI06G Cylindrical Ball End 21 mm T4 2 27 4 6 R216.42-04030-AI04G Cylindrical Ball End 20 mm N: Number of Flutes/Inserts, H: Helix Angle, D: Tool Diameter, SD: Shank Diameter Table 3.3 : List of Operations. 67 Chapter 3. Process Simulation In an effort to measure cutting forces during machining, casting was mounted on a measurement platform: Kistler 9281B Dynamometer, which was then mounted on the machine tool’s table. MAL’s 4-Channel I/O Box with MALDAQ measurement software was used to collect cutting force data. In addition to forces, the sound data from a high quality ICP microphone placed inside the machine was also collected with the same sampling rate as for cutting forces. Table 3.4 : Average cutting force coefficients for Uddeholm Carmo - Carmo Cast. T63 T20 T20F T16 T8 Ktc [N/mm²] 2312.7 2294.2 2024.1 2068.1 2156 Krc [N/mm²] Kac [N/mm²] 1443.1 -167.5 1137.5 -536 1132.9 42.2 1120.9 -99.54 924.1 172.9 Kte [N/mm] 62.82 63.53 53.73 41.4 29.44 Kre [N/mm] 112.74 100.24 64.2 72.63 19.75 Kae [N/mm] 20.15 5.45 1.54 19.71 1.24 Measured and simulated cutting forces in three principal directions of the machine tool are shown in Figure 3.25 for semi-finishing of the walls with T63. Note that, although complete time history of cutting forces was captured by the dynamometer, the simulation contained only minimum and maximum values of cutting forces at each engagement map. In general, a close match between measured and simulated cutting forces is evident; however, there are sections where cutting forces deviate from simulation drastically. For example around 36th [second], although the measured axial (z) force has a minimum of -200 [N], predicted value is still at zero. In order to determine the source of such deviation, a detailed analysis was performed starting with the engagement map shown in Figure 3.26.a. The intersection between the tool and workpiece reveals that cutter is currently in down milling with about 25 [degree] immersion. Considering that T63 has 4 flutes with 90 [degree] pitch angle, cutter must stay out of cut 90-25=65 [degrees] within each tooth period. With the tool rotating at 1011 [rpm], the time spent out of cut ( t out ) and spindle period ( t period ) 68 Chapter 3. Process Simulation 60 ( 65 ) 60 ( 360 ) are calculated as t out = -------------------------- = 10.7 [ μs ] , t period = -------------------------- = 59.3 [ μs ] . Figure 360 ( 1011 ) 360 ( 1011 ) 3.26.b shows cutting forces for a time window of two spindle periods, i.e. approximately 0.12 [seconds]. Four peaks due to four flutes are clear in z cutting forces within one spindle period, however, they never drop to zero when all flutes are out of cut. Moreover, x and y cutting forces have more than four peaks in one spindle period. This kind of measurement distortion is attributed to the reduced bandwidth of the dynamometer due to the heavy workpiece (~150 [kg]); therefore, it can be concluded that simulation results are still in good agreement with the measured cutting forces when dynamometer distortion is discarded. Simulated cutting forces are compared against measured cutting forces at different sections of the part as well. Results are presented in Figures 3.27-3.30. In finishing operations, although cutting forces were recorded, the dynamometer signal was dominated by noise due to extremely light cuts. Reducing the scale setting on the charge amplifier did not help to reduce noise levels, therefore, these measurements were not left out of the analysis. In addition to measurement errors, there are other sources of errors in predicting cutting forces, some of which can be listed as follows: - Cutting Force Coefficients: One set of average cutting force coefficients was used for each cutter which might have led to deviation in extremely low and high depths of cut (see Figure 3.23). Since identification tests were run for multiple depths and chip thickness values, this type of error can be reduced by using a more complex cutting force coefficient model taking depth variation and size effect into account. 69 Cutting Force Coefficient Chapter 3. Process Simulation Average Axial Depth of Cut Figure 3.23 : Cutting force coefficient variation. - Real Feedrate: Since the CNC characteristics were estimated, simulated feedrates did not exactly match the values during real machining (see Figure 3.24). Any error on feedrate estimation was directly reflected on cutting forces. Feedrate Commanded Actual Along Tool Path Figure 3.24 : Actual versus commanded feedrate. - Engagement: At some sections of the workpiece, cutter workpiece intersection was either noisy or the map resolution was not small enough to capture true intersection, which might have led to under or over estimation of cutting forces. -Flute Geometry: Most of the cutters had large run-outs resulting in one flute cutting more than the other, which in turn led to conservative cutting force prediction. Moreover, periphery and centre inserts had different geometries although they were assumed identical for simulations. Lastly, tool wear over time affecting the edge geometry hence 70 Chapter 3. Process Simulation resulting in increase in cutting force coefficients was not taken into account. In all simulations, coefficients identified with sharp tools were used. 71 Chapter 3. Process Simulation Measured Forces Fx [N] 1000 0 -1000 0 20 40 60 80 100 120 140 160 180 0 20 40 60 80 100 120 140 160 180 0 20 40 60 80 100 time [s] 120 140 160 180 Fy [N] 1000 0 -1000 Fz [N] 500 0 -500 VMS Simulated Forces Maximum Fx [N] 1000 0 Minimum -1000 0 20 40 60 80 100 120 140 160 180 0 20 40 60 80 100 120 140 160 180 0 20 40 60 80 100 time [s] 120 140 160 180 Fy [N] 1000 0 -1000 Fz [N] 500 0 -500 Figure 3.25 : T63 - Operation 3, finishing of the walls at z = 105.26 [mm]. 72 Depth of Cut [mm] Chapter 3. Process Simulation Immersion Angle [deg] (a) Fx [N] Measured Forces 800 600 400 200 0 -200 35.94 35.96 35.98 36 36.02 36.04 36.06 35.98 36 36.02 36.04 36.06 36 36.02 time [s] 36.04 36.06 Fy [N] 1000 500 0 -500 35.96 Fz [N] 400 200 0 -200 Spindle Period tout 35.96 35.98 (b) Figure 3.26 : T63 - Operation 3, (a) engagement map at one point, (b) measured cutting forces for two spindle revolutions at and around the engagement map in (a) at z = 105.26 [mm]. 73 Chapter 3. Process Simulation Measured Forces Fx [N] 2000 0 -2000 0 20 40 60 80 100 120 140 160 180 0 20 40 60 80 100 120 140 160 180 0 20 40 60 80 100 time [s] 120 140 160 180 Fy [N] 2000 0 -2000 3000 Fz [N] 2000 1000 0 -1000 VMS Simulated Forces Fx [N] 2000 0 -2000 0 20 40 60 80 100 120 140 160 180 0 20 40 60 80 100 120 140 160 180 0 20 40 60 80 100 time [s] 120 140 160 180 Fy [N] 2000 0 -2000 3000 Fz [N] 2000 1000 0 -1000 Figure 3.27 : T63-Operation 1, roughing of the top part at z=163.5 [mm]. 74 Chapter 3. Process Simulation Measured Forces Fx [N] 2000 0 -2000 0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 time [s] 60 70 80 90 100 Fy [N] 2000 0 -2000 3000 Fz [N] 2000 1000 0 -1000 VMS Simulated Forces Fx [N] 2000 0 -2000 0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 time [s] 60 70 80 90 100 Fy [N] 2000 0 -2000 3000 Fz [N] 2000 1000 0 -1000 Figure 3.28 : T63 - Operation 1, roughing at z = 88.51 [mm]. 75 Chapter 3. Process Simulation Measured Forces Fx [N] 2000 1000 0 -1000 3 4 5 6 7 8 9 10 11 12 13 3 4 5 6 7 8 9 10 11 12 13 3 4 5 6 7 9 10 11 12 13 Fy [N] 2000 1000 0 -1000 Fz [N] 2000 1000 0 8 time [s] VMS Simulated Forces Fx [N] 2000 1000 0 -1000 5 10 15 20 25 10 15 20 25 10 15 time [s] 30 Fy [N] 2000 1000 0 -1000 5 Fz [N] 2000 1000 0 5 20 25 30 Figure 3.29 : T20 - Operation 2, semi-finishing of the top. 76 Chapter 3. Process Simulation Measured Forces Fx [N] 1000 0 -1000 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 time [s] 14 16 18 20 Fy [N] 2000 0 -2000 Fz [N] 2000 1000 0 VMS Simulated Forces Fx [N] 1000 0 -1000 265 270 275 280 285 290 265 270 275 280 285 290 265 270 275 280 time [s] 285 290 Fy [N] 2000 0 -2000 Fz [N] 2000 1000 0 Figure 3.30 : T20 - Operation 2, semi - finishing of door handle. 77 Chapter 4 Dynamics of Milling Operations 4.1. Introduction During a milling operation, cutting forces acting on the tool excite the tool - tool holder - spindle structure, leading to a variation in cut chip thickness, which in-turn changes cutting forces exciting the structure. In the literature, such self-excited vibrations during metal cutting are called regenerative chatter vibrations. If regeneration is not controlled, the cutting process may go to an unstable zone where exponentially increasing vibrations can potentially damage the tool, the workpiece and even the whole machine tool structure. Dynamics of milling processes has been studied by many researchers after pioneering works of Tlusty [62] and Tobias [63]. A considerable amount of work has focused on the stability of a milling process with an objective to predict a stability chart that would allow selection of chatter free cutting conditions. The dynamics of milling operations are developed and presented for 3-axis machining applications in this chapter. Stability of milling is analyzed with different methods including analytical model of chatter vibrations presented by Altintas and Budak [67], the Semi Discretization method proposed by Stepan et al [73], and the Multi Frequency solution by Budak and Altintas [67], and Merdol and Altintas [75]. Moreover, the analytical zero order solution is generalized to model stability of end mills with complex geometries. Finally, an alternative eigenvalue search method is proposed to accelerate iterative calculation of stability lobes when the Multi Frequency solution is used. 78 Chapter 4. Dynamics of Milling Operations 4.2. Dynamics of Milling A typical machine tool consists of the column, the housing, the spindle and the table. Tool holder is the physical interface between the cutting tool and the machine tool, and depending on the application, different tool holder - cutting tool assemblies are used in production. Since all of these structural components are physically connected to each other, together, they define a multidegree-of-freedom structural dynamic system. Some of the important sources of flexibilities are shown in Figure 4.1. In low speed machining, typically large and heavy components with low natural frequencies such as the spindle housing, the column and the fixture are excited; whereas high speed machining always excites the spindle, the tool holder and the tool that have high natural frequencies. Once excited, structural flexibilities cause the tool to vibrate, which undulates the chip thickness generated along the cutting flute, and the tool starts to leave a wavy surface behind as seen in Figure 4.2. This wavy surface is then removed by the successive flute, which also vibrates; therefore the difference between two surfaces becomes dynamically varying. Since prediction of cutting forces, torque, power, and vibrations are closely related to the instantaneous chip geometry, dynamic chip thickness distribution along the cutting edge has to be accurately modeled. Figure 4.1 : Main sources of flexibilities in milling. 79 Chapter 4. Dynamics of Milling Operations In general, the cutter geometry can be defined by a unit vector perpendicular to the cutter body (Figure 4.2) for flute j as: n j ( z ) = sin κ ( z ) ⋅ sin φ j ( z ) ⋅ i + sin κ ( z ) ⋅ cos φ j ( z ) ⋅ j – cos κ ( z ).k , (4.1) where i, j, and k are the unit vectors in the principal directions x, y and z, respectively. According to the regeneration phenomenon, the difference between the current vibrations and the vibrations that happened one tooth-pass period before is important. The vibration vector Δ ( t ) is defined as: Δ ( t ) = Δx ( t ) ⋅ i + Δy ( t ) ⋅ j + Δz ( t ) ⋅ k , (4.2) where Δx ( t ) = x ( t ) – x ( t – τ ) , Δy ( t ) = y ( t ) – y ( t – τ ) , Δz ( t ) = z ( t ) – z ( t – τ ) , and τ is the tooth period defined as: τ = 2π ⁄ ω T , (4.3) where ω T is the tooth passing frequency: ωT = N ⋅ Ω , (4.4) N is the number of flutes, and Ω [rad/s] is the rotational speed of the spindle. Vibration marks are directly imprinted on the chip because of the combined rigid and vibration motions of the cutter. The dynamic chip thickness can be obtained by taking the projection of vibrations on the unit vector perpendicular to the cutter geometry: h d, j ( z, t ) = n j ( z ) • Δ ( t ) . 80 (4.5) Chapter 4. Dynamics of Milling Operations y φst Present dynamic displacement φj φex [x(t), y(t)] x Previous dynamic displacement Cross section at elevation z [x(t-T), y(t-T)] z Dynamic Chip Chip(t,φ,κ) κ(z) Δ(t) nj(z) dz dFt, j dFr, j dFa, j dS x hd, j Figure 4.2 : Undulations on the chip, and dynamic chip thickness. For two different end mill geometries, the vibration vector, the unit geometry vector and the resulting dynamic chip thickness are presented in Figure 4.3. For a cylindrical cutter, vibration in the axial direction of the cutter (z-axis) does not contribute to dynamic chip thickness; whereas, for a ball end mill, three orthogonal vibrations influence chip thickness. Differential cutting forces in tangential, radial and axial directions can be described using the linear edge force model [27] as: ⎧ ⎪ dF r, j ( t, z ) ⎪ ⎨ dF t, j ( t, z ) ⎪ ⎪ dF a, j ( t, z ) ⎩ ⎧ ⎫ ⎪ K re ⎪ ⎪ ⎪ ⎬ = ⎨ K te ⎪ ⎪ ⎪ K ae ⎪ ⎩ ⎭ ⎫ ⎧ ⎪ ⎪ K rc ⎪ ⎪ ⎬ ⋅ dS ( z ) + ⎨ K tc ⎪ ⎪ ⎪ ⎪ K ac ⎭ ⎩ 81 ⎫ ⎪ ⎪ ⎬ ⋅ h d, j ( t, z ) ⋅ dS ( z ) , ⎪ ⎪ ⎭ (4.6) Chapter 4. Dynamics of Milling Operations where dS(z) is the differential contact length expressed as: dS ( z ) = dz ⁄ sin κ ( z ) , (4.7) where dz is the differential axial depth of cut, K rc, K tc, K ac are the cutting force coefficients due to the shear and K re, K te, K ae are the edge cutting force coefficients due to the rubbing of the tool flank with the workpiece, in radial, tangential and axial directions, respectively [86]. The helix angle is taken as zero since its effect on the stability of a milling operation [64] is negligible for most of the engagement conditions. Ball End Mill Cylindrical End Mill Cutter Type κ(z) nj(z) z x Dynamic Chip Thickness nj(z) x Δ = Δx . i + Δy . j + Δz . k Vibration Vector Unit Geometry Vector z nj(z) = sin φj(t,z) . i + cos φj(t,z) . j Δx . sin φj(t,z) + Δy . cos φj(t,z) nj(z) = sin κ(z) . sin φj(t,z) . i + sin κ(z) . cos φj(t,z) . j - cos κ(z) . k Δx . sin κ(z) . sin φj(t,z) + Δy . sin κ(z) . cos φj(t,z) - Δz . cos κ(z) Figure 4.3 : Dynamic chip thickness definition for general end mills. Cutting forces acting on flute j are then calculated as: ⎧ ⎫ ⎧ ⎪ F r, j ( t ) ⎪ ⎪ K re ⎪ ⎪ ⎪ ⎨ F t, j ( t ) ⎬ = ⎨ K te ⎪ ⎪ ⎪ ⎪ F a, j ( t ) ⎪ ⎪ K ae ⎩ ⎭ ⎩ ⎫ ⎧ ⎪ ⎪ K rc ⎪ a ⎪ - + ⎨ K tc ⎬ ⋅ ---------⎪ sin κ ⎪ ⎪ ⎪ K ac ⎭ ⎩ 82 ⎫ ⎪ ⎪ a, ⎬ ⋅ h d, j ( t ) ⋅ ---------sin κ ⎪ ⎪ ⎭ (4.8) Chapter 4. Dynamics of Milling Operations where a is the axial depth of cut, κ is the average axial immersion angle, and h d, j ( t ) is the average dynamic chip thickness obtained from Eq (4.5) after replacing κ with κ in Eq (4.1). For complex cutters such as ball end and bull nose, axial immersion angle varies along the z-axis; therefore, it can be assumed to be acting in the middle of the maximum axial depth of cut [96]. For instance, a ball end mill has an average axial immersion angle of κ = π ⁄ 4 . Note that cutting forces become independent of the z variable with the introduction of the average axial immersion angle. Cutting forces in rotating coordinates are projected onto Cartesian coordinates through a geometric transformation matrix: ⎧ ⎫ ⎪ Fx ( t ) ⎪ – sin κ sin φ j ( t ) ⎪ ⎪ ⎨ F y ( t ) ⎬ = – sin κ cos φ j ( t ) ⎪ ⎪ cos κ ⎪ Fz ( t ) ⎪ ⎩ ⎭ – cos φ j ( t ) sin φ j ( t ) 0 ⎧ – cos κ sin φ j ( t ) ⎪ F r ( t ) ⎪ – cos κ cos φ j ( t ) ⎨ F t ( t ) ⎪ – sin κ ⎪ Fa ( t ) ⎩ ⎫ ⎪ ⎪ ⎬ , ⎪ ⎪ ⎭ (4.9) where total rotating cutting forces are calculated by adding the contribution of each flute as: ⎧ ⎫ ⎪ Fr ( t ) ⎪ ⎪ ⎪ ⎨ Ft ( t ) ⎬ = ⎪ ⎪ ⎪ Fa ( t ) ⎪ ⎩ ⎭ ⎧ ⎫ ⎪ F r, j ( t ) ⎪ ⎪ ⎪ ( t ) F ⎨ ∑ ⎪ t, j ⎬⎪ . j = 1 ⎪ F (t) ⎪ a, j ⎩ ⎭ N (4.10) Substituting Eq (4.8) into Eq (4.10), and then into Eq (4.9), relation between cutting forces and vibrations is derived as: 83 Chapter 4. Dynamics of Milling Operations ⎧ ⎫ ⎪ Fx ( t ) ⎪ a xx ( t ) a xy ( t ) a xz ( t ) ⎧ Δx ( t ) ⎫ ⎪ ⎪ ⎪ ⎪ 1--⎨ F y ( t ) ⎬ = 2 ⋅ a ⋅ K tc ⋅ a yx ( t ) a yy ( t ) a yz ( t ) ⎨ Δy ( t ) ⎬ , ⎪ ⎪ ⎪ ⎪ a zx ( t ) a zy ( t ) a zz ( t ) ⎩ Δz ( t ) ⎭ ⎪ Fz ( t ) ⎪ ⎩ ⎭ (4.11) where time-varying directional coefficients in each principal direction are obtained as: ∑ axx, j = ∑ –gj ⋅ [ sin 2φj + ( 1 – cos 2φj ) ( Kr sin κ + Ka cos κ ) ] a xx ( t ) = j=1 j=1 N N ∑ axy, j = ∑ –gj ⋅ [ ( 1 + cos 2φj ) + sin 2φj ( Kr sin κ + Ka cos κ ) ] a xy ( t ) = j=1 j=1 N N j 1 N a yy ( t ) = a yz ( t ) = (4.12) ∑ axz, j = ∑ 2gj ⋅ [ cos φj ⁄ tan κ + sin φj cos κ ( Kr + Ka ⁄ tan κ ) ] a xz ( t ) = a yx ( t ) = , j 1 N ∑ ayx, j = ∑ gj ⋅ [ ( 1 – cos 2φj ) – sin 2φj ( Kr . sin κ + Ka . cos κ ) ] j=1 j=1 N N ∑ a yy, j = ∑ j=1 j=1 N N g j ⋅ [ sin 2φ j – ( 1 + cos 2φ j ). ( K r . sin κ + K a . cos κ ) ] ∑ ayz, j = ∑ 2gj ⋅ [ – sin φj ⁄ tan κ + cos φj . cos κ . ( Kr + Ka ⁄ tan κ ) ] 84 , (4.13) Chapter 4. Dynamics of Milling Operations a zx ( t ) = a zy ( t ) = ∑ azx, j = ∑ 2gj ⋅ [ sin φj . cos κ. Kr – sin φj . sin κ.Ka ] j=1 j=1 N N ∑ azy, j = ∑ 2gj . [ cos φj . cos κ. Kr – cos φj . sin κ.Ka ] j=1 , (4.14) j=1 N N ∑ azz, j = ∑ 2gj . [ – cos κ ⁄ tan κ .Kr + cos κ.Ka ] a zz ( t ) = j 1 j 1 where g j = g ( φ j ) is the switching function which is equal to unity while cutting and zero otherwise, K r = K rc ⁄ K tc and K a = K ac ⁄ K tc . One important observation about directional coefficients is that they significantly change the direction of excitation as the tool rotates and this constitutes the fundamental difference between the dynamics of turning and milling. Eq (4.11) can be expressed in the time domain as: 1 f ( t ) = --- ⋅ a ⋅ K tc ⋅ A ( t ) ⋅ Δ ( t ), 2 (4.15) T where the vibration difference vector is Δ ( t ) = { Δx, Δy, Δz } , the dynamic cutting forces T vector is f ( t ) = { F x, F y, F z } and the directional coefficient matrix is N A(t) = ∑ Aj ( t ) j=1 N = ∑ j=1 85 a xx, j a xy, j a xz, j a yx, j a yy, j a yz, j . a zx, j a zy, j a zz, j (4.16) Chapter 4. Dynamics of Milling Operations 4.3. Stability of Milling xn fxn x2 fx2 x1 fx1 . . . Figure 4.4 : Multi degree-of-freedom milling dynamics. In order to achieve high material removal rates in milling, a careful dynamic stability analysis has to be carried out and cutting parameters must be selected accordingly to avoid unstable (chatter) vibrations. Eq (4.15) describes how dynamic cutting forces in milling are generated. A dynamical system can be defined by expressing the equation of motion at the discrete points along the cutting tool and applying cutting forces as external excitations: M ⋅ x·· ( t ) + C ⋅ x· ( t ) + K ⋅ x ( t ) = f x , T (4.17) where x ( t ) = { x 1, x 2, …, x n } represent physical degrees of freedom, f x = { f x1, f x2, …, f xn } T is the excitation force vector, and n is the number of degrees of freedom. M, C, K are n-dimen- 86 Chapter 4. Dynamics of Milling Operations sional mass, damping and stiffness matrices, respectively. Physical and modal displacements are related by mode shapes as: x = Px . μ , (4.18) T where μ = { μ 1, μ 2, …, μ mx } are modal displacements corresponding to each mode, and m x is the number of modes, respectively. The mode shape matrix (n-by- m x ) is defined as: Px = p x, 11 p x, 12 . . p x, 1mx p x, 21 p x, 22 . . p x, 2mx . . . . . . . . . . . . p x, n1 p x, n2 . (4.19) p x, nmx Each column of the mode shape matrix ( P x ) represents a ratio of structural displacement amplitudes at each node in response to each mode, i.e. eigenvectors. When Eq (4.18) is substituted into Eq (4.17), the equation of motion in modal coordinates is obtained as: ·· ( t ) + C ⋅ μ· ( t ) + K ⋅ μ ( t ) = f , Mμ ⋅ μ μ μ xμ T T (4.20) where M μ = P x .M.P x , and K μ = P x .K.P x are modal mass and stiffness matrices, and f xμ is the modal force vector: 87 Chapter 4. Dynamics of Milling Operations ⎧ p x, 11 f x1 + p x, 21 f x2 + … + p x, n1 f xn ⎪ ⎪ p x, 12 f x1 + P x, 22 f x2 + … + p x, n2 f nx T f xμ = P x f x = ⎨ … ⎪ ⎪ ⎩ p x, 1mx f x1 + p x, 2mx f 2x + … + p x, nmx f xn ⎫ ⎪ ⎪ ⎬. ⎪ ⎪ ⎭ (4.21) When the system has a proportional damping, the modal damping matrix is expressed as T C μ = P x .C.P x . A similar analysis can be performed for the rest of the orthogonal directions (y and z): T ·· ( t ) + C ⋅ η· ( t ) + K ⋅ η ( t ) = f , Mη ⋅ η η η yη (4.22) M ν ⋅ ν·· ( t ) + C ν ⋅ ν· ( t ) + K ν ⋅ ν ( t ) = f zν , (4.23) T where f yη = P y f y and f zν = P z f z are modal force vectors of length m y and m z ; m y and m z are the number of modes; P y and P z are mode shape matrices of size (n-by- m y ) and (n-by- m z ); and f y and f z are excitation force vectors in y and z directions, respectively. Using mode shapes, physical and modal displacements are related in y and z directions: y = Py ⋅ η , (4.24) z = Pz ⋅ ν , (4.25) where η and ν are modal displacement vectors in y and z directions, respectively. Eqs (4.35)(4.23) can be assembled into a single matrix equation of the form: 88 Chapter 4. Dynamics of Milling Operations ·· ( t ) + C ⋅ Γ· ( t ) + K ⋅ Γ ( t ) = T ⋅ f ( t ) , MΓ ⋅ Γ Γ Γ full all (4.26) where T full is the transformation matrix between physical and modal forces and defined as: Px T full = T 0 0 Py 0 0 T , 0 0 Pz (4.27) T mt × mt T ⎧ ⎫ f all ( t ) = ⎨ f x ( t ) f y ( t ) f z ( t ) ⎬ , ⎩ ⎭ (4.28) T Γ ( t ) = { μ 1, μ 2, …, μ mx, η 1, η 2, …, η my, ν 1, ν 2, …, ν mz } , Mμ 0 MΓ = 0 , 0 Mη 0 0 (4.29) 0 Mν (4.30) mt × mt where m t = m x + m y + m z is the total number of modes, and C Γ and K Γ are defined similar to Eq (4.30). Without loss of generality, dynamics can be analyzed at the tool tip. Let all modes be normalized to unity at the tool tip (Point #1): ⎫ p x, 11 = p x, 12 = … = p x, 1mx = 1 ⎪ ⎪ p y, 11 = p y, 12 = … = p y, 1my = 1 ⎬ ; ⎪ p z, 11 = p z, 12 = … = p z, 1mz = 1 ⎪ ⎭ 89 (4.31) Chapter 4. Dynamics of Milling Operations and forces be applied only at the tool tip: ⎫ f x1 = F x ( t ), f x2 = 0, …, f xn = 0 ⎪ ⎪ f y1 = F y ( t ), f 2y = 0, …, f yn = 0 ⎬ . ⎪ f z1 = F z ( t ), f z2 = 0, …, f zn = 0 ⎪ ⎭ (4.32) The equation of motion now reduces to: ·· ( t ) + C ⋅ Γ· ( t ) + K ⋅ Γ ( t ) = T ⋅ f ( t ) , MΓ ⋅ Γ Γ Γ (4.33) where T is the reduced transformation matrix: T = { 1 } mx 0 0 0 { 1 } my 0 0 0 { 1 }mz , (4.34) ⎧ ⎫ f ( t ) = ⎨ Fx ( t ) Fy ( t ) Fz ( t ) ⎬ . ⎩ ⎭ (4.35) mt × 3 and f(t) is the cutting force vector: T Note that varying dynamics along the depth of cut can be considered by simply considering the distributed forces ( f xn, f yn, f zn ≠ 0 ) and displacements along the cutter axis. Cutting forces in milling are previously defined in Eq (4.15), and they can be re-expressed using modal coordinates and the transformation matrix T as: 90 Chapter 4. Dynamics of Milling Operations 1 T f ( t ) = --- ⋅ a ⋅ K tc ⋅ A ( t ) ⋅ T ⋅ { Γ ( t ) – Γ ( t – τ ) } . 2 (4.36) Substituting Eq (4.36) into Eq (4.33), and reorganizing the terms: ·· ( t ) + C ⋅ Γ· ( t ) + K ⋅ Γ ( t ) = T ⋅ 1--- ⋅ a ⋅ K ⋅ A ( t ) ⋅ T T ⋅ { Γ ( t ) – Γ ( t – τ ) } , MΓ ⋅ Γ Γ Γ tc 2 ·· ( t ) + C ⋅ Γ· ( t ) + ( K – L ( t ) ) ⋅ Γ ( t ) = – L ( t ) ⋅ Γ ( t – τ ) , MΓ ⋅ Γ Γ Γ Γ Γ (4.37) 1 T L Γ ( t ) = --- ⋅ a ⋅ K tc ⋅ T ⋅ A ( t ) ⋅ T , 2 (4.38) where the general equation of motion in Eq (4.37) can be expressed in a state space form as: · Θ(t) = U(t) ⋅ Θ(t) + V(t) ⋅ Θ(t – τ) , (4.39) ⎧ Γ( t) ⎫ Θ(t) = ⎨ ⎬, ⎩ Γ· ( t ) ⎭ (4.40) where U(t) = 0 I –1 MΓ ( LΓ ( t ) – KΓ ) 91 –1 –MΓ CΓ , (4.41) Chapter 4. Dynamics of Milling Operations V(t) = 0 –1 0 –MΓ LΓ ( t ) 0 . (4.42) Modal parameters (mass, damping and stiffness) of the structure are identified from a set of impact tests by attaching an accelerometer at the tool tip and tapping it with an instrumented hammer in the direction of interest. Measured acceleration and force data are then used to obtain the frequency response function (FRF), from which modal parameters are extracted using a modal analysis software [93]. 4.3.1. Time Domain Based Stability Solution The equation of motion in milling - represented by Eq (4.39) - is a time periodic delay differential equation (DDE). The regenerative effect due to the wavy surface on both sides of the chip appears as the delay term; whereas, rotating cutting forces represented by directional coefficients introduces periodicity into the equation of motion. In milling, the delay and the periodicity of the system are both equal to the tooth period ( τ ). Stability of linear DDE depends on infinite number of characteristic roots making stability properties of delayed systems extremely complex [73]. Stepan et al. [73] proposed a semi-analytical solution called the Semi-Discretization (SD) method for the stability of linear delayed systems. The SD method in reference [73] was applied to two degree-of-freedom milling with a single mode in each direction. In this section, the SD method is applied to investigate stability characteristics of dynamics of generalized milling applications. General form of three degree-of-freedom milling dynamics with multiple modes is derived in the the previous section and the resulting state space form is expressed in Eq (4.39). 92 Chapter 4. Dynamics of Milling Operations The idea behind the SD method is to approximate the time delayed and time-dependant terms with piecewise constant functions so that the DDE can be reduced into series of ordinary differential equations (ODEs). Since the actual time domain terms are not discretized and left unchanged, an exact solution for each ODE can be obtained. The SD starts with construction of the time interval division as shown in Figure 4.5. Delay time (tooth period) is divided into finite number of elements of length Δt such that Δt = T ⁄ k , (4.43) where k is an integer and the time at each interval is expressed as: t i = i ⋅ Δt , i = 0, 1, …, k . (4.44) Note that time step or size of a discrete element must be set much smaller than the period of the highest structural mode so that dynamics can be properly captured within the solution. In this study, time step is set equal to or smaller than one tenth of the period of the largest mode. Once the discretization is completed, the time-dependant matrices, which contain dynamic parameters and directional coefficients, are approximated with constant matrices for each discretization interval [t i, t i + 1 ) : ti + 1 1 U i = ------Δt ∫ ti + 1 1 U ( t ) dt , V i = ------Δt ti ∫ ti 93 V ( t ) dt , (4.45) Chapter 4. Dynamics of Milling Operations and the delayed term is approximated as a weighted linear combination of the delayed discrete values Θ i – k and Θ i – k + 1 : Θ ( t – τ ) ≈ Θ ( t i – kΔt ) ≈ w b ⋅ Θ i – k + w a ⋅ Θ i – k + 1 , (4.46) where w a = w b = 1 ⁄ 2 when the time delay is equal to the time period of the system, which is always the case for a milling system. When both approximations, Eqs (4.45)-(4.46), are substituted into non-autonomous DDE expressed in Eq (4.39), the system becomes a linear autonomous ODE within a time interval: · Θ ( t ) = Ui ⋅ Θ ( t ) + Vi ⋅ ( wb ⋅ Θi – k + wa ⋅ Θi – k + 1 ) , ( ti ≤ t < ti + 1 ) . Approximation of the delayed term Q T=k*Δt t ti-k ti-k+1 ti ti+1 Discretization of the delay (tooth period) 0 1 2 Counter k 0 Δt 2Δt Time kΔt Figure 4.5 : Approximation of the delayed term using the semi-discretization (SD) by Stepan and Insperger [73]. 94 (4.47) Chapter 4. Dynamics of Milling Operations Being treated as an initial value problem, i.e. Θ ( t i ) = Θ i , Eq (4.47) has an exact solution of the form: Θ(t) = e Ui ( t – ti ) ⋅ Θi + ( e Ui ( t – t i ) – I ) ⋅ Ui –1 ⋅ Vi ⋅ ( wb ⋅ Θi – k + wa ⋅ Θi – k + 1 ) . (4.48) Substituting t = t i + 1 in Eq (4.48), first, the state at the exit of the ith time interval is found as: Θi + 1 = Si ⋅ Θi + Ri ⋅ ( wb ⋅ Θi – k + wa ⋅ Θi – k + 1 ) , ( ti ≤ t < ti + 1 ) , (4.49) and then a discrete map is defined as: full full ϒi + 1 = Ψi full ϒi , (4.50) where Si = e U i Δt Ri = ( Si – I ) ⋅ Ui are ( 2m t -by- 2m t ) matrices; 95 , –1 (4.51) ⋅ Vi , (4.52) Chapter 4. Dynamics of Milling Operations Si 0 0 … 0 wa ⋅ Ri wb ⋅ Ri Ψi full = I 0 0 … 0 0 I 0 … 0 0 0 0 0 . . . . . . . . . . . . 0 0 0 … 0 0 0 0 0 I 0 0 0 0 0 0 0 I 0 0 0 . . . . . . . . . (4.53) is a [ 2m t ( k + 1 ) -by- 2m t ( k + 1 ) ] matrix; and T ϒi full ⎧ ⎫ · · = ⎨ Γ Γ· i Γ Γ … Γ Γ i–1 i–k ⎬ i–1 i–k ⎩ i ⎭ (4.54) is a vector of length ( 2m t ( k + 1 ) ). It is clear from Eq (4.42) that Γ· ( t – τ ) does not play any role in Eq (4.39); therefore, Θ i + 1 in Eq (4.49) does not depend on Γ· i – k or Γ· i – k + 1 . Using this observation, size of the problem can be reduced and a new discrete map is defined as: ϒi + 1 = Ψi ⋅ ϒi , where Ψ i is a square matrix of size m t .(k+2): 96 (4.55) Chapter 4. Dynamics of Milling Operations S i ( 1…m t , 1…m t ) S i ( m +1…2m t Ψi = t , 1…m t ) S i ( 1…m , t m t +1…2m t ) S i ( m +1…2m t t , m t +1…2m t ) 0 0 … … w a R i ( 1…m t , 1…m t ) 0 0 … … w a R i ( m +1…2m t … … … … I 0 … … 0 I 0 … I 0 0 0 0 0 . . . . . . . . . 0 0 0 0 0 … . . . . . . . . . t , 1…m t ) w b R i ( 1…m t , 1…m t ) w b R i ( m +1…2m t 0 0 0 0 0 0 . . . . . . I 0 t , 1…m t ) (4.56) and ϒ i is a vector of length m t .(k+2): T ⎧ ⎫ ϒ i = ⎨ Γ Γ· i Γ . … Γ Γ i–1 i–k–1 i–k ⎬ ⎩ i ⎭ (4.57) Note that additional subscripts of S i and R i in Eq (4.56) are used to indicate sub matrices defined within the specified row and column ranges. The transition matrix over one period can now be defined by coupling all the discrete maps: Φ = Ψ k – 1 ⋅ Ψ k – 2 …Ψ 0 , (4.58) where ϒ k = Φ ⋅ ϒ 0 . According to the Floquet theory, the stability of a linear periodic system depends on the characteristic multipliers (eigenvalues) of the so-called principal or transition matrix and if all eigenvalues of Φ are in modulus less than one, then the system is asymptotically stable [73]. 4.3.2. Frequency Domain Based Stability Solutions An alternative method to the time domain based solution of chatter stability is a frequency domain based method. Similar to cutting forces, the directional coefficient matrix, A ( t ) , is also periodic 97 Chapter 4. Dynamics of Milling Operations at the tooth passing frequency ( ω T ) or at the tooth period ( τ ), i.e., A ( t ) = A ( t + τ ) . Harmonic functions can be expressed by Fourier series [67] as: ∞ ∑ A(t) = Ar ⋅ e irω T t , (4.59) r = –∞ where i is the imaginary unit, and the Fourier coefficients are defined as: τ 1 A r = --τ ∫ A(t) ⋅ e – irω T t dt . (4.60) 0 In order to obtain Fourier coefficients, first, Eq (4.16) is substituted into Eq (4.60): τ 1 A r = --τ ∫ 0 ⎛ N ⎞ ⎜ ⎟ ⋅ e – irω T t dt = --1A ( t ) j ∑ ⎜ ⎟ τ ⎝j = 1 ⎠ τ ∫ 0 ⎛ N ⎞ ⎜ ⎟ –irNΩt dt . A ( t ) ⎜∑ j ⎟ ⋅e ⎝j = 1 ⎠ (4.61) Since φ j = φ + ( j – 1 )φ p = Ω ( t + ( j – 1 ) ⋅ τ ) , where φ p = 2π ⁄ N is the pitch angle for an even-pitch cutter, the summation in the above relation can be taken out by changing integration limits: N 1 A r = --- ⋅ τ ∑ j=1 ⎛ ⎜ ⎜ ⎝ ( j + 1 )τ ∫ Aj ( γ ) ⋅ e jτ 98 – irNΩγ ⎞ dγ ⎟ , ⎟ ⎠ (4.62) Chapter 4. Dynamics of Milling Operations where γ is still a time variable like t; however, acting between different limits. Using the relation between the rotation angle, and the rotational speed of the tool, i.e., θ = Ωτ , integration can be redefined as: N 1 A r = ------- ⋅ τΩ ∑ j=1 ⎛ ⎜ ⎜ ⎝ ( j + 1 )φ p ∫ Aj ( θ ) ⋅ e – irNθ jφ p ⎞ dθ ⎟⎟ . ⎠ (4.63) Since A j ( θ ) contains the switching function g j , the above expression can be simplified by taking all flutes into account within the immersion zone: φ ex N A r = -----2π ∫ A(θ) ⋅ e – irNθ dθ , (4.64) φ st where A ( θ ) is defined equal to A j ( t ) in Eq (4.16) without the dependency on flute number j, i.e., g j being unity and φ j = θ . The right-hand-side of Eq (4.15) is a function of not only the tooth passing frequencies due to directional coefficients A ( t ) , but also the vibration (chatter) frequency due to vibration terms, Δ ( t ) . Thus, as a reaction, dynamic milling forces ( f ( t ) ) should contain both the tooth passing frequencies ( k ⋅ ω T ) and the chatter frequency, ω c when the system is marginally stable. In general, Floquet’s theorem [66][98][99] states that for a second order differential equation with periodic and piece-wise continuous coefficients, the solution has the following form: f(t) = e sλ t ⋅ p(t) , 99 (4.65) Chapter 4. Dynamics of Milling Operations where p ( t ) is a periodic function of the tooth passing frequency ( ω T ). It is evident from Eq (4.65) that the milling system is stable if the real part of exponent s λ is negative. Since the cutter and the workpiece vibrate at the chatter frequency on the stability border, s λ can be replaced by s λ = iω c to ensure marginal stability [67]. When periodic function p ( t ) is also expanded by using Fourier series, the following expression is obtained for the dynamic milling forces: ∞ f(t) = e iω c t ∑ ∞ pk e ikω T t = k = –∞ ∑ pk e i ( ω c + kω T )t . (4.66) k = –∞ Spindle is a multi degree of freedom dynamical system whose dynamics is reflected at the tool tip during machining process. Considering each principal direction, the displacement response to an acting force f q ( t ) can be expressed as: q ( D ) = G ( D ) ⋅ f q ( t ) , (q = x,y,z), (4.67) where D is the differential operator d ( … ) ⁄ dt , G(D) is the transfer function expressed in terms of modal parameters such as natural frequency ( ω n ), modal mass (m) and damping ratio ( ζ ): M G(D) = 1 ⁄ mj -, ∑ --------------------------------------------------------------2 2 D +2⋅ζ ⋅ω ⋅D+ω j=1 j n, j (4.68) n, j where M is the total number of modes. Considering all directions with direct (diagonal elements) and cross (off-diagonal elements) transfer functions, a matrix form of displacements can be written as: 100 Chapter 4. Dynamics of Milling Operations ⎧ x(t) ⎫ ⎪ ⎪ ⎨ y(t) ⎬ = G(D) ⋅ f(t) , ⎪ ⎪ ⎩ z(t) ⎭ (4.69) where G ( D ) is the overall transfer function matrix: G xx ( D ) G xy ( D ) G xz ( D ) G ( D ) = G yx ( D ) G yy ( D ) G yz ( D ) . (4.70) G zx ( D ) G zy ( D ) G zz ( D ) Note that if there is no coupling between orthogonal axes, the dynamic displacement for each axis is obtained by multiplying the transfer function of each axis with the milling force, i.e., the transfer function matrix G ( D ) becomes a diagonal matrix. The vibration difference vector is expressed by using the transport delay (tooth period, τ ) as: Δ(t) = (1 – e – Dτ ⎧ x( t) ⎫ ⎪ ⎪ ) ⋅ ⎨ y( t) ⎬ . ⎪ ⎪ ⎩ z(t) ⎭ (4.71) Substituting Eqs (4.66), (4.69), and (4.71) into Eq (4.15) yields ∞ 1 – Dτ f ( t ) = --- ⋅ a ⋅ K tc ⋅ A ( t ) ⋅ ( 1 – e ) 2 ∑ G ( D ) ⋅ pk ⋅ e i ( ω c + kω T )t . (4.72) k = –∞ Using the definition of parametric transfer functions of linear periodic systems given in reference [100], Minis and Yanushevsky [66] expressed Eq (4.72) as: 101 Chapter 4. Dynamics of Milling Operations ∞ 1 f ( t ) = --- ⋅ a ⋅ K tc ⋅ A ( t ) ⋅ 2 ∑ e i ( ω c + kω T )t ⋅ (1 – e – i ( ω c + kω T )τ ) ⋅ G ( i ( ω c + kω T ) ) ⋅ p k . (4.73) k = –∞ Eq (4.65) is substituted into the left-hand-side of Eq (4.73), and both sides are divided by e iω c t to obtain: ∞ – iω T 1 p ( t ) = --- ⋅ a ⋅ K tc ⋅ ( 1 – e c ) 2 where e – i ( ω c + kω T )τ equation by 1 ⁄ τ ⋅ e is replaced by e – ipω T t ∑ e ikω T t ⋅ A ( t ) ⋅ G ( i ( ω c + kω T ) ) ⋅ p k , (4.74) k = –∞ – iω c τ as ω T ⋅ τ = 2π . Multiplying both sides of the above , integrating from 0 to τ , and using the orthogonality principle [67], the following is obtained: ∞ τ --1τ ∫ p ( t )e – ipω T t – iω τ 1 dt = --- a ⋅ K tc ⋅ ( 1 – e c ) 2 0 ⎛ 1 ∑ ⎜⎜ --τk = –∞ ⎝ τ ∫ 0 A(t) ⋅ e – i ( p – k )ω T t ⎞ dt⎟ ⋅ G ( i ( ω c + kω T ) ) ⋅ p k ⎟ ⎠ which can be simplified to ∞ – iω τ 1 p p = --- a ⋅ K tc ⋅ ( 1 – e c ) 2 ∑ A p – k ⋅ G ( iω k ) ⋅ p k , ( p = 0, ± 1, ± 2, ... ) , (4.75) k = –∞ where ω k = ω c + k ⋅ ω T . The infinite limits of the summation must be truncated in order to solve the system of equations. If the maximum number of harmonics is limited to h r , and the summation in Eq (4.75) is expressed in matrix form, then the problem is reduced to an eigenvalue problem: 102 Chapter 4. Dynamics of Milling Operations ( I – λ ⋅ Φ ( ω c, ω T ) ) ⋅ p = 0 , where λ = 0.5 ⋅ a ⋅ K tc ⋅ ( 1 – e – iω c τ (4.76) ) is the eigenvalue, p is the eigenvector defined as: T ⎧ ⎫ p = ⎨ p 0 p 1 p –1 … p h p –h ⎬ , r r ⎩ ⎭ (4.77) Φ is a 3(2 h r +1)-by-3(2 h r +1) matrix called oriented transfer function: ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ ⎛ k = 0 1 –1 . . h –h ⎞ r r ⎠ ⎝ ( 0, 0 ) Go ( 1, 0 ) Φ ( ω c, ω T ) = Go ( – 1, 0 ) Go ( – 1, 1 ) Go ( 0, 1 ) Go ( – 2, 1 ) Go ( 1, – 1 ) Go ( 2, – 1 ) Go ( 0, – 1 ) Go . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Go ( ( p – k ), k ) ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ p= 0 1 –1 , . . hr (4.78) –hr where ( ( p – k ), k ) Go = A p – k ⋅ G ( iω k ) . (4.79) Directional coefficients for the multi frequency case are obtained by carrying out the integration in Eq (4.64): 103 Chapter 4. Dynamics of Milling Operations (r) (r) (r) α xx α xy α xz N A r = ------ α ( r ) α ( r ) α ( r ) , 2π yx yy yz (r) (r) (4.80) (r) α zx α zy α zz where matrix elements are φ ex (r) α xx = [ – c 1 ( θ ) – ( c 0 ( θ ) – c 2 ( θ ) ) ( K r ⋅ sin κ + K a ⋅ cos κ ) ] (r) α xy = [ – c 0 ( θ ) – c 2 ( θ ) – c 1 ( θ ) ( K r ⋅ sin κ + K a ⋅ cos κ ) ] φ ex φ st φ ex (r) α xz = 2 [ c 4 ( θ ) ⁄ tan κ + c 3 ( θ ) ⋅ cos κ ⋅ ( K r + K a ⁄ tan κ ) ] (r) α yx = [ c 0 ( θ ) – c 2 ( θ ) – c 1 ( θ ) ⋅ ( K r ⋅ sin κ + K a ⋅ cos κ ) ] φ st φ st φ ex φ st (r) α yy = [ c 1 ( θ ) – ( ( c 0 ( θ ) + c 2 ( θ ) ) ⋅ ( K r ⋅ sin κ + K a ⋅ cos κ ) ) ] (r) α yz = 2 [ – c 3 ( θ ) ⁄ tan κ + c 4 ( θ ) ⋅ cos κ ⋅ ( K r + K a ⁄ tan κ ) ] (r) φ ex (r) φ ex α zx = 2 [ c 3 ( θ ) ( K r ⋅ cos κ – K a ⋅ sin κ ) ] α zy = 2 [ c 4 ( θ ) ( K r ⋅ cos κ – K a ⋅ sin κ ) ] (r) φ st φ st α zz = 2 [ c 0 ( θ ) ⋅ cos κ ⋅ ( – K r / tan κ + K a ) ] 104 φ ex φ st φ ex φ st φ ex φ st ⎫ ⎪ ⎪ ⎪ ⎬, ⎪ ⎪ ⎪ ⎭ ⎫ ⎪ ⎪ ⎪ ⎬, ⎪ ⎪ ⎪ ⎭ ⎫ ⎪ ⎪ ⎪ ⎬, ⎪ ⎪ ⎪ ⎭ (4.81) (4.82) (4.83) Chapter 4. Dynamics of Milling Operations ⎧ ⎪ θ if r = 0 c0 ( θ ) = ⎨ – irNθ ⎪ – i ⁄ ( rN ) ⋅ e otherwise ⎩ 2 2 c 1 ( θ ) = ( ( 2 cos 2θ + irN sin 2θ ) ⁄ ( r N – 4 ) ) ⋅ e 2 – irNθ 2 c 2 ( θ ) = ( ( – 2 sin 2θ + irN cos 2θ ) ⁄ ( r N – 4 ) ) ⋅ e 2 2 c 3 ( θ ) = ( ( cos θ + irN sin θ ) ⁄ ( r N – 1 ) ) ⋅ e 2 2 – irNθ – irNθ c 4 ( θ ) = ( ( – sin θ + irN cos θ ) ⁄ ( r N – 1 ) ) ⋅ e – irNθ ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬. ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ (4.84) Note that some coefficients in Eq (4.84) are undefined when rN = ± 1 , ± 2 , therefore, undefined elements of the directional coefficient matrix, A r , in Eq (4.80) are separately calculated as • when rN = 2, (r>0): (r) α xx 1 – i4θ – i2θ – i4θ = --- [ i ( K r ⋅ sin κ + K a ⋅ cos κ ) ( e – 4e – i4θ + 1 ) + ( e + i4θ – 1 ) ] 8 1 (r) – i4θ – i4θ – i2θ α xy = --- [ ( K r ⋅ sin κ + K a ⋅ cos κ ) ( e + i4θ – 1 ) + ( – ie – i 4e – 4θ + i ) ] 8 φ ex φ st φ ex φ st φ ex (r) α yx 1 – i4θ – i4θ – i2θ = --- [ ( K r ⋅ sin κ + K a ⋅ cos κ ) ( e + i4θ – 1 ) + ( – ie + i4e – 4θ + i ) ] 8 (r) α yy 1 – i4θ – i2θ – i4θ = --- [ i ( K r ⋅ sin κ + K a ⋅ cos κ ) ( – e – 4e + i4θ + 1 ) + ( – e – i4θ + 1 ) ] 8 • when rN = 1, (r>0): 105 φ st φ ex φ st ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ , (4.85) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ Chapter 4. Dynamics of Milling Operations (r) α xz = ( – cosκ ⁄ 2 ) [ ( e (r) α yz = ( cosκ ⁄ 2 ) [ i ( e – i2θ – i2θ + i2θ – 1 ) ( K r + K a ⁄ tan κ ) – ( i ⁄ sin κ ) ( e – i 2θ ) ( K r + K a ⁄ tan κ ) + ( 1 ⁄ sin κ ) ( e (r) α zx = { – ( K r ⋅ cos κ – K a ⋅ sin κ ) ⋅ [ iθ – sin2 θ – 0.5 sin 2θ ] } (r) α zy = { – ( K r ⋅ cos κ – K a ⋅ sin κ ) ⋅ [ i sin2 θ – 0.5 sin 2θ – θ ] } – i2θ – i2θ – i 2θ – 1 ) ] + i2θ – 1 ) ] φ ex φ st φ ex φ st φ ex φ st φ ex φ st ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ . (4.86) ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ Once positive coefficients (r>0) are calculated, negative coefficients (r<0) simply become equal to complex conjugates of the positive ones. Finally, a non-trivial solution can be found when the determinant vanishes: det ( I – λ ⋅ Φ ( ω c, ω T ) ) = 0 , (4.87) where ω c is the vibration frequency, and ω T is the tooth passing frequency. Eq (4.87) describes the characteristic equation of the closed loop milling system. When the oriented transfer function Φ ( ω c, ω T ) is known, all eigenvalues can be determined, from which critical depths of cut are obtained as: 2λ q a lim, q = --------------------------------------= a R, q + i ⋅ a I, q , – iω T K tc ⋅ ( 1 – e c ) (4.88) where λ q = λ R, q + i ⋅ λ I, q , q = 1…3 ( 2h r + 1 ) , i is the imaginary unit and subscripts R and I represent real and imaginary parts, respectively. Substituting e – iω c τ = cos ( ω c τ ) – i ⋅ sin ( ω c τ ) into Eq (4.88), real and imaginary parts of depth of cut are expressed as: 106 Chapter 4. Dynamics of Milling Operations 1 λ R, q ⋅ ( 1 – cos ( ω c τ ) ) + λ I, q sin ( ω c τ ) a R, q = ------- ------------------------------------------------------------------------------------------- , K tc 1 – cos ( ω c τ ) (4.89) 1 λ R, q ⋅ ( 1 – cos ( ω c τ ) ) – λ I, q sin ( ω c τ ) a I, q = ------- ---------------------------------------------------------------------------------------- . K tc 1 – cos ( ω c τ ) (4.90) In order to obtain a physical value, the imaginary part of the limiting depth of cut must be zero, which results in the following root condition: a I, q = 0 ⇒ λ I, q sin ( ω c τ ) cos ( ω c τ ⁄ 2 ) ---------- = ------------------------------- = -----------------------------. λ R, q 1 – cos ( ω c τ ) sin ( ω c τ ⁄ 2 ) (4.91) Substituting Eq (4.91) into Eq (4.89), real depth of cut is solved as: λ R, q λ I, q 2 a R, q = ----------- 1 + ⎛ -----------⎞ . ⎝ λ R, q⎠ K tc (4.92) The root condition in Eq (4.91) can be expressed alternatively as: λ R, q ⋅ cos ( ω c τ ⁄ 2 ) – λ I, q ⋅ sin ( ω c τ ⁄ 2 ) = 0 . (4.93) Dividing both sides by the magnitude of the eigenvalue: λ R, q ⋅ cos ( ω c τ ⁄ 2 ) – λ I, q ⋅ sin ( ω c τ ⁄ 2 ) ----------------------------------------------------------------------------------------------- = 0 , λq where λ q = 2 2 λ R, q + λ I, q , the root condition for each eigenvalue can be simplified as: 107 (4.94) Chapter 4. Dynamics of Milling Operations cos ( ω c τ ⁄ 2 + ϕ q ) = 0 , (4.95) λ I, q ϕ q = atan ⎛ -----------⎞ , ⎝ λ R, q⎠ (4.96) where is the phase shift of the eigenvalue. Finally, the general root function can be determined by multiplying individual root conditions for each eigenvalue: 3 ( 2h r + 1 ) f ( ωc ) = ∏ cos ( ω c τ ⁄ 2 + ϕ q ) = 0 , (4.97) q=1 where h r is the maximum number of harmonics. Since roots (chatter frequencies) of Eq (4.97) guarantee real-valued depths of cut, this equation will be useful when stability lobes are constructed iteratively. Note that, in the most general case, eigenvalues are functions of both excitation and tooth passing frequencies; therefore, ϕ q = ϕ q ( ω c, ω T ) . 4.3.2.1. Zero Order (ZO) Solution The most simplistic approximation to directional coefficient matrix can be done by considering the average component of the Fourier series expansion, i.e., when r = 0. In this case, general characteristic equation in Eq (4.87) reduces to: ( 0, 0 ) det ( I – λ ⋅ G o ) = det ( I – λ ⋅ A 0 ⋅ G ( iω c ) ) = 0 , 108 (4.98) Chapter 4. Dynamics of Milling Operations and the overall transfer function matrix, G ( iω c ) , becomes independent of spindle speed or ω T . This simplifies the problem as eigenvalues of Eq (4.98) can be solved by just specifying a chatter frequency, ω c . Having found the eigenvalues, allowable depths of cut are calculated from Eq (4.92) only when real parts of the eigenvalues are positive. Among all positive depths of cut, the lowest one is selected to mark the stability boundary. The next step is to calculate spindle speeds for this depth of cut. Using the root condition in Eq (4.95), the phase shift can be related to the chatter frequency and the spindle speed in a general form as: π ω c τ ⁄ 2 + ϕ q = --- ( 2k + 1 ) , 2 (4.99) where k = 0, 1, 2… represents "lobe" numbers. From Eq (4.99): ω c τ = 2πk + ε , (4.100) where ε = π – 2ϕ q is the phase angle between current and previous vibration marks. Since eigenvalues hence phase angle can be calculated for a given ω c , a set of spindle speeds, n [rev/min], is obtained from solving for tooth period τ as: 60ω c n = ---------------------------- . N ( 2πk + ε ) (4.101) 4.3.2.2. Generalized Zero Order (GZO) Solution Cutters used in 3-axis milling operations are generally complex cutters with varying geometry along the axis such as ball end and bull nose types. Rotating dynamic cutting forces acting on the 109 Chapter 4. Dynamics of Milling Operations tool must be projected onto fixed machine tool axes using the local geometry of the cutter. It is important how forces are distributed between orthogonal directions along which structural flexibilities are defined for stability consideration. When the stability equation in Eq (4.98) is analyzed, it becomes evident that this equation is limited to an end mill with a single geometry because it contains only one directional coefficient matrix, which is earlier defined in Eq (4.11) to relate dynamic displacements to dynamic cutting forces. Moreover, Eq (4.98) contains only one cutting force coefficient, which is unlikely for general end mills as change of cutting force coefficient along the axis can be significant. In conclusion, zero order stability equation for a general end mill has to be expressed in a more general form. For general end mills, the envelope of the cutter varies along the cutter axis as shown in Figure 4.6.a. Although the cutter can be divided into more axial elements, the geometry in this figure is divided into three parts for simplicity and the stability analysis can still be done without loss of generality. For each section, a relative depth of cut term is assigned, each of which is limited between zero and an upper boundary as shown in Figure 4.6.a. Separating the cutting force coefficient from the depth of cut and considering the influence of each geometry with a different directional coefficient matrix, stability for a general end mill can now be described as: det ( I – ( E 1 ⋅ a 1 ⋅ A 0, 1 + E 2 ⋅ a 2 ⋅ A 0, 2 + E 3 ⋅ a 3 ⋅ A 0, 3 ) ⋅ G ( iω c ) ) = 0 , where E i = 0.5 ⋅ ( 1 – e – iω c τ (4.102) ).K tc, i is the delay term containing the cutting force coefficient K tc, i , a i is the relative depth of cut, and A 0, i is the directional coefficient matrix of each section (i=1,2,3). 110 Chapter 4. Dynamics of Milling Operations z z z a3 ! M @ (a) Sz a2 Nz a1 Mz O # # # N S S S ! r M N @ a1 O ! r M @ a2 a1=Mz O (b) N r (c) Figure 4.6 : a) General end mill with various sections; b) Solution of stability for the bottom most part; c) Switching to second relative depth of cut for stability solution. The characteristic equation of general end milling dynamics given in Eq (4.102) must be analyzed in multiple steps as it varies along the axis. The solution starts with the lowest relative depth of cut, a 1 , having the following characteristic equation: det ( I – E 1 ⋅ a 1 ⋅ A 0, 1 ⋅ G ( iω c ) ) = 0 , (4.103) which can alternatively be represented as a standard eigenvalue problem as: det ( I – λ 1 ⋅ A 0, 1 ⋅ G ( iω c ) ) = 0 , ( 0 < a 1 < M z ) , (4.104) where λ 1 = E 1 ⋅ a 1 . Since Eq (4.104) is identical to Eq (4.98), allowable depth of cut and spindle speeds are solved analytically as described in Section 4.3.2.1 (ZO solution). One important difference is that only depths of cut within the range of characteristic equation, i.e. 0 < a 1 < M z , are accepted to construct stability lobes. For out of range values, i.e. when the critical depth of cut 111 Chapter 4. Dynamics of Milling Operations exceeds the upper boundary of the corresponding geometry, the stability problem is reformulated by considering the effect of successive geometry resulting in a new characteristic equation: det ( I – ( E 1 ⋅ M z ⋅ A 0, 1 + E 2 ⋅ a 2 ⋅ A 0, 2 ) G ( iω c ) ) = 0 . (4.105) The above equation can be represented in an equivalent general eigenvalue problem as: det ( C ( ω c, τ ) – λ 2 ⋅ A 0, 2 ⋅ G ( iω c ) ) = 0 , (4.106) where C ( ω c, τ ) = I – E 1 ⋅ M z ⋅ A 0, 1 ⋅ G ( iω c ) and λ 2 = E 2 ⋅ a 2 . There are three unknowns in this equation: the tooth period ( τ ) or the spindle speed, the depth of cut ( a 2 ) and the excitation frequency ( ω c ). Since chatter vibrations occur at frequencies close to natural frequencies of the structure, the excitation frequency can be scanned around the modes of the structure; however, this will not be sufficient to solve for eigenvalues due to presence of the C matrix. As stated earlier, a physical depth of cut is possible only for eigenvalues satisfying the general root function in Eq (4.97); therefore, first, a frequency band is established in between the smallest ( ω start ) and largest ( ω end ) natural frequencies of the structure. This band is then divided into a finite number of elements at which eigenvalues and root function, f ( ω ) , are calculated at a fixed spindle speed (or similarly tooth period τ ) using standard eigenvalue solver. Brent’s Algorithm [92] is used to solve for roots (chatter frequencies) between frequencies for which root function changes sign. Graphical representation of root search is shown in Figure 4.7. Once critical relative depths of cut, a 2 , are calculated at chatter frequencies, the smallest depth of cut is selected and a range check, i.e. 0 < a 2 < ( N z – M z ) , is performed before accepting it as a valid solution. Advancing through the remaining sections of the tool, the complete stability chart is constructed for a general end 112 Chapter 4. Dynamics of Milling Operations mill. Although GZO solution is an iterative approach, it provides a better estimate for stability as characteristic equation is constructed more realistically. ωstart + +/ Approximate solution sought in between _ + _ ωend _ _ Sign of the root function Freq. Step Size - Δω General Search Fine Tuning Large Eigenvalue Calculation Method Standard Small Standard or Approximate Figure 4.7 : Graphical representation of iterative root search. 4.3.2.3. Multi Frequency (MF) Solution Cutter, workpiece and machine tool structures are subjected to periodic and transient vibrations due to the intermittent engagement of the cutter teeth and periodically varying milling forces. Wave-forms of milling forces have harmonic components of the tooth passing frequency depending on the width of cut and/or the number of teeth. This dependency is such that the smaller the width of cut and/or the less the number of teeth on the cutter, the higher the number of harmonics is involved in cutting forces; therefore, the frequency domain solution should cover wider spectrum of frequencies involved in the cutting process. In this section, multi frequency (MF) solution to dynamic milling equation, which includes multiple harmonics of tooth passing frequency, is analyzed. In Section 4.3.2, marginal stability of milling process is described as a standard eigenvalue problem (Eq (4.76)) including effects of multiple harmonics due to directional coefficient matrix and dynamic cutting forces. Similar to GZO solution, stability lobes for MF case are constructed via iterative solution of Eq (4.76) by assigning spindle speed (or ω T ), and searching for positive-real depths of cut by iterating chatter frequency ( ω c ). Such analysis for a single mode can be found in 113 Chapter 4. Dynamics of Milling Operations the author’s previous work [75]. Unlike GZO solution where the smallest of all depths of cut is selected to represent the stability boundary at a fixed spindle speed, extra care must be given in the case of MF. Although some of these solutions (depths of cut) are related to neighboring branches or lobes, some exist as a result of multiple harmonics. In order to distinguish acceptable solutions from the unacceptable, the physical validity of each solution has to be checked, and in this thesis, it is proposed to do this through physical interpretation of eigenvalue solutions. In Eq (4.76), an eigenvector p is defined for each eigenvalue. Elements of eigenvector represent relative strength of each harmonic on dynamic cutting forces, which are already expanded into Fourier series as shown in Eq (4.66). Since dynamic cutting forces are directly related to regeneration or vibration difference vector dominated by excitation frequency ω c , the first component of eigenvector, p 0 , has to be the biggest in order for the solution to be physically viable. In other words, the root condition in Eq (4.91) is a necessary but not a sufficient condition, and any solution (depth of cut) resulting in an eigenvector whose first component is not being the largest is mathematically possible but physically impossible to exist, therefore, it needs to be eliminated from stability lobes. Properties of different stability methods presented in this chapter are compared in Table 4.1. "Variable" directional coefficient matrix means that the solution method allows the cutter to be divided into a number of axial elements at which directional coefficient matrix can be defined separately. In "scanning" type solution, stability of milling operation is checked by point by point investigation of the parameter domain, i.e., spindle speed, depth of cut, and width of cut; whereas, "analytical" and "iterative" solution types seek for a specific solution defining the stability boundary. 114 Chapter 4. Dynamics of Milling Operations Table 4.1 : Summary of different stability solution methods. SD: Semi-discretization (Section 4.3.1); ZO: Zero order (Section 4.3.2.1); GZO: Generalized zero order (Section 4.3.2.2); MF: Multi frequency (Section 4.3.2.3); Dir. Coef.: Directional coefficient matrix A in Eq (4.16); Sln. Type: Solution type, Scanning means that stability is solved (chatter or stable) for a given set of cutting conditions; Cmp. Time: Computation time, lower number indicates shorter computation time. Method Input Output Dir. Coef. Sln. Type Cmp. Time SD Spindle speed Stability state: Variable Scanning 3 Depth of cut Stable/Unstable ZO Chatter frequency ω c Spindle speeds Fixed Analytical 1 Depth of cut GZO Spindle speed Depth of cut Variable Iterative 2 Chatter frequency ω c MF Spindle speed Depth of cut Fixed Iterative 3 Chatter frequency ω c 4.3.2.4. Iterative Stability Solution and Eigenvalue re-analysis The iterative MF solution becomes computationally demanding because the step size associated with the chatter frequency needs to be small enough to account for all possible roots in the vicinity of a particular structural mode. Structural modes for machine tools are usually in the order of hundreds and even thousands of Hertz; whereas, the step size sometimes becomes as low as 10 –3 Hz, even if only three digit accuracy is demanded, which is not enough for convergence in many cases (two harmonics and up). Another computational difficulty is due to eigenvalue solution. If h r is the total number of harmonics of the tooth passing frequency taken into account, 3(2 h r +1) complex eigenvalues are obtained from Eq (4.87), which means that at each iteration step corresponding to a very small change in chatter frequency estimate, a significant number of eigenvalues have to be recalculated. A better and more intelligent technique is to estimate eigenvalues based on previously calculated ones, hence a search algorithm with a high convergence rate and an eigenvalue estimation method are indispensable. In this section, algorithms presented by Lu et 115 Chapter 4. Dynamics of Milling Operations al. [83] and Lou et al. [85] for eigenvalue re-analysis are applied to increase computational efficiency of the solution of the machine tool stability problem. By eliminating the need for recalculation of eigenvalues, stability lobes are solved and constructed more intelligently. Consider the generalized eigenvalue problem of n-by-n complex matrix pairs A and B: ⎫ A xi = λi B xi ⎪ ⎬, * * yi A = yi λi B ⎪ ⎭ (4.107) where λ i ’s are the eigenvalues, x i ’s and y i ’s are the associated right and left eigenvectors, respectively, i = 1, …, n , and "*" is the complex conjugate transpose operator. Suppose that the original eigenvalue problem is represented when A= A (0) , B= B (0) (0) (0) , x i = x i , y i = y i , and λ i = λ 0i . For a finite change of matrices A and B, the modified eigenvalue problem becomes: A (1) (1) xi ( 1 )* ( 1 ) yi A where A (1) =A (0) + ΔA , B (1) =B (0) = λ 1i B = (1) (1) xi ( 1 )* y i λ 1i (1) B (1) (0) ⎫ ⎪ ⎬, ⎪ ⎭ (4.108) (1) (0) + ΔB , x i = x i +Δx i , y i = y i +Δy i , λ 1i = λ 0i + Δλ i . Note that Δ terms represent perturbation due to a parameter change in matrices A and B, which is the frequency in this case. Since the procedure is identical for both right and left eigenvalues, in the following section, derivations for only right eigenvalue problem is presented. Substituting the perturbed values in right eigenvalue problem in Eq (4.108) and omitting a small second order component [83]: 116 Chapter 4. Dynamics of Milling Operations (A (0) – λ 0i B (0) (0) ) ⋅ Δx i + ( ΔA – λ 0i ΔB ) ⋅ x i = (B (1) (0) xi +B (0) Δx i ) ⋅ Δλ i . (4.109) When it is assumed that the system has distinct eigenvalues, the right and left eigenvectors form a complete set of vectors, i.e., they become linearly independent. Thus, perturbations Δx i and Δy i (0) can be expressed as complex linear combinations of all eigenvectors, x i n Δx i = (0) g ji x j ∑ j = 1( ≠ i) n Δy i = ∑ * (0) h ji y j j = 1( ≠ i) (1) (0) Note that, in ( x i = x i + Δx i ), the effect of i (1) because the mode vector x i hence the i th th (0) and y i ⎫ ⎪ ⎪ ⎪ ⎬. ⎪ ⎪ ⎪ ⎭ (0) vector x i [85]: (4.110) can be combined into the first term only describes the proportional relation among its elements [85], term is excluded in Eq (4.110). This statement is also valid for the left eigenvector. ( 0 )* When Eq (4.109) is multiplied by y i and the orthogonality condition of the eigenvectors, being * * y i A x j = 0 , y i B x j = 0, i ≠ j , (4.111) is applied, the variation in eigenvalues can be easily calculated as: ( 0 )* (0) ( 0 )* (1) (1) (0) ( 0 )* (1) (0) y i ( ΔA – λ 0i ΔB ) x i y i ( A – λ 0i B ) x i yi A xi - = -------------------------------------------------------------- = ------------------------------ – λ 0i . (4.112) Δλ i = ------------------------------------------------------------( 0 )* ( 1 ) ( 0 ) ( 0 )* ( 1 ) ( 0 ) ( 0 )* ( 1 ) ( 0 ) yi B xi yi B xi yi B xi Finally, eigenvalue estimate for the modified system becomes: 117 Chapter 4. Dynamics of Milling Operations ( 0 )* (1) (0) yi A xi -. = -----------------------------( 0 )* ( 1 ) ( 0 ) yi B xi λ 1i (4.113) This is just a specific form of the Rayleigh quotient for the eigenvalue problem in Eq (4.107). If ( 0 )* Eq (4.110) is substituted into Eq (4.109) and the result is multiplied by y r ( r ≠ i ), the follow- ing form is obtained: ( 0 )* yj (A (0) – λ 0i B ( 0 )* = yj (B (1) (0) –B (0) (0) (0) Δλ i ) ( x 0 g 0i + x 1 g 1i + … ) , (0) Δλ i – ΔA + λ 0i ΔB ) x i . (4.114) Using the orthogonality condition in Eq (4.111), the above equation is simplified as follows: ( 0 )* yj (A (0) ( 0 )* yj – λ 0i B (A (0) –B (0) (0) (0) ( 0 )* Δλ i ) x j g ji = y j (0) ( 0 )* ( λ 0j – λ 1i )B (0) (0) x j g ji = yj (1) (1) – λ 1i B ( 0 )* yj (0) ( 0 )* ) x j g ji = y j (B (B (1) ( 0 )* (1) (0) Δλ i – ΔA + λ 0i ΔB ) x i Δλ i – A ( λ 1i B (1) (1) + λ 0i B –A (1) (1) (0) ) xi (0) ) xi (0) y j ( A – λ 1i B ) x i - ( j = 1, …, n ), ( j ≠ i ). g ji = --------------------------------------------------------------( 0 )* ( 0 ) ( 0 ) y j B x j ( λ 1i – λ 0j ) (4.115) Coefficients for the left eigenvalue problem in Eq (4.108) are similarly obtained as: ( 0 )* (1) (1) (0) y i ( A – λ 1i B ) x j - , ( j = 1, …, n ), ( j ≠ i ), h ji = --------------------------------------------------------------( 0 )* ( 0 ) ( 0 ) y j B x j ( λ 1i – λ 0j ) and finally eigenvectors of the perturbed system are easily found by: 118 (4.116) Chapter 4. Dynamics of Milling Operations (1) xi (1) yi ⎫ + Δx i ⎪ ⎬. (0) = y i + Δy i ⎪ ⎭ (0) = xi (4.117) Since both eigenvalue and eigenvector estimates are expressed in terms of modified matrices, it is possible to use them in an iterative sense. Iteration starts with the exact solution of unperturbed (0) (0) system: λ 0i , x i , and y i and the eigenvalue solution of the modified system is estimated using (1) (1) Eqs (4.113) and (4.117): λ 1i , x i , and y i . Convergence at each step is checked by comparing the rate of change of eigenvalue estimate: ( λ 1i – λ 0i ) ⁄ λ 1i ≤ tol where tol is user-defined tolerance. If the iteration is not successful, a new set of eigenvalues are estimated by replacing the (0) (1) original eigenvalues and eigenvectors with the new estimates, i.e. λ 0i = λ 1i , x i = x i (0) (1) and 3 y i = y i . Although the perturbation method described above involves operation count of n , which is identical to a standard n-by-n eigenvalue problem operation count, computation time is still reduced by utilizing available information (eigenvalues and eigenvectors of the unperturbed problem) as initial estimates in the solution of the modified eigenvalue problem. Comparing Eqs (4.87) and (4.107), matrix pairs A and B for chatter stability are simply equal to: A = I, B = Φ ( ω, ω T ) . (4.118) It is essential for the eigenvalue estimation to trace the same eigenvalue when stepping from one frequency to another. At the beginning, eigenvalues are calculated for a set of frequencies; however, their orders are independent of each other at each frequency. Figure 4.8.a shows an example of variation in eigenvalue indices when frequency is slightly changed. Note that index of each eigenvalue corresponds to the element number of the eigenvalue vector obtained from solution of 119 Chapter 4. Dynamics of Milling Operations Eq (4.87). An appropriate sorting criterion to obtain consistent indices is to check the relative displacement of eigenvalues. Eigenvalues obtained at the starting frequency, ω start , are taken as the reference for indexing, and when the frequency is varied slightly, the distance of each new eigenvalue from the previous one is calculated. By sorting the distances, proximity of each eigenvalue is determined, and proper indexing is achieved. Figure 4.8.b shows the effect of sorting on indexing eigenvalues and how eigenvalues are successfully traced. Eigenvalue estimation is incorporated into the stability solution to increase computational efficiency. The only difference from the previously presented method is that the fine tuning stage of the root search (Figure 4.7) is now achieved by estimating eigenvalues instead of solving them using a standard eigenvalue solution. λ2 λ5 λ6 λI λ2 ω0 ω1 λ5 λ3 λ4 λ1 λ1 λ4 λ6 λ5 λ6 λI λ6 λ3 ω0 ω1 λ2 λ5 λ2 λ4 λ4 λ1 λ1 λ3 λ3 λR λR (b) (a) Figure 4.8 : Shift of eigenvalues when the frequency is changed from ω 0 to ω 1 : (a) unsorted eigenvalues; (b) sorted eigenvalues. 4.4. Comparison of Stability Methods In an effort to compare various stability methods presented above, the stability of a milling operation is investigated. Dynamic parameters of a 76.2 [mm] face milling cutter with four 25.4 [mm] 120 Chapter 4. Dynamics of Milling Operations circular inserts (N=4) mounted on a three axis horizontal machining centre are given in Table 4.2 [101]. Stability is investigated for a half immersion down milling ( φ st =0 and φ ex = π ⁄ 2 ) of a P20 steel with the following constant tangential, radial and axial cutting force coefficients, respec2 2 2 tively: K tc =3146 [N/ mm ], K rc =1680 [N/ mm ], K ac =419.4 [N/ mm ]. In addition to time (SD) and frequency (ZO, GZO, MF) domain based stability methods, a comprehensive time domain simulation is also used to construct stability lobes for comparison. Direction X Y Table 4.2 : Dynamic parameters [101]. Natural Frequency Stiffness [Hz] [N/m] 28 2.54 e7 55 4.30 e7 28 2.15 e8 55 6.18 e8 Damping [-] 0.17 0.12 0.10 0.06 In order to account for the variation of directional coefficients along tool axis, the cutter is divided into a number of axial elements and individual directional coefficients are calculated. Since depth of cut is scanned in the semi-discretization (SD) and generalized-zero-order (GZO) solutions, overall directional coefficients at a given depth are calculated by aggregation of the axial elements. On the other hand, for methods seeking for a depth of cut such as zero-order (ZO) and multi-frequency (MF), average directional coefficients are used at mid point, i.e., when κ = π ⁄ 4 . Number of discrete time elements for the SD is varied depending on the spindle speed to accommodate either minimum of 40 elements, or a time step that is one tenth of the lowest period (inverse of the highest natural frequency) in the system to prevent aliasing. The sampling time of the time-domain (TD) solution is also varied with a similar approach to capture dynamics accurately. 121 Chapter 4. Dynamics of Milling Operations A complete set of stability lobes are plotted in Figure 4.9 and chatter frequencies obtained from frequency domain based solutions are given in Figure 4.10. The ZO solution is improved with the proposed GZO by considering the variation of the directional coefficient matrix along the cutter, and it has become closer to the TD solution. The SD and the GZO are in good agreement below 425 rpm; however, the minimum stability limit, which is expected to be constant, has increased from 6 mm to 8.7 mm above 425 rpm when the SD solution is used. Similar increase is also observed in the TD solution. In an effort to investigate this behavior, frequencies are checked in this region. Chatter frequency at 475 rpm, where tooth passing frequency is ω T ⁄ 2π =31.7 Hz, is found to be ω c ⁄ 2π =60.2 Hz. As mentioned before, when chatter occurs, dynamic cutting forces Eq (4.66) - have components of not only chatter frequency but also frequencies that are a tooth passing frequency away from chatter frequency. Unless radial immersion is very low, these harmonics do not influence the stability [75]. In this case-study, although radial immersion is considerably high, it is found out that harmonics play an important role around a specific spindle speed region. For example at 475 rpm, ( ω c - ω T )/ 2π =28.5 Hz corresponds to one of the structural modes resulting in the distribution of vibration energy between two neighboring modes and increased stability. In order to observe the effect of harmonics, stability lobes are generated by using MF solution. Low ( h r =1) and high ( h r =5) numbers of harmonics are used for comparison. Although, one harmonic is enough to capture the stability increase, better accuracy is obtained with the addition of more harmonics. It is observed that the stability solution has not changed when more than five harmonics are used. Eigenvectors are analyzed at different spindle speeds, and magnitudes of them are normalized with respect to the magnitude of eigenvectors at the chatter vibration frequency, i.e., 122 Chapter 4. Dynamics of Milling Operations T ⎧ ⎫ ⎪ ⎪ p p p p h – h 1 –1 p = ⎨ { 1 } T ------- ---------… --------r- -----------r- ⎬ . p0 p0 p0 p0 ⎪ ⎪ ⎩ ⎭ (4.119) Figure 4.11 displays eigenvector magnitudes in the x direction at two spindle speeds. At 475 [rpm], where increased stability occurs, eigenvectors corresponding to harmonics have larger magnitudes, and the strength of the first harmonic (1,-1) is almost 80% of the magnitude of eigenvector corresponding to chatter vibration frequency (see Figure 4.11.a). Since the magnitude of each eigenvector component represents the relative vibration strength at the corresponding frequency, the increase in stability can be explained by the fact that vibration energy is split mostly between two modes of the structure. This observation is similar to the working principle of a tuned mass damper (harmonic absorber). At 310 [rpm], where the MF solution becomes almost identical to the ZO solution, the magnitudes of harmonics are generally smaller and strength of the first harmonic has dropped to less than 40% of that of the one at chatter vibration frequency. When general stability trend is concerned, it is observed that the MF solution is capable of capturing the increase in stability when a mode is located at a frequency that is one or a couple of tooth passing frequencies away from the dominant chatter frequency. 123 Chapter 4. Dynamics of Milling Operations 18 14 12 MF-1 10 SD MF-5 TD 8 6 GZO ZO 4 200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 1100 1150 Spindle Speed [rpm] Figure 4.9 : Stability lobes obtained using SD: Semi-discretization, ZO: Zero order solution, GZO: Generalized zero order solution, MF-1: Multi frequency solution with one harmonic ( h r =1), and MF-5: Multi frequency solution with five harmonics ( h r =5). 80 Chatter Frequency [Hz] Depth of Cut [mm] 16 75 70 MF-1 MF-5 65 60 ZO 55 200 250 300 350 400 450 500 550 600 Spindle Speed [rpm] Figure 4.10 : Chatter frequencies obtained from frequency domain based solutions. 124 Normalized Eigenvector Magnitues Chapter 4. Dynamics of Milling Operations 1 ωc ωc (ωc−ωT) 0.8 0.6 (ωc+ωT) (ωc−ωT) 0.4 (ωc+ωT) 0.2 0 -5 -4 -3 -2 -1 0 1 2 3 4 5 -5 -4 -3 -2 -1 0 1 2 3 4 5 (a) (b) Figure 4.11 : Normalized X-direction eigenvector magnitudes of MF-5 at (a) n=475 [rpm]; ω T / 2π =31.7 [Hz], ω c / 2π =60.2 [Hz]; (b) n=310 [rpm], ω T / 2π =20.7 [Hz], ω c / 2π =58.5 [Hz]. The same spindle speeds are analyzed using the SD solution by selecting depths of cut right on the stability border, where modulus of the largest eigenvalue is unity. Eigenvalues of the transition matrix - Eq (4.58) - are solved for spindle speeds of 475 [rpm] and 310 [rpm], and the results are shown in Figure 4.12. At 475 [rpm] and depth of cut of 12.8 [mm], more eigenvalues appear closer to the unit circle compared to that of the transition matrix with 310 [rpm] and depth of cut of 7.6 [mm]. Although the marginally stable vibration is dominated by the eigenvalue pair equal to unity while others diminishing at a steady state, another eigenvalue pair being very close to unity (modulus of 0.94) at 475 [rpm] is suspected of causing the unusual increase in stability. This observation is parallel to the results obtained from the multi-frequency (MF) solution, however, no analytical proof is available. 125 Chapter 4. Dynamics of Milling Operations 1 Unit Circle 4 0 1.0 Imaginary 0.9 Imaginary 0.5 Unit Circle 0 -0.5 -1 -1 -0.5 0 0.5 1 -1 Real -0.5 0 0.5 1 Real (a) (b) Figure 4.12 : Eigenvalues of the transition matrix - Eq (4.58) - of the semidiscretization method (k=41) when cutting conditions are selected on stability border at (a) 475 [rpm], 12.8 [mm]; (b) 310 [rpm], 7.6 [mm]. The performance of the eigenvalue estimation technique is also tested using the same milling system. Since the root function is critical to solve for chatter frequencies, its estimation is used for performance evaluation. The root function considering one harmonic is calculated with exact eigenvalues obtained at 1000 sample frequency points evenly spaced between ω start / 2π =14 [Hz] and ω end / 2π =100 [Hz]. The same root function is approximated at sample points using eigenvalue estimation in between limited number of exactly calculated eigenvalues, and comparisons are presented in Figure 4.13. The approximation in Figure 4.13.a is based on 10 exact eigenvalues and it is clear that the estimation suffers. On the other hand, increasing the number of exact eigenvalues to 50 improved estimation significantly as seen in Figure 4.13.b. In summary, approximating the eigenvalue problem with 50 exact points requires less computational load while maintaining acceptable accuracy. 126 Chapter 4. Dynamics of Milling Operations Root Function, f (ω) [-] 1 Exact Approximate Exact Approximate 0.5 0 -0.5 10 20 30 40 50 60 70 80 90 100 10 Frequency [Hz] 20 30 40 50 60 70 80 90 100 Frequency [Hz] (a) (b) Figure 4.13 : Comparison of exact and approximate root functions of one harmonic case using (a) 10 ( Δω / 2π =8.6 [Hz]); (b) 50 ( Δω / 2π =1.7 [Hz]) exact frequency elements. 1000 ( Δω / 2π =0.086 [Hz]) sampling points are used and eigenvalues are estimated between exact values obtained from a standard eigenvalue solver. Spindle speed is 475 [rpm]. 4.5. Summary Several frequency and time domain chatter stability techniques are presented in this chapter. Each method has its advantages, disadvantages and limitations. The zero order (ZO) solution is the only and fastest approach that solves for chatter free cutting conditions directly; however, it is limited to single immersion cases with constant cutting force coefficients. The proposed Generalized Zero Order (GZO) solution on the other hand, is an iterative approach that includes the effect of varying cutter geometry (i.e. varying directional coefficient matrix) and cutting force coefficients along the axis of the cutter. Cutting conditions are specified as inputs; therefore, complex cutter engagements (multiple depths and widths of cut) can be successfully handled with this solution, which is a requirement for virtual milling of parts. Semi-discretization (SD) and multi frequency (MF) methods are capable of capturing the effect of higher harmonics mostly due to multiple mode excitation or highly interrupted cutting, which happens when radial immersion is lower than 127 Chapter 4. Dynamics of Milling Operations 10% of the cutter diameter in general. Both solutions are indirect as stability lobes are generated iteratively or by checking stability for specified cutting conditions. Although the problem size is increased drastically compared to the zero order solution, longer computation time is justified with the increased accuracy in stability prediction for finish milling operations. Stability lobes provide guidelines for chatter free selection of spindle speed, axial depth and radial depth of cut during process planning, as well as checking whether the process will experience chatter along the tool path. The optimization of milling process in the Virtual Machining environment is presented in the following chapter, where the chatter is used as a fundamental constraint. Chatter free cutting conditions are used as a safe starting point and further optimized with respect to other physical constraints. Since the biggest improvement in cycle time reduction can be achieved during roughing, and semi finishing operations, which have simple and large cutter engagements, the zero order solution is frequently used in the next chapter to solve for stability lobes. 128 Chapter 5 Process Optimization 5.1. Introduction In the manufacturing industry, significant inefficiency occurs due to improper selection of machining parameters such as speeds, feedrates, and depths of cut. Traditionally, industry uses machinability data handbooks as a source for selection of parameters, which generally represent conservative machining conditions, serving as a safe operating point. These parameters are then tuned for individual operation by costly trial and error type experiments, which lead to machine down times and scrap parts. In this chapter, the optimization of milling processes are based on process physics, and linear and nonlinear programming. The physics based planning, provides a useful tool to process engineers to tailor machining parameters in a computer environment during the Computer Aided Manufacturing (CAM) stage and simulate the performance for different cutting conditions for efficiency comparison. This chapter presents a novel generalized optimization strategy that combines different machining constraints including cutting forces, chip thickness, spindle torque-power, form errors on the workpiece and even stability of the system to determine the most efficient machining parameters. The nature of optimization varies with the availability of the design variables and constraints; therefore, both linear and non-linear programming algorithms are used to obtain optimum machining parameters. 129 Chapter 5. Process Optimization 5.2. Process Optimization Manufacturing of a part follows a well-defined methodology (Figure 5.1) that starts with the conceptual design of a part generated in a solid modeler by a design engineer. This step is known as the Computer-Aided Design (CAD). Later, in the CAM stage, a process engineer takes over the solid model to generate necessary numerical control (NC) code that contains cutting tool motions in workpiece space and process parameters such as feedrate and spindle speed values. The CAM stage requires the most expertise due to its direct impact on the process efficiency and the quality of finished product. Selection and even sequence of cutting operations make a significant difference in cycle time of a part. After part programming is completed, the part is taken over by a shop floor engineer who tests the machinability of the part under specified cutting conditions and makes changes to the part program if necessary. At the same time, costly cutting trials are conducted in order to increase the productivity. A much preferred alternative, however, is to be able to simulate and optimize the process in virtual environment. This will not only reduce the set-up time drastically, but it will also cut down unnecessary expenditure due to increased machine down times and scrapped parts. Some of the critical cutting variables involved in process planning are classified as geometric parameters such as depth and width of cut, and performance parameters like feedrate and spindle speed. Proper selection of all will lead to satisfactory surface finish tolerance with optimum power consumption, longer tool life, and higher productivity. Necessary relations to simulate machining operation of a part have been discussed in an earlier chapter. In this chapter, machining state variables (process outputs) will be used to obtain optimum machining conditions. 130 Chapter 5. Process Optimization Concept Solid model CAD STAGE Process Engineer Shop Floor Engineer NC Code generation (process parameters) Application of NC code on real machine Determination of machining failures CAM STAGE Failure Success Design Engineer Experiments Production CURRENT PROPOSED Virtual Process Simulation & Optimization Module Optimum Production Figure 5.1 : Manufacturing flow chart, current and proposed methodologies. The optimization problem is separated into two independent parts. The first part is called Pre-Process and the second part is referred as Post-Process Optimization. Although these optimization strategies, which will be presented in the following sections, are different, they both share the same objective which is to maximize productivity of milling operations. Highly efficient manufacturing can only be achieved by maximizing the Material Removal Rate (MRR); therefore, the MRR will be the objective function throughout this chapter. In manufacturing, the MRR is defined as: MRR = a ⋅ B ⋅ c ⋅ n ⋅ N , (5.1) where a is the axial depth of cut, B is the radial width of cut, c is the feed per tooth, n is the spindle speed and N is the number of flutes on the milling cutter, see Figure 5.2. In an optimization problem, parameters a, B, c, n, and N are called design variables. Objective function in Eq (5.1) is subject to two types of constraints named linear and nonlinear constraints. Linear constraints are bounds on the design variables and can be simply expressed as: q lb ≤ q ≤ q ub , 131 (5.2) Chapter 5. Process Optimization where q is a design variable limited by lower bound q lb and upper bound q ub . Nonlinear constraints; on the other hand, are rather expressed in a function form such as: f nc ( a, B, n, N ) ≤ 0 , (5.3) where f nc is a nonlinear constraint that is expressed as a function of design variables. c n B a Figure 5.2 : Design variables. In milling operations, sources of constraints can be classified under two categories. First is the limitations due to machine tool characteristics, which involve machine tool torque/power curves, structural dynamic properties of cutter-machine tool assembly, minimum/maximum spindle speed, and axis limitations (maximum feedrate). The second source of constraints is due to the machining properties. Form error due to static deflection of the tool, minimum/maximum chip load, maximum surface speed and tool life, are the most important limitations in this category. Some of the constraints listed above directly become upper and lower bounds on one of the design variables, i.e., they are linear constraints. For example, minimum/maximum spindle speed and feedrate constraints can be simply expressed as n min ≤ n ≤ n max and c min ≤ c ≤ c max , respectively. The cutting speed influences tool wear, and the chip thickness below a threshold will cause the tool to rub the workpiece material, whereas, the chip thickness above a maximum value will lead to cutting edge chipping and eventually tool breakage. On the other hand, machining process 132 Chapter 5. Process Optimization outputs such as form error, torque/power demands, and dynamic characteristics of the cutter/ machine tool indirectly limit design variables; therefore, they impose nonlinear constraints. In the following sections, linear or nonlinear relationship between process outputs and design variables will be derived depending on the type of the optimization problem (Pre or Post Process Optimization). 5.3. Pre-Process Optimization The ultimate goal of the virtual machining concept is to be able to simulate and optimize milling operations at the process planning stage. In a sense, the Pre-Process Optimization serves as a tool given to a process engineer to obtain the most optimum cutting parameters based on constraints during the process planning stage. In other words, it provides the CAM software with the upper bound of cutting parameters to obtain an efficient milling operation (NC program). Figure 5.3 shows at which stage of manufacturing planning, the Pre-Process Optimization is used. Given in Eq (5.1), the objective function has five design variables in total. The cutter is usually chosen before selecting cutting conditions and generating the tool path; therefore, one of the design variables, number of teeth (N) is known prior to the optimization. Spindle speed (n), on the other hand, is selected based on either surface finish and tool life requirements specified in handbooks, or stability of cutting due to dynamic of milling, details of which will be presented later in this section. Finally, objective function for the Pre-Process Optimization reduces down to the following form: f obj = a ⋅ B ⋅ c . 133 (5.4) Chapter 5. Process Optimization 5.3.1. Stability Constraint Chatter stability is one of the limiting factors for the selection of cutting parameters. Stability of a milling process is dependant on three parameters: depth of cut, width of cut, and spindle speed. Improper selection of any of these leads to an unacceptable surface finish due to chatter vibration marks as depicted in Figure 5.4. The mathematical relationship between design variables and dynamics of a milling process has already been comprehensively presented in the Chapter 4 of this thesis; therefore, stability of a milling process is taken as input in this section. Process Planning (CAM) CAD Model Milling Cutters / Operations Cutting Conditions Depth of Cut (a) Width of Cut (B) Feed Rate (c) Spindle Speed (n) ... Facing Pocketing Contouring Plunging Ramping ... CAM Tool Path Generator Optimized NC Code Pre-Process Optimization Constraints: Machining Machine Tool Figure 5.3 : Optimized NC code generation flow chart with the Pre-Process Optimization module. Stability boundaries represented as a function of cutting parameters are called stability lobes. For a constant width of cut (B), stability can be solved to obtain chatter stability diagram presented in Figure 5.5.a. Operation points that fall under the graph lead to chatter free machining. Spindle speeds at and around that the peaks of each lobe such as n lobe, 1 and n lobe, 2 have advantage because depth of cut can be increased while maintaining the stability of the process. Alternatively, stability lobes can be obtained for various widths of cut while keeping axial depth of cut constant. 134 Chapter 5. Process Optimization For convenience, instead of actual width of cut, normalized width of cut, b, is preferred as it varies between 0 and 1 corresponding to no and full immersion cases, respectively: ⎧ ⎪ ( 1 – cos φ ex ) ⁄ 2 for up milling b = B⁄D = ⎨ ⎪ ( 1 + cos φ st ) ⁄ 2 for down milling ⎩ ⎫ ⎪ ⎬, ⎪ ⎭ (5.5) where D is the cutter diameter, φ st and φ ex are the cutter entry and exit angles, respectively. n Vibr atio nM arks B a Figure 5.4 : Design variables that are important to suppress chatter vibrations. Stability lobes for constant width of cut are shown in Figure 5.5.b. When both stability lobes are compared, it can be seen that locations of the lobes remain same. This is due to the fact that locations of the lobes are not determined by the type of solution but rather by the dynamic parameters of the system, which are functions of structural elements such as tool, tool holder, spindle and machine tool. Having located optimum spindle speeds, design space - shown in Figure 5.6 - can now be formed by extracting maximum stable depth and width of cut at a constant spindle speed from various stability charts similar to the ones in Figure 5.5. but obtained with different widths and depths of cut in parallel to the proposed method in reference [61]. The use of the stability design space will become more significant when case studies are presented at the end of this section. 135 Chapter 5. Process Optimization 1.0 a = const b - Width of Cut a - Depth of Cut B = const Stable Region Stable Region nlobe,2 0 n lobe,1 n - Spindle Speed n lobe,2 nlobe,1 n - Spindle Speed (a) (b) Figure 5.5 : Chatter stability lobes for constant (a) radial width of cut, (b) axial depth of cut. 1.0 b - Width of Cut n = const Feasible Region (safe to operate) 0 a - Depth of Cut Figure 5.6 : Design space bounded by chatter stability. 136 Chapter 5. Process Optimization 5.3.2. Machine Tool Torque/Power Constraints Apart from stability, torque/power limit of the machine tool shall not be exceeded at any time. Both torque and power are dependant on the tangential component of the cutting force, F t ( φ ) , and vary as the cutter rotates: Torque ( φ ) = F t ( φ ) ⋅ D ⁄ 2 [Nm], (5.6) D.π.n –3 Power ( φ ) = F t ( φ ) ⋅ ⎛ ---------------⎞ ⋅ 1.34102 × 10 [hp], ⎝ 60 ⎠ (5.7) where D is the cutter diameter in [m], n is the spindle speed in [rpm], and F t is in [N]. In this section, an analytical solution for cylindrical end mills is derived to determine the maximum torque and power rapidly as design variables (cutting conditions) vary so that design space can be constructed for torque/power constraint. Consider a Cutter Engagement Feature (CEF) for a cylindrical end mill composed of depth of cut "a" and immersion angles " φ st " and " φ ex " as shown in Figure 5.7. Angular positions of the cutter at which the first flute completely enters into and exits the CEF are called φ cr, 0 and φ cr, 1 , respectively and they can be expressed as: φ cr, 0 = φ st + k i0 a , (5.8) φ cr, 1 = φ ex + k i0 a , (5.9) where k i 0 = ( 2 tan i 0 ) ⁄ D and i 0 is the helix angle. Two different cases occur when φ cr, 0 < φ ex and φ cr, 0 ≥ φ ex , and they are named as "Case I" (Figure 5.7.a) and "Case II" (Figure 5.7.b), respectively. 137 Chapter 5. Process Optimization a φst φcr,0 φex φcr,1 e ut te 0 a Fl Depth of Cut z F lu Depth of Cut z φ 0 Immersion Angle φst φex φcr,0 φ φcr,1 Immersion Angle (a) (b) Figure 5.7 : Important parameters for analytical solution: (a) Case I ( φ cr, 0 < φ ex ), (b) Case II ( φ cr, 0 ≥ φ ex ). For a cylindrical end mill, tangential cutting force for flute j can be obtained by integrating differential cutting forces along: t Fj ( φj ) = z j, 1 ∫z ( K te + K tc .c. sin φ ) dz , (5.10) j, 0 where φ = φ j – k i0 z is the angular position of differential element dz, φ j = φ + ( j – 1 )φ p is the current angular position of flute j ( j = 1, …, N ), N is the number of teeth, φ p is the pitch angle and defined as φ p = 2π ⁄ N for a uniform pitch cutter. Integration can be defined between alternative boundaries using dφ = – k i0 dz conversion: – 1 φ j, 1 F tj ( φ j ) = ------ ∫ ( K te + K tc .c. sin φ ) dφ , k i 0 φ j, 0 (5.11) which can be further simplified as: 1 F tj ( φ j ) = ----- { K tc ⋅ c ⋅ ( cos φ j, 1 – cos φ j, 0 ) + K te ( φ j, 0 – φ j, 1 ) } , ki0 (5.12) where lower ( φ j, 0 ) and upper ( φ j, 1 ) integration boundaries are the angles at which flute enters and exits the cutter workpiece engagement area (see Figure 5.8). 138 Chapter 5. Process Optimization Depth of Cut (a) Immersion Angle,( φ ) y y Flute j φj,1 φj,0 x x Figure 5.8 : Cutter workpiece engagement, integration boundaries. Total tangential cutting force acting on the tool can be obtained by summing individual contributions of all flutes: N Ft = ∑ Ftj . (5.13) j=1 Variation of lower ( φ j, 0 ) and upper ( φ j, 1 ) integration limits are graphically represented in Figure 5.9, where grey area represents the "in-cut" zone. Since the effect of helix is already embedded in the integration limit, the cutting flute can be represented as a vertical line moving to the left of the figure as time advances. The cutting zone is divided into three distinct regions and the definition of each integration limit is tabulated in Table 5.1 for two different cases. 139 Chapter 5. Process Optimization Table 5.1 : Integration boundaries for different regions. CASE I Region 1 Region 2 Region 3 φj,1 φst φj - kio*a φj - kio*a φj,0 φj φj φex CASE II Region 1 Region 2 Region 3 φj,1 φst φst φj - kio*a φj,0 φj φex φex where Region 1: ( φ st ≤ φ j ≤ φ b1 ) , Region 2: ( φ b1 < φ j ≤ φ b2 ) , Region 3: ( φ b2 < φ j ≤ φ cr, 1 ) , and φj,0 3 2 φj,1 φst 1 φst φcr,0 φj φex φcr,1 (5.14) CASE II φex φ0 φ1 2 3 1 φst φst φex φcr,0 φex φ0 φ1 Integration Boundaries (φ0,φ1) CASE I FLUTE "j" Integration Boundaries (φ0,φ1) φ b1 ⎧ ⎫ ⎧ ⎫ ⎪ φ cr, 0 Case I ⎪ ⎪ φ ex Case I ⎪ = ⎨ ⎬ , φ b2 = ⎨ ⎬ ⎪ φ ex Case II ⎪ ⎪ φ cr, 0 Case II ⎪ ⎩ ⎭ ⎩ ⎭ φcr,1 Immersion Angle (φ) Immersion Angle (φ) (a) (b) Figure 5.9 : (a) Case I ( φ cr, 0 < φ ex ), (b) Case II ( φ cr, 0 ≥ φ ex ). Total tangential force given in Eq (5.13) can be expressed in a general form when integration boundaries in Table 5.1 are substituted. In order to isolate the dependency of integration bound- 140 Chapter 5. Process Optimization aries on the immersion angle, sum/difference formulae for sine and cosine are utilized, and analytic tangential force is reorganized as: N Ft = ∑ ( Ktc ⋅ c ⁄ ki ) ( Cc, j ⋅ cos φ + Cs, j ⋅ sin φ + Cp, j ⋅ φ + const ) , (5.15) 0 j=1 where const is the collection of all terms that does not contain any angle ( φ ) or time dependent term next to it. Coefficients, C, are derived for each different case and region, and given in Table 5.2. In order to locate the extremum, the first derivative of tangential force is equated to zero: dF t -------- = ( K tc ⋅ c ⁄ k i ) ( – C c ⋅ sin φ + C s ⋅ cos φ + C p ) = 0 , 0 dφ where C q = ∑ Cq, j , (5.16) ( q = c, s , p ) . j Table 5.2 : Generalized coefficients for tangential force. CASE I/II Region 1 CASE I Region 3 CASE II Region 2 φ st ≤ φ j ≤ φ b1 φ b2 < φ j ≤ φ cr, 1 φ b1 < φ j ≤ φ b2 Cc,j - cos((j-1)φp) cos((j-1)φp- kio*a) cos((j-1)φp- kio*a) - cos((j-1)φp) 0 Cs,j sin((j-1)φp) - sin((j-1)φp- kio*a) - sin((j-1)φp- kio*a) + sin((j-1)φp) 0 Cp,j Kte /(Ktc.c) 0 - Kte /(Ktc.c) Depending on the angular location of the cutter, flutes that are in contact with the workpiece vary as shown in Figure 5.10. At one point in time, only flutes 1, 2, and 3 are in cut at regions 1, 2, and 3; however, another moment same flutes cut regions 2, 2, and 3. In order to determine all possible intersection configurations, pitch angle is divided into a couple of pieces and tool position is swept in that range. Since cutting forces in milling are periodic at tooth period, it is guaranteed 141 Chapter 5. Process Optimization that no combination will be left out. Moreover, such sweeping process is very fast because it does not require any heavy mathematical operation (only boolean operations are used). φ0,φ1 Integration Boundaries j=1 1 2 2 3 3 4 4 3 shift of flutes with time 2 1 φj=1 = φ(t) φj=2 = φ+φp φp φ Figure 5.10 : Possible intersection scenarios. Once all possible scenarios are identified, critical angular positions of the cutter at which minimum and maximum forces appear can be calculated by solving Eq (5.16). Using the trigonometric identity, sin φ = 2 1 – cos φ , Eq (5.16) can be reduced into quadratic equation and two roots, if they exist, are calculated as: φ roots ⎛ –C C ± C C2 + C2 – C2 ⎞ s p c c s p⎟ = acos ⎜⎜ ---------------------------------------------------------------2 2 ⎟ , ( 0 ≤ φ roots ≤ 2π ) . C + C c s ⎝ ⎠ (5.17) In addition to roots obtained above, angular positions φ cr, 0 , φ cr, 1 , φ st , and φ ex , which are located at the boundaries of piece-wise integration limits, are also taken into account and stored as roots. At final stage, all possible roots are simply substituted into Eq (5.15) to obtain tangential forces. Roots giving the global minimum and maximum forces are selected as the critical angular t t position, φ min and φ max , respectively. 142 Chapter 5. Process Optimization Once tangential forces are known, torque and power demands denoted by Torque ( φ ) and Power ( φ ) are easily calculated. In order to prevent violation of torque/power constraints, cumulative average torque and power shall not exceed machine tool torque/power limitations at spindle speed n: Torque cav ≤ T mo ( n ) ⎫ ⎬, Power cav ≤ P mo ( n ) ⎭ (5.18) T mo and P mo are machine tool torque and power curves, respectively, and they are provided by the machine tool manufacturer (see Figure 5.11). T1 T2 Scale n1 n2 Figure 5.11 : Typical torque-power characteristics of a machine tool. Subscript "cav" indicates cumulative average which is mathematically defined as: q cav = t 2 t 2 [ q ( φ min ) + q ( φ max ) ] ⁄ 2 , (5.19) where q={Torque,Power}. Standard piecewise logarithmic relation of motor torque/power with spindle speed can be expressed as: q mo ( n ) = q 2 ⋅ e – log ( q 2 ⁄ q 1 ) log ( n 2 ⁄ n ) ------------------------------------------------------log ( n 2 ⁄ n 1 ) 143 , (5.20) Chapter 5. Process Optimization where q={T,P}, and {n: n 1 < n < n 2 } is the spindle speed. Note that since torque/power curves are logarithmic, the constraints associated with these two quantities become non-linear inequality constraints. 5.3.3. Chip Thinning Constraint When the width of cut is smaller than the radius of the tool (b < 0.5), maximum chip thickness will never reach the commanded feed per tooth value in up/down milling. Moreover, for operations with ball end or face mills, chip starts to thin as depth of cut decreases as shown in Figure 5.12. Reduced chip thickness not only results in reduced material removal rate but it also makes cutting edge rub on workpiece rather than cut it. In order to avoid these problems, the feedrate has to be modified as cutting conditions, depth and width of cut, change. z R c κ a x hex Figure 5.12 : Chip thickness for a face mill with a round insert. Earlier in Eq (3.27), the chip thickness for 3-axis milling was defined as: h ex = c. ( f xy ⋅ sin κ ⋅ sin φ – f z ⋅ cos κ ) , (5.21) where c is the feed per tooth, f xy and f z are the components of a unit vector in the direction of resultant feed and tool axis, κ is the axial immersion angle measured from x axis in a clock-wise direction and φ is the radial immersion angle measured from y axis also in a clock-wise direction. Depending on the depth and width of cut, maximum chip thickness is reached at certain φ and κ 144 Chapter 5. Process Optimization chip chip values denoted by φ max and κ max , respectively. There are two distinct cases for 3-axis milling operations: (a) Ramping Up or Planar Milling ( f xy ≥ 0, f z ≥ 0 ) : Only the front of the cutter ( 0 ≤ φ ≤ π ) can chip engage in a cut. In this case, maximum chip thickness is obtained when φ max = π ⁄ 2 , chip κ max = π ⁄ 2 , or at values closest to these. This rule can be generalized as follows: Down Milling ( φ ex = π ) If (φ st > π ⁄ 2 ), Up Milling ( φ st = 0 ) chip chip φ max =φ st If (φ ex < π ⁄ 2 ), φ max =φ ex chip chip If (φ st ≤ π ⁄ 2 ), φ max =π ⁄ 2 If (φ ex ≥ π ⁄ 2 ), φ max =π ⁄ 2 chip and axial immersion angle at κ max = κ ( z = a ) if ( a < R ) . (b) Ramping Down ( f xy ≥ 0, f z < 0 ) : Both front and back of the cutter can engage in a cut. At the chip front of the cutter ( 0 ≤ φ ≤ π ), maximum chip thickness is reached when φ max = π ⁄ 2 , and the axial immersion angle is found as follows: dh ex ----------- = c ⋅ ( f xy ⋅ cos κ ⋅ sin φ chip max + f z ⋅ sin κ ) = 0 , dκ chip chip κ max ⎛ – f xy ⋅ sin φ max ⎞ = atan ⎜ ----------------------------------⎟ . fz ⎝ ⎠ (5.22) At the back of the cutter (π ≤ φ ≤ 2π ); on the other hand, maximum chip thickness is reached chip chip chip when φ max = π or φ max = 2π , and κ max = 0 , or at values closest to these. If the desired maximum chip load is denoted by h max and the necessary feedrate to achieve this maximum chip load is denoted by c max , feedrate of the tool is updated by the following relation: h max c max = --------------- , h factor 145 (5.23) Chapter 5. Process Optimization chip chip chip where h factor = f xy sin ( κ max ) sin ( φ max ) – f z cos ( κ max ) is the chip thinning factor. Note that the chip thinning factor must be limited with an upper value in order to avoid extreme magnification of the feedrate when very small depth and width of cut are present. 5.4. Verification for Pre-Process Optimization 5.4.1. Case Study I An example machining operation was simulated using the proposed analytical algorithm. 20 [mm] diameter, 4-flute cylindrical end mill was used to go through an immersion zone of φ st =0, φ ex =30 [degrees] with a feedrate of 0.05 [mm/tooth.rev] at 5000 [rpm]. Al 7075 was selected as workpiece with the following tangential cutting and edge force coefficients: K tc = 796 [N/ 2 mm ] and K te = 27.7 [N/mm]. Integration boundaries at different depths of cut are presented in Figure 5.13.a along with simulated tangential cutting forces. The variation in maximum tangential forces with change in depth of cut is obtained using the efficient analytical algorithm and the result is presented in Figure 5.13.b. Maximum tangential cutting force can be effectively used to determine optimum radial width and axial depth of cut. One constraint of the process is the torque and power limits of the machine tool, which is given in Figure 5.14. By varying both axial depth of cut (a) and normalized radial depth of cut (b) for a fixed feedrate and spindle speed, and calculating maximum tangential force which leads to torque and power demands, the safe operation region can be identified as shown in Figure 5.15. Any point picked in this feasible region will not violate the machine tool’s torque and 2 power characteristics. Cutting force coefficients used in this example are K tc = 600 [N/ mm ] and K te = 42 [N/mm]. 146 Ft [N] 30 25 20 15 10 5 -100 -100 80 -300 0 -100 -200 -300 -400 a = 30 mm 60 -200 -400 a = 15 mm 40 -300 0 φ0 φ1 20 -200 -400 a = 10 mm 0 -300 0 φ0 φ1 φ0 φ1 -200 -400 a = 5 mm Ft [N] φ0,φ1 [deg] 30 25 20 15 10 5 0 -100 Ft [N] φ0,φ1 [deg] 30 25 20 15 10 5 φ0 φ1 Ft [N] φ0,φ1 [deg] 30 25 20 15 10 5 φ0,φ1 [deg] Chapter 5. Process Optimization 100 120 140 160 180 0 50 100 φ(t) [deg] 150 200 250 300 350 φ(t) [deg] (a) 700 |Ft,max| [N] 600 500 400 300 200 100 0 0 5 10 15 20 25 30 35 40 a [mm] - depth of cut (b) Figure 5.13 : (a) Integration boundaries with various depths of cut, simulated tangential forces, (b) analytically obtained maximum tangential forces. 147 Chapter 5. Process Optimization 24.12 20.12 35.3 Power [hp] Torque [N-m] 95.5 11 100 1500 3500 5000 1.34 13000 Spindle Speed [rpm] Figure 5.14 : Torque - Power curves of the machine tool. Similar analysis can also be conducted from stability point of view. The stability of the process is dependent on both radial width of cut and axial depth of cut as explained before. Using stability lobes, combination of radial width of cut and axial depth of cut for each spindle speed is obtained. In this example, dynamic parameters in two orthogonal directions are taken from reference [61]: natural frequency ω nx = 600 [Hz], ω ny = 660 [Hz]; stiffness k x = k y = 5.6 × 10 6 [N/m], damping ratio ζ x = ζ y = 0.035 . Using these modal parameters, sample stability diagram for full immersion is obtained and presented in Figure 5.16. It is evident that cutting conditions that fall under stability pockets provide efficient machining; therefore, spindle speed that has the highest possible value yet that is still in the range of machine limits (Figure 5.14) is selected for analysis. For this example case, this value is 12600 [rpm]. The stability limit along with torque and power limits at 12600 [rpm] are plotted in Figure 5.15. From this figure, it is clear that not only stability of the process limits the selection of optimum process parameters, but also torque and power demands put additional boundaries on the feasibility region. Since the main objective for the selection of parameters is to increase production by maximizing MRR, the MRR is calculated 148 Chapter 5. Process Optimization using Eq (5.1) and plotted in Figure 5.17. The depth of cut at which maximum MRR can be achieved is around 5 [mm] and for this particular case, it is limited mainly by the torque limit of the machine. Note that stability of the process is not function of feedrate when linear force model is used; whereas, torque and power are. Although the stability limit in Figures 5.15 and 5.17 will remain the same regardless of the feedrate, both the torque and power limits will vary. If the feedrate is low, then the MRR is limited mainly by the stability; however, as the feedrate gets higher, torque and power demands will increase, hence the torque and power limits will become more conservative in determining the suitable depth of cut - radial width of cut pair. Stability Limit Torque Limit Power Limit 1 c = 0.25 mm/tooth.rev n = 12600 rpm D = 20 mm N=3 0.8 b [-] 0.6 Feasible Region (safe to operate) 0.4 0.2 0 1 2 4 6 8 10 12 a [mm] Figure 5.15 : Radial and axial depths of cut curves for upmilling, feasible region (12600 [rpm]). 149 Chapter 5. Process Optimization 4 b = 1 (slotting) N=3 a [mm] 3 2 1 Stable Region 0 2 4 6 8 10 12 14 16 18 20 3 n x 10 [rpm] Figure 5.16 : Stability lobes for full immersion milling case. 1200 Stability Limit Torque Limit Power Limit MRR [cm3/min] 1000 c = 0.25 mm/tooth.rev n = 12600 rpm D = 20 mm N=3 800 Max. MRR 600 400 Feasible Region (safe to operate) 200 0 1 2 4 6 8 10 12 a [mm] Figure 5.17 : MRR curves of up-milling based on stability, torque and power limits. Note that the second stability lobe located at 6300 [rpm] is not selected as a design point because stability lobes significantly decrease as spindle speed is lowered. Although the torque and power 150 Chapter 5. Process Optimization demands are lower, process parameters are limited by the stability of the system, which can be seen clearly in Figure 5.18. Stability Limit Torque Limit Power Limit 1 c = 0.25 mm/tooth.rev n = 6300 rpm D = 20 mm N=3 0.8 b [-] 0.6 0.4 Feasible Region (safe to operate) 0.2 0 0 2 4 6 8 10 12 14 16 18 20 a [mm] Figure 5.18 : Radial and axial depths of cut curves for upmilling, feasible region (6300 [rpm]). 5.4.2. Case Study II Proposed optimization algorithm was used for process planning of an industrial part, the front face of a helicopter gear box cover. A two fluted helical end mill with 10 [mm] diameter and a 25 [mm] indexable cutter with 2 inserts were used. The original two-sided part, solid model of front face during machining, and a picture of machined part are shown in Figure 5.19. In order to prevent tool holder - workpiece collision and allow enough room for chip evacuation, the tool over hang was selected as 45 [mm]. Dynamic characteristics of both tools were identified by impulse model testing, and measured frequency response functions (FRFs) are presented in Figures 5.20 and 5.21. The work material was Al 6061 and the experimentally identified cutting coefficients are given in Table 5.3. 151 Chapter 5. Process Optimization (a) (b) (c) Figure 5.19 : Gear box cover: (a) complete solid model (double sided machining), (b) solid model for machining front features, (c) real cut part with front features. 2.5 x 1e-6 2.5 Model Experimental X - Direction 2 Magnitude [m/N] Magnitude [m/N] Model Experimental Y - Direction 2 1.5 1 1.5 1 0.5 0.5 0 1000 x 1e-6 1500 2000 2500 3000 3500 Frequency [Hz] 4000 4500 5000 0 1000 1500 2000 2500 3000 3500 Frequency [Hz] 4000 4500 5000 Figure 5.20 : Transfer functions of 10 [mm] endmill in x and y directions (HSK63A-Shrink Fit). 152 Chapter 5. Process Optimization 2.5 x 1e-7 3 Model Experimental X - Direction Model Experimental Y - Direction 2.5 Magnitude [m/N] 2 Magnitude [m/N] x 1e-7 1.5 1 2 1.5 1 0.5 0.5 0 0 500 1000 1500 2000 2500 3000 Frequency [Hz] 3500 4000 4500 5000 500 1000 1500 2000 2500 3000 Frequency [Hz] 3500 4000 4500 5000 Figure 5.21 : Transfer functions of 25 [mm] indexable cutter in x and y directions (HSK63ACorogrip). Tool Type 10 [mm] 25 [mm] Table 5.3 : Cutting force coefficients of Al 6061 K te K re K tc 2 [N/mm] [N/mm] [N ⁄ mm ] 38.38 -9.792 403.396 21.34 23.141 662.296 K rc 2 [N ⁄ mm ] 98.911 72.86 Based on material data and dynamic characteristics of tool-tool holder-machine tool assembly, stability lobes were calculated (see Figure 5.22). For 10 [mm] end mill, the maximum spindle speed of the machine tool, 20000 [rpm], was selected, whereas, for the 25 [mm] indexable cutter, wide enough stability pocket around 17000 [rpm] was selected. In order to give a general idea of how stability lobes vary as a function of immersion conditions, Figure 5.23 is provided. It can be observed that the location of a high stability pocket around 20000 [rpm] propagates as immersion (or width of cut) is varied. Since the objective is to maximize the production, the MRR variation as a function of immersion and spindle speed is plotted in 153 Chapter 5. Process Optimization Figure 5.24, which justifies the argument that high productivity is achieved around pockets of stability with the highest possible spindle speed. 15 2.5 Slotting Slotting a - Depth of Cut [mm] a - Depth of Cut [mm] 2 1.5 1 10 5 0.5 0 0 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 n - Spindle Speed [rpm] x 1e4 n - Spindle Speed [rpm] x 1e4 (a) (b) Figure 5.22 : Stability lobes for full immersion milling (b=1) of (a) 10 [mm] end mill, (b) 25 [mm] indexable cutter. Cutting coefficients from Table 5.3 are used. 154 2.1 Chapter 5. Process Optimization 16 14 12 a [mm] 10 8 6 4 2 0 0.5 0 10000 12000 14000 16000 18000 20000 22000 b [-] 1 n [rpm] Figure 5.23 : Stability lobes for various immersion cases, 10 [mm] end mill. Cutting coefficients from Table 5.3 are used. 160 140 MRR [cm3/min] 120 100 80 60 40 0 0.5 20 10000 12000 14000 16000 18000 20000 22000 1 b [-] n [rpm] Figure 5.24 : MRR characteristics of 10 [mm] end mill based only on stability in Figure 5.23 with feed variation in Figure 5.26.d where h max =0.1 [mm/tooth.rev], and N=2. 155 Chapter 5. Process Optimization Process planning charts were constructed at selected spindle speeds based on analysis summarized in Case Study I. Both stability and torque/power curves of the machine tool (Figure 5.25) were used as constraints. Design spaces with supplementary charts are given in Figures 5.26 and 5.27. Charts 5.26-5.27.a present the design spaces that were used to select proper width and depth of cut combination to determine machining strategy for optimum production, i.e. maximum material removal rate (see Charts 5.26-5.27.b). For 10 [mm] tool, deep depths of cut with reduced width of cut were used; whereas, for 25 [mm] tool shallow depths of cut with increased widths of cut were preferred for part programming to maximize the MRR. Charts 5.26-5.27.c provide alternative representations for the design space, start angle of cut versus depth of cut. Note that since down milling was used, exit angle φ ex was always equal to 180 [degrees]. 24.81 20.12 30.0 11.67 Power [HP] Torque [N-m] 95.4 17.6 8.8 200 1500 4500 10000 2.48 20000 Spindle Speed [rpm] Figure 5.25 : Mori Seiki SH403 - Torque/Power characteristics. 156 Chapter 5. Process Optimization 180 150 1 φs t [deg] 0.8 b 0.6 0.4 0.2 0 2 4 6 8 10 12 14 100 50 0 15 2 4 6 8 a [mm] 120 100 80 2 4 6 12 14 15 (c) 140 Feed Rate [mm/th.rev] MRR [cm3/min] (a) 60 10 a [mm] 8 10 12 14 15 0.14 0.13 0.12 0.11 0.1 0.1 0.2 0.3 0.4 0.5 a [mm] 0.6 0.7 0.8 0.9 1 b (b) (d) Figure 5.26 : Process planning charts for 10 [mm] end mill, n=20000 [rpm], N=2: (a) stability limited design space, (b) objective function, (c) start angle for down milling, (d) feedrate 1 φs t [deg] b 0.8 0.6 0.4 180 150 100 50 0.2 0 7 8 9 10 11 12 13 14 0 15 7 8 9 10 a [mm] (a) 800 600 400 8 9 10 11 12 13 14 15 Feed Rate [mm/th.rev] MRR [cm3/min] 1000 7 12 13 14 15 (c) 1200 200 11 a [mm] 0.4 0.36 0.32 0.28 0.24 0.2 0 0.1 0.2 0.3 0.4 0.5 a [mm] b (b) (d) 0.6 0.7 0.8 0.9 1 Figure 5.27 : Process planning charts for 25 [mm] indexable cutter, n=17000 [rpm], N=2: (a) torque limited design space, (b) objective function, (c) start angle for down milling, (d) feedrate variation. Finally, Charts 5.26-5.27.d show the variation of feedrate as the width of cut changes. This chart ensures that the maximum chip load specified by the process planner is always maintained even 157 Chapter 5. Process Optimization when the width of cut changes (for details, refer to the "Chip Thinning Constraint" section of this chapter). The gear box cover was successfully cut without presence of chatter vibrations and spindle stall. The workpiece was mounted on a table dynamometer to measure cutting forces acting on the cutter in machine coordinates. For further comparison, process was also simulated in virtual environment using developed closed form algorithms. 5.5. Post-Process Optimization Post-process optimization is an effective tool to optimize cutting parameters of an existing NC program; therefore, parameters such as depth and width of cut, which will effect the existing planned tool path, are not changed whereas performance parameters such as feedrate and spindle speed are optimized. CAD Model Process Planning (CAM) Post-Process Optimization Output NC Code Constraints: Machining Machine Tool Optimized NC Code Figure 5.28 : Optimized NC code generation flow chart with Post-Process Optimization module. The flow chart for process optimization is given in Figure 5.29. The idea here is to maximize feedrate and spindle speed by virtually simulating the process after planning stage and recommending new feedrates and spindle speeds by checking the process outputs against both machining and machine tool constraints. Some important machining constraints are form errors left on the part, 158 Chapter 5. Process Optimization maximum chip load and maximum resultant force on the tool, where as important machine tool constraints are the torque and power characteristics. START END Figure 5.29 : Process optimization flow chart. In milling, sinusoidal variation of the chip thickness as the cutter rotates causes cutting forces, shown in Figure 5.30, to oscillate between maximum and minimum points. In terms of optimization, only extremum points limit the range of cutting parameters during process planning. In the previous chapter, it has been shown that an analytical solution for cutting forces is not always possible. In the same sense, the maximum and minimum of cutting forces cannot always be obtained analytically and a numerical method is needed. Note that the tasks of maximization and minimization are trivially related to each other, since the maximum of function g is same as the minimum of -g. In this section, both maximization and minimization problems will be expressed as a minimization problem, i.e. if F ( φ ) represents varying cutting forces, then Minimum Force ⇒ min ( F ( φ ) ) . Maximum Force ⇒ min ( – F ( φ ) ) 159 Chapter 5. Process Optimization F(φ) φmax φmin No of flutes = 3 φp=120 240 360 φ [deg] Figure 5.30 : Sample force variation for a uniform pitch cutter. There are numerous numerical solutions to tackle the minimization problem. A commonly used method for one-dimensional minimization (to minimize a function of one variable) without calculation of the derivative is called Brent’s Method [92]. This method takes three initial immersion 0 0 0 points ( φ 1, φ 2, φ 3 ) that potentially bracket the minimum, i.e., satisfy the following relations: 0 0 0 0 F ( φ 2 ) < F ( φ 1 ) & F ( φ 2 ) < F ( φ 3 ) , then it estimates the next minimum point by inverse parabolic interpolation, i.e., finding the location of extremum of the parabola defined by these three points. Inverse parabolic interpolation is derived as follows. First, a parabolic curve is drawn passing 0 0 0 through three points, namely end points φ a = φ 1 and φ c = φ 3 , and a midpoint φ b = φ 2 . Note that there are two possibilities, one with a minimum and one with a maximum as seen in Figure 5.31. 160 Chapter 5. Process Optimization F(φ) min parabola max parabola φb φa φc φ Figure 5.31 : Two possible quadratic functions. Either parabola can be described as: 2 F ( φ ) = α2 φ + α1 φ + α0 . (5.24) Substituting given immersion angles, a set of equations are obtained: 2 F ( φa ) = α2 φa + α1 φa + α0 2 F ( φb ) = α2 φb + α1 φb + α0 . (5.25) 2 F ( φc ) = α2 φc + α1 φc + α0 The extremum of the parabola can be found by equating the first derivative of function F to zero: d F ( φ ) = 0 ⇒ 2α φ + α = 0 , 2 1 dφ (5.26) and from Eq (5.26), φ is obtained as: –α1 φ = --------- . 2α 2 161 (5.27) Chapter 5. Process Optimization Combining Eqs (5.25) and (5.27), and solving system of linear equations for α ’s, the formula for the abscissa φ for the extremum of the parabola is obtained as: 2 2 1 ( φb – φa ) [ F ( φb ) – F ( φc ) ] – ( φb – φc ) [ F ( φb ) – F ( φa ) ] φ = φ b – --- ------------------------------------------------------------------------------------------------------------------------------------------ . 2 ( φb – φa ) [ F ( φb ) – F ( φc ) ] – ( φb – φc ) [ F ( φb ) – F ( φa ) ] (5.28) Since F( φ ) is not known to be minimum or maximum, minimization scheme cannot solely depend on Eq (5.28). In this case, the signum of α 2 can be used for final decision such that if it is positive, the extremum is minimum, if negative, one has a maximum. Inverse parabolic interpolation for two iteration steps is depicted in Figure 5.32. F(φ) first parabola second parabola F(φ) φ Iteration Steps 0 φ01 1 φ11 φ21 2 φ02 φ12 φ22 φ03 φ13 φ23 Figure 5.32 : Inverse parabolic interpolation. When the function is not co-operative, for example, when three points are co-linear ( φ a = φ b = φ c ) so that denominator of Eq (5.28) becomes zero or there are sharp changes in the function value then the method switches to sure but slow technique, like golden search [92]. Golden search is analogous to the bisection method and used for finding the minimum of a function using a brute 162 Chapter 5. Process Optimization force approach. Bearing in mind that golden ratio is equal to G R = 0.61803, given a bracketing triplet of immersions, the next point to be tried is placed at a fraction (1- G R ) into the larger of the two intervals measuring from the middle point φ b : φ new ⎧ ⎪ ( φ b – Δ φ ) if ( ( φ b – φ a ) > ( φ c – φ b ) ) = ⎨ ⎪ ( φ b + Δφ ) if ( ( φ c – φ b ) > ( φ b – φ a ) ) ⎩ ⎫ ⎪ ⎬, ⎪ ⎭ (5.29) where Δφ = ( 1 – G R ) × max { ( φ b – φ a ), ( φ c – φ b ) } . Then the update decision is again made based on bracketing the minimum, i.e., if F ( φ new ) < F ( φ b ) then φ new replaces the midpoint and φ b becomes an end point, else if F ( φ b ) < F ( φ new ) then φ b remains midpoint with φ new replacing one of the end points, φ a or φ c . Either way the width of the bracket, ( φ c - φ a ), reduces until a desired tolerance is achieved. In order to account for all local extremums, one tooth period (or one spindle period for non-uniform pitch cutter) is divided into three sections, extremums for each range are determined, and finally the global solution is selected as the minimum of all. The bracketing condition for the ini0 0 0 tial points (φ 1, φ 2, φ 3 ) is satisfied through an iterative algorithm which starts from any set of points and moves in the downhill direction until the downhill trend of the function stops so that the minimum is bracketed. 5.5.1. Feedrate Optimization In the previous chapter, it is shown that majority of the process outputs becomes linearly dependent on the feedrate in the following form provided that cutting force coefficients are independent of the feedrate: q q F q ( φ ) = A 0 ( φ ) + A 1 ( φ ) .c , ( q ≡ x, y, z, t, r, trq, pwr), 163 (5.30) Chapter 5. Process Optimization where subscripts x, y, z, t, r are cutting forces in the feed, perpendicular to feed, axial, tangential, and radial directions, trq and pwr are spindle torque and power, respectively. Force coefficients q q A 0 and A 1 are functions of process type, tool geometry, and time (or angular position). In Chapter 3, a numerical search technique has been presented to obtain the angular position at which a q process output becomes maximum, which is denoted as φ max in this section. At the most critical instance, process output -Eq (5.30)- reaches the user-defined maximum, F q, max : q F q ( φ max ) = F q, max , (5.31) for which maximum allowable feedrate can be solved as: q q c max q F q, max – A 0 ( φ max ) -. = --------------------------------------------q q A 1 ( φ max ) (5.32) There are two other important process outputs that do not have as straight forward a relationship with feedrate: 1) The resultant force causing bending moment; 2) the form error due to static deflection of the tool. The resultant force in xy-plane is previously defined as: F res ( φ ) = res res 2 res A0 ( φ ) + A1 ( φ ) c + A2 ( φ ) c . (5.33) Similar to Eq (5.31), maximum resultant force is equated to user-defined maximum F res, max , i.e. res F res ( φ max ) = F res, max , and the maximum feedrate is obtained by solving quadratic equation as: res res c max res res 2 res – A 1 ( φ max ) ± Δ = ---------------------------------------------, res res 2A 2 ( φ max ) res 2 res (5.34) res res where Δ = [ A 1 ( φ max ) ] – 4 ⋅ [ A 0 ( φ max ) – ( F res, max ) ] ⋅ A 2 ( φ max ) . From two possible solutions in Eq (5.34), smaller positive real answer is selected. 164 Chapter 5. Process Optimization Form errors in milling result from the static tool deflection at the entry ( φ = 0 ) and exit ( φ = π ) points of the cutter due to forces only in y-direction. Generating the surface becomes complex when the end mill has helical flutes and there are more than one engagement zone (CEF). At any tool position, tool deflection in y-direction is calculated using cantilever beam approximation that has a cylindrical cross section with an effective radius R eff = 0.8R tool , where 0.8 is the approximate scale factor due to flutes [94]. An example case is given in Figure 5.33. In order to find the actual form error, tool deflection must be calculated not necessarily at the tool tip but at the axial point where a helical flute generates the surface. The axial location of the surface generation point is referred as z gp . Total tool deflection reflected at the surface generation point can be found by adding up the deflections at z gp caused by y direction forces of each CEF: # of CEFs δ = ∑ δi , (5.35) i=1 where deflection is defined by [95] as: ⎛ ⎜ ⎜ δi = ⎜ ⎜ ⎜ ⎝ 2 F y, i .z i ---------------- ( 3z i – z gp ), if 0<z gp < z i 6EI , 2 F y, i .z i ---------------- ( 3z gp – z i ), if z i < z gp 6EI (5.36) where E is the Young’s Modulus and I is the area moment of inertia of the tool. Cutting forces acting on the tool are concentrated at the middle point of their respective CEFs. When Eq (5.36) is analyzed, it can be seen that the deflection of the tool is directly proportional to the magnitude of force in y-direction; therefore linear relationship between feedrate and deflection is still applicable: 165 Chapter 5. Process Optimization dfl dfl z - Axial Depth of Cut δ ( φ ) = A0 ( φ ) + A1 ( φ ) ⋅ c . FL FL UT E UT E -2 -3 FL UT E (5.37) CEF 3 -1 CEF 2 CEF 1 0 δ Fy,3 Fy,1 φ - Immersion Angle z1 π z3 z2 zgp Fy,2 Figure 5.33 : CEFs and static tool deflection at surface generation point, z gp . dfl The tool position, φ max , at which maximum surface deflection is generated can be calculated; dfl dfl however, it is not trivial to come up with deflection coefficients A 0 and A 1 explicitly. A more practical approach to obtain these two coefficients is using linear interpolation. Given two arbitrary feedrates, c 1 and c 2 , the corresponding deflections δ 1 and δ 2 are calculated. Utilizing those two points, maximum allowable feedrate for a user-defined deflection constraint δ max can be found as: δ max – δ 1 dfl c max = c 1 + ( c 2 – c 1 ) ---------------------- . δ2 – δ1 (5.38) dfl Note that if c max is solved as a negative number, it means that edge forces are already enough to create a deflection more than the specified deflection constraint δ max ; therefore, under no circumstances, can such a constraint be satisfied. 166 Chapter 5. Process Optimization 5.5.2. Feedrate and Spindle Speed Optimization As mentioned earlier, highly efficient manufacturing can only be achieved by the maximizing material removal rate (MRR) - Eq (5.1). Cutter locations are generated when process planning is completed at the CAM stage; therefore, it is required that any optimization taking place afterwards shall not change any physical cutter location information. Regarding axial depth of cut, width of cut, and number of flutes as system input parameters, i.e., constants, the number of free parameters (design variables) in the MRR finally reduces from five to two: feedrate and spindle speed. Dropping the constants from the MRR, a simplified objective function called "reduced MRR" denoted by rMRR is defined for optimization problem as: f obj = rMRR = c ⋅ n . (5.39) In addition to constraints presented previously, which are functions of only feed, two more constraints that are a function of spindle speed can be added. The first constraint is torque-power characteristics of the machine, which are given in Eqs (5.18) and (5.20), and the second constraint is the chatter stability. Stability lobes can be generated for each process using dynamic characteristics of the tool/workpiece. As mentioned before, the depth of cut (a) is already provided in the NC code. The fixed depth of cut corresponds to a horizontal line in stability lobes shown in Figure 5.34. Regions that fall under the curve ensure that the process is stable, therefore, their respective lower and upper spindle speeds are noted to be used as a constraint on spindle speed (design variable): n ≤ n max, i ⎫ ⎬, – n ≤ – n min, i ⎭ 167 (5.40) Chapter 5. Process Optimization where i = 1,2,...,M, and M is the number of stable regions at depth of cut a. There are two special cases for the selection of the spindle speed range. If the depth of cut in NC code is less than the minimum of the stability lobe denoted by a lim, min , then full spindle speed range of torque/power curve is to be used. However, if the depth of cut (a) is more than the peak point of the stability border, a lim, max , then the process is already unstable, therefore, optimization will not be applicable. For optimization to be successful, every effort must be made to select stable depths of cut. In the complete virtual machining package, necessary stability lobes are provided prior to process planning. High Torque Region Low Torque Region a Stable Region nmax,2 nmin,1 alim,min nmin,2 Unstable Region nmax,1 Depth of Cut (a) alim,max nmax nmin Spindle Speed (n) Figure 5.34 : Sample stability lobes. The described non-linear optimization problem with linear and non-linear inequality constraints is solved by routines of MATLAB’s optimization tool-box [97]. The constrained nonlinear optimization (nonlinear programming) is employed by sequential quadratic programming (SQP). Two main stages of implementation, i.e. the Hessian matrix update, and quadratic programming (QP) 168 Chapter 5. Process Optimization have the following properties: the Hessian matrix is updated using the BFSG method and quadratic programming method is an active set strategy [97]. Graphical representation of the design space is obtained for a case that has only torque and power constraints and plotted in Figure 5.35. Any cutting condition selected within the feasible region will respect torque and power characteristics of the machine; however, the most optimum solution (maximum MRR) is obtained only at the limits of the constraints which is marked with a star in the figure. The coordinates for this point are c = 0.2088 [mm/tooth.rev], and n = 5250 [rpm]. Sudden changes of the constraints like the torque constraint around 5500 [rpm] increases the number of iterations before the algorithm converges, or even sometimes results in non-convergence. Nonetheless, convergence is ensured by scanning each spindle speed range where the torque and power curves change behavior and the biggest of all is selected as the final optimum solution. 5.6. Filtering "Unrealizable" Feed Commands Once optimization is completed, a new set of feedrates is generated. Optimized feedrates might fluctuate sharply depending on the rate of change of intersection between the cutter and workpiece, i.e., how much the workpiece geometry varies along the tool path. For example, if the tool goes through an existing hole on the workpiece, the engagement changes fast leading to rapid feed adjustment within a small travel distance. A sample output of Post-Process Optimization engine is given in Figure 5.36. Sudden changes in the feedrate within short tool motion, marked in the figure, cannot be achieved due to limited acceleration and deceleration of the motors. In this section, such feed commands are described as "unrealizable". Every machine tool has acceleration and deceleration characteristics based on its motor torque capability. NC unit of the machine tool ensures that variations in commanded feedrate do not overload motors with excessive accelera- 169 Chapter 5. Process Optimization tion or jerk. Since feed motion planning is performed by machine tool’s NC unit, the objective of this section is not to duplicate feed motion planning but to come up with a set of rules to filter out sharp feedrate changes obtained during constraint-based optimization. 12 -2 T co ns tra int 50 0 ) D th irec e ob tion je of ct ive inc fu rea nc se tio fo n r (e M RR 00 15 um tim 00 n - Spindle Speed x 1e3 [rpm] bj, j ,- Op 10 int 0 9 8 fo f ob f obj, tra ns f obj,-50 10 f obj,- P co 11 fo bj, -20 00 Infeasible Region 7 6 5 4 Feasible Region 3 2 1 Pconstraint =P(φmax)-Pmo(n)=0 Tconstraint =T(φmax)-Tmo(n)=0 0.05 0.1 0.15 c - Feed rate [mm/tooth.rev] 0..2 0.25 Figure 5.35 : Graphical representation of design space. Star is located at c=0.2088 [mm/tooth.rev] and n=5250 [rpm]. 5.6.1. Classification of Feed Blocks In order to automate the filtering process, feed steps are classified under four main groups depending on start ( f s ), commanded ( f c ), and end ( f e ) feedrates. All of the groups are shown in Figure 5.37. Black dots represent nodes at the start and end points of each block. Previous and next feed blocks can only be attached at these nodal points. Dashed lines at both ends indicate the 170 Chapter 5. Process Optimization allowable trends for previous and next feed blocks. Due to such restrictions, there are a limited number of feed block combinations as presented in Figure 5.38, where circles represent different feed types. Optimized Feedrate 15 ! 14 ! f - Feedrate x 1e3 [mm/min] 13 12 11 10 ! 9 ! 8 7 ! 6 5 4 0 20 40 60 80 100 120 140 160 180 s - Path Length [mm] Figure 5.36 : An example of optimized feedrate with sharp variations. Feed type stated in the middle circle can only be preceded by the types shown on its left; and succeed by the ones on its right. Each type has also its unique acceleration and deceleration properties: TYPE 1 has both acceleration and deceleration, TYPE 2 has only acceleration, TYPE 3 is constant velocity block, whereas TYPE 4 has only deceleration. TYPE 1 is further divided into three subgroups by considering the relationship between start and end feed commands. Four types and their sub-groups will later be used to build rules for filtering unwanted feed commands. 171 Chapter 5. Process Optimization TYPE 1 fc fs fe (fe>fs) fs fc fc (fs>fe) fe fs fe= fs (1.a) (1.b) (1.c) TYPE 2 TYPE 3 TYPE 4 fc fs fe fc fs fc fe fs fe Figure 5.37 : Types of feed command used for filtering. 2 3 2 1 3 1 4 3 1 1 2 3 4 3 4 1 2 2 4 4 Figure 5.38 : Allowable combination of different feed block types (TYPE numbers are indicated in circles). Feedrates passed from optimization module are first scanned and classified according to four main types described above. Total travel length is also calculated using tool positions and stored for each feed block. Figure 5.39 shows necessary parameters for the ith feed block. 172 Chapter 5. Process Optimization Optimized Feedrate f - Feedrate fci fsi ith block fi 0 s1 si-1 sN si s - Path Length Figure 5.39 : Parameter definitions for a single feed block i. 5.6.2. Filtering Rules For a given step feed command between initial and final feed values, machine tool uses certain amount of time to accelerate to a commanded feed value and decelerate to a final feed while cruising at a constant velocity in between. Kinematic profiles for a constant acceleration and deceleration case are shown in Figure 5.40. The objective of motion planning is to achieve commanded feedrate and stay at constant velocity only if possible. The minimum distance of travel required for such motion is obtained by calculating the area under the velocity curve when T1 = i i fc – fs ⁄ A , T2 = 0 , T3 = i i f c – f e ⁄ D as: i 2 i 2 i 2 i 2 ( fc ) – ( fs ) ( fc ) – ( fe ) L i = -------------------------------- + -------------------------------- , 2A 2D 173 (5.41) Chapter 5. Process Optimization where A and D are acceleration and deceleration limits, respectively. If an acceleration or deceleration stage does not exist for the feed block, then its corresponding term should not be considered in the above equation. In order for a feed block i where i = 1, …, N to be realizable, the actual path length has to be long enough to accommodate minimum length requirement; ( si – si – 1 ) ≥ Li . (5.42) The equality case of the above inequality represents a motion with no constant velocity section, i.e., T 2 = 0 . si Position [mm] si-1 Time [sec] fci Velocity [mm/sec] i fs T1 T2 T3 fei Time [sec] Acceleration [mm/sec2] A Time [sec] -D Figure 5.40 : Constant acceleration kinematic profiles. When Eq (5.42) is not satisfied, feed block needs to be updated by reducing commanded feedrate to either preceding or succeeding feedrate depending on the feed block type. This update will not only affect the current block but it will also change the type and limiting feed values of neighboring blocks; therefore, a set of update rules has to be established. 174 Chapter 5. Process Optimization Consider that all feedrate commands are stored in an array called the Feed-array, and i represents the position of current feed block to be processed in that array. The Feed-array is first scanned to classify each feed block using the convention presented earlier. Later, each feed block is checked to determine whether the minimum length requirement given in Eq (5.42) is satisfied. In Table 5.4, all possible feed block configurations are shown and rules of update are given for each feed block combinations. In the "Update Rule" section, circles represent types of feed blocks. Current feed block with index i is represented with a shaded circle. When the current feed block cannot be realized due to the minimum length requirement, the update rule is initiated and the current feed block is converted into a different feed block by merging neighboring block or blocks depending on the formula of update rule. The type of the new feed block can be the same as one of the blending feed blocks or it can be a different type. Update rules are classified under three major categories: FORWARD, BACKWARD, and FULL UPDATEs depending on the involvement of succeeding, preceding, or both blocks. FORWARD UPDATE reduces the current feedrate to the value of succeeding block. Since feed blocks prior to the current one are not affected from this update, the current block’s index, i, remains as it is, however, the total length of the Feed-array reduces by one. The BACKWARD UPDATE is required when the current feedrate is reduced to the commanded feedrate of the preceding block. In the Feed-array, this update has to be reflected onto the prior block. Once update is completed, current block’s index, i, and the length of the Feed-array is reduced by one. FULL UPDATE is conducted only when a feed block of TYPE 1.c is not realizable. The current feedrate is dropped down to the value of neighbouring feed commands, and in the Feed-array, both preceding and successive blocks have to be updated due to 175 Chapter 5. Process Optimization i i symmetric nature of TYPE 1.c, i.e. f s = f e . To complete the update, the current block’s index, i, is reduced by one whereas the total length of the Feed-array is reduced by two. Table 5.4 : All possible update configurations with proposed update rules, shaded circle represents current step, i.e., the step being processed. Update Category: FORWARD UPDATE i i fc i fc i+1 FWD I 1.a i+1 fe FWD II +3= 1.a 2 fc i+1 fe i+1 fe fs Update Code i+1 fc i i fs Update Rule i+1 fc i+1 fe fc i fs i+1 fc Configuration i fc i fs FWD III +4= 4 1 FWD IV +3= 3 4 +4= 4 Update Category: BACKWARD UPDATE i i i-1 fc Configuration i fc fc fc i-1 fc i-1 fs i-1 fc i-1 i-1 fs BKWD I 2 i-1 fe fs Update Rule i fe fc i-1 fs i i fe Update Code i fc i fe + 1.b = BKWD II 3 1 BKWD III + 1.b = 2 4 +2= BKWD IV 2 3 +2= 3 Update Category: FULL UPDATE i fc fc i +1 i-1 fc Configuration i i fc fc i +1 i-1 fc i+1 fe i-1 fc i-1 fs i+1 fc fc i+1 fe i+1 fe i-1 fs Update Code Update Rule i-1 fs FULL I 2 + 1.c + 3 = FULL II 2 2 + 1.c + 4 = FULL III 1 3 + 1.c + 3 = 3 A simplified flow diagram of filtering algorithm is given in Figure 5.41. The Feed_Update() is a sub-routine that organizes update events. It takes the Feed-array and current block’s index, i, as 176 Chapter 5. Process Optimization Figure 5.41 : Update algorithm flow chart. 177 Chapter 5. Process Optimization input, and it returns updated the Feed-array with updated index i. In some cases, when very narrow feed blocks appear next to each other, one update will not be sufficient to obtain a realizable feed block. Moreover, two of the update options, FWD III and IV, affect the previous feed block and possibly change it from a realizable block to an unrealizable one. To make filtering more robust to handle mentioned special cases, the Feed_Update() is used as a recursive function. Using recursion, feed block is continuously updated until minimum length requirement is satisfied before program moves to the next feed block. An example of feed filtering with recursive updating is given in Figure 5.42. Optimized Optimized & Filtered Feedrate x 1e3 [mm/min] STEP I 12 STEP II STEP III STEP IV 10 8 6 B1 B1 B0 B0 B0 B2 B3 2 0 B1 B2 4 1 2 3 4 5 6 7 8 9 10 Path Length x 10 [mm] Objective: Filter feed blocks B0-B3 Path Length Path Length B1 is realizable: (i) Move to next feed block, B2 B2 is not realizable: (i) Update B2 using FWD III (ii) Update block numbers (iii) Check B1 by initiating recursion B1 is no more realizable: (i) Update B1 using BKWD I (ii) Update block numbers B0 B1 Path Length Filtering is completed! Figure 5.42 : Feed filtering example with recursive updating. When a feed block is not realizable, instead of reducing feedrate to one of the neighboring feed values, an alternative method can be developed by using minimum path length requirement given in Eq (5.41) to solve for maximum possible feedrate. In this thesis, this method is named as Maximum Feed Correction (MFC). For a known path length of feed block i, L i = ( s i – s i – 1 ) , maximum feedrate for each feed block type is calculated as 178 Chapter 5. Process Optimization i f max ⎛ ⎜ ⎜ ⎜ ⎜ = ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ i 2 i 2 ( fs ) ⋅ D + ( fe ) ⋅ A + 2 ⋅ A ⋅ D ( si – si – 1 ) --------------------------------------------------------------------------------------------------- TYPE 1 A+D i 2 ( fs ) + 2 ⋅ A ⋅ ( si – si – 1 ) ---------------------------------------------------------2 TYPE 2 . (5.43) i 2 ( fe ) + 2 ⋅ D ⋅ ( si – si – 1 ) ----------------------------------------------------------2 TYPE 4 The previous feed filtering example is solved once more by using MFC and steps are shown in Figure 5.43. The result of filtering with MFC shows that obtained feedrates tend to follow optimized feedrates more closely, which results in increased productivity compared to filtering without MFC. On the other hand, filtering without MFC generates more uniform feedrates leading to decreased feed fluctuations along the tool path especially when many more feed blocks are considered. Bearing in mind that feed fluctuations along tool path will cause feed marks on the finished product, filtering without MFC becomes the method of choice during finish milling, whereas filtering with MFC is rather preferred for roughing operations due to decreased cycle time. 179 Chapter 5. Process Optimization Optimized Optimized & Filtered Feedrate x 1e3 [mm/min] STEP I 12 STEP II STEP III STEP IV 10 8 6 B1 B1 B1 B1 B2 4 B0 B0 B0 B0 0 1 2 3 4 5 6 7 8 9 B2 B2 B3 2 10 Path Length x 10 [mm] Objective: Filter feed blocks B0-B3 Path Length Path Length B1 is realizable: (i) Move to next feed block, B2 Path Length B2 is not realizable: (i) Update B2 using FWD III (ii) Update block numbers (iii) Check B1 by initiating recursion B1 is no more realizable: (i) Update B1 using maximum feed correction Filtering is completed! Figure 5.43 : Feed filtering example with recursive updating & MFC. Optimized feedrates given in Figure 5.36 are filtered out using the proposed algorithm. 1000 2 [mm/ s ] ( ≅ 1 [g]) is used for both machine acceleration and deceleration. Filtered feedrates along with original optimized (raw) feedrates are presented in Figure 5.44 for comparison. 15 15 Optimized Optimized & Filtered 14 13 Feedrate x 1e3 [mm/min] Feedrate x 1e3 [mm/min] 14 12 11 10 9 8 7 13 12 11 10 9 8 7 6 6 5 5 4 0 Optimized Optimized & Filtered 4 20 40 60 80 100 120 140 Path Length [mm] 160 180 0 20 40 60 80 100 120 140 Path Length [mm] 160 180 (b) (a) Figure 5.44 : Filtered optimum-feed-commands for 1 [g] acc/dec; filtering (a) without MFC, (b) with MFC. 180 Chapter 5. Process Optimization 5.7. Updating the Original CL File A complete flow diagram of a CL file optimization algorithm is given in Figure 5.45. The first step is to read the original CL file and obtain cutter-workpiece engagement along the tool path. CL file provides target tool coordinates with a type of motion desired in between i.e. linear or circular. Since geometry of the workpiece might change during this point to point motion, tool path is divided into many more segments in between target points, and the cutter-workpiece engagement is obtained at this higher sampling rate. The increment for sampling is determined based on the stock size. In the flow diagram, these additional points are numbered by appending the previous CL file point number. For example, three additional points after the first coordinate but before the second coordinate in the CL file are numbered as 1.1, 1.2, and 1.3. Once the intersection between cutter and workpiece is successfully obtained, the next step is to adjust the feedrate and spindle speed based on user defined constraints. Critical (minimum/maximum) process outputs are calculated using cutter engagement features (CEFs) and necessary feedrate and spindle speed adjustments are determined. Optimization might result in highly fluctuating feedrates due to continuously varying workpiece geometry. Sharp feedrate changes may not only saturate axes motors, but they will also result in undesired feed-marks on the finished surface; therefore, optimized feedrates are first filtered. The final step is to update the original CL file with optimized and filtered feedrate and spindle speed commands. As described before due to higher sampling, there will be intermediate points in between original CL file points. If these additional points fall in between a linear tool command (G01), then all of them are added into a final optimized CL file along their optimized feeds and speeds. On the other hand if they are located in between a circular tool command (G02-G03), then intermediate points are not added between original CL file target 181 Chapter 5. Process Optimization points, but the feedrate will be updated with minimum feedrate of all intermediate points. This update strategy is depicted in Figure 5.46 for both linear and circular tool commands. GOTO / X1, Y1, Z1 FEDRAT/ F1,MMPM SPINDL/ S1,RPM,CLW GOTO / X2, Y2, Z2 GOTO / X3, Y3, Z3 FEDRAT/ F2,MMPM SPINDL/ S2,RPM,CLW GOTO / X4, Y4, Z4 GOTO / X5 , Y5, Z5 GOTO / X6 , Y6, Z6 . . . F2, S2 GOTO / X1, Y1, Z1 FEDRAT/ OFF1,MMPM SPINDL/ OFS1,RPM,CLW GOTO / X1.1, Y1.1, Z1.1 FEDRAT/ OFF2,MMPM SPINDL/ OFS2,RPM,CLW GOTO / X1.2, Y1.2, Z1.2 FEDRAT/ OFF3,MMPM SPINDL/ OFS3,RPM,CLW GOTO / X1.3, Y1.3, Z1.3 FEDRAT/ OFF4,MMPM SPINDL/ OFS4,RPM,CLW GOTO / X2 , Y2, Z2 . . . FEDRAT/ OFF12,MMPM SPINDL/ OFS12,RPM,CLW GOTO / X6 , Y6, Z6 . . . F: Feedrate S: Spindle Speed OF : Optimized F OS: Optimized S OFF : Optimized & Filtered F OFS: Optimized & Filtered S Figure 5.45 : CL file optimization routine. 182 Chapter 5. Process Optimization G01: Linear Segment (X,Y,Z)1 (X,Y,Z)1.1 (X,Y,Z)1.2 (X,Y,Z)1.3 (X,Y,Z)2 Feedrate can be varied along the path OFF1 OFF2 OFF3 OFF4 G02-G03: Circular Segment Uniform feedrate is required along the path (X,Y,Z)1.3 (X,Y,Z)1.2 (X,Y,Z)1.1 OFF3 OFF4 (X,Y,Z)2 (X,Y,Z)2 OFF1 = min(OFF1,...,OFF4) OFF2 OFF1 (X,Y,Z)1 (X,Y,Z)1 : Tool coordinate from Original CL File : Intermediate point due to sampling for optimization Figure 5.46 : Feedrate update strategy for linear versus circular segments. 5.8. Verification for Post-Process Optimization The existing NC program for the die part described in the Chapter 3 was optimized using the proposed method. The machine tool used for producing the die was a 3-axis horizontal machining centre with maximum torque and power of 98 [Nm] and 15 [kW], respectively. Torque and power curves of the spindle shown in Figure 5.47 were obtained from the documentation of the machine. Figure 5.47 : Estimated Torque and Power Curves. 183 Chapter 5. Process Optimization Controller parameters tabulated in Table 5.5 were also required for making reasonable estimations of the cycle time and feedrate variation due to acceleration and deceleration of the table. The rest of the machine tool related data are tabulated in Table 5.6 for documentation purposes. Table 5.5 : Assumed Characteristics of the Heidenhain Controller. Table 5.6 : Other parameters of the machine tool used for verification. Torque and power curves of the machine were used as constraints for all cutters although they were most critical only for roughing operations, i.e., large depth and width of cut. Moreover, maximum chip thickness was also specified for each cutter based on Sandvik Coromant’s Metal Cutting Technical Guide and Main Catalogue. Lastly, maximum spindle speed was set for some of 184 Chapter 5. Process Optimization the operations considering dynamic characteristics, speed limitations of the spindle, and recommended surface speed. A summary of user-defined constraints is given in Table 5.7. Note that the optimized feed per tooth was capped at approximately 3.7 times the maximum chip thickness for very low immersions and depths of cut. This constant value corresponds to a cutting condition with 15% of the diameter of radial and axial immersions. Table 5.7 : Optimization Criteria for Different Operations and Cutting Tools. Opt. Maximum Chip Maximum Spindle Comments Code Thickness [mm] Speed [rpm] Opt 63 0.25 n/a For roughing with T63 Opt 63 0.056 n/a For finishing with T63 Opt 20 0.15 3000 For semi-finishing with T20 Opt 16 0.10 n/a For semi-finishing with T16 Opt 20F1 0.15 3000 For semi-finishing with T20F Opt 20F2 0.05 3000 For finishing with T20F Opt 20F3 0.10 3000 For finishing of outer walls with T20F Opt 12 0.05 5300 For rest milling with T12 Opt Small 0.025 6000 For rest milling with T8, T6, T4 Second stage optimization was run, necessary feedrate adjustments were automatically calculated and placed inside the original CL file creating the optimized CL file. The optimized CL file contained 940,000 lines as opposed to the original file’s 315,699 lines due to newly added optimized feedrate commands. Cycle times of operations were recorded during machining. There were some data logging difficulties during machining with the original CL file; therefore, time data could be recorded only for a few operations. Moreover, since spindle speed and feedrate were changed by the programmer for operations 6-9 in the original file just before machining at the site, these operations were not 185 Chapter 5. Process Optimization included in the comparison. Leaving operations 6-9 out of analysis do not affect results significantly because they had light or no cutting content. Cycle times were also predicted and results are tabulated in Table 5.8. The discrepancy between real and predicted cycle times was expected since CNC characteristics were assumed for the controller on the machine (see Table 5.5). Table 5.8 : Cycle time results of operations. Cycle Times [hr:min:second] Operation 1 Operation 2 Operation 3 Original Real Predicted 1:58:00 1:57:00 0:03:10 2:58:33 Operation 4 Operation 5 (Semi-Finishing) Operation 5 (Finishing) Operation 6 Operation 7 Operation 8 Operation 9 0:10:50 0:08:04 0:46:12 % Error -0.8% -25.0% 2:30:14 VMS Real Predicted 1:28:00 1:17:24 0:03:13 0:02:29 2:05:00 2:28:00 %Error -12.0% -22.7% 18.4% Cycle Time Difference -25.4% -21.5% -17.1% 0:05:00 1:01:00 0:04:13 0:58:26 -15.6% -4.2% -47.7% 26.5% 3:30:00 3:55:00 11.9% 37.7% 0:11:00 0:19:27 0:14:06 0:18:50 Cycle time is one parameter to evaluate the optimization, however, it is not the only one when other aspects of machining operation are considered. Indeed the goal of optimization is to maximize material removal rate and reduce cycle time, however, this must be done in a controlled way by respecting physical constraints as explained before. The cycle time reduction of Operation 1 (roughing) is reduced by 25.4% compared to the original program. More important than the cycle time was the experience and observation made during this operation. When the part was machined with the original program, there were some problems with the tool. At first, tool temperature increased to such a high level that inserts turned into red and finally chipped due to reduced chip thickness. At small feedrates, the tool was more in a grinding state than in a cutting state resulting in increased temperature and wear of the inserts (see 186 Chapter 5. Process Optimization Figure 5.48.a for chip formation). Inserts had to be replaced or rotated many times before the operation ended. In addition to tool wear and temperature, there were serious chatter problems at some sections of the program leading up to chipping of inserts. The power consumption read from operator’s panel was around 20% during this operation. Since it was not possible for the programmer to foresee how the cutting conditions would change during machining, the feedrate was conservatively selected. Optimized CL file of Operation 1; on the other hand, contained feedrates that were constantly changing along tool path. This resulted in increased chip thickness hence, a decreased temperature (see Figure 5.48.b for chip formation). The cutting operation was smoother with increased efficiency (power consumption was around 40%) Chatter problems also disappeared due to a decrease in cutting force coefficients when chip thickness became larger. Comparison of the tool T63 after both original and optimized programs are visually made in Figure 5.49. Parts of the damage to T63 (with original program) were due to sudden insert failures. (a) (b) Figure 5.48 : Chip formation with (a) original, (b) optimized NC programs. 187 OPTIMIZED ORIGINAL Chapter 5. Process Optimization Figure 5.49 : Pictures of T63 from front and bottom after execution of original and optimized CL files. An opposite case, a case with increased cycle time after optimization is also analyzed. The first half of Operation 5 was semi-finishing of the top section of the die. The tool path of this operation showed a steep gradient near the beginning and end resulting in sudden loading and unloading of the tool (see Figure 5.50 for the tool path). During optimization, not only maximum chip thick- Figure 5.50 : T20F - Part of the tool path from Operation 5, Semi-Finishing. 188 Chapter 5. Process Optimization ness was lowered, but the feedrate was also controlled during sudden geometry changes so that loading on the tool was kept uniform throughout the tool path. Inserts used during original and optimized programs were analyzed after completion of Operation 5, and pictures of them are shown in Figure 5.51 for comparison. The insert used in the uncontrolled original program chipped at the tip section whereas the insert used in the optimized program showed only regular wear marks without any apparent damage to itself. Virtual optimization is a tool to control cutting conditions to meet user-defined constraints so that part program can be executed with reduced or no problems. The user is expected to define what the most critical machining parameters to be controlled are through constraints such as maximum chip thickness, torque/power limitations, and cutting loads on the tool. For cases when an original program is too conservative, the optimization will result in a faster process with increased feedrate and spindle speed thus shortening machining time. On the other, machining operation is slowed down if an aggressive feedrate scheduling is detected. In summary, it can be concluded that cycle time will reduce or increase depending on the priorities set by the user-defined constraints. (a) (b) Figure 5.51 : Pictures of T20F inserts after the (a) original, (b) optimized machining. 189 Chapter 6 Conclusions 6.1. Conclusions A novel, physics based 3-axis virtual milling system is introduced in this thesis. The system receives the NC part program, and the geometry of blank and part from an industry standard Computer Aided Manufacturing (CAM) environment. The part-cutter engagement conditions are identified along the tool path, and used as boundary conditions to simulate the physics of milling operations. The system is capable of simulating cutting forces, torque, power, vibrations, dimensional form errors, chip loads, and bending load on the spindle bearings along the tool path through computationally efficient, generalized mathematical models of milling processes. The system allows prediction of cutting performance and possible damage to the part and the machine in a virtual environment so that remedial action can be taken before actual production takes place. In parallel, a series of physical constraints can be imposed and cutting conditions can be automatically adjusted along the tool path to achieve highest material removal rates without damaging the part and the machine. The virtual milling system is composed of three components: 1) Cutter-part engagement boundaries; 2) computationally efficient mathematical models of the mechanics and dynamics of generalized milling operations; 3) a graphical user interface that allows the user to set the physical parameters and visualize the process. The academic contributions of the thesis belong to the second component, which governs the physics of the machining 190 Chapter 6. Conclusions processes. The cutter-part engagement conditions are extracted from a Z-buffer based commercial software that provides a grip map of cutter-part intersection. The grip map is then interpolated, smoothed and projected on the cutter body to form a three dimensional envelope of the cutter-part intersection. The thesis contributions can be summarized as follows: The generalized geometric model of arbitrary cutter geometries with helical or indexable cutters is introduced. The intersection of general helical flute with a three-dimensional cutter-part engagement domain is divided into regions and closed form integrations are obtained to solve the mechanics and dynamics of milling for arbitrary cutters. Analytical and numerical models are developed to locate cutter positions when minimum and maximum values of process states are attained along the tool path and the process simulation time of a complete part is reduced significantly. This work is published in an archival journal [Merdol, S.D., Altintas, Y., 2006, Virtual Simulation and Optimization of Milling Operations Part I: Process Simulation, ASME Journal of Manufacturing Science and Engineering, accepted]. The generalized mechanics and dynamics of milling are developed, which can handle cutters having arbitrary geometries. The model has two stages; pre-process planning and post-process simulation. First, chatter stability of the machine tool - workpiece material is solved. Several mathematical approaches are introduced and their performances are compared. Analytical zero order solution is generalized to model stability of tools with complex geometries. Multi-frequency solution is analyzed in detail and a new eigenvector based chatter decision criteria is proposed. It is shown that the frequency domain stability solutions with zero and multi frequency components are computationally most feasible in identifying chatter free spindle speed, depth and width of 191 Chapter 6. Conclusions cut. The generalized mathematical model allows dividing the engagement into features and assembling their equations of motion matrixes analytically in order to solve the stability problem efficiently in the frequency domain. The chatter stability lobes are combined with torque and power limits of the machine to obtain the most efficient chatter free spindle speed, depth and width of cut combinations, which are provided to the planners ahead of NC program generation. Machining of a part in a virtual environment as the cutter moves along the tool path is simulated after the NC program is generated. The process mechanics and dynamics are decoupled in two parts: geometric orientation depending on the geometry of cutter-part engagement conditions, and the mechanics and dynamics of the process. Since the minimum and maximum values of the process states, i.e. force, torque, power, tool deflection, and chip thickness, are used as process measures, the geometric parts of the equations are reformulated as closed form analytical equations whenever possible or represented by discrete numerical models when cutters with arbitrary geometries are used. The proposed algorithms led to a rapid prediction of process states as the cutter moves along the tool path with varying engagement conditions. A comprehensive optimization formulation to maximize material removal rate is introduced. A large number of process constraints including dynamics of the process, machine tool torque and power limits, and cutting forces acting on the tool are taken into account in optimal selection of cutting parameters: spindle speed, feedrate, depth and width of cut. The material removal rates are optimized by automatically adjusting feeds and speeds along the tool path. The novel algorithms are published in archival journals [Merdol, S.D., Altintas, Y., 2006, Virtual Simulation and Optimization of Milling Operations Part II: Optimization and Feedrate Scheduling, ASME Journal of Manufacturing Sci- 192 Chapter 6. Conclusions ence and Engineering, accepted] - [Merdol, S.D., Altintas, Y., 2007, Virtual High Performance Milling, Annals of the CIRP Manufacturing Technology, Vol. 56-1,pp. 81-84]. The process is interrogated at any engagement station by first checking the stability of the process using discrete time solution of the delayed differential equation in modal coordinates. The time varying periodic directional factors of the milling are considered by checking the stability of the system in linear time-domain solution. On the other hand, the complete time history of process states such as vibrations, chip load, and dimensional surface finish are simulated numerically, allowing nonlinear cutting force coefficients, tool jumping out of cut conditions, tool run-out, serrated and variable pitch cutters. The algorithms of this work are finally integrated into a prototype industrial product: Virtual Milling System. 6.2. Future Research Directions There are remaining key challenges before achieving a fully virtual machining of parts using digital models of the machine tool and the cutting process. The present Z-buffer based cutter-part engagement identification algorithms do not have high accuracy at small grid sizes. The computational cost increases dramatically if the grid size is reduced. It is necessary to develop computationally efficient yet accurate cutter-part engagement identification algorithms. The process mechanics and dynamics models must be able to cover 5-axis milling, turning, boring, drilling, reaming and thread cutting in order to simulate the machining of a complex part on a machining centre. The stability of turning, drilling, boring, reaming and low speed machining with sufficiently practical accuracy remain unsolved. The contact between the flank face of the tool and wavy finish 193 Chapter 6. Conclusions surface add additional stiffness and damping but their physics are yet to be modeled with an acceptable accuracy. CNC dynamics, i.e. feed and spindle servo systems, change the effective chip loads that directly affects all process states along the tool path. The dynamics of the CNC system must be integrated to the process physics in order to improve the accuracy of the virtual simulation of the machining process and machine tool motion. 194 Bibliography [1] Taylor, F. W., 1907, On the art of cutting metals, Transactions of the ASME, Vol. 28. [2] Airey, J., and Oxford, C. J., 1921, On the art of milling, Transactions of the ASME, Vol. 43, pp. 549. [3] Sawin, N. N., 1926, Theory of milling cutters, Mechanical Engineering, Vol. 48, pp. 12031209. [4] Salomon, C., 1926, Die frasarheit, Werkstattstechnik, Vol. 20, pp. 469-474. [5] Parsons, F., 1923, Power required for cutting metal, Transactions of the ASME, Vol. 49, pp. 193-227. [6] Boston, O.W., and Kraus, C.E., 1932, Elements of milling, Transactions of the ASME, Vol. 54, pp. 71-104. [7] Boston, O.W., and Kraus, C.E., 1932, Elements of milling – Part II, Transactions of the ASME, Vol. 56, pp. 355-371. [8] Martelotti, M.E., 1941, An analysis of the milling process, Transactions of the ASME, pp. 677-700. [9] Fussell, B.K., Jerard, R.B., Hemmett, J.G., 2001, Robust feedrate selection for 3-Axis NC machining using discrete models, Transactions of the ASME, Vol. 123, pp. 214-224. [10] Erdim, H., Lazoglu, I., Ozturk, B., 2006, Feedrate scheduling strategies for free-form surfaces, International Journal of Machine Tools and Manufacture, Vol. 46, pp. 747-757. [11] Available from http://www.dptechnology.com. [12] Fu, H.J., Devor, R.E., Kapoor, S.G., 1984, A mechanistic model for the prediction of the force system in face milling operation, ASME Journal of Engineering for Industry, Vol.106-1, pp. 81-88. [13] Sutherland, J.W., Devor, R.E., 1986, An improved method for cutting force and surface error prediction in flexible end milling system, ASME Journal of Engineering for Industry, Vol. 108, pp. 269-279. [14] Bayoumi, A.E., Yucesan, G., Kendall, L.A., 1994, An analytical mechanistic cutting force model for milling operations: a theory and methodology, Transactions of ASME, Vol. 116, pp. 324-330. 195 Bibliography [15] Spiewak, S.A., 1994, Analytical modeling of cutting point trajectories in milling, ASME Journal of Engineering for Industry, Vol. 116-4, pp. 440-448. [16] Lazoglu, I., Liang, S.Y., 1997, Analytical modeling of force system in ball end milling, Journal of Machining Science and Technology, Vol. 1-2, pp. 219-234. [17] Altintas, Y., Lee, P., 1998, Mechanics and dynamics of ball end milling, ASME Journal of Manufacturing Science and Engineering, Vol. 120, pp. 684-692. [18] Yucesan, G., Altintas Y., Prediction of ball end milling forces, ASME Journal of Engineering for Industry, Vol. 1-1, pp. 95-103. [19] Ramaraj, T.C., Eleftheriou, E., 1994, Analysis of the mechanics of machining with tapered end milling cutters, Transactions of ASME, Vol. 1216, pp. 398-404. [20] Engin, S., Altintas, Y., 2001, Mechanics and dynamics of general milling cutters. Part I: helical end mills, International Journal of Machine Tools and Manufacture, Vol. 41, pp. 2195-2212. [21] Salomon. C., 1926, Die Frasarheit, Werkststtstechnik, Vol. 20, pp. 469-474. [22] Sabberwal, A.J.P., 1961, Chip section and cutting force during the milling operation, Annals of the CIRP, Vol. 10, pp. 197-203. [23] Koenigsberger, F., Sabberwal, A.J.P, 1961, An investigation into the cutting force pulsations during milling operations, International Journal of Machine Tool Design and Research, Vol. 1, pp. 15-33. [24] Tlusty, J., McNeil P., 1970, Dynamics of Cutting Forces in End Milling, Annals of the CIRP, Vol.24, pp.21-25. [25] Kline, W.A., DeVor, R.E., Zdeblick, W.J., 1980, A mechanistic model for the force system in end milling with application to machining airframe structures, Manufacturing Engineering Transactions, pp. 297. [26] Altintas, Y., and Spence, A., 1991, End Milling Force Algorithms for CAD Systems, Manufacturing Technology CIRP Annals, Vol. 40, pp. 31-34. [27] Armarego, E.J.A, Epp, C.J., 1970, An Investigation of Zero Helix Peripheral Up-Milling, Intenational Journal of Machine Tool Design and Research, Vol. 10, pp. 273-291. [28] Armarego, E.J.A., Whitfield R.C., 1985, Computer Based Modeling of Popular Machining Operations for Force and Power Predictions, Annals of the CIRP, Vol. 34, pp. 65-69. [29] Brown, R.H., Armarego, E.J.A., 1964, Oblique Machining with a Single Cutting Edge, Intenational Journal of Machine Tool Design and Research, Vol. 4, pp. 9-25. 196 Bibliography [30] Budak, E., Altintas, Y., Armarego, E.J.A., 1996, Prediction of Milling Force Coefficients from Orthogonal Cutting Data, Transactions of ASME, Vol. 118, pp. 216-224. [31] Altintas, Y., Spence, A. D., 1994, A Solid modeller based milling process simulation and planning system, ASME Journal of Engineering for Industry, Vol. 116, pp. 61-69. [32] Smith, S., Tlusty, J., 1991, An overview of modeling and simulation of the milling process, ASME Journal of Engineering for Industry, Vol. 113-2, pp.169-175. [33] Armarego, E.J.A., Deshpande, N.P., 1991, Computerized end-milling force predictions with cutting models allowing eccentricity and cutter deflections, Annals of the CIRP, Vol. 40-1, pp. 25-29. [34] Tlusty, J., Ismail, F., 1981, Basic nonlinearity in machining chatter, Annals of the CIRP, Vol. 30, pp. 21–25. [35] Sutherland, J. W., and DeVor, R. E., 1988, A dynamic method for the cutting force system in the end milling process, Sensors and Controls for Manufacturing, ASME, New York, PED-33, pp. 269–279. [36] Montgomery, D., Altintas, Y., 1991, Mechanism of cutting force and surface generation in dynamic milling, ASME Journal of Engineering for Industry, Vol. 113, pp. 160-168. [37] Altintas, Y., Lee, P., 1998, Mechanics and dynamics of ball end milling, ASME Journal of Manufacturing Science and Engineering, Vol. 120, pp. 684-692. [38] Smith, S., and Tlusty, J., 1993, Efficient simulation programs for chatter in milling, Annals of the CIRP, Vol. 42, pp. 463–466. [39] Bryne, G., Dornfeld, D., and Denkena, B., 2003, Advancing Cutting Technology, Annals of the CIRP, Vol. 52-2. [40] Fussell, B.K., Ersoy, C., and Jerard, R.B., 1992, Computer Generated CNC Machining Feedrates, ASME Japan/USA Symposium on Flexible Automation, Vol. 1, pp. 377-384. [41] Hemmett, J.G., Fussell, B.K., and Jerard, R.B., 2000, A Robust and Efficient Approach to Feedrate Selection for 3-Axis Milling, Proceedings of IMECE Symposium on Dynamics and Control of Material Removal Processes, DSC-Vol. 2, pp. 729-736. [42] Jerard, R.B., Drysdale, R.L., Hauck, K., Schaudt, B., and Magewick, J., 1989, Methods for Detecting Errors in Sculptured Surface Machining, IEEE Computer Graphics and Applications, Vol. 9-1, pp. 26-39. [43] El-Mounayri, H., Spence, A.D. and Elbestawi, M.A., 1998, Milling Process Simulation - A Generic Solid Modeller Based Paradigm, Journal of Manufacturing Science and Engineering, Vol. 120, No. 2, p. 213-221. 197 Bibliography [44] Altintas, Y., and Spence, A. D., 1994, A Solid modeller based milling process simulation and planning system, ASME Journal of Engineering for Industry, Vol. 116, pp. 61-69. [45] Spatial Technology Inc., Boulder, CO, USA, Available from http://www.spatial.com/products/acis.html. [46] Saturley, P.V. and Spence, A.D., 2000, Integration of Milling Process Simulation with OnLine Monitoring and Control, International Journal of Advanced Manufacturing Technologies, Vol.16, pp. 92-99. [47] Yip-Hoi, D., and Huang, X., 2006, Cutter/Workpiece Engagement Feature Extraction from Solid Models for End Milling, Journal of Manufacturing Science and Engineering, Vol. 128, pp. 249-260. [48] Weinert, K., Du, S., Damm, P., and Stautner, M., 2004, Swept Volume Generation for the Simulation of Machining, International Journal of Machine Tools and Manufacture, Vol. 44, pp. 617-628. [49] Gradisek J., Kalveram, M., Weinert, K., 2004, Mechanistic identification of specific force coefficients for a general end mill, International Journal of Machine Tools and Manufacture, Vol. 44, pp. 401-414. [50] Abrari, F., Elbestawi, M. A., 1996, Closed form formulation of cutting forces for ball and flat end mills, International Journal of Machine Tools and Manufacture, Volume 37, No. 1, pp.17-27. [51] CGTech, Irvine, CA, USA, Available from http://www.cgtech.com/usa/optipath%c2%a9/. [52] Kline, W.A., DeVor, R.E., Lindberg, J.R., 1982, Prediction of cutting forces in end milling with application to cornering cuts, International Journal of Machine Tool Design Research, Vol. 22. [53] Wang, W.P., 1988, Solid modeling for optimizing metal removal three-dimensional NC end milling, ASME Journal of Engineering for Industry, Vol. 113-2, pp. 169-175. [54] Ip, R.W.L., Lau, H.C.W., Chan, F.T.S., 2003, An economical sculptured surface machining approach using fuzzy models and ball-nosed cutters, Journal of Materials Processing Technology, Vol. 138, pp .579–585. [55] Available from http://www.mastercam.com. [56] Yazar, Z., Koch, K.F., Merrick, T., Altan, T., 1994, Feed rate optimization based on cutting force calculations in 3-axis milling of dies and molds with sculptured surfaces, International Journal of Machine Tool and Manufacture, Vol. 34, pp. 365–377. 198 Bibliography [57] Lim, M., Hsiang, M.C., 1997, Integrated planning for precision machining of complex surfaces, Part 1: cutting-path and feedrate optimization, International Journal of Machine Tools and Manufacture, Vol. 37, pp. 61–75. [58] Bailey, T., Elbestawi, M.A., El-Wardany, T.I., Fitzpatrick, P., 2002, Generic simulation approach for multi-axis machining, part 2: model calibration and feed rate scheduling, ASME Journal of Manufacturing Science and Engineering, Vol. 124, pp. 634–642. [59] Erdim, H., Lazoglu, I., Guzel, B., 2006, Feedrate scheduling strategies for free-form surfaces, International Journal of Machine Tools and Manufacture, Vol. 46, n 7-8, pp. 747757. [60] Fussell, B.K., Jerard, R.B., Hemmett, J.G., 2001, Robust feedrate selection for 3-axis NC machining using discrete models, ASME Journal of Manufacturing Science and Engineering, Vol. 123, pp. 214-224. [61] Budak, E., Tekeli, A., 2005, Maximizing chatter free material removal rate in milling through optimal selection of axial and radial depth of cut pair, Annals of the CIRP, Vol. 54/ 1, pp. 353-356. [62] Tlusty, J., 1965, A Method of Analysis of Machine Tool Stability, Proceeding MTDR, pp.514. [63] Tobias, S.A., 1965, Machine Tool Vibration, Blackie and Sons Ltd. [64] Altintas, Y., 2000, Manufacturing Automation, Cambridge University Press. [65] Merritt, H.E., 1965, Theory of Self-Excited Machine Tool Chatter, ASME Journal of Engineering for Industry, Vol. 87, pp.447-454. [66] Minis I., Yanushevsky T., 1993, A new theoretical approach for the prediction of machine tool chatter in milling, ASME Journal of Engineering for Industry, Vol. 115, pp. 1-8. [67] Altintas Y., Budak E., 1995, Analytical prediction of stability lobes in milling'', Annals of CIRP, Vol. 44/1, pp. 357-362. [68] Insperger T., Stepan G., Bayly P.V., Mann B.P., 2003, Multiple chatter frequencies in milling processes, Journal of Sound and Vibration, Vol 262, pp. 333-345. [69] Davies, M.A., Pratt J.R., Dutterer B., Burns T.J., 2002, Stability prediction for radial immerision milling, ASME Journal of Manufacturing Science and Engineering, Vol. 124, pp. 217-225. [70] Bayly P.V., Halley J. E., Mann B.P., Davies M.A., 2003, Stability of interrupted cutting by temporal finite element analysis, ASME Journal of Manufacturing Science and Engineering, Vol. 125, pp. 220-225. 199 Bibliography [71] Bayly, P.V., Schmitz, T.L., Gabor, S., Mann, B.P., Peters D.A., Insperger T., 2002, Effects of radial immersion and cutting direction on chatter instability in end-milling, Proceedings of the ASME Manufacturing Engineering Division, Vol. 13, pp. 351-363. [72] Stepan, G., Insperger, T., 2002, Semi-discretization method for delayed systems, International Journal for Numerical Methods in Engineering, Vol. 55-5, pp. 503-518. [73] Insperger, T., Stepan, G., 2004, Updated semi-discretization method for periodic delay-differential equations with discrete delay, International Journal for Numerical Methods in Engineering, Vol. 61-1, pp. 117-141. [74] Olgac, N., Sipahi, R., 2003, The direct method for stability analysis of time delayed LTI systems, Proceedings of the American Control Conference, Vol. 1, pp 869-874. [75] Merdol, S.D., Altintas, Y., 2004, Multi frequency solution of chatter stability for low immersion milling, ASME Journal of Manufacturing Science and Engineering, Vol. 126, pp. 459-466. [76] Chan, S.H., Yang, X.W., Lian, H.D., 2001, Comparison of several eigenvalue re-analysis methods for modified structures, Structural and Multidisciplinary Optimization, Vol. 21-1, pp. 84. [77] Courant, R., Hilbert, D., 1945, Methods of mathematical physics, Interscience Publishers Inc., Vol. 1., New York. [78] Nayfeh, AH., 1973, Perturbation methods, Wiley, New York. [79] Chen J.C., Wada, B.K., 1977, Matrix perturbation for structural dynamics, AIAA Journal, Vol. 15, pp.1095–1100. [80] Murthy D.V., Haftka R.T., 1988, Approximations to eigenvalues of modified general matrices, Computers and Structures, Vol. 29, pp. 903-917. [81] Chen S-H., Song D-T., Ma A-J., 1994, Eigensolution re-analysis of modified structures using perturbations and Rayleigh quotients, Communications in Numerical Methods in Engineering, vol. 10, pp. 111–119. [82] Bickford WB., 1987, An improved computational technique for perturbations of the generalized symmetric linear algebraic eigenvalue problem, International Journal for Numerical Methods in Engineering, Vol. 24, pp. 529 –541. [83] Lu, Z.H., Shao, C., Feng, Z.D., 1993, High accuracy step-by-step perturbation method of the generalized eigenvalue problem and its application to the parametric history analysis of structural vibration characteristics, Journal of Sound and Vibration, Vol. 164-3, pp. 459469. 200 Bibliography [84] Liu, X.L., Oliveira, C.S., 2003, Iterative modal perturbation and re-analysis of eigenvalue problem, Communications in Numerical Methods in Engineering, Vol. 19-4, pp.263-274. [85] Lou, M., Chen, G, 2003, Modal perturbation method and its applications in structural systems, Journal of Engineering Mechanics, Vol. 129-8, pp. 935-943. [86] Armarego, E.J.A., 1993, Material Removal Process: An Intermediate Course, The University of Melbourne, pp. 48-49. [87] Merchant, M.E., 1994, Basic Mechanics of the Metal Cutting Process, Transactions of ASME Journal of Applied Mechanics, pp. 168-175. [88] Stabler, G.V., 1964, The Chip Flow Law and Its Consequences, Advances in Machine Tool Design and Research, pp. 243-251. [89] Takata, S., 1989, A Cutting Simulation System for Machinability Evaluation Using a Workpiece Model, Annals of the CIRP, 38-1, pp. 417-420. [90] Hemmett, J.G., Fussell, B.K., Jerard, R.B., 2000, A Robust and Efficient Approach to Feedrate Selection for 3-Axis Milling, Proc. IMECE Symp. on Dynamics and Control of Material Removal Processes, DSC-Vol. 2, pp. 729-736. [91] Jerard, R.B., Drysdale, R.L., Hauck, K., Schaudt, B., Magewick, J., 1989, Methods for Detecting Errors in Sculptured Surface Machining, IEEE Computer Graphics and Applications, Vol. 9-1, pp. 26-39. [92] W. H. Press et al., 1988, Numerical Recipes in C: The Art of Scientific Computing. Cambridge University Press. [93] Campomanes, M., 1998, Dynamics of Milling Flexible Structures, The University of British Columbia - Master of Applied Science Thesis. [94] Kops, L., Vo, D.T., 1990, Determination of equivalent diameter of an end mill based on its conpliance, Annals of the CIRP, Vol. 39/1, p. 93. [95] Budak, E., Altintas, Y., 1994, Peripheral Milling Conditions for Improved Dimensional Accuracy, International Journal of Machine Tool Design and Research, pp.907-918. [96] Shamoto, E., Altintas, Y., 1999, Prediction of shear angle in oblique cutting with maximum shear stress and minimum energy principles, ASME Journal of Manufacturing Science and Engineering, Vol. 121-3, pp. 399-407. [97] Coleman, T., Branch, M.A., and Grace, A., Optimization toolbox user’s guide, Matlab The Mathworks Inc., version 2. [98] Coddington, E.A., Levinson, N., 1955, Theory of ordinary differential equations, Mc-Graw Hill, N.Y. 201 Bibliography [99] Yakubovitch, V.A., 1975, Strazhinskii, Linear differential equations with periodic coefficients, John Wiley and Sons, N.Y. [100] Rozenvasser, E.N., 1972, Computation and transformation of transfer functions of linear periodic systems, Automation and Remote Control, Vol. 33-2.1, pp. 220-227. [101] Analytical prediction of three dimensional chatter stability in milling, 2001, JSME International Journal Series C: Mechanical Systems, Machine Elements and Manufacturing, Vol. 44-3-September, pp. 717-723. 202
Featured Collection
UBC Theses and Dissertations
Virtual three-axis milling process simulation and optimization Merdol, Doruk Sūkrū 2008
Item Metadata
Download
- Media
- 24-ubc_2008_spring_merdol_doruk.pdf [ 4.35MB ]
- Metadata
- JSON: 24-1.0066338.json
- JSON-LD: 24-1.0066338-ld.json
- RDF/XML (Pretty): 24-1.0066338-rdf.xml
- RDF/JSON: 24-1.0066338-rdf.json
- Turtle: 24-1.0066338-turtle.txt
- N-Triples: 24-1.0066338-rdf-ntriples.txt
- Original Record: 24-1.0066338-source.json
- Full Text
- 24-1.0066338-fulltext.txt
- Citation
- 24-1.0066338.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url: http://iiif.library.ubc.ca/presentation/dsp.24.1-0066338/manifest