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 model- ing of the machining processes is a fundamental requirement in identifying optimal cutting condi- tions 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 opti- mize 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 com- plexity. The dynamics of generalized milling operations are modeled and the stability of the pro- cess 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 optimiza- tion 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. Opti- mized 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 ............................................................................38iv 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 ......................................................................................144v 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 ...............................................................................................................................195vi vii 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 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 (2- fluted) 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.........................................59viii 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..................124ix 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) .153x 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 .....................................................................183xi 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 ...............189xii Nomenclature a axial depth of cut limiting axial depth of cut for chatter stability A closed form coefficients 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 feed vector in the direction of resultant tool motion commanded feedrate end feedrate start feedrate 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 alim A t( ) f fc fe fs fuxiii uncut dynamic chip thickness maximum number of harmonics i, j, k unit vectors in the principal directions x, y and z helix angle I area moment of inertia I identity matrix j flute number tangential, radial, and axial milling-cutting force coefficients tangential, radial, and axial milling-edge force coefficients m modal mass n spindle speed unit vector perpendicular to the cutter body N number of teeth on the cutter p eigenvector P mode shape matrix 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 lower and upper axial integration limits of process outputs friction angle normal rake angle vibration vector hd hr i0 Ktc Krc Kac, , Kte Kre Kae, , nj R zlim 0, zlim 1,, βa αn Δxiv tool deflection reflected on the surface phase shift between current and previous vibration marks damping ratio chip flow angle axial immersion angle eigenvalue imaginary part of the eigenvalue real part of the eigenvalue tooth period, time delay shear stress oriented transfer function matrix rotation angle of the cutter shear angle pitch angle of the cutting tool radial start and exit angles in milling phase shift of the complex eigenvalue lag angle chatter frequency natural frequency tooth passing frequency transpose operator ACIS Andy Charles Ian’s System (solid modeling kernel) δ ε ζ η κ λ λI λR τ τs Φ φ φc φp φst φex, ϕ ψ ωc ωn ωT . Txv 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-Orderxvi xvii 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 sup- port 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 diffi- culties 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. Chapter 1 Introduction Milling operations are most widely used to more accurately produce net shapes in the manufactur- ing 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 indus- try. 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. Produc- tivity 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. Neverthe- less, 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 condi- tions are modeled. First, the chatter vibration free spindle speeds, depth and width of cut are pre- dicted 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 mate- rial removal rate while respecting the physical limits of the machine and cutting tool. The devel- oped 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. Computa- tionally efficient, closed-form formulations are derived for cutters with analytically defined cut- ting 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 avail- able 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 dis- cussed. 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 pro- gramming. 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 Reviewtool 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 fol- lowed 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: , (2.1) where h is the instantaneous uncut chip thickness, c is the feed per tooth and is the tooth immer- sion angle. h φ( ) c φsin⋅= φ Y X X φ st= 0 Feed φ φ ex= π Feed Y (a) (b) Figure 2.1 : Methods of milling operations: (a) up milling; (b) down milling.6 Chapter 2. Literature ReviewSinusoidal 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 cal- culate the uncut chip thickness. Lazoglu et al. [10]; on the other hand, separated conventional hor- izontal 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. A great selection of cutting tools are used in the industry to machine different parts. In the aero- space industry, cylindrical end mills are heavily involved in roughing and semi-finishing of alu- minum 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 devel- oped 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]. Figure 2.2 : 3-axis milling examples (source: Esprit [11]).7 Chapter 2. Literature ReviewIn 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 vibra- tions. 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 Koenigs- berger [22][23] used similar exponential specific cutting pressures to model milling forces analyt- ically. 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 pres- sure (K): , (2.2) where is the specific cutting pressure, C and m are constants depending on the work- piece material and the milling cutter, is the chip area, a is the depth of cut and h is the instan- taneous 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 FCutting Force K a h⋅ ⋅= K C hm⋅= a h⋅8 Chapter 2. Literature Reviewmaterial on the flank face of the cutting edge. In this thesis, the linear edge force model by Armar- ego 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 thick- ness as: , (2.3) where 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. FCutting Force Ke a⋅ Kc a h⋅ ⋅+= Ke, and Kc9 Chapter 2. Literature ReviewIn 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 out- puts. (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 Reviewemphasize that the effect of previous time steps are taken into account. Tlusty and Ismail [34] pre- sented 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 ele- ments 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 influ- ence 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 com- plete 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 Review2.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. Altin- tas and Spence [26] used constructive solid geometry to identify cutter-part intersection, and re- Y X Exact Chip Thickness Workpiece End Mill Previous Surface Combined Dynamics of Cutter, Workpiece, and Machine Tool Current Surface (being generated) Figure 2.3 : Exact kinematics of a milling tool for calculation of the chip thickness.12 Chapter 2. Literature Reviewformulated 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. Cut- ting forces were then calculated using the removed material volume and other geometric informa- tion. 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 repre- sentation (B-Rep) modeller [46]. Yip-Hoi et al. [47] presented feature based cutter-part intersec- tion model for cylindrical end mills using ACIS solid modeler. Although more accurate cutter- part 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 cut- ter-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 project- ing the chip area onto a set of reference coordinate planes; however, this method cannot be gener- alized 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 geo- metric 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 ReviewIn 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 aggre- gation 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 Reviewproposed 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 analy- sis 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 compre- hensive machining simulation and optimization scheme that considers not only the geometric fac- tors 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 sys- tematic 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 Reviewvarious rough milling strategies to reduce machining time and cost. Lim et al. [57] proposed cut- ting 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 vol- ume removal. Hence, there is still a demand for a generalized optimization strategy that combines 16 Chapter 2. Literature Reviewdifferent 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. 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 condi- tions, 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 Strategy analysis simulation results improved path planning Spindle dynamics (FRF) Feed drive FRF and control Cutter-material cutting constants Kinematics of the machine Volumetric error model of the machine MACHINING DATABASE NC Tool Path Cutter Geometry CAD MODEL Optimized speed, feed, depth, width, error compensation FINAL PROCESS PLAN Peak force, torque, power, tracking error, modal frequencies MONITORING AND CONTROL DATA CL File PATH PLANNER M A C H I N E T O O L Cutter-part intersection geometry calculations Machining process simulation CAD BASED VIRTUAL MACHINING PROCESS SIMULATION SYSTEM Figure 2.4 : A comprehensive process simulation and optimization flow chart by Altintas [26].17 Chapter 2. Literature Reviewleave 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 dif- ference 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 work- piece, 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 opera- tions. 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: , (2.4) where 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 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 coeffi- alim 1– 2Kt b d -- N 2 --- Gmin α⋅ ⋅ ⋅ ⋅ -------------------------------------------------= Kt Gmin α⋅18 Chapter 2. Literature Reviewcient 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). Minis and Yanushevsky [66] applied the theory of periodic differential equations and obtained the resulting characteristic equation of infinite order: , (2.5) where is the Kronecker delta, I is a (2m+1)-by-(2m+1) identity matrix, a is the depth of cut, is the cutting force coefficient in tangential direction, (r,k) are harmonic numbers, i.e , is the tooth passing frequency and is the oriented transfer function, which consists of directional milling coefficients and two dimensional transfer func- h(t) y(t) y(t-T) h0 εh0 Tool Workpiece Tool Feed Direction a Chuck n Workpiece Figure 2.5 : Wave generation in turning. det δrk I⋅ 12-- a Kt 1 e λT––( ) Wr k– λ ikωT+( )⋅ ⋅ ⋅ ⋅– 0= δrk Kt r k, 0 1 2 … m±,±,±,= ωT Wr k–19 Chapter 2. Literature Reviewtions. 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 at the limit of stability (marginal stability) and obtained chatter free- axial depth of cut as: , (2.6) where N is the number of teeth. and are the real and imaginary values of , which is obtained as an eigenvalue of the following equation: , (2.7) where 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 indi- vidually. 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: , (2.8) λ i± ωc= alim 2– πΛR NKt ----------------- 1 ΛI ΛR ------⎝ ⎠⎛ ⎞ 2 += ΛR ΛI Λ det I Λ A G⋅ ⋅+( ) 0= A G⋅ ωH ωc k ωT⋅±=( ) k, 0 1 2…, ,={ }20 Chapter 2. Literature Reviewwhere is the dominant chatter frequency and is the tooth passing frequency. The second set referred to the Period Doubling Bifurcation leading to . (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 condi- tions - 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 fre- quency 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-free- dom. 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 ωc ωT 2π T⁄= ωPD ωT 2⁄( ) k ωT⋅±=( ) k, 0 1 2…, ,={ }21 Chapter 2. Literature Reviewdetermine stability. Bayly et al. [71], later, extended their stability analysis to cover two-degree- of-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 Reviewtaken into account in characterizing the directional coefficients; therefore, the standard calcula- tion of eigenvalues between iteration steps becomes computationally expensive. Since the system defined by the characteristic equation varies slightly at each iteration step, eigenvalues and eigen- vectors 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; effi- ciency 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. Clas- sical 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-by- step as suggested by Lu et al. [83]. The step-by-step perturbation method starts to lose its accu- racy 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 Reviewapplied 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 sim- plified 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 geome- try covering all possible helical end mills to model mechanics of 3-axis milling operations. Gener- alizations 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 indepen- dent 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 pro- file 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 SimulationIn 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 vari- able. In a general milling operation, as illustrated in Figure 3.2, there are two different coordinate sys- tems 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 cut- ting 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, (rep- resented by ) 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 rota- tion matrix between two coordinate frames can be defined as: . (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 vari- ables of an end mill have to be analyzed in detail before proceeding with the modeling of cutting fxy X Y Z⎩ ⎭⎪ ⎪ ⎨ ⎬⎪ ⎪ ⎧ ⎫ αcos αsin– 0 αsin αcos 0 0 0 1 x y z⎩ ⎭⎪ ⎪ ⎨ ⎬⎪ ⎪ ⎧ ⎫ =27 Chapter 3. Process Simulationforces. 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 cut- ter 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. When a cutter is rotated at a certain spindle speed, n [rpm], the angular position of the cutter var- ies over time as: . (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: , , (3.3) x y α X Y MTCS CCS fxy Workpiece z Z XY Figure 3.2 : Coordinate systems in a milling operation. φ t( ) 2πn60--------- t⋅= φj t z,( ) φj φ t( ) φpj ψ z( )–+= = j 1 … N, ,=28 Chapter 3. Process Simulationwhere N is the number of teeth, , 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, is the pitch angle of flute j and is expressed as: , (3.4) and , (3.5) 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 geome- try is defined for each portion. Once all the geometric relations are known, the unit vector perpen- dicular to cutter body at axial elevation z and angular position can be defined as: , (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 ( ) 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 ( ): ψ z( ) φpj φpj j 1–( )φp= φp 2π N⁄= φj nj z( ) κ z( ) φj z( ) i⋅ κ z( ) φj z( ) j⋅ κ z( ) k⋅cos–cos⋅sin+sin⋅sin= ψ z( ) 0 z Mz≤<29 Chapter 3. Process Simulation. (3.7) To prevent singularity at the tool tip, z value must start from a finite but very small number, for example from z = as stated in reference [20]. • Arc Section ( ): . (3.8) • Tapered Section ( ), (3.9) where is the helix angle, and , the nominal helix for a constant lead cutter, is defined as: , (3.10) where is the lead of the helical flute. Boundaries of each zone, TIP, ARC, and TAP are given as follows: , (3.11) , (3.12) , (3.13) , (3.14) ψTIP z( ) z αcot⋅( )ln i0tan⋅ αcos---------------------------------------------= Mr 20⁄ Mz z Nz≤ ≤ ψARC z( ) z Mz–( ) R Rz–+( ) i0tan⋅ R --------------------------------------------------------------- ψTIP Mz( )+= Nz z h≤ ≤ ψTAP z( ) 1 z N– z( ) Nr⁄( ) βtan⋅+( )ln i0tan⋅ βsin------------------------------------------------------------------------------------- ψARC Nz( )+ if β 0≠( ) & Constant Helix z Nz–( ) Nr⁄( ) i0tan⋅ ψARC Nz( )+ if β 0=( ) & Constant Helix z Nz–( ) Nr⁄( ) istan⋅ ψARC Nz( )+ Constant Lead⎩⎪ ⎪⎨ ⎪⎪ ⎧ = i0 is is 2πNr( ) Llead. βcos( )⁄[ ]atan= Llead Mr Rz αtan Rr R2 Rr2–( ) αtan2 2RzRr αtan R2 Rz2–( )+ ++ + αtan2 1+---------------------------------------------------------------------------------------------------------------------------------------------= Mz Mr αtan= Nr Rr u–( ) βtan Rz R2 Rz2–( ) βtan2 2Rz Rr u–( ) βtan R2 Rr u–( )2–+ +–+ βtan2 1+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------= Nz u Nr βtan+=30 Chapter 3. Process Simulationwhere , , and . 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= = = =0, =D/2) with helix angle and diameter D, , (3.15) where . • Tapered Helical Ball End Mill ( = =0, =R= ) with taper angle , Ball End Section with ball radius : for , (3.16) Tapered Section with constant helix ( ): for , (3.17) Tapered Section with constant lead ( ): for , (3.18) where , , , , . 3.3.2. Radius ( ) and Axial Immersion Angle ( ) 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, , 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: u 1 α βtantan–( )D 2⁄= 0 α π 2⁄<≤ β π 2⁄< Rz α β Rr i0 ψ z( ) ki0z= ki0 2 i0tan( ) D⁄= Rr α Rz Rb0 β Rb0 ψb z( ) i0tan Rb0⁄( ) z⋅= z z1≤ i0 ψ z( ) c1 1 c2 z z1–( )⋅+[ ]ln⋅ ψb z1( )+= z z1≥ Llead ψ z( ) c3 z z1–( )⋅ ψb z1( )+= z z1≥ c1 i0tan βsin⁄= c2 βtan R0⁄= c3 istan R0⁄= is 2πR0 Llead. βcos( )⁄( )atan= R0 Nr= z1 Mz= R z( ) κ z( ) κ z( )31 Chapter 3. Process Simulation, (3.19) . (3.20) • Arc Section: , (3.21) . (3.22) • Tapered Section , (3.23) . (3.24) R z( ) z αtan⁄= κ z( ) α= R z( ) R2 Rz z–( )2– Rr+= κ z( ) R z( ) Rr–R---------------------⎝ ⎠ ⎛ ⎞asin= R z( ) u z βtan+= κ z( ) π 2⁄ β–= M N M O P C R L M N h ST z r r z NRz Rr D α TIP ARC TAP ψ(z) φj(z) φpj φ+ φ(t) z xy x y κ(z) ! @ # !: Reference Cutter Position @: Tip Position for Flute j #: Position of Element Located at Axial Height z on Flute j β Figure 3.3 : (a) Generalized tool geometry, (b) angle convention. (a) (b)32 Chapter 3. Process Simulation3.4. Mechanics of Milling - Generalized Kinematics Process simulation of a milling operation requires extensive and accurate modeling of the mate- rial removed by the cutter at each time step. The most fundamental parameter in cutting mechan- ics 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: , (3.25) where c [mm/rev.tooth] is the feed per tooth and is the unit feed vector in the direction of resultant tool motion expressed as: , . (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 . 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 onto the geometry vector via a dot product of two vectors: . (3.27) As material is removed from the workpiece, forces are exerted on the tool. Three orthogonal dif- ferential 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: f c fu⋅= fu fu fxy i⋅ fz k⋅+= fu 1= fxy fx 2 fy 2+= f nj z( ) ( )• hj φ z,( ) nj z( ) f• c fxy κ z( ) φj z( )sin⋅sin⋅ fz κ z( )cos⋅–( )⋅= =33 Chapter 3. Process Simulation, (3.28) where dS(z) is the differential contact length expressed as: , (3.29) and dz is the axial height of the differential element. are cutting force coefficients due to the shear and are edge cutting force coefficients due to the rubbing of the tool flank with the workpiece, in radial, tangential and axial directions, respectively. The oblique cutting coefficients, [N/ ], are either obtained from mechanistic calibration tests or transformed from an orthogonal cutting database that contains shear stress ( ), shear angle ( ), friction angle ( ) and the edge cutting force coefficients [86]. Although, in dFrj φ z,( ) dFtj φ z,( ) dFaj φ z,( )⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ Kre Kte Kae⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ dS z( )⋅ Krc Ktc Kac⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ hj φ z,( ) dS z( )⋅ ⋅+= dS z( ) dz κ z( )sin⁄= Krc Ktc Kac, , Kre Kte Kae, , z x y Cutting Edge( j) n zj,0 zj,1 . . Cutting Zone y x φj κ(z) z x dFrj dFtj Chip(ψ,φ,κ) dz dS Differential Chip dFaj hj dz nj(z) f Figure 3.4 : Mechanics and kinematics of 3-axis milling. Krc Ktc Kac, , mm2 τ φc βa34 Chapter 3. Process Simulationgeneral, 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: . (3.30) By substituting Eq (3.28) into Eq (3.30), the differential cutting forces can be represented as , (3.31) where , , (3.32) . (3.33) dFxj φ z,( ) dFyj φ z,( ) dFzj φ z,( )⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ κ z( ) φj z( )sinsin– φj z( )cos– κ z( ) φj z( )sincos– κ z( ) φj z( )cossin– φj z( )sin κ z( ) φj z( )coscos– κ z( )cos 0 κ z( )sin– dFrj φ z,( ) dFtj φ z,( ) dFaj φ z,( )⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ = dFj xyz φ z,( ) Aexyz j, dz⋅ Acxyz j, c⋅+ dz⋅= dFj xyz φ z,( ) dFxj φ z,( ) dFyj φ z,( ) dFzj φ z,( )⎩ ⎭⎨ ⎬ ⎧ ⎫T = Ae xyz j, φjKresin– φjcos κsin-------------Kte φjsin κtan------------Kae–– φjcos Kre– φjsin κsin------------Kte φjcos κtan-------------Kae–+ 1 κtan---------- Kre Kae–⎩ ⎭⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎧ ⎫ = Ac xyz j, φjsin κKrcsin– φjcos Ktc φjsin– κKaccos– φjcos κsin Krc– φjsin Ktc φjcos κcos Kac–+ κcos Krc κsin Kac–⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ fxy φjsin fz κtan----------–⎝ ⎠⎛ ⎞⋅=35 Chapter 3. Process SimulationDifferential 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 and 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: , (3.34) where . The integrations are carried out to reduce cut- ting force experienced by each flute into the following form: , (3.35) where . (3.36) And finally, total cutting forces at time t ( ) are calculated by simply adding contribution of each flute: . (3.37) 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 coefficients , and need to be calculated once for a constant cutter-workpiece zj 0, zj 1, Fj xyz φ( ) Aexyz j, Acxyz j, c+( ) zd zj 0, zj 1, ∫= Fj xyz φ( ) Fxj φ( ) Fyj φ( ) Fzj φ( )⎩ ⎭⎨ ⎬ ⎧ ⎫T = Fj xyz φ( ) A0xyz j, φ( ) A1xyz j, φ( ) c⋅+= A0 xyz j, A0 x j, A0 y j, A0 z j, ⎩ ⎭⎨ ⎬ ⎧ ⎫T Ae xyz j, zd zj 0, zj 1,∫= = A1 xyz j, A1 x j, A1 y j, A1 z j, ⎩ ⎭⎨ ⎬ ⎧ ⎫T Ac xyz j, zd zj 0, zj 1,∫= = ⎭⎪ ⎪⎪ ⎬⎪ ⎪⎪ ⎫ φ φ t( )= Fxyz φ( ) Fjxyz φ( ) j 1= N ∑ A0xyz j, j 1= N ∑ A1xyz j, j 1= N ∑ c⋅+= = A0 xyz j, A1 xyz j,36 Chapter 3. Process Simulationengagement 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: , (3.38) where .dz and , from Eq (3.28). Similar to cutting forces, the overall torque is obtained by carrying out the inte- gration between engagement boundaries and adding the contribution of all flutes: (3.39) Since power is equal to torque times spindle speed (n), it can be trivially calculated as: , (3.40) where 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 . (3.41) dTorquej φ z,( ) R z( ) dFtj φ z,( )⋅= dFtj φ z,( ) Aet j, Act j, c⋅+( )= Aet j, Kte κsin⁄= Act j, Ktc fxy φjsin fz κtan⁄–( )⋅= Torque φ( ) Torquej φ z,( ) j 1= N ∑ R z( ) Aet j, Act j, c⋅+( ) z d⋅ zj 0, zj 1,∫ j 1= N ∑ A0 trq φ( ) A1trq+ φ( ) c⋅ = = = Power φ( ) n T⋅ orque φ( ) A0pwr φ( ) A1pwr φ( ) c⋅+= = Ai pwr φ( ) n Aitrq φ( )⋅= Fres φ( ) Fx φ( )( )2 Fy φ( )( )2+=37 Chapter 3. Process SimulationTaking the square of both sides and substituting Eq (3.37) for the x and y forces, the following form is obtained: , (3.42) where coefficients are derived as: , (3.43) and (p = 0,1; q = x,y). All of these additional process outputs will be analyzed more closely when they are used as constraints for feedrate optimization in the optimization chap- ter. 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-depen- dent 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 con- ducted 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 meth- ods are given in the following sub-sections. Fres φ( ) A0res φ( ) A1res φ( ) c⋅ A2res φ( ) c2⋅+ += A0 res A0 x( )2 A0y( ) 2 += A1 res 2 A0 x A1 x⋅ A0y A1y⋅+( )= A2 res A1 x( )2 A1y( ) 2 += ⎭⎪ ⎪⎬ ⎪⎪ ⎫ Ap q Ap q j, j ∑=38 Chapter 3. Process Simulation3.5.1. Orthogonal to Oblique Cutting Transformation Helical end milling is an oblique cutting technique where the cutting edge is no longer perpendic- ular to the cutting speed but inclined with an inclination angle ( ), 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 ( ), the shear angle ( ) and the friction angle ( ). These parameters can be deter- mined from orthogonal cutting tests for different rake angles ( ) 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. Using orthogonal cutting forces, namely feed/radial ( ) and tangential ( ), material cutting parameters, the shear stress, the shear angle and the friction angle can be obtained as: , (3.44) i0 τs φc βa αn η Workpiece Chip Tool V,F Vc Fr t Fa Chip-flow angle i0 Cutting edge inclination angle Tool Workpiece Chip Vc Flank face Rake face b h Orthogonal Cutting Oblique Cutting i0 Fr V,Ft Figure 3.5 : Orthogonal and oblique cutting geometries [64]. Frc Ftc τs Ftc φccos Frc φcsin– b hφcsin ------------- -------------------------------------------------=39 Chapter 3. Process Simulation, (3.45) , (3.46) where h is the uncut chip thickness, b is the width of cut, is the rake angle and is the chip compression ratio, which is defined as a function of uncut (h) and cut chip thickness ( ) as, . 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.e. , as pro- posed 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: , (3.47) where . In the oblique cutting transformation, the orthogonal shear angle is assumed to be equal to the normal shear angle in oblique cutting ( ); the normal rake angle is equal to the rake angle in orthogonal cutting ( ) and the chip flow angle is equal to the oblique angle ( = ). The edge cutting force coefficients in tangential and radial directions, i.e. , are calculated from the orthogonal cutting data by simply finding the force axis intercept φc rc αrcos 1 rc αrsin– --------------------------⎝ ⎠⎛ ⎞atan= βa αr Ffc Ftc ------⎝ ⎠⎛ ⎞atan+= αr rc hc rc h hc⁄= η i0 η i0= Ktc τs φnsin ------------- βn αn–( )cos i0 η βnsintantan+ φn βn αn–+( )cos 2 ηtan 2 βnsin 2+ --------------------------------------------------------------------------------------= Krc τs φnsin i0cos -------------------------- βn αn–( )sin φn βn αn–+( )cos 2 ηtan 2 βnsin 2+ --------------------------------------------------------------------------------------= Kac τs φnsin ------------- βn αn–( )cos i0tan η βnsintan– φn βn αn–+( )cos 2 ηtan 2 βnsin 2+ --------------------------------------------------------------------------------------= βntan βatan ηcos= φn φc≡ αn αr≡ η i0 Kte Kre,40 Chapter 3. Process Simulationof 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 ; 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 vari- ety of milling cutter geometries using three orthogonal parameters (shear stress, shear angle, fric- tion 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, . Since the total material removed per tooth is constant, the helix angle can be ignored for calculat- ing 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 . Finally, analytical average milling forces can be described as: , (3.48) Kae FMxyz φp Fxyz 1φp ----- g φj( )Fjxyz φ( ) j 1= N ∑⎝ ⎠⎜ ⎟ ⎜ ⎟⎛ ⎞ φd 0 φp∫ 1φp----- Fj=1 xyz φ( ) φd⋅ φst φex∫= =41 Chapter 3. Process Simulationwhere and are the cutter entry and exit angles, respectively and is the switching function defined as: . (3.49) Considering a cylindrical end mill with 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): . (3.50) Moreover, the immersion angle, , in and reduces to the following form when the helix angle is not in effect; . Using Eqs (3.32) and (3.33) with and planar milling feed conditions ( =1, =0), the average cutting forces in Eq (3.48) are calculated as: , (3.51) where . The above equation can be expressed in an alternative form, which isolates unknown cutting force coefficients: , (3.52) where , and the transformation matrix T(c) represents the effect of both cutting geometry and feedrate (c) and is given by: φst φex g φj( ) g φj( ) 1 if φst φj φex≤ ≤( ) 0 otherwise ⎩ ⎭⎨ ⎬ ⎧ ⎫ = κ z( ) π 2⁄= Fj xyz φ( ) Aexyz j, Acxyz j, c⋅+( ) a⋅= φj z( ) Aexyz j, Acxyz j, φj z( ) φj φ j 1–( )φp+= = φj 1= z( ) φ= fxy fz Fxyz Na2π------ 0.25 Ktc 2φcos Krc 2φ 2φsin–( )–( ) c⋅ Kte φsin– Kre φcos+( )+ 0.25 Krc 2φcos Ktc 2φ 2φsin–( )+( ) c⋅ Kre φsin Kte φcos+( )– Kac c φcos⋅ ⋅ Kae– φ φst φex = Fxyz Fx Fy Fz⎩ ⎭⎨ ⎬ ⎧ ⎫T = Fxyz c( ) T c( ) K⋅= K Kre Kte Kae Krc Ktc Kac⎩ ⎭⎨ ⎬ ⎧ ⎫T =42 Chapter 3. Process Simulation, (3.53) where , , , and . Averages of the measured cutting forces at m different feedrates are equated to analytically obtained average cutting forces: , (3.54) where is a 3m-by-6 matrix, and ’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: . (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, . For example, for a ball end mill: , (3.56) and . (3.57) T c( ) Na8π------ 4Δc1 4Δs1– 0 c 2Δp1 Δs2–( )– cΔc2 0 4Δs1– 4– Δc1 0 cΔc2 c 2Δp1 Δs2–( ) 0 0 0 4– Δp1 0 0 4cΔc1 = Δc1 φexcos φstcos–= Δs1 φexsin φstsin–= Δc2 2φexcos 2φstcos–= Δs2 2φexsin 2φstsin–= Δp1 φex φst–= FMxyz T K⋅= T T c1( ) T c2( ) . . T cm( )⎩ ⎭⎨ ⎬ ⎧ ⎫T = ci K T TT( ) 1– T TFMxyz= κ z( ) 1 κ z( )sin------------------ zd∫ Rb0 z Rb0–z 2Rb0 z–( )-------------------------------⎝ ⎠⎜ ⎟ ⎛ ⎞ atan Rb0 κ z( )( )cot–( )atan= = 1 κtan z( )------------------ zd∫ z 2Rb0 z–( )⋅=43 Chapter 3. Process SimulationNote that for general end mills, full or half immersion cases have to be used for cutting force coef- ficient 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 actu- ally removes material. The shape of the intersection zone varies as a function of tool and work- piece 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 visual- ization 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 horizon- tal 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 map- ping 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 Fig- ure 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 Simulationpreferred due to the smaller computational burden, deriving a generalized analytical relation is not always possible. 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 pre- φR z 0 πR Feed Direction φ R xy φ[rad] π0 z CEF CEF CEF CEP MAP Flute UNF OLD Figure 3.6 : Cutter - workpiece intersection, CEF mapping. 0 A xi al D ep th o f C ut Immersion Angle π Flute CEF Figure 3.7 : Generalized CEF.45 Chapter 3. Process Simulationdiction. In many other studies, the Z-buffer method has been extensively used for machining pro- cess 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 intersec- tion 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. Figure 3.8.a shows five distinct cases for flute-CEF intersection, one of which is detailed in Fig- ure 3.8.b. Parameters bounding the CEF, i.e , , , are provided by a CAD sys- tem hence known. For a given tool position , angular limits of the flute are determined Figure 3.8 : Rectangular CEF; (a) possible intersection cases, (b) parameter definition for CEFs. (a) (b) φ π0 z φst φex 1 0 3 2 4 φ π0 z zlim,0 (zj,0) zlim,1 zj,1 φj,0φj,1 zlim 0, zlim 1, φst φex φ φ t( )=46 Chapter 3. Process Simulationboth at the bottom and the top of the CEF using Eq (3.3), i.e., and . 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: In Figure 3.9, 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, , in Eq (3.3). For a cylindrical end mill, the lag angle, , is substituted and z is simply solved as: , (3.58) where 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 can only vary along either or line; however, such a useful property does not hold when more complex cutters like ball end mills are considered. φj 0, φj zlim 0,( )= φj 1, φj zlim 1,( )= φj,0>φst &φj,0<=φex YES ( ) NO ( )YES ( ) φj,1>φst &φj,1<=φex zj,1= z(φst) NO ( ) φj,0>φex &φj,1<φex YES ( ) NO ( ) NO ( )YES ( ) φj,1>φst zj,0 = z(φj,0) zj,1= z(φst)zj,0= z(φj,1) zj,1= z(φj,1) 10 2 3 4 zj,0 = 0.0 zj,1 = 0.0 zj,0 = z(φex) Figure 3.9 : Calculation of integration boundaries for general end mills. z φj( ) φj ψ z( ) z i0tan D 2⁄( )⁄= z φj( ) φ t( ) φpj φj–+( ) D⋅ 2 i0tan ------------------------------------------------ (Cylindrical End Mill)= i0 zj 0, zj 1,, zlim 0, zlim 1,,[ ]∈ φj φst= φj φex=47 Chapter 3. Process SimulationConsider 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 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. and ; 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, = , can be calculated using local cutter radius as follows: y x wyosR(z) + 0 60 120 180 1 2 3 4 5 6 7 8 9 φ[deg] z [m m ] D = 19.05 mm yos = 2 mm w = 8 mm φbnd(z) = φex(z) CEF D φbnd(z) = φst(z) x y z workpiece φst(z) φex(z) Figure 3.10 : (a) Planar machining with a ball end mill, (b) example CEF for ball end milling. (a) (b) yos φst φst z( )= φex φex z( )= φbnd z( ) φst z( ) φex z( ),{ }48 Chapter 3. Process SimulationIf 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 , and the variation of the integration boundaries as tool rotates is depicted in Figure 3.11. yval >= 0 YES ( ) R(z) <= yval φbnd (z) = 0 φbnd (z) = cos-1(yval / R(z)) R(z) <= |yval| φbnd (z) =π φbnd (z) = π/2 + sin-1(|yval| / R(z)) NO ( ) Start Angle: yval = yos + w/2, bnd = st Exit Angle: yval = yos - w/2, bnd = ex NO ( ) YES ( ) YES ( ) NO ( ) φbnd z( )49 Chapter 3. Process SimulationAssume that at a certain time step, flute j intersects with the CEF, curve, at axial height , which leads to the following relation at the intersection point: . (3.59) This equation might not be solved for analytically due to the presence of non-linear terms both in and . Redefining Eq (3.59) in a residue form: , (3.60) roots satisfying can be sought numerically. Brent’s numerical root finding algo- rithm [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 and when and 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, algo- rithm fine tunes the roots which lie between brackets. 0 π Immersion Angle Immersion Angle Advancement of the flute as cutter rotates Immersion Angle φbnd (z) Integration B oundaries z zmin zmax Figure 3.11 : Change of integration boundaries over time for the most general case. φbnd z( ) zroot φbnd zroot( ) φj zroot( ) φ t( ) j 1–( )φp ψ zroot( )–+= = zroot φbnd z( ) ψ z( ) g z( ) φbnd z( ) φ– t( ) j 1–( )φp– ψ z( )+= g zroot( ) 0= z0 z2 g z0( ) g z2( )50 Chapter 3. Process SimulationBrent’s method utilizes two methods for its fine tuning: the bisection method and inverse qua- dratic interpolation. The bisection method evaluates the function value at the midpoint of the bracket, checks its sign against and , and replaces the new point on the side that has the same sign as or . 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 func- tion (z = g( )) and determines the next root estimate at y=0 [92] as: , (3.61) where , i = 0, 1, 2, and is selected such that . Although a quadratic method has a higher convergence rate, its success depends solely on the smoothness of the func- tion. In conclusion, Brent’s method successfully combines the fast convergence of inverse qua- dratic 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 z0 z g(z) z1 z2 zroot,1 zroot,2 Figure 3.12 : Root bracketing. g z0( ) g z2( ) g z0( ) g z2( ) y2 z g0g1z2 g2 g0–( ) g2 g1–( ) ------------------------------------------- g1g2z0 g0 g1–( ) g0 g2–( ) ------------------------------------------- g2g0z1 g1 g2–( ) g1 g0–( ) -------------------------------------------+ += gi g zi( )= z1 z0 z1 z2< <51 Chapter 3. Process SimulationBrent’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 func- tion of axial height z: axial immersion angle , and lag angle . Except for cylindrical end mills, one or both of these terms bring non-linearity into the integrand. Moreover, cutting force coefficients, , 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: , (3.62) where g(z) is expressed by Eqs (3.32) and (3.33) depending on type of the machining state vari- able. Exact solution to Eq (3.62) can be approximated numerically using Trapezium rule as: , (3.63) where represent the Trapezoidal approximation with panels, i.e., step size of = , (3.64) and error estimate : κ z( ) ψ z( ) Krc Ktc Kac, , I g z( ) zd zroot 1, zroot 2,∫= g z( ) zd zroot 1, zroot 2,∫ Ik0≈ Ik 0 2k δk zroot 2, zroot 1,–( ) 2k⁄ O δ2( )52 Chapter 3. Process Simulation, (3.65) where I(0) is the approximate integration value, and are constants. A better approximation can be obtained by exploiting the error expansion. For example, consider two trapezoidal approx- imations with one has a half step size with respect to the other: , (3.66) . (3.67) The first error terms are eliminated by multiplying Eq (3.66) by 4 and subtracting it from Eq (3.67), an improved estimate is obtained with error estimate O( ): . (3.68) The method of error elimination is called Richardson extrapolation and is called the first Rich- ardson Extrapolant. In general, is the mth Richardson Extrapolant obtained by halving the step size, , at each stage, and the successive Richardson Extrapolants can be calculated by the follow- ing recursion: , (3.69) where , ,and M is the maximum number of refinements. The error term for this approximation becomes O( ). The recursive formula can also be defined in a matrix form as follows: Ik 0 I 0( ) c2δ2 c4δ4 …+ + += c2 c4 I1 0 I 0( ) c2δ2 c4δ4 …+ + += I2 0 I 0( ) 4c2δ2 16c4δ4 …+ + += δ4 I1 1 4I1 0 I2 0– 3 ----------------- I 0( ) 4 c4– δ4 …+= = I1 1 Ik m δ Ik m 1+ Ik 1+ m Ik 1+ m Ik m– 4m 1– ----------------------+= m 0 … M, ,= k 1 … m 1+, ,= δ2 m 1+53 Chapter 3. Process Simulation. 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-lineari- ties is required. In simulations, the maximum number of refinements, M, was limited to 20. Adap- tive 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 ele- ments 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 or . The singularity can appear only when since is not physically possible as seen in Figure (3.3). From Eq (3.22), angle of a ball end tool reduces to , (3.70) where , is the ball radius, and from the geometry of ball end: . (3.71) I1 0 I2 0 I1 1 I3 0 I2 1 I1 2 I4 0 I3 1 I2 2 I1 3 . . . . . . . . . . . Ik 0 Ik 1– 1 . . . . I1 m δ 1 κ z( )sin⁄ 1 κ z( )tan⁄ κ z( ) 0= κ z( ) π= κ κ z( )sin R z( ) Rb0⁄= R z( ) z 2Rb0 z–( )= Rb0 κcos z( ) Rb0 z–( ) Rb0⁄=54 Chapter 3. Process SimulationIf 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 nei- ther smooth nor non-singular at or near the tool tip. The following change of variable is proposed to relax singularity condition: , (3.72) and it is substituted into radius: . (3.73) Using the above conversion, new integrands are now obtained as: , (3.74) . (3.75) Note that at the tool tip, the relation 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. z Rb0 1 ξsin+( )= R ξ( ) Rb0 1 ξsin+( ) 2Rb0 Rb0 1 ξsin+( )–( ) Rb0 1 ξsin 2–( ) Rb0 ξcos= = = 1 κ z( )sin------------------dz 1 ξcos-----------Rb0 ξdξcos→ Rb0dξ= 1 κtan z( )------------------dz ξsin– ξcos-------------Rb0 ξdξcos→ Rb0 ξdsin– ξ= z 0=( ) ξ 3π 2⁄=( )→55 Chapter 3. Process Simulation3.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 loca- tions at a fixed tool position. This special case, shown in Figure 3.14.b, has to be taken into con- sideration during process simulation for accurate calculation of cutting forces. 0 2 4 6 -1.5 -1.0 -0.5 0 0.5 z [mm] -2.0 -1.5 -1 -0.5 0 -0.2 -0.1 0 0.1 0.2 ξ [rad] A e [N ] ( x 1e 3 ) y, j Figure 3.13 : Variation of edge integrand in y direction for a 12 mm ball end mill. Depth of Cut (a) Immersion Angle,( )φ #1 #2 #1 #2 Fl ut e N um be rs Chip Geometry zmax zmax 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. (a) (b)56 Chapter 3. Process SimulationConsider a case when cutter-workpiece engagement is composed of multi CEFs as shown in Fig- ure 3.15. The Cutter Engagement Plane (CEP) bounded between and defines the possible Cutting Region at the front of the cutter, where is the high- est 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, , 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. During process simulation, flute re-cut condition can be detected by checking the number of occurrence: , (3.76) where , (3.77) z zmin 0= zmax,[ ]∈ φ 0 π,[ ]∈ zmax ψ π zmax φ[rad] π0 zmax z Cutting Region φ[rad] π0 CEF-1 CEF-2 CEF-3 z CEF-1 CEF-2 CEF-3 Figure 3.15 : Cutting region at the front of the cutter. nocc npi 2 ------ 1 if npi is even+ npi 1–( ) 2 -------------------- 1 if npi is odd+⎩⎪ ⎪⎨ ⎪⎪ ⎧ = npi ceil ψ zmax( ) π----------------------⎝ ⎠ ⎛ ⎞=57 Chapter 3. Process Simulationand 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 >1. Figure 3.16 presents different cutting configura- tions 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., . 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 coefficients: = 796 [N/ ], = 169 [N/ ], = 222 [N/ ], = 28 [N/mm], = 31 [N/mm], and = 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. nocc φ φ π 2π,[ ]∈ Ktc mm 2 Krc mm 2 Kac mm 2 Kte Kre Kae58 Chapter 3. Process SimulationReal Flute-Position Equivalent Flute-Positions 34 2 1 π 2π-6π -5π -4π -3π -2π -π 0 π 2π π 2π π 2π 34 2 1 π 2π-6π -5π -4π -3π -2π -π π 2π π 2π 34 2 1 π 2π-6π -5π -4π -3π -2π -π 4ππ 2π π 2π π 2π 34 2 1 π 2π-6π -5π -4π -3π -2π -π φ φ φ φ φ φ φ φ φ φ φ φ z z z z z z z z z z z z (a) (b) (c) (d) npi noccSection 1 1 2 2 3 2 4 3 a b c d + + + + = = = = Flute idF!Fflute idF@ idF#+ += flute flute flute+... Figure 3.16 : Flute re-cut example, same depth, different helix angles.59 Chapter 3. Process SimulationA 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. Cut- ting 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- φ z Immersion Angle Ax ia l D ep th o f C ut [m m ] Tool Rotation [deg] -300 180 3600 -200 0 -100 100 -200 0 200 400 600 -500 500 1000 0 Fz [N ] Fy [N ] Fx [N ] Measured Simulated Cutter 3 7 11 feed direction FO R C E MAP Figure 3.17 : Multiple immersion milling, engagement map, simulated versus measured cutting forces.60 Chapter 3. Process Simulationlog, 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. 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 cut- ting 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 mag- nitudes between simulation and actual cut is due to the trajectory generation, and the limited max- imum feed speed of the machine, which was set to 5000 [mm/min]. Table 3.1 : Cutting force coefficients of Al 6061 Tool Type [N/mm] [N/mm] [N/ ] [N/ ] 10 [mm] 38.38 -9.792 403.396 98.911 25 [mm] 21.34 23.141 662.296 72.86 Kte Kre Ktc mm2 Krc mm261 Chapter 3. Process SimulationFigure 3.18 : (a) Face milling, (b) right side profiling, (c) top hole helical milling, (d) middle hole helical milling.62 Chapter 3. Process Simulation0 5 10 15 20 25 30 33 0 5 10 15 20 25 30 35 40 45 -800 -400 0 400 800 time [s] (a) -800 -400 0 400 800 Fx [N ] Fy [N ] time [s] (b) -800 -400 0 400 800 Fx [N ] -800 -400 0 400 800 Fy [N ] Figure 3.19 : Face milling a) measured, (b) simulated cutting forces.63 Chapter 3. Process Simulation0 1.5 3.5 5.5 7.5 9.5 11.5 12 time [s] time [s] (b) (a) -1500 -1000 -500 0 500 1000 1500 Fx [N ] 0 2 4 6 8 10 12 14 16 18-1500 -1000 -500 0 500 1000 1500 Fy [N ] -1500 -1000 -500 0 500 1000 1500 Fx [N ] -1500 -1000 -500 0 500 1000 1500 Fy [N ] Figure 3.20 : Profile milling a) measured, (b) simulated cutting forces.64 Chapter 3. Process SimulationIn co-operation with Sandvik Coromant, Sweden, a die was machined to verify the proposed algo- rithms. 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 pro- vided by Uddeholm. Different views of both stock (casting) and final part are given in Figure 3.22.a-d. -2000 -1000 0 1000 2000 Fx [N ] 0.1 0.3 0.5 0.7 0.9 1.1 1.3 1.5 1.7 -2000 -1000 0 1000 2000 time [s] Fy [N ] -2000 -1000 0 1000 2000 Fx [N ] 0 0.5 1 1.5 2 2.5 3 3.5 -2000 -1000 0 1000 2000 time [s] Fy [N ] (a) (b) Figure 3.21 : Top hole helical milling a) measured, (b) simulated cutting forces.65 Chapter 3. Process SimulationThe 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]. Cut- ting 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 (a) (b) (c) (d) Figure 3.22 : Test part stock and final shape (a) isometric; (b) left side; (c) front; (d) top views.66 Chapter 3. Process Simulationwas 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 cut- ting 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 SimulationIn 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. 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 cur- rently 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 ( ) and spindle period ( ) Ktc [N/mm²] Krc [N/mm²] Kac [N/mm²] Kte [N/mm] Kre [N/mm] Kae [N/mm] T63 2312.7 1443.1 -167.5 62.82 112.74 20.15 T20 2294.2 1137.5 -536 63.53 100.24 5.45 T20F 2024.1 1132.9 42.2 53.73 64.2 1.54 T16 2068.1 1120.9 -99.54 41.4 72.63 19.71 T8 2156 924.1 172.9 29.44 19.75 1.24 Table 3.4 : Average cutting force coefficients for Uddeholm Carmo - Carmo Cast. tout tperiod68 Chapter 3. Process Simulationare calculated as , . Figure 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, there- fore, 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 thick- ness 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. tout 60 65( ) 360 1011( )------------------------- 10.7 μs[ ]= = tperiod 60 360( ) 360 1011( )------------------------- 59.3 μs[ ]= =69 Chapter 3. Process Simulation- 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 fee- drate estimation was directly reflected on cutting forces. - 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 iden- tical for simulations. Lastly, tool wear over time affecting the edge geometry hence C u tt in g F o rc e C o e ff ic ie n t Axial Depth of Cut Average Figure 3.23 : Cutting force coefficient variation. F e e d ra te Along Tool Path Commanded Actual Figure 3.24 : Actual versus commanded feedrate.70 Chapter 3. Process Simulationresulting in increase in cutting force coefficients was not taken into account. In all simu- lations, coefficients identified with sharp tools were used.71 Chapter 3. Process Simulation 0 20 40 60 80 100 120 140 160 180 -1000 0 1000 Fx [N ] VMS Simulated Forces 0 20 40 60 80 100 120 140 160 180 -1000 0 1000 Fy [N ] 0 20 40 60 80 100 120 140 160 180 -500 0 500 time [s] Fz [N ] 0 20 40 60 80 100 120 140 160 180 -1000 0 1000 Measured Forces Fx [N ] 0 20 40 60 80 100 120 140 160 180 -1000 0 1000 Fy [N ] 0 20 40 60 80 100 120 140 160 180 -500 0 500 Fz [N ] time [s] Figure 3.25 : T63 - Operation 3, finishing of the walls at z = 105.26 [mm]. Maximum Minimum72 Chapter 3. Process Simulation 35.94 35.96 35.98 36 36.02 36.04 36.06 -200 0 200 400 600 800 Measured Forces Fx [N ] 35.96 35.98 36 36.02 36.04 36.06 -500 0 500 1000 Fy [N ] 35.96 35.98 36 36.02 36.04 36.06 -200 0 200 400 Fz [N ] time [s] D ep th o f C ut [m m ] (a) Immersion Angle [deg] (b) Spindle Period 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]. tout73 Chapter 3. Process Simulation 0 20 40 60 80 100 120 140 160 180 -2000 0 2000 Fx [N ] VMS Simulated Forces 0 20 40 60 80 100 120 140 160 180 -2000 0 2000 Fy [N ] 0 20 40 60 80 100 120 140 160 180 -1000 0 1000 2000 3000 time [s] Fz [N ] 0 20 40 60 80 100 120 140 160 180 -2000 0 2000 Measured Forces Fx [N ] 0 20 40 60 80 100 120 140 160 180 -2000 0 2000 Fy [N ] 0 20 40 60 80 100 120 140 160 180 -1000 0 1000 2000 3000 Fz [N ] time [s] Figure 3.27 : T63-Operation 1, roughing of the top part at z=163.5 [mm].74 Chapter 3. Process Simulation0 10 20 30 40 50 60 70 80 90 100 -2000 0 2000 Measured Forces Fx [N ] 0 10 20 30 40 50 60 70 80 90 100 -2000 0 2000 Fy [N ] 0 10 20 30 40 50 60 70 80 90 100 -1000 0 1000 2000 3000 Fz [N ] time [s] 0 10 20 30 40 50 60 70 80 90 100 -2000 0 2000 Fx [N ] VMS Simulated Forces 0 10 20 30 40 50 60 70 80 90 100 -2000 0 2000 Fy [N ] 0 10 20 30 40 50 60 70 80 90 100 -1000 0 1000 2000 3000 time [s] Fz [N ] Figure 3.28 : T63 - Operation 1, roughing at z = 88.51 [mm].75 Chapter 3. Process Simulation 3 4 5 6 7 8 9 10 11 12 13 -1000 0 1000 2000 Measured Forces Fx [N ] 3 4 5 6 7 8 9 10 11 12 13 -1000 0 1000 2000 Fy [N ] 3 4 5 6 7 8 9 10 11 12 13 0 1000 2000 Fz [N ] time [s] 5 10 15 20 25 30 -1000 0 1000 2000 Fx [N ] VMS Simulated Forces 5 10 15 20 25 -1000 0 1000 2000 Fy [N ] 5 10 15 20 25 30 0 1000 2000 time [s] Fz [N ] Figure 3.29 : T20 - Operation 2, semi-finishing of the top.76 Chapter 3. Process Simulation 0 2 4 6 8 10 12 14 16 18 20 -1000 0 1000 Measured Forces Fx [N ] 0 2 4 6 8 10 12 14 16 18 20 -2000 0 2000 Fy [N ] 0 2 4 6 8 10 12 14 16 18 20 0 1000 2000 Fz [N ] time [s] 265 270 275 280 285 290 -1000 0 1000 Fx [N ] VMS Simulated Forces 265 270 275 280 285 290 -2000 0 2000 Fy [N ] 265 270 275 280 285 290 0 1000 2000 time [s] Fz [N ] 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 excit- ing 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 mill- ing 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 eigen- value search method is proposed to accelerate iterative calculation of stability lobes when the Multi Frequency solution is used.78 Chapter 4. Dynamics of Milling Operations4.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 multi- degree-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 nat- ural 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 fre- quencies. 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 geome- try, 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 OperationsIn general, the cutter geometry can be defined by a unit vector perpendicular to the cutter body (Figure 4.2) for flute j as: , (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 is defined as: , (4.2) where , , , and is the tooth period defined as: , (4.3) where is the tooth passing frequency: , (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 vec- tor perpendicular to the cutter geometry: . (4.5) nj z( ) κ z( ) φj z( ) i⋅ κ z( ) φj z( ) j⋅ κ z( ).kcos–cos⋅sin+sin⋅sin= Δ t( ) Δ t( ) Δx t( ) i⋅ Δy t( ) j⋅ Δz t( ) k⋅+ += Δx t( ) x t( ) x t τ–( )–= Δy t( ) y t( ) y t τ–( )–= Δz t( ) z t( ) z t τ–( )–= τ τ 2π ωT⁄= ωT ωT N Ω⋅= Ω hd j, z t,( ) nj z( ) Δ t( )•=80 Chapter 4. Dynamics of Milling OperationsFor 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 lin- ear edge force model [27] as: , (4.6) κ(z) z x x y dF r, j dFt, j Chip(t,φ,κ) dz dS Dynamic Chip dFa, j φ Present dynamic displacement Previous dynamic displacement [x(t-T), y(t-T)] [x(t), y(t)] Cross section at elevation z st φex φj Δ(t) nj(z) hd, j Figure 4.2 : Undulations on the chip, and dynamic chip thickness. dFr j, t z,( ) dFt j, t z,( ) dFa j, t z,( )⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ Kre Kte Kae⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ dS z( )⋅ Krc Ktc Kac⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ hd j, t z,( ) dS z( )⋅ ⋅+=81 Chapter 4. Dynamics of Milling Operationswhere dS(z) is the differential contact length expressed as: , (4.7) where dz is the differential axial depth of cut, are the cutting force coefficients due to the shear and 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. Cutting forces acting on flute j are then calculated as: , (4.8) dS z( ) dz κ z( )sin⁄= Krc Ktc Kac, , Kre Kte Kae, , κ(z) Ball End MillCylindrical End Mill nj(z) nj(z) = sin φj(t,z) . i + cos φj(t,z) . j nj(z) = sin κ(z) . sin φj(t,z) . i + sin κ(z) . cos φj(t,z) . j - cos κ(z) . k nj(z)x z x z Vibration Vector Cutter Type Unit Geometry Vector Dynamic Chip Thickness Δ= Δx . i + Δy . j + Δz . k Δx . sin φj(t,z) + Δy . cos φj(t,z) Δ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. Fr j, t( ) Ft j, t( ) Fa j, t( )⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ Kre Kte Kae⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ a κsin----------⋅ Krc Ktc Kac⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ hd j, t( ) aκsin----------⋅ ⋅+=82 Chapter 4. Dynamics of Milling Operationswhere a is the axial depth of cut, is the average axial immersion angle, and is the aver- age dynamic chip thickness obtained from Eq (4.5) after replacing with in Eq (4.1). For com- plex 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 . 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: , (4.9) where total rotating cutting forces are calculated by adding the contribution of each flute as: . (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: κ hd j, t( ) κ κ κ π 4⁄= Fx t( ) Fy t( ) Fz t( )⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ κ φj t( )sinsin– φj t( )cos– κ φj t( )sincos– κ φj t( )cossin– φj t( )sin κ φj t( )coscos– κcos 0 κsin– Fr t( ) Ft t( ) Fa t( )⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ = Fr t( ) Ft t( ) Fa t( )⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ Fr j, t( ) Ft j, t( ) Fa j, t( )⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ j 1= N ∑=83 Chapter 4. Dynamics of Milling Operations, (4.11) where time-varying directional coefficients in each principal direction are obtained as: , (4.12) , (4.13) Fx t( ) Fy t( ) Fz t( )⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ 1 2 -- a Ktc axx t( ) axy t( ) axz t( ) ayx t( ) ayy t( ) ayz t( ) azx t( ) azy t( ) azz t( ) Δx t( ) Δy t( ) Δz t( )⎩ ⎭⎪ ⎪ ⎨ ⎬⎪ ⎪ ⎧ ⎫ ⋅ ⋅ ⋅= axx t( ) axx j, j 1= ∑ gj 2φjsin 1 2φjcos–( ) Kr κsin Ka κcos+( )+[ ]⋅– j 1= ∑= = axy t( ) axy j, j 1= N ∑ gj 1 2φjcos+( ) 2φj Kr κsin Ka κcos+( )sin+[ ]⋅– j 1= N ∑= = axz t( ) axz j, j 1 N ∑ 2gj φjcos κtan⁄ φj κcossin Kr Ka κtan⁄+( )+[ ]⋅ j 1 N ∑= = ayx t( ) ayx j, j 1= N ∑ gj 1 2φjcos–( ) 2φj Kr. κsin Ka. κcos+( )sin–[ ]⋅ j 1= N ∑= = ayy t( ) ayy j, j 1= N ∑ gj 2φjsin 1 2φjcos+( ). Kr. κsin Ka. κcos+( )–[ ]⋅ j 1= N ∑= = ayz t( ) ayz j, N ∑ 2gj φjsin– κtan⁄ φj. κcoscos . Kr Ka κtan⁄+( )+[ ]⋅ N ∑= =84 Chapter 4. Dynamics of Milling Operations, (4.14) where is the switching function which is equal to unity while cutting and zero other- wise, and . One important observation about directional coeffi- cients 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: , (4.15) where the vibration difference vector is , the dynamic cutting forces vector is and the directional coefficient matrix is . (4.16) azx t( ) azx j, j 1= ∑ 2gj φj. κ.cossin Kr φj. κ.Kasinsin–[ ]⋅ j 1= ∑= = azy t( ) azy j, j 1= N ∑ 2gj. φjcos . κ.cos Kr φjcos . κ.Kasin–[ ] j 1= N ∑= = azz t( ) azz j, j 1 N ∑ 2gj. κcos– κtan⁄ .Kr κ.Kacos+[ ] j 1 N ∑= = gj g φj( )= Kr Krc Ktc⁄= Ka Kac Ktc⁄= f t( ) 12-- a Ktc A t( ) Δ t( )⋅ ⋅ ⋅ ⋅= Δ t( ) Δx Δy Δz, ,{ }T= f t( ) Fx Fy Fz, ,{ }T= A t( ) Aj t( ) j 1= N ∑ axx j, axy j, axz j, ayx j, ayy j, ayz j, azx j, azy j, azz j,j 1= N ∑= =85 Chapter 4. Dynamics of Milling Operations4.3. Stability of Milling 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 cut- ting tool and applying cutting forces as external excitations: , (4.17) where represent physical degrees of freedom, is the excitation force vector, and n is the number of degrees of freedom. M, C, K are n-dimen- fx1 fx2 fxn . . . x1 x2 xn Figure 4.4 : Multi degree-of-freedom milling dynamics. M x·· t( )⋅ C x· t( )⋅ K x t( )⋅+ + fx= x t( ) x1 x2 … xn, , ,{ }T= fx fx1 fx2 … fxn, , ,{ }T=86 Chapter 4. Dynamics of Milling Operationssional mass, damping and stiffness matrices, respectively. Physical and modal displacements are related by mode shapes as: , (4.18) where are modal displacements corresponding to each mode, and is the number of modes, respectively. The mode shape matrix (n-by- ) is defined as: . (4.19) Each column of the mode shape matrix ( ) represents a ratio of structural displacement ampli- tudes 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: , (4.20) where , and are modal mass and stiffness matrices, and is the modal force vector: x Px . μ= μ μ1 μ2 … μmx, , ,{ } T= mx mx Px px 11, px 12, . . px 1mx, px 21, px 22, . . px 2mx, . . . . . . . . . . px n1, px n2, . . px nmx, = Px Mμ μ·· t( )⋅ Cμ μ· t( )⋅ Kμ μ t( )⋅+ + fxμ= Mμ Px T.M.Px= Kμ Px T.K.Px= fxμ87 Chapter 4. Dynamics of Milling Operations= = . (4.21) When the system has a proportional damping, the modal damping matrix is expressed as . A similar analysis can be performed for the rest of the orthogonal directions (y and z): , (4.22) , (4.23) where = and = are modal force vectors of length and ; and are the number of modes; and are mode shape matrices of size (n-by- ) and (n-by- ); and and are excitation force vectors in y and z directions, respectively. Using mode shapes, physical and modal displacements are related in y and z directions: , (4.24) , (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: fxμ Px Tfx px 11, fx1 px 21, fx2 … px n1, fxn+ + + px 12, fx1 Px 22, fx2 … px n2, fnx+ + + … px 1mx, fx1 px 2mx, f2x … px nmx, fxn+ + +⎩ ⎭⎪ ⎪ ⎪ ⎪⎨ ⎬ ⎪ ⎪⎪ ⎪ ⎧ ⎫ Cμ Px T.C.Px= Mη η·· t( )⋅ Cη η· t( )⋅ Kη η t( )⋅+ + fyη= Mν ν·· t( )⋅ Cν ν· t( )⋅ Kν ν t( )⋅+ + fzν= fyη Py Tfy fzν Pz Tfz my mz my mz Py Pz my mz fy fz y Py η⋅= z Pz ν⋅= η ν88 Chapter 4. Dynamics of Milling Operations, (4.26) where is the transformation matrix between physical and modal forces and defined as: , (4.27) , (4.28) , (4.29) , (4.30) where is the total number of modes, and and 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): ; (4.31) MΓ Γ·· t( )⋅ CΓ Γ· t( )⋅ KΓ Γ t( )⋅+ + Tfull fall t( )⋅= Tfull Tfull Px T 0 0 0 Py T 0 0 0 Pz T mt mt× = fall t( ) fx t( ) fy t( ) fz t( )⎩ ⎭⎨ ⎬ ⎧ ⎫T = Γ t( ) μ1 μ2 … μmx η1 η2 … ηmy ν1 ν2 … νmz, , , , , , , , , , ,{ } T= MΓ Mμ 0 0 0 Mη 0 0 0 Mν mt mt× = mt mx my mz+ += CΓ KΓ px 11, px 12, … px 1mx, 1= = = = py 11, py 12, … py 1my, 1= = = = pz 11, pz 12, … pz 1mz, 1= = = = ⎭⎪ ⎪⎬ ⎪⎪ ⎫89 Chapter 4. Dynamics of Milling Operationsand forces be applied only at the tool tip: . (4.32) The equation of motion now reduces to: , (4.33) where T is the reduced transformation matrix: , (4.34) and f(t) is the cutting force vector: . (4.35) Note that varying dynamics along the depth of cut can be considered by simply considering the distributed forces ( ) 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: fx1 Fx t( ) fx2, 0 … fxn, , 0= = = fy1 Fy t( ) f2y, 0 … fyn, , 0= = = fz1 Fz t( ) fz2, 0 … fzn, , 0= = = ⎭⎪ ⎪⎬ ⎪⎪ ⎫ MΓ Γ·· t( )⋅ CΓ Γ· t( )⋅ KΓ Γ t( )⋅+ + T f t( )⋅= T 1{ }mx 0 0 0 1{ }my 0 0 0 1{ }mz mt 3× = f t( ) Fx t( ) Fy t( ) Fz t( )⎩ ⎭⎨ ⎬ ⎧ ⎫T = fxn fyn fzn, , 0≠90 Chapter 4. Dynamics of Milling Operations. (4.36) Substituting Eq (4.36) into Eq (4.33), and reorganizing the terms: , , (4.37) where , (4.38) the general equation of motion in Eq (4.37) can be expressed in a state space form as: , (4.39) where , (4.40) , (4.41) f t( ) 12-- a Ktc A t( ) T T Γ t( ) Γ t τ–( )–{ }⋅ ⋅ ⋅ ⋅ ⋅= MΓ Γ·· t( )⋅ CΓ Γ· t( )⋅ KΓ Γ t( )⋅+ + T 12-- a Ktc A t( ) T T Γ t( ) Γ t τ–( )–{ }⋅ ⋅ ⋅ ⋅ ⋅ ⋅= MΓ Γ·· t( )⋅ CΓ Γ· t( )⋅ KΓ LΓ t( )–( ) Γ t( )⋅+ + LΓ t( ) Γ t τ–( )⋅–= LΓ t( ) 12-- a Ktc T A t( ) T T⋅ ⋅ ⋅ ⋅ ⋅= Θ· t( ) U t( ) Θ t( )⋅ V t( ) Θ t τ–( )⋅+= Θ t( ) Γ t( ) Γ· t( )⎩ ⎭⎨ ⎬ ⎧ ⎫ = U t( ) 0 I MΓ 1– LΓ t( ) KΓ–( ) MΓ– 1– CΓ =91 Chapter 4. Dynamics of Milling Operations. (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 ham- mer 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 intro- duces periodicity into the equation of motion. In milling, the delay and the periodicity of the sys- tem 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). V t( ) 0 0 MΓ– 1– LΓ t( ) 0 = τ92 Chapter 4. Dynamics of Milling OperationsThe 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 differen- tial 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 inter- val division as shown in Figure 4.5. Delay time (tooth period) is divided into finite number of ele- ments of length such that , (4.43) where k is an integer and the time at each interval is expressed as: , . (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 discret- ization interval : , , (4.45) Δt Δt T k⁄= ti i Δt⋅= i 0 1 … k, , ,= [ti ti 1+ ), Ui 1 Δt------ U t( ) td ti ti 1+ ∫= Vi 1Δt------ V t( ) td ti ti 1+ ∫=93 Chapter 4. Dynamics of Milling Operationsand the delayed term is approximated as a weighted linear combination of the delayed discrete values and : , (4.46) where 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 substi- tuted into non-autonomous DDE expressed in Eq (4.39), the system becomes a linear autonomous ODE within a time interval: , . (4.47) Θi k– Θi k– 1+ Θ t τ–( ) Θ ti kΔt–( ) wb Θi k–⋅ wa Θi k– 1+⋅+≈ ≈ wa wb 1 2⁄= = Θ· t( ) Ui Θ t( )⋅ Vi wb Θi k–⋅ wa Θi k– 1+⋅+( )⋅+= ti t ti 1+<≤( ) ti t Q ti+1ti-k ti-k+1 T=k*Δt Δt 2Δt kΔt 0 1 2 k 0 Counter Time Discretization of the delay (tooth period) Approximation of the delayed term Figure 4.5 : Approximation of the delayed term using the semi-discretization (SD) by Stepan and Insperger [73].94 Chapter 4. Dynamics of Milling OperationsBeing treated as an initial value problem, i.e. , Eq (4.47) has an exact solution of the form: . (4.48) Substituting in Eq (4.48), first, the state at the exit of the ith time interval is found as: , , (4.49) and then a discrete map is defined as: , (4.50) where , (4.51) , (4.52) are ( -by- ) matrices; Θ ti( ) Θi= Θ t( ) eUi t ti–( ) Θi⋅ eUi t ti–( ) I–( ) Ui 1– Vi wb Θi k–⋅ wa Θi k– 1+⋅+( )⋅ ⋅ ⋅+= t ti 1+= Θi 1+ Si Θi⋅ Ri wb Θi k–⋅ wa Θi k– 1+⋅+( )⋅+= ti t ti 1+<≤( ) ϒi 1+ full Ψifull ϒifull= Si e Ui Δt= Ri Si I–( ) Ui 1– Vi⋅ ⋅= 2mt 2mt95 Chapter 4. Dynamics of Milling Operations(4.53) is a [ -by- ] matrix; and (4.54) is a vector of length ( ). It is clear from Eq (4.42) that does not play any role in Eq (4.39); therefore, in Eq (4.49) does not depend on or . Using this observation, size of the problem can be reduced and a new discrete map is defined as: , (4.55) where is a square matrix of size .(k+2): Ψi full Si 0 0 … 0 wa Ri⋅ wb Ri⋅ I 0 0 … 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 = 2mt k 1+( ) 2mt k 1+( ) ϒi full Γi Γ· i Γi 1– Γ· i 1– … Γi k– Γ· i k–⎩ ⎭⎨ ⎬ ⎧ ⎫T = 2mt k 1+( ) Γ· t τ–( ) Θi 1+ Γ· i k– Γ· i k– 1+ ϒi 1+ Ψi ϒi⋅= Ψi mt96 Chapter 4. Dynamics of Milling Operations(4.56) and is a vector of length .(k+2): . (4.57) Note that additional subscripts of and 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: , (4.58) where . 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, , is also periodic Ψi Si 1…mt 1…mt,( ) Si 1…mt mt+1…2mt,( ) 0 0 … … waRi 1…mt 1…mt,( ) wbRi 1…mt 1…mt,( ) Si mt+1…2mt 1…mt,( ) Si mt+1…2mt mt+1…2mt,( ) 0 0 … … waRi mt+1…2mt 1…mt,( ) wbRi mt+1…2mt 1…mt,( ) I 0 … … … … 0 0 0 0 I 0 … … 0 0 0 0 0 I 0 … 0 0 . . . . . . . . . . . . . . . . . . . . . . . . 0 0 0 0 0 … I 0 = ϒi mt ϒi Γi Γ· i Γi 1– … Γi k 1–– Γi k–⎩ ⎭⎨ ⎬ ⎧ ⎫T = Si Ri Φ Ψk 1– Ψk 2– …Ψ0⋅= ϒk Φ ϒ0⋅= Φ A t( )97 Chapter 4. Dynamics of Milling Operationsat the tooth passing frequency ( ) or at the tooth period ( ), i.e., = . Harmonic functions can be expressed by Fourier series [67] as: , (4.59) where i is the imaginary unit, and the Fourier coefficients are defined as: . (4.60) In order to obtain Fourier coefficients, first, Eq (4.16) is substituted into Eq (4.60): . (4.61) Since , where is the pitch angle for an even-pitch cutter, the summation in the above relation can be taken out by changing integration limits: , (4.62) ωT τ A t( ) A t τ+( ) A t( ) Ar eirωTt⋅ r ∞–= ∞ ∑= Ar 1 τ-- A t( ) e irωTt– td⋅ 0 τ ∫= Ar 1 τ-- Aj t( ) j 1= N ∑⎝ ⎠⎜ ⎟ ⎜ ⎟⎛ ⎞ e irωTt– td⋅ 0 τ ∫ 1τ-- Aj t( ) j 1= N ∑⎝ ⎠⎜ ⎟ ⎜ ⎟⎛ ⎞ e irNΩt– td⋅ 0 τ ∫= = φj φ j 1–( )φp+ Ω t j 1–( ) τ⋅+( )= = φp 2π N⁄= Ar 1 τ-- Aj γ( ) e irNΩγ– γ d⋅ jτ j 1+( )τ ∫⎝ ⎠⎜ ⎟ ⎜ ⎟⎛ ⎞ j 1= N ∑⋅=98 Chapter 4. Dynamics of Milling Operationswhere 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: . (4.63) Since contains the switching function , the above expression can be simplified by taking all flutes into account within the immersion zone: , (4.64) where is defined equal to in Eq (4.16) without the dependency on flute number j, i.e., being unity and . The right-hand-side of Eq (4.15) is a function of not only the tooth passing frequencies due to directional coefficients , but also the vibration (chatter) frequency due to vibration terms, . Thus, as a reaction, dynamic milling forces ( ) should contain both the tooth passing frequencies ( ) and the chatter frequency, when the system is marginally stable. In gen- eral, Floquet’s theorem [66][98][99] states that for a second order differential equation with peri- odic and piece-wise continuous coefficients, the solution has the following form: , (4.65) γ θ Ωτ= Ar 1 τΩ------- Aj θ( ) e irNθ– θ d⋅ jφp j 1+( )φp ∫⎝ ⎠⎜ ⎟ ⎜ ⎟⎛ ⎞ j 1= N ∑⋅= Aj θ( ) gj Ar N 2π------ A θ( ) e irNθ– θd⋅ φst φex ∫= A θ( ) Aj t( ) gj φj θ= A t( ) Δ t( ) f t( ) k ωT⋅ ωc f t( ) esλt p t( )⋅=99 Chapter 4. Dynamics of Milling Operationswhere is a periodic function of the tooth passing frequency ( ). It is evident from Eq (4.65) that the milling system is stable if the real part of exponent is negative. Since the cutter and the workpiece vibrate at the chatter frequency on the stability border, can be replaced by to ensure marginal stability [67]. When periodic function is also expanded by using Fourier series, the following expression is obtained for the dynamic milling forces: . (4.66) 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 can be expressed as: , (q = x,y,z), (4.67) where D is the differential operator , G(D) is the transfer function expressed in terms of modal parameters such as natural frequency ( ), modal mass (m) and damping ratio ( ): , (4.68) 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 writ- ten as: p t( ) ωT sλ sλ sλ iωc= p t( ) f t( ) eiωct pk eikωTt k ∞–= ∞ ∑ pk ei ωc kωT+( )t k ∞–= ∞ ∑= = fq t( ) q D( ) G D( ) fq t( )⋅= d …( ) dt⁄ ωn ζ G D( ) 1 mj⁄ D2 2 ζj ωn j, D⋅ ⋅ ⋅ ωn j,2+ + --------------------------------------------------------------- j 1= M ∑=100 Chapter 4. Dynamics of Milling Operations, (4.69) where is the overall transfer function matrix: . (4.70) 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 trans- fer function matrix becomes a diagonal matrix. The vibration difference vector is expressed by using the transport delay (tooth period, ) as: . (4.71) Substituting Eqs (4.66), (4.69), and (4.71) into Eq (4.15) yields . (4.72) Using the definition of parametric transfer functions of linear periodic systems given in reference [100], Minis and Yanushevsky [66] expressed Eq (4.72) as: x t( ) y t( ) z t( )⎩ ⎭⎪ ⎪ ⎨ ⎬⎪ ⎪ ⎧ ⎫ G D( ) f t( )⋅= G D( ) G D( ) Gxx D( ) Gxy D( ) Gxz D( ) Gyx D( ) Gyy D( ) Gyz D( ) Gzx D( ) Gzy D( ) Gzz D( ) = G D( ) τ Δ t( ) 1 e Dτ––( ) x t( ) y t( ) z t( )⎩ ⎭⎪ ⎪ ⎨ ⎬⎪ ⎪ ⎧ ⎫ ⋅= f t( ) 12-- a Ktc A t( ) 1 e Dτ––( ) G D( ) pk ei ωc kωT+( )t⋅ ⋅ k ∞–= ∞ ∑⋅ ⋅ ⋅ ⋅=101 Chapter 4. Dynamics of Milling Operations. (4.73) Eq (4.65) is substituted into the left-hand-side of Eq (4.73), and both sides are divided by to obtain: , (4.74) where is replaced by as . Multiplying both sides of the above equation by , integrating from 0 to , and using the orthogonality principle [67], the following is obtained: which can be simplified to , , (4.75) where . 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 , and the summa- tion in Eq (4.75) is expressed in matrix form, then the problem is reduced to an eigenvalue prob- lem: f t( ) 12-- a Ktc A t( ) e i ωc kωT+( )t 1 e i ωc kωT+( )τ––( ) G i ωc kωT+( )( ) pk ⋅ ⋅ ⋅ k ∞–= ∞ ∑⋅ ⋅ ⋅ ⋅= e iωct p t( ) 12-- a Ktc 1 e iωcT––( ) eikωTt A t( ) G i ωc kωT+( )( ) pk⋅ ⋅ ⋅ k ∞–= ∞ ∑⋅ ⋅ ⋅= e i ωc kωT+( )τ– e iωcτ– ωT τ⋅ 2π= 1 τ e ipωTt–⋅⁄ τ 1 τ-- p t( )e ipωTt– td 0 τ ∫ 12--a Ktc 1 e iωcτ––( ) 1τ-- A t( ) e i– p k–( )ωTt td⋅ 0 τ ∫⎝ ⎠⎜ ⎟ ⎜ ⎟⎛ ⎞ G i ωc kωT+( )( ) pk⋅ ⋅ k ∞–= ∞ ∑⋅ ⋅= pp 1 2 --a Ktc 1 e iωcτ––( ) Ap k– G iωk( ) pk⋅ ⋅ k ∞–= ∞ ∑⋅ ⋅= p 0 1 2 ...,±,±,=( ) ωk ωc k ωT⋅+= hr102 Chapter 4. Dynamics of Milling Operations, (4.76) where is the eigenvalue, p is the eigenvector defined as: , (4.77) is a 3(2 +1)-by-3(2 +1) matrix called oriented transfer function: , (4.78) where . (4.79) Directional coefficients for the multi frequency case are obtained by carrying out the integration in Eq (4.64): I λ Φ ωc ωT,( )⋅–( ) p⋅ 0= λ 0.5 a Ktc 1 e iωcτ––( )⋅ ⋅ ⋅= p p0 p1 p 1– … phr p hr–⎩ ⎭⎨ ⎬ ⎧ ⎫T = Φ hr hr Φ ωc ωT,( ) k = 0 1 1– . . hr hr–⎝ ⎠⎛ ⎞ Go 0 0,( ) Go 1 1,–( ) Go 1 1–,( ) . . . . Go 1 0,( ) Go 0 1,( ) Go 2 1–,( ) . . . . Go 1 0,–( ) Go 2– 1,( ) Go 0 1–,( ) . . . . . . . . . . . . . . . . . . . . . . . . Go p k–( ) k,( ) ⎭⎪ ⎪⎪ ⎪⎪ ⎬⎪ ⎪⎪ ⎪⎪ ⎫ p = 0 1 1– . . hr hr– = ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ Go p k–( ) k,( ) Ap k– G iωk( )⋅=103 Chapter 4. Dynamics of Milling Operations, (4.80) where matrix elements are , (4.81) , (4.82) , (4.83) Ar N 2π------ αxxr( ) αxyr( ) αxzr( ) αyxr( ) αyyr( ) αyzr( ) αzxr( ) αzyr( ) αzzr( ) = αxxr( ) c1 θ( ) c0 θ( ) c2 θ( )–( ) Kr κsin⋅ Ka κcos⋅+( )––[ ] φst φex= αxyr( ) c0 θ( ) c2 θ( ) c1 θ( ) Kr κsin⋅ Ka κcos⋅+( )–––[ ] φst φex= αxzr( ) 2 c4 θ( ) κtan⁄ c3 θ( ) κ Kr Ka κtan⁄+( )⋅cos⋅+[ ] φst φex= ⎭⎪ ⎪⎪ ⎬⎪ ⎪⎪ ⎫ αyxr( ) c0 θ( ) c2 θ( ) c1 θ( ) Kr κsin⋅ Ka κcos⋅+( )⋅––[ ] φst φex= αyyr( ) c1 θ( ) c0 θ( ) c2 θ( )+( ) Kr κsin⋅ Ka κcos⋅+( )⋅( )–[ ] φst φex= αyzr( ) 2 c3– θ( ) κtan⁄ c4 θ( ) κ Kr Ka κtan⁄+( )⋅cos⋅+[ ] φst φex= ⎭⎪ ⎪⎪ ⎬⎪ ⎪⎪ ⎫ αzxr( ) 2 c3 θ( ) Kr κcos Ka– κsin⋅ ⋅( )[ ] φst φex= αzyr( ) 2 c4 θ( ) Kr κcos Ka– κsin⋅ ⋅( )[ ] φst φex= αzzr( ) 2 c0 θ( ) κcos Kr – / κtan Ka+( )⋅ ⋅[ ] φst φex= ⎭⎪ ⎪⎪ ⎬⎪ ⎪⎪ ⎫104 Chapter 4. Dynamics of Milling Operations. (4.84) Note that some coefficients in Eq (4.84) are undefined when rN = , , therefore, undefined elements of the directional coefficient matrix, , in Eq (4.80) are separately calculated as • when rN = 2, (r>0): , (4.85) • when rN = 1, (r>0): c0 θ( ) θ if r 0= i– rN( )⁄ e irNθ–⋅ otherwise⎩⎪ ⎨⎪ ⎧ = c1 θ( ) 2 2θ irN 2θsin+cos( ) r2N2 4–( )⁄( ) e irNθ–⋅= c2 θ( ) 2– 2θ irN 2θcos+sin( ) r2N2 4–( )⁄( ) e irNθ–⋅= c3 θ( ) θ irN θsin+cos( ) r2N2 1–( )⁄( ) e irNθ–⋅= c4 θ( ) θ irN θcos+sin–( ) r2N2 1–( )⁄( ) e irNθ–⋅= ⎭⎪ ⎪⎪ ⎪⎪ ⎪⎪ ⎬⎪ ⎪⎪ ⎪⎪ ⎪⎪ ⎫ 1± 2± Ar αxxr( ) 18-- i Kr κsin⋅ Ka κcos⋅+( ) e i4θ– 4e i2θ–– i4θ– 1+( ) e i4θ– i4θ 1–+( )+[ ] φst φex = αxyr( ) 18-- Kr κsin⋅ Ka κcos⋅+( ) e i4θ– i4θ 1–+( ) ie– i4θ– i– 4e i2θ– 4θ– i+( )+[ ] φst φex = αyxr( ) 18-- Kr κsin⋅ Ka κcos⋅+( ) e i4θ– i4θ 1–+( ) ie– i4θ– i4e i2θ– 4θ– i+ +( )+[ ] φst φex = αyyr( ) 18-- i Kr κsin⋅ Ka κcos⋅+( ) e i4θ–– 4e i2θ–– i4θ 1+ +( ) e i4θ–– i4θ 1+–( )+[ ] φst φex = ⎭⎪ ⎪⎪ ⎪⎪ ⎪⎪ ⎬⎪ ⎪⎪ ⎪⎪ ⎪⎪ ⎫105 Chapter 4. Dynamics of Milling Operations. (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: , (4.87) where is the vibration frequency, and is the tooth passing frequency. Eq (4.87) describes the characteristic equation of the closed loop milling system. When the oriented transfer function is known, all eigenvalues can be determined, from which critical depths of cut are obtained as: , (4.88) where , , i is the imaginary unit and subscripts R and I represent real and imaginary parts, respectively. Substituting = into Eq (4.88), real and imaginary parts of depth of cut are expressed as: αxzr( ) cosκ 2⁄–( ) e i2θ– i2θ 1–+( ) Kr Ka κtan⁄+( ) i κsin⁄( ) e i2θ– i– 2θ 1–( )–[ ] φst φex = αyzr( ) cosκ 2⁄( ) i e i2θ– i– 2θ( ) Kr Ka κtan⁄+( ) 1 κsin⁄( ) e i2θ– i2θ 1–+( )+[ ] φst φex = αzxr( ) Kr κcos Ka– κsin⋅ ⋅( ) iθ θsin2 0.5 2θsin––[ ]⋅–{ } φst φex= αzyr( ) Kr κcos Ka– κsin⋅ ⋅( ) i θsin2 0.5 2θsin– θ–[ ]⋅–{ } φst φex= ⎭⎪ ⎪⎪ ⎪⎪ ⎬⎪ ⎪⎪ ⎪⎪ ⎫ det I λ Φ ωc ωT,( )⋅–( ) 0= ωc ωT Φ ωc ωT,( ) alim q, 2λq Ktc 1 e iωcT––( )⋅ --------------------------------------- aR q, i aI q,⋅+= = λq λR q, i λI q,⋅+= q 1…3 2hr 1+( )= e iωcτ– ωcτ( ) i ωcτ( )sin⋅–cos106 Chapter 4. Dynamics of Milling Operations, (4.89) . (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: . (4.91) Substituting Eq (4.91) into Eq (4.89), real depth of cut is solved as: . (4.92) The root condition in Eq (4.91) can be expressed alternatively as: . (4.93) Dividing both sides by the magnitude of the eigenvalue: , (4.94) where , the root condition for each eigenvalue can be simplified as: aR q, 1 Ktc ------- λR q, 1 ωcτ( )cos–( )⋅ λI q, ωcτ( )sin+ 1 ωcτ( )cos– ------------------------------------------------------------------------------------------= aI q, 1 Ktc ------- λR q, 1 ωcτ( )cos–( ) λI q,– ωcτ( )sin⋅ 1 ωcτ( )cos– ---------------------------------------------------------------------------------------= aI q, 0= ⇒ λI q, λR q, ---------- ωcτ( )sin 1 ωcτ( )cos– ------------------------------- ωcτ 2⁄( )cos ωcτ 2⁄( )sin -----------------------------= = aR q, λR q, Ktc ---------- 1 λI q, λR q, ----------⎝ ⎠⎛ ⎞ 2 += λR q, ωcτ 2⁄( )cos⋅ λI q, ωcτ 2⁄( )sin⋅– 0= λR q, ωcτ 2⁄( )cos⋅ λI q, ωcτ 2⁄( )sin⋅– λq ---------------------------------------------------------------------------------------------- 0= λq λR q,2 λI q,2+=107 Chapter 4. Dynamics of Milling Operations, (4.95) where , (4.96) is the phase shift of the eigenvalue. Finally, the general root function can be determined by multi- plying individual root conditions for each eigenvalue: , (4.97) where 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 con- structed iteratively. Note that, in the most general case, eigenvalues are functions of both excita- tion and tooth passing frequencies; therefore, . 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 char- acteristic equation in Eq (4.87) reduces to: , (4.98) ωcτ 2⁄ ϕq+( )cos 0= ϕq λI q, λR q, ----------⎝ ⎠⎛ ⎞atan= f ωc( ) ωcτ 2⁄ ϕq+( )cos q 1= 3 2hr 1+( ) ∏ 0= = hr ϕq ϕq ωc ωT,( )= det I λ Go0 0,( )⋅–( ) det I λ A0 G iωc( )⋅ ⋅–( ) 0= =108 Chapter 4. Dynamics of Milling Operationsand the overall transfer function matrix, , becomes independent of spindle speed or . This simplifies the problem as eigenvalues of Eq (4.98) can be solved by just specifying a chatter frequency, . 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: , (4.99) where represents "lobe" numbers. From Eq (4.99): , (4.100) where is the phase angle between current and previous vibration marks. Since eigen- values hence phase angle can be calculated for a given , a set of spindle speeds, n [rev/min], is obtained from solving for tooth period as: . (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 G iωc( ) ωT ωc ωcτ 2⁄ ϕq+ π2-- 2k 1+( )= k 0 1 2…, ,= ωcτ 2πk ε+= ε π 2ϕq–= ωc τ n 60ωc N 2πk ε+( )---------------------------=109 Chapter 4. Dynamics of Milling Operationstool 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 flexi- bilities are defined for stability consideration. When the stability equation in Eq (4.98) is ana- lyzed, 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 coeffi- cient 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 coeffi- cient from the depth of cut and considering the influence of each geometry with a different direc- tional coefficient matrix, stability for a general end mill can now be described as: , (4.102) where is the delay term containing the cutting force coefficient , is the relative depth of cut, and is the directional coefficient matrix of each section (i=1,2,3). det I E1 a1 A0 1,⋅ ⋅ E2 a⋅ 2 A0 2,⋅ E3 a⋅ 3 A0 3,⋅+ +( )– G iωc( )⋅( ) 0= Ei 0.5 1 e iωcτ––( ).Ktc i,⋅= Ktc i, ai A0 i,110 Chapter 4. Dynamics of Milling OperationsThe 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, , having the following characteristic equation: , (4.103) which can alternatively be represented as a standard eigenvalue problem as: , , (4.104) where . 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 differ- ence is that only depths of cut within the range of characteristic equation, i.e. , are accepted to construct stability lobes. For out of range values, i.e. when the critical depth of cut M O S S z Mz N z r @ # ! Nza1 a2 a3 M O S N z r @ # ! a1 M O S N z r @ # ! a1=Mz a2 (a) (b) (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. a1 det I E1 a1 A0 1,⋅ ⋅ – G iωc( )⋅( ) 0= det I λ1– A0 1, G iωc( )⋅ ⋅( ) 0= 0 a1 Mz< <( ) λ1 E1 a1⋅= 0 a1 Mz< <111 Chapter 4. Dynamics of Milling Operationsexceeds 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: . (4.105) The above equation can be represented in an equivalent general eigenvalue problem as: , (4.106) where and . There are three unknowns in this equation: the tooth period ( ) or the spindle speed, the depth of cut ( ) and the excitation frequency ( ). 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 ear- lier, 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 ( ) and largest ( ) natural frequencies of the structure. This band is then divided into a finite number of elements at which eigenvalues and root function, , 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, , are calculated at chatter frequencies, the smallest depth of cut is selected and a range check, i.e. , 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 det I E1 M⋅ z A0 1,⋅ E2 a2 A0 2,⋅ ⋅+( ) – G iωc( )( ) 0= det C ωc τ,( ) λ2– A0 2, G iωc( )⋅ ⋅( ) 0= C ωc τ,( ) I E1 M⋅ z– A0 1, G iωc( )⋅ ⋅= λ2 E2 a2⋅= τ a2 ωc ωstart ωend f ω( ) τ a2 0 a2 Nz Mz–( )< <112 Chapter 4. Dynamics of Milling Operationsmill. Although GZO solution is an iterative approach, it provides a better estimate for stability as characteristic equation is constructed more realistically. 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 depend- ing 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 spec- trum 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 prob- lem (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 ), and searching for positive-real depths of cut by iterating chatter frequency ( ). Such analysis for a single mode can be found in startω endω Approximate solution sought in between General Search Sign of the root function Fine Tuning + + _ +/ _ _ _ ΔωFreq. Step Size - Eigenvalue Calculation Method StandardLarge Small Standard or Approximate Figure 4.7 : Graphical representation of iterative root search. ωT ωc113 Chapter 4. Dynamics of Milling Operationsthe 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 , the first component of eigen- vector, , 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 mathe- matically possible but physically impossible to exist, therefore, it needs to be eliminated from sta- bility 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 bound- ary. ωc p0114 Chapter 4. Dynamics of Milling Operations4.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 vicin- ity 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 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 is the total number of harmonics of the tooth passing frequency taken into account, 3(2 +1) complex eigenvalues are obtained from Eq (4.87), which means that at each iteration step corre- sponding to a very small change in chatter frequency estimate, a significant number of eigenval- ues 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 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 Depth of cut Stability state: Stable/Unstable Variable Scanning 3 ZO Chatter frequency Spindle speeds Depth of cut Fixed Analytical 1 GZO Spindle speed Chatter frequency Depth of cut Variable Iterative 2 MF Spindle speed Chatter frequency Depth of cut Fixed Iterative 3 ωc ωc ωc 10 3– hr hr115 Chapter 4. Dynamics of Milling Operationsal. [83] and Lou et al. [85] for eigenvalue re-analysis are applied to increase computational effi- ciency of the solution of the machine tool stability problem. By eliminating the need for recalcu- lation 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: , (4.107) where ’s are the eigenvalues, ’s and ’s are the associated right and left eigenvectors, respectively, , and "*" is the complex conjugate transpose operator. Suppose that the original eigenvalue problem is represented when A= , B= , , , and . For a finite change of matrices A and B, the modified eigenvalue problem becomes: , (4.108) where = + , = + , = + , = + , = + . 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]: A xi λi B xi= yi * A yi * λi B= ⎭⎪ ⎬⎪ ⎫ λi xi yi i 1 … n, ,= A 0( ) B 0( ) xi xi 0( )= yi yi 0( )= λi λ0i= A 1( )xi 1( ) λ1i B 1( )xi 1( )= yi 1( )*A 1( ) yi 1( )*λ1i B 1( )= ⎭⎪ ⎬⎪ ⎫ A 1( ) A 0( ) ΔA B 1( ) B 0( ) ΔB xi 1( ) xi 0( ) Δxi yi 1( ) yi0( ) Δyi λ1i λ0i Δλi Δ116 Chapter 4. Dynamics of Milling Operations. (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 and can be expressed as complex linear combinations of all eigenvectors, and [85]: . (4.110) Note that, in ( = + ), the effect of vector can be combined into the first term because the mode vector only describes the proportional relation among its elements [85], hence the term is excluded in Eq (4.110). This statement is also valid for the left eigenvector. When Eq (4.109) is multiplied by and the orthogonality condition of the eigenvectors, being , (4.111) is applied, the variation in eigenvalues can be easily calculated as: . (4.112) Finally, eigenvalue estimate for the modified system becomes: A 0( ) λ0i B 0( )–( ) Δxi⋅ ΔA λ0i ΔB–( ) xi 0( )⋅+ B 1( )xi 0( ) B 0( )Δxi+( ) Δλi⋅= Δxi Δyi xi 0( ) yi 0( ) Δxi gjixj 0( ) j 1 i≠( )= n ∑= Δyi hji* yj 0( ) j 1 i≠( )= n ∑= ⎭⎪ ⎪⎪ ⎬⎪ ⎪⎪ ⎫ xi 1( ) xi 0( ) Δxi ith xi 0( ) xi 1( ) ith yi 0( )* yi * A xj 0= , yi * B xj 0 i j≠,= Δλi yi 0( )* ΔA λ0i ΔB–( ) xi0( ) yi 0( )*B 1( )xi 0( )------------------------------------------------------------- yi 0( )* A 1( ) λ0i B 1( )–( ) xi0( ) yi 0( )*B 1( )xi 0( )--------------------------------------------------------------- yi 0( )*A 1( )xi 0( ) yi 0( )*B 1( )xi 0( )------------------------------ λ0i–= = =117 Chapter 4. Dynamics of Milling Operations. (4.113) This is just a specific form of the Rayleigh quotient for the eigenvalue problem in Eq (4.107). If Eq (4.110) is substituted into Eq (4.109) and the result is multiplied by ( ), the follow- ing form is obtained: , = . (4.114) Using the orthogonality condition in Eq (4.111), the above equation is simplified as follows: ( ), ( ). (4.115) Coefficients for the left eigenvalue problem in Eq (4.108) are similarly obtained as: , ( ), ( ), (4.116) and finally eigenvectors of the perturbed system are easily found by: λ1i yi 0( )*A 1( )xi 0( ) yi 0( )*B 1( )xi 0( )------------------------------= yr 0( )* r i≠ yj 0( )* A 0( ) λ0iB 0( )– B 0( )Δλi–( ) x00( )g0i x10( )g1i …+ +( ) yj 0( )* B 1( ) Δλi ΔA– λ0i ΔB+( ) xi0( ) yj 0( )* A 0( ) λ0iB 0( )– B 0( )Δλi–( ) xj 0( )gji yj0( )* B 1( )Δλi ΔA– λ0i ΔB+( ) xi 0( )= yj 0( )* A 0( ) λ1i B 0( )–( ) xj 0( )gji yj0( )* B 1( )Δλi A 1( )– λ0i B 1( )+( ) xi0( )= yj 0( )* λ0j λ1i–( )B 0( )xj 0( )gji yj 0( )* λ1i B 1( ) A 1( )–( ) xi 0( )= gji yj 0( )* A 1( ) λ1i B 1( )–( ) xi0( ) yj 0( )*B 0( )xj 0( ) λ1i λ0j–( ) ---------------------------------------------------------------= j 1 … n, ,= j i≠ hji yi 0( )* A 1( ) λ1i B 1( )–( ) xj0( ) yj 0( )*B 0( )xj 0( ) λ1i λ0j–( ) ---------------------------------------------------------------= j 1 … n, ,= j i≠118 Chapter 4. Dynamics of Milling Operations. (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 system: , , and and the eigenvalue solution of the modified system is estimated using Eqs (4.113) and (4.117): , , and . Convergence at each step is checked by comparing the rate of change of eigenvalue estimate: where tol is user-defined toler- ance. If the iteration is not successful, a new set of eigenvalues are estimated by replacing the original eigenvalues and eigenvectors with the new estimates, i.e. = , = and = . Although the perturbation method described above involves operation count of , 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: . (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; how- ever, 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 xi 1( ) xi 0( ) Δxi+= yi 1( ) yi 0( ) Δyi+= ⎭⎪ ⎬⎪ ⎫ λ0i xi 0( ) yi 0( ) λ1i xi1( ) yi 1( ) λ1i λ0i–( ) λ1i⁄ tol≤ λ0i λ1i xi0( ) xi1( ) yi 0( ) yi 1( ) n3 A I, = B Φ ω ωT,( )=119 Chapter 4. Dynamics of Milling OperationsEq (4.87). An appropriate sorting criterion to obtain consistent indices is to check the relative dis- placement of eigenvalues. Eigenvalues obtained at the starting frequency, , are taken as the reference for indexing, and when the frequency is varied slightly, the distance of each new eigen- value 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 index- ing eigenvalues and how eigenvalues are successfully traced. Eigenvalue estimation is incorporated into the stability solution to increase computational effi- ciency. 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. 4.4. Comparison of Stability Methods In an effort to compare various stability methods presented above, the stability of a milling opera- tion is investigated. Dynamic parameters of a 76.2 [mm] face milling cutter with four 25.4 [mm] ωstart λ5 λ6 λ2 λ1λ1 λ3 λ4 λR λR ω0ω1 ω0ω1 λ3 λ6 λ2 λ5 λ4 λI λI λ5 λ6 λ2 λ1λ1 λ3 λ4λ2 λ4 λ6 λ5 λ3 (a) (b) Figure 4.8 : Shift of eigenvalues when the frequency is changed from to : (a) unsorted eigenvalues; (b) sorted eigenvalues. ω0ω1120 Chapter 4. Dynamics of Milling Operationscircular 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 ( =0 and = ) of a P20 steel with the following constant tangential, radial and axial cutting force coefficients, respec- tively: =3146 [N/ ], =1680 [N/ ], =419.4 [N/ ]. 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. 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 ele- ments. 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 = . Number of discrete time elements for the SD is varied depending on the spindle speed to accom- modate 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 accu- rately. Table 4.2 : Dynamic parameters [101]. Direction Natural Frequency [Hz] Stiffness [N/m] Damping [-] X 28 2.54 e7 0.17 55 4.30 e7 0.12 Y 28 2.15 e8 0.10 55 6.18 e8 0.06 φst φex π 2⁄ Ktc mm 2 Krc mm 2 Kac mm 2 κ π 4⁄121 Chapter 4. Dynamics of Milling OperationsA 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 =31.7 Hz, is found to be =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 har- monics do not influence the stability [75]. In this case-study, although radial immersion is consid- erably high, it is found out that harmonics play an important role around a specific spindle speed region. For example at 475 rpm, ( - )/ =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 ( =1) and high ( =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 chat- ter vibration frequency, i.e., ωT 2π⁄ ωc 2π⁄ ωc ωT 2π hr hr122 Chapter 4. Dynamics of Milling Operations. (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 eigen- vector 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 fre- quency, 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 captur- ing 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. p 1{ }T p1p0 ------- p 1– p0 ---------- … phrp0 -------- p hr– p0 ----------- ⎩ ⎭⎪ ⎪ ⎨ ⎬⎪ ⎪ ⎧ ⎫T =123 Chapter 4. Dynamics of Milling Operations200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 1100 1150 4 6 8 10 12 14 16 18 Spindle Speed [rpm] D e p th o f C u t [m m ] GZO ZO SD MF-5 MF-1 TD 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 ( =1), and MF-5: Multi frequency solution with five harmonics ( =5).hr hr 200 250 300 350 400 450 500 550 600 55 60 65 70 75 80 Spindle Speed [rpm] C h a tt e r F re q u e n c y [ H z ] ZO MF-5 MF-1 Figure 4.10 : Chatter frequencies obtained from frequency domain based solutions.124 Chapter 4. Dynamics of Milling OperationsThe 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. -5 -4 -3 -2 -1 0 1 2 3 4 5 0 0.2 0.4 0.6 0.8 1 ωc (ωc+ωT) (ωc−ωT) ωc (ωc+ωT) (ωc−ωT) -5 -4 -3 -2 -1 0 1 2 3 4 5 N o rm a li z e d E ig e n v e c to r M a g n it u e s Figure 4.11 : Normalized X-direction eigenvector magnitudes of MF-5 at (a) n=475 [rpm]; / =31.7 [Hz], / =60.2 [Hz]; (b) n=310 [rpm], / =20.7 [Hz], / =58.5 [Hz]. ωT 2π ωc 2π ωT 2π ωc 2π (a) (b)125 Chapter 4. Dynamics of Milling OperationsThe performance of the eigenvalue estimation technique is also tested using the same milling sys- tem. 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 / =14 [Hz] and / =100 [Hz]. The same root function is approximated at sample points using eigen- value 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 eigen- values to 50 improved estimation significantly as seen in Figure 4.13.b. In summary, approximat- ing the eigenvalue problem with 50 exact points requires less computational load while maintaining acceptable accuracy. -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 Im a g in a ry Im a g in a ry Real Real Unit Circle Unit Circle 0. 94 1.0 0 Figure 4.12 : Eigenvalues of the transition matrix - Eq (4.58) - of the semi- discretization method (k=41) when cutting conditions are selected on stability border at (a) 475 [rpm], 12.8 [mm]; (b) 310 [rpm], 7.6 [mm]. (a) (b) ωstart 2π ωend 2π126 Chapter 4. Dynamics of Milling Operations4.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 10 20 30 40 50 60 70 80 90 10010 20 30 40 50 60 70 80 90 100 -0.5 0 0.5 1 Exact Approximate Exact Approximate R o o t F u n c ti o n , f (ω ) [- ] Frequency [Hz]Frequency [Hz] Figure 4.13 : Comparison of exact and approximate root functions of one harmonic case using (a) 10 ( / =8.6 [Hz]); (b) 50 ( / =1.7 [Hz]) exact frequency elements. 1000 ( / =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]. Δω 2π Δω 2π Δω 2π (a) (b)127 Chapter 4. Dynamics of Milling Operations10% 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 environ- ment 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 cut- ting 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 Optimization5.2. Process Optimization Manufacturing of a part follows a well-defined methodology (Figure 5.1) that starts with the con- ceptual 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 differ- ence 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 con- ducted 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 OptimizationThe optimization problem is separated into two independent parts. The first part is called Pre-Pro- cess 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 manu- facturing 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: , (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 spin- dle 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: , (5.2) Design Engineer Concept Solid model Process Engineer NC Code generation (process parameters) Shop Floor Engineer Application of NC code on real machine Determination of machining failures Failure Su cc es s Production Experiments Optimum Production PROPOSED Virtual Process Simulation & Optimization Module CURRENT CAD STAGE CAM STAGE Figure 5.1 : Manufacturing flow chart, current and proposed methodologies. MRR a B c n N⋅ ⋅ ⋅ ⋅= qlb q qub≤ ≤131 Chapter 5. Process Optimizationwhere q is a design variable limited by lower bound and upper bound . Nonlinear con- straints; on the other hand, are rather expressed in a function form such as: , (5.3) where is a nonlinear constraint that is expressed as a function of 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 and , respec- tively. 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 qlb qub fnc a B n N, , ,( ) 0≤ fnc c B a n Figure 5.2 : Design variables. nmin n nmax≤ ≤ cmin c cmax≤ ≤132 Chapter 5. Process Optimizationoutputs 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 Optimiza- tion). 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 vari- ables, 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 follow- ing form: . (5.4)fobj a B c⋅ ⋅=133 Chapter 5. Process Optimization5.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. 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 and 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. CAD Model Cutting Conditions Depth of Cut (a) Width of Cut (B) Feed Rate (c) Spindle Speed (n) . . . Facing Pocketing Contouring Plunging Ramping . . . Milling Cutters / Operations Pre-Process Optimization Constraints: Machining Machine Tool CAM Tool Path Generator Optimized NC Code Process Planning (CAM) Figure 5.3 : Optimized NC code generation flow chart with the Pre-Process Optimization module. nlobe 1, nlobe 2,134 Chapter 5. Process OptimizationFor convenience, instead of actual width of cut, normalized width of cut, b, is preferred as it var- ies between 0 and 1 corresponding to no and full immersion cases, respectively: , (5.5) where D is the cutter diameter, and are the cutter entry and exit angles, respectively. 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 loca- tions 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 sec- tion. b B D⁄ 1 φexcos–( ) 2 for up milling⁄ 1 φstcos+( ) 2 for down milling⁄⎩ ⎭⎪ ⎪ ⎨ ⎬⎪ ⎪ ⎧ ⎫ = = φst φex B a n Vibration Marks Figure 5.4 : Design variables that are important to suppress chatter vibrations.135 Chapter 5. Process Optimizationn - Spindle Speed a - D ep th o f C ut Stable Region nlobe,1nlobe,2 B = const n - Spindle Speed b - W id th o f C ut nlobe,1n lobe,2 a = const Stable Region 0 1.0 Figure 5.5 : Chatter stability lobes for constant (a) radial width of cut, (b) axial depth of cut. (a) (b) Feasible Region (safe to operate) a - Depth of Cut b - W id th o f C ut n = const 0 1.0 Figure 5.6 : Design space bounded by chatter stability.136 Chapter 5. Process Optimization5.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, , and vary as the cutter rotates: [Nm], (5.6) [hp], (5.7) where D is the cutter diameter in [m], n is the spindle speed in [rpm], and is in [N]. In this sec- tion, 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 con- structed for torque/power constraint. Consider a Cutter Engagement Feature (CEF) for a cylindrical end mill composed of depth of cut "a" and immersion angles " " and " " 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 and , respectively and they can be expressed as: , (5.8) , (5.9) where and is the helix angle. Two different cases occur when and , and they are named as "Case I" (Figure 5.7.a) and "Case II" (Figure 5.7.b), respectively. Ft φ( ) Torque φ( ) Ft φ( ) D 2⁄⋅= Power φ( ) Ft φ( ) D.π.n60--------------⎝ ⎠⎛ ⎞ 1.34102 10 3–×⋅ ⋅= Ft φst φex φcr 0, φcr 1, φcr 0, φst ki0a+= φcr 1, φex ki0a+= ki0 2 i0tan( ) D⁄= i0 φcr 0, φex< φcr 0, φex≥137 Chapter 5. Process OptimizationFor a cylindrical end mill, tangential cutting force for flute j can be obtained by integrating differ- ential cutting forces along: , (5.10) where is the angular position of differential element dz, is the current angular position of flute j ( ), N is the number of teeth, is the pitch angle and defined as for a uniform pitch cutter. Integration can be defined between alter- native boundaries using conversion: , (5.11) which can be further simplified as: , (5.12) where lower ( ) and upper ( ) integration boundaries are the angles at which flute enters and exits the cutter workpiece engagement area (see Figure 5.8). φst φex φ0 z φcr,0 φcr,1 φst φex φ0 z φcr,0 φcr,1 Immersion Angle Immersion Angle D ep th o f C ut D ep th o f C ut Flute Flute a a Figure 5.7 : Important parameters for analytical solution: (a) Case I ( ), (b) Case II ( ).φcr 0, φex< φcr 0, φex≥ (a) (b) Fj t φj( ) Kte Ktc.c. φsin+( ) zd zj 0, zj 1,∫= φ φj ki0z–= φj φ j 1–( )φp+= j 1 … N, ,= φp φp 2π N⁄= dφ ki0– dz= Ftj φj( ) 1–ki0 ----- Kte Ktc.c. φsin+( ) φdφj 0, φj 1,∫= Ftj φj( ) 1ki0 ----- Ktc c φj 1,cos φj 0,cos–( )⋅ ⋅ Kte φj 0, φj 1,–( )+{ }= φj 0, φj 1,138 Chapter 5. Process OptimizationTotal tangential cutting force acting on the tool can be obtained by summing individual contribu- tions of all flutes: . (5.13) Variation of lower ( ) and upper ( ) 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. φj,1 φj,0 Depth of Cut (a) Immersion Angle,( )φ x y y x Flute j Figure 5.8 : Cutter workpiece engagement, integration boundaries. Ft Ftj j 1= N ∑= φj 0, φj 1,139 Chapter 5. Process OptimizationTable 5.1 : Integration boundaries for different regions. where , , , and , (5.14) 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- Region 1 Region 2 Region 3 φj,1 φj,0 φst φj - kio*a φj φj φex CASE I Region 1 Region 2 Region 3 φj,1 φj,0 φst φj - kio*a φj φex CASE II φst φex φj - kio*a Region 1: φst φj φb1≤ ≤( ) Region 2: φb1 φj< φb2≤( ) Region 3: φb2 φj< φcr 1,≤( ) φb1 φcr 0, Case I φex Case II⎩ ⎭⎪ ⎪ ⎨ ⎬⎪ ⎪ ⎧ ⎫ = φb2 φex Case I φcr 0, Case II⎩ ⎭⎪ ⎪ ⎨ ⎬⎪ ⎪ ⎧ ⎫ = φImmersion Angle ( ) In te gr at io n B ou nd ar ie s ( ) φst φ 0 ,φ 1 φexφcr,0 φcr,1 φ0φ1 φImmersion Angle ( ) φ0φ1 1 φst φ ex φ cr ,0 φcr,1 2 3 CASE I CASE II 1 2 3FL U TE " j" φj φj,0 φj,1 φst φex φst φex In te gr at io n B ou nd ar ie s ( ) φ 0 ,φ 1 Figure 5.9 : (a) Case I ( ), (b) Case II ( ).φcr 0, φex< φcr 0, φex≥ (a) (b)140 Chapter 5. Process Optimizationaries on the immersion angle, sum/difference formulae for sine and cosine are utilized, and ana- lytic tangential force is reorganized as: , (5.15) 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: , (5.16) where . Table 5.2 : Generalized coefficients for tangential force. 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 Ft Ktc c⋅ ki0⁄( ) Cc j, φcos⋅ Cs j, φsin⋅+ Cp j, φ⋅ const+ +( ) j 1= N ∑= φ dFt dφ------- Ktc c⋅ ki0⁄( ) Cc– φsin⋅ Cs φcos⋅ Cp+ +( ) 0= = Cq Cq j, j ∑ q = c s p, ,( ),= Region 1 Region 2Region 3 Cc,j Cs,j - cos((j-1)φp) CASE I cos((j-1)φp- kio*a) - cos((j-1)φp) cos((j-1)φp- kio*a) - sin((j-1)φp- kio*a) + sin((j-1)φp)- sin((j-1)φp- kio*a) sin((j-1)φp) CASE II 0 0 CASE I/II Cp,j Kte /(Ktc.c) - Kte /(Ktc.c) 0 φst φj φb1≤ ≤ φb1 φj< φb2≤φb2 φj< φcr 1,≤141 Chapter 5. Process Optimizationthat 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). Once all possible scenarios are identified, critical angular positions of the cutter at which mini- mum and maximum forces appear can be calculated by solving Eq (5.16). Using the trigonometric identity, , Eq (5.16) can be reduced into quadratic equation and two roots, if they exist, are calculated as: , . (5.17) In addition to roots obtained above, angular positions , , , and , 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 position, and , respectively. φ In te gr at io n B ou nd ar ie s 1 φ0,φ1 2 3 j = 1 φj=2 =φ+φpφj=1 =φ(t) φp shift of flutes with time 2 3 41 2 3 4 Figure 5.10 : Possible intersection scenarios. φsin 1 φcos 2–= φroots C– sCp Cc Cc 2 Cs 2 Cp 2–+± Cc 2 Cs 2+ ---------------------------------------------------------------- ⎝ ⎠⎜ ⎟ ⎜ ⎟⎛ ⎞acos= 0 φroots 2π≤ ≤( ) φcr 0, φcr 1, φst φex φmint φmaxt142 Chapter 5. Process OptimizationOnce tangential forces are known, torque and power demands denoted by and are easily calculated. In order to prevent violation of torque/power constraints, cumu- lative average torque and power shall not exceed machine tool torque/power limitations at spindle speed n: , (5.18) and are machine tool torque and power curves, respectively, and they are provided by the machine tool manufacturer (see Figure 5.11). Subscript "cav" indicates cumulative average which is mathematically defined as: , (5.19) where q={Torque,Power}. Standard piecewise logarithmic relation of motor torque/power with spindle speed can be expressed as: , (5.20) Torque φ( ) Power φ( ) Torquecav Tmo n( )≤ Powercav Pmo n( )≤ ⎭⎬ ⎫ Tmo Pmo n2 T1 T2 Scale n1 Figure 5.11 : Typical torque-power characteristics of a machine tool. qcav q φmint( ) 2 q φmaxt( ) 2 +[ ] 2⁄= qmo n( ) q2 e q2 q1⁄( )log– n2 n⁄( )log n2 n1⁄( )log-------------------------------------------------------⋅=143 Chapter 5. Process Optimizationwhere q={T,P}, and {n: } 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 opera- tions 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. Earlier in Eq (3.27), the chip thickness for 3-axis milling was defined as: , (5.21) where c is the feed per tooth, and 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 n1 n n2< < hex a R κ z c x Figure 5.12 : Chip thickness for a face mill with a round insert. hex c. fxy κ φsin⋅sin⋅ fz κcos⋅–( )= fxy fz κ φ φ κ144 Chapter 5. Process Optimizationvalues denoted by and , respectively. There are two distinct cases for 3-axis milling operations: (a) Ramping Up or Planar Milling : Only the front of the cutter ( ) can engage in a cut. In this case, maximum chip thickness is obtained when , , or at values closest to these. This rule can be generalized as follows: and axial immersion angle at if . (b) Ramping Down : Both front and back of the cutter can engage in a cut. At the front of the cutter ( ), maximum chip thickness is reached when , and the axial immersion angle is found as follows: , . (5.22) At the back of the cutter ( ); on the other hand, maximum chip thickness is reached when or , and , or at values closest to these. If the desired maximum chip load is denoted by and the necessary feedrate to achieve this maximum chip load is denoted by , feedrate of the tool is updated by the following relation: , (5.23) Down Milling ( ) Up Milling ( ) φmaxchip κmaxchip fxy 0≥ fz 0≥,( ) 0 φ π≤ ≤ φmaxchip π 2⁄= κmaxchip π 2⁄= φex π= φst 0= If (φst π 2⁄ ), φmaxchip=φst> If (φex π 2⁄ ), φmaxchip=φex< If (φst π 2⁄ ), φmaxchip=π 2⁄≤ If (φex π 2⁄ ), φmaxchip=π 2⁄≥ κmaxchip κ z a=( )= a R<( ) fxy 0≥ fz 0<,( ) 0 φ π≤ ≤ φmaxchip π 2⁄= dhex dκ---------- c fxy κcos φmax chipsin⋅ ⋅ fz κsin⋅+( )⋅ 0= = κmaxchip fxy φmaxchipsin⋅– fz ----------------------------------⎝ ⎠⎜ ⎟ ⎛ ⎞ atan= π φ 2π≤ ≤ φmaxchip π= φmaxchip 2π= κmaxchip 0= hmax cmax cmax hmax hfactor --------------=145 Chapter 5. Process Optimizationwhere is the chip thinning factor. Note that the chip thinning factor must be limited with an upper value in order to avoid extreme magnifica- tion 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 =0, =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: [N/ ] and [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 power characteristics. Cutting force coefficients used in this example are [N/ ] and [N/mm]. hfactor fxy κmaxchip( ) φmaxchip( )sinsin fz κmaxchip( )cos–= φst φex Ktc 796 = mm2 Kte 27.7= Ktc 600= mm 2 Kte 42=146 Chapter 5. Process Optimization0 20 40 60 80 100 120 140 160 180 5 10 15 20 25 30 (t) [deg]φ φ 0 ,φ 1 [d eg ] 5 10 15 20 25 30 φ 0 ,φ 1 [d eg ] 5 10 15 20 25 30 φ 0 ,φ 1 [d eg ] 5 10 15 20 25 30 φ 0 ,φ 1 [d eg ] a = 5 mm a = 10 mm a = 15 mm a = 30 mm 0 50 100 150 200 250 300 350 -400 -300 -200 -100 0 (t) [deg]φ F t [ N ] -400 -300 -200 -100 0 F t [ N ] -400 -300 -200 -100 0 F t [N ] -400 -300 -200 -100 0 F t [ N ] 0 5 10 15 20 25 30 35 40 0 100 200 300 400 500 600 700 |F t,m ax | [ N ] a [mm] - depth of cut (a) (b) φ0φ1 φ0φ1 φ0φ1 φ0φ1 Figure 5.13 : (a) Integration boundaries with various depths of cut, simulated tangential forces, (b) analytically obtained maximum tangential forces.147 Chapter 5. Process OptimizationSimilar 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 [Hz], [Hz]; stiffness [N/m], damping ratio . 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 P ow er [h p] To rq ue [N -m ] Spindle Speed [rpm] 100 1500 3500 5000 13000 95.5 35.3 11 1.34 20.12 24.12 Figure 5.14 : Torque - Power curves of the machine tool. ωnx 600= ωny 660= kx ky 5.6 106×= = ζx ζy 0.035= =148 Chapter 5. Process Optimizationusing 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 fee- drate 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. Feasible Region (safe to operate) 2 4 6 8 10 12 0 0.2 0.4 0.6 0.8 1 a [mm] b [ -] Stability Limit Torque Limit Power Limit c = 0.25 mm/tooth.rev n = 12600 rpm D = 20 mm N = 3 1 Figure 5.15 : Radial and axial depths of cut curves for up- milling, feasible region (12600 [rpm]).149 Chapter 5. Process Optimization 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 2 4 6 8 10 12 14 16 18 20 x 103 0 1 2 3 4 n [rpm] a [ m m ] b = 1 (slotting) N = 3 Stable Region Figure 5.16 : Stability lobes for full immersion milling case. Feasible Region (safe to operate) 0 200 400 600 800 1000 1200 Stability Limit Torque Limit Power Limit c = 0.25 mm/tooth.rev n = 12600 rpm D = 20 mm N = 3 M R R [ cm 3 / m in ] Max. MRR 1 2 4 6 8 10 12 a [mm] Figure 5.17 : MRR curves of up-milling based on stability, torque and power limits.150 Chapter 5. Process Optimizationdemands are lower, process parameters are limited by the stability of the system, which can be seen clearly in Figure 5.18. 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 pre- vent 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. Feasible Region (safe to operate) 0 0.2 0.4 0.6 0.8 1 a [mm] b [-] Stability Limit Torque Limit Power Limit c = 0.25 mm/tooth.rev n = 6300 rpm D = 20 mm N = 3 0 2 4 6 8 10 12 14 16 18 20 Figure 5.18 : Radial and axial depths of cut curves for up- milling, feasible region (6300 [rpm]).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. 1000 1500 2000 2500 3000 3500 4000 4500 5000 0 0.5 1 1.5 2 2.5 x 1e-6 Frequency [Hz] M ag ni tu de [m /N ] Model Experimental X - Direction 1000 1500 2000 2500 3000 3500 4000 4500 5000 0 0.5 1 1.5 2 2.5 x 1e-6 Frequency [Hz] M ag ni tu de [m /N ] Model Experimental Y - Direction Figure 5.20 : Transfer functions of 10 [mm] endmill in x and y directions (HSK63A-Shrink Fit).152 Chapter 5. Process OptimizationBased 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 maxi- mize the production, the MRR variation as a function of immersion and spindle speed is plotted in Table 5.3 : Cutting force coefficients of Al 6061 Tool Type [N/mm] [N/mm] [ ] [ ] 10 [mm] 38.38 -9.792 403.396 98.911 25 [mm] 21.34 23.141 662.296 72.86 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0 0.5 1 1.5 2 2.5 x 1e-7 Frequency [Hz] M ag ni tu de [m /N ] Model Experimental X - Direction 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0 0.5 1 1.5 2 2.5 3 x 1e-7 Frequency [Hz] M ag ni tu de [m /N ] Model Experimental Y - Direction Figure 5.21 : Transfer functions of 25 [mm] indexable cutter in x and y directions (HSK63A- Corogrip). Kte Kre Ktc N mm⁄ 2 Krc N mm⁄ 2153 Chapter 5. Process OptimizationFigure 5.24, which justifies the argument that high productivity is achieved around pockets of sta- bility with the highest possible spindle speed. 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 x 1e4 0 0.5 1 1.5 2 2.5 n - Spindle Speed [rpm] a - D ep th o f C ut [m m ] Slotting 0 5 10 15 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 x 1e4n - Spindle Speed [rpm] a - D ep th o f C ut [m m ] Slotting 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. (a) (b)154 Chapter 5. Process Optimization0 2 4 6 8 10 12 14 16 0 0.5 110000 12000 14000 16000 18000 20000 22000 b [-] n [rpm] a [m m ] Figure 5.23 : Stability lobes for various immersion cases, 10 [mm] end mill. Cutting coefficients from Table 5.3 are used. 0 0.5 110000 12000 14000 16000 18000 20000 22000 20 40 60 80 100 120 140 160 b [-] n [rpm] M R R [c m 3 / m in ] 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.hmax155 Chapter 5. Process OptimizationProcess planning charts were constructed at selected spindle speeds based on analysis summa- rized 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 mate- rial 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 alter- native representations for the design space, start angle of cut versus depth of cut. Note that since down milling was used, exit angle was always equal to 180 [degrees].φex P ow er [H P ] To rq ue [N -m ] Spindle Speed [rpm] 200 1500 4500 10000 20000 95.4 30.0 8.8 2.48 20.12 24.81 17.6 11.67 Figure 5.25 : Mori Seiki SH403 - Torque/Power characteristics.156 Chapter 5. Process OptimizationFinally, 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 0 0.2 0.4 0.6 0.8 1 b 2 4 6 8 10 12 14 15 60 80 100 120 140 M R R [c m 3 / m in ] a [mm] 0 50 100 150 180 φ s t 2 4 6 8 10 12 14 15 a [mm] [d eg ] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.11 0.12 0.13 0.14 Fe ed R at e [m m /th .re v] b 2 4 6 8 10 12 14 15 a [mm] (a) (b) (c) (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 0 0.2 0.4 0.6 0.8 1 7 8 9 10 11 12 13 14 15200 400 600 800 1000 1200 b 0 50 100 150 180 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.24 0.28 0.32 0.36 0.4 M R R [c m 3 / m in ] a [mm] φ s t a [mm] [d eg ] Fe ed R at e [m m /th .re v] b (a) (b) (c) (d) 7 8 9 10 11 12 13 14 157 8 9 10 11 12 13 14 15 a [mm] 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.157 Chapter 5. Process Optimizationwhen 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 cut- ter in machine coordinates. For further comparison, process was also simulated in virtual environ- ment 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. The flow chart for process optimization is given in Figure 5.29. The idea here is to maximize fee- drate and spindle speed by virtually simulating the process after planning stage and recommend- ing 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, CAD Model Process Planning (CAM) Post-Process Optimization Constraints: Machining Machine Tool Optimized NC CodeNC Code Output Figure 5.28 : Optimized NC code generation flow chart with Post-Process Optimization module.158 Chapter 5. Process Optimizationmaximum chip load and maximum resultant force on the tool, where as important machine tool constraints are the torque and power characteristics. 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 optimiza- tion, 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 pos- sible. 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 minimiza- tion 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 mini- mization problem, i.e. if represents varying cutting forces, then . START END Figure 5.29 : Process optimization flow chart. F φ( ) Minimum Force min F φ( )( )⇒ Maximum Force min F φ( )–( )⇒159 Chapter 5. Process OptimizationThere 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 calcu- lation of the derivative is called Brent’s Method [92]. This method takes three initial immersion points ( ) that potentially bracket the minimum, i.e., satisfy the following relations: & , 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 through three points, namely end points = and = , and a midpoint = . Note that there are two possibilities, one with a minimum and one with a maximum as seen in Figure 5.31. F(φ) φ[deg] 360240 No of flutes = 3 φp=120 φ m in φ m ax Figure 5.30 : Sample force variation for a uniform pitch cutter. φ10 φ20 φ30, , F φ20( ) F φ10( )< F φ20( ) F φ30( )< φa φ10 φc φ30 φb φ20160 Chapter 5. Process OptimizationEither parabola can be described as: . (5.24) Substituting given immersion angles, a set of equations are obtained: . (5.25) The extremum of the parabola can be found by equating the first derivative of function F to zero: , (5.26) and from Eq (5.26), is obtained as: . (5.27) min parabola max parabola F(φ) φφa φb φc Figure 5.31 : Two possible quadratic functions. F φ( ) α2φ2 α1φ α0+ += F φa( ) α2φa2 α1φa α0+ += F φb( ) α2φb2 α1φb α0+ += F φc( ) α2φc2 α1φc α0+ += φd d F φ( ) 0 2α2φ α1+⇒ 0= = φ φ α1–2α2 --------=161 Chapter 5. Process OptimizationCombining 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: . (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 can be used for final decision such that if it is positive, the extremum is minimum, if negative, one has a maximum. Inverse parabolic interpola- tion for two iteration steps is depicted in Figure 5.32. When the function is not co-operative, for example, when three points are co-linear ( = = ) 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 α φ φ φb 12--– φb φa–( )2 F φb( ) F φc( )–[ ] φb φc–( )2 F φb( ) F φa( )–[ ]– φb φa–( ) F φb( ) F φc( )–[ ] φb φc–( ) F φb( ) F φa( )–[ ]– -----------------------------------------------------------------------------------------------------------------------------------------= φ α2 F(φ) first parabola second parabola φ10 φ30φ20 φ21 φ31φ11 φ22 φ32φ12 Iteration Steps 0 1 2 F(φ) φ Figure 5.32 : Inverse parabolic interpolation. φa φb φc162 Chapter 5. Process Optimizationforce approach. Bearing in mind that golden ratio is equal to = 0.61803, given a bracketing triplet of immersions, the next point to be tried is placed at a fraction (1- ) into the larger of the two intervals measuring from the middle point : , (5.29) where . Then the update decision is again made based on bracketing the minimum, i.e., if < then replaces the midpoint and becomes an end point, else if < then remains midpoint with replac- ing one of the end points, or . Either way the width of the bracket, ( - ), reduces until a desired tolerance is achieved. In order to account for all local extremums, one tooth period (or one spindle period for non-uni- form 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 ini- tial points ( ) 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 depen- dent on the feedrate in the following form provided that cutting force coefficients are independent of the feedrate: , ( x, y, z, t, r, trq, pwr), (5.30) GR GR φb φnew φb Δ– φ( ) if φb φa–( ) φc φb–( )>( ) φb Δφ+( ) if φc φb–( ) φb φa–( )>( )⎩ ⎭⎪ ⎪ ⎨ ⎬⎪ ⎪ ⎧ ⎫ = Δφ 1 GR–( ) max φb φa–( ) φc φb–( ),{ }×= F φnew( ) F φb( ) φnew φb F φb( ) F φnew( ) φb φnew φa φc φc φa φ10 φ20 φ30, , Fq φ( ) A0q φ( ) A1q φ( ) .c+= q ≡163 Chapter 5. Process Optimizationwhere 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 and are functions of process type, tool geometry, and time (or angular position). In Chap- ter 3, a numerical search technique has been presented to obtain the angular position at which a process output becomes maximum, which is denoted as in this section. At the most critical instance, process output -Eq (5.30)- reaches the user-defined maximum, : , (5.31) for which maximum allowable feedrate can be solved as: . (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: . (5.33) Similar to Eq (5.31), maximum resultant force is equated to user-defined maximum , i.e. , and the maximum feedrate is obtained by solving quadratic equation as: , (5.34) where . From two possible solutions in Eq (5.34), smaller positive real answer is selected. A0 q A1 q φmaxq Fq max, Fq φmaxq( ) Fq max,= cmax q Fq max, A0 q φmaxq( )– A1 q φmaxq( ) ---------------------------------------------= Fres φ( ) A0res φ( ) A1res φ( ) c A2res φ( ) c2+ += Fres max, Fres φmaxres( ) Fres max,= cmax res A1 res φmaxres( )– Δ± 2A2 res φmaxres( ) ----------------------------------------------= Δ A1res φmaxres( )[ ] 2 4 A0 res φmaxres( ) Fres max,( )2–[ ] A2res φmaxres( )⋅ ⋅–=164 Chapter 5. Process OptimizationForm errors in milling result from the static tool deflection at the entry ( ) 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 , where 0.8 is the approxi- mate 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 . Total tool deflection reflected at the surface generation point can be found by adding up the deflections at caused by y direction forces of each CEF: , (5.35) where deflection is defined by [95] as: , (5.36) where E is the Young’s Modulus and I is the area moment of inertia of the tool. Cutting forces act- ing 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 applica- ble: φ 0= φ π= Reff 0.8Rtool= zgp zgp δ δi i 1= # of CEFs ∑= δi Fy i, .zi 2 6EI --------------- 3zi zgp–( ), if 0<zgp zi< Fy i, .zi 2 6EI --------------- 3zgp zi–( ), if zi zgp<⎝⎜ ⎜⎜ ⎜⎜ ⎛ =165 Chapter 5. Process Optimization. (5.37) The tool position, , at which maximum surface deflection is generated can be calculated; however, it is not trivial to come up with deflection coefficients and explicitly. A more practical approach to obtain these two coefficients is using linear interpolation. Given two arbi- trary feedrates, and , the corresponding deflections and are calculated. Utilizing those two points, maximum allowable feedrate for a user-defined deflection constraint can be found as: . (5.38) Note that if is solved as a negative number, it means that edge forces are already enough to create a deflection more than the specified deflection constraint ; therefore, under no circum- stances, can such a constraint be satisfied. δ φ( ) A0dfl φ( ) A1dfl φ( ) c⋅+= φ- Immersion Angle z - A xi al D ep th o f C ut FLUTE - 1 1 3 Fy,2 Fy,1 Fy,3CEF CEF δ zgp 0 π FLUTE - 2 FLUTE - 3 2 CEF z1 z2 z3 Figure 5.33 : CEFs and static tool deflection at surface generation point, .zgp φmaxdfl A0 dfl A1 dfl c1 c2 δ1 δ2 δmax cmax dfl c1 c2 c1–( ) δmax δ1– δ2 δ1– ---------------------+= cmax dfl δmax166 Chapter 5. Process Optimization5.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 after- wards 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: . (5.39) In addition to constraints presented previously, which are functions of only feed, two more con- straints 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/work- piece. 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): , (5.40) fobj rMRR c n⋅= = n nmax i,≤ n– nmin i,–≤ ⎭⎬ ⎫167 Chapter 5. Process Optimizationwhere 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 , 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, , then the process is already unstable, therefore, optimization will not be applica- ble. 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. 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 optimi- zation (nonlinear programming) is employed by sequential quadratic programming (SQP). Two main stages of implementation, i.e. the Hessian matrix update, and quadratic programming (QP) alim min, alim max, Spindle Speed (n) D ep th o f C ut ( a) alim,min alim,max nmin nmax a n m in ,1 n m ax ,1 n m in ,2 n m ax ,2 High Torque Region Low Torque Region Stable Region Unstable Region Figure 5.34 : Sample stability lobes.168 Chapter 5. Process Optimizationhave the following properties: the Hessian matrix is updated using the BFSG method and qua- dratic 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]. Sud- den 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 work- piece, 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 fig- ure, 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 Optimizationtion 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. 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 ( ), commanded ( ), and end ( ) 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 Feasible Region 0.05 0.1 0.15 0.2 0.25 1 2 3 4 5 6 7 8 9 10 11 12 c - Feed rate [mm/tooth.rev] n - S pi nd le S pe ed x 1 e3 [r pm ] P constraint f obj,-2000 f obj,-2500 f obj,-1500 f obj,-1000 f obj,-500 Pconstraint =P(φmax)-Pmo(n)=0 Tconstraint =T(φmax)-Tmo(n)=0 T constraint f obj,Optim um Di rec tio n o f in cre as e f or the ob jec tiv e f un cti on (e MR R) Infeasible Region Figure 5.35 : Graphical representation of design space. Star is located at c=0.2088 [mm/tooth.rev] and n=5250 [rpm]. fs fc fe170 Chapter 5. Process Optimizationallowable 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. Feed type stated in the middle circle can only be preceded by the types shown on its left; and suc- ceed by the ones on its right. Each type has also its unique acceleration and deceleration proper- ties: 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. 0 20 40 60 80 100 120 140 160 180 4 5 6 7 8 9 10 11 12 13 14 15 f - F ee dr at e x 1 e3 [m m /m in ] s - Path Length [mm] Optimized Feedrate ! ! ! ! ! Figure 5.36 : An example of optimized feedrate with sharp variations.171 Chapter 5. Process OptimizationFeedrates 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. TYPE 1 (1.a) (1.b) (1.c) TYPE 2 TYPE 3 TYPE 4 fc fs fe fc fc fc fc fc fs fs fs fs fs fe fe fe fe=fe fs (fe>fs) (fs>fe) Figure 5.37 : Types of feed command used for filtering. 1 2 3 3 4 2 2 3 1 2 3 1 4 1 2 4 1 4 3 4 Figure 5.38 : Allowable combination of different feed block types (TYPE numbers are indicated in circles).172 Chapter 5. Process Optimization5.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 cruis- ing at a constant velocity in between. Kinematic profiles for a constant acceleration and decelera- tion 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 , , as: , (5.41) f - F ee dr at e s - Path Length Optimized Feedrate fc fs f i i i sisi-1 ith block s10 sN Figure 5.39 : Parameter definitions for a single feed block i. T1 fc i fs i – A⁄= T2 0= T3 fc i fe i – D⁄= Li fc i( )2 fs i( ) 2 – 2A -------------------------------- fc i( )2 fe i( ) 2 – 2D --------------------------------+=173 Chapter 5. Process Optimizationwhere A and D are acceleration and deceleration limits, respectively. If an acceleration or deceler- ation 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 = to be realizable, the actual path length has to be long enough to accommodate minimum length requirement; . (5.42) The equality case of the above inequality represents a motion with no constant velocity section, i.e., . 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 neighbor- ing blocks; therefore, a set of update rules has to be established. 1 … N, , si si 1––( ) Li≥ T2 0= A Position [mm] Velocity [mm/sec] Acceleration [mm/sec2] Time [sec] Time [sec] Time [sec] fc i si si-1 fs i fe i -D T 1 T 2 T 3 Figure 5.40 : Constant acceleration kinematic profiles.174 Chapter 5. Process OptimizationConsider 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 blend- ing feed blocks or it can be a different type. Update rules are classified under three major catego- ries: 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 com- mands, and in the Feed-array, both preceding and successive blocks have to be updated due to 175 Chapter 5. Process Optimizationsymmetric nature of TYPE 1.c, i.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. 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 fs i fe i= FORWARD UPDATE FULL UPDATE fc fc 1.a 2+ 3 i i+1 = 1.a 1+ 4 = 4 3+ 3 = 4 4+ 4 = 2 1+ 1.b = 3 4+ 1.b = 2 2+ 2 = 3 3+ 2 = 2 2+ 1.c =+ 3 2 1+ 1.c =+ 4 3 3+ 1.c =+ 3 FWD IUpdate Code Update Rule Configuration Update Category: FWD II FWD III FWD IV BACKWARD UPDATE Update Code Update Rule Configuration Update Category: BKWD I BKWD II BKWD III BKWD IV Update Code Update Rule Configuration Update Category: FULL I FULL II FULL III fc i fc i fc i fc i fc i fc i fc i fc ifc ifc i fc i+1 fc i+1 fc i+1 fc i-1 fc i-1 fc i-1 fc i-1 fc i-1fc i-1 fc i +1 fc i+1fc i-1 fc i +1 fs i fs i fs i fs i fe i+1 fe i+1 fe i+1 fe i+1 fe i+1 fe i+1 fe i+1 fs i-1 fs i-1 fs i-1 fs i-1 fs i-1 fs i-1 fs i-1 fe i fe i fe ife i176 Chapter 5. Process OptimizationFigure 5.41 : Update algorithm flow chart.177 Chapter 5. Process Optimizationinput, and it returns updated the Feed-array with updated index i. In some cases, when very nar- row 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 satis- fied before program moves to the next feed block. An example of feed filtering with recursive updating is given in Figure 5.42. 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 Max- imum Feed Correction (MFC). For a known path length of feed block i, , maxi- mum feedrate for each feed block type is calculated as Optimized Optimized & Filtered 2 4 6 8 10 12 0 1 2 3 4 5 6 7 8 9 10 Path Length x 10 [mm] Fe ed ra te x 1 e3 [m m /m in ] B1 is realizable: (i) Move to next feed block, B2 B0 B1 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 B3 Path Length Path Length Path Length Filtering is completed! STEP I STEP II STEP III STEP IV Objective: Filter feed blocks B0-B3 B0 B1B1 B2 B0 B1 B0 Figure 5.42 : Feed filtering example with recursive updating. Li si si 1––( )=178 Chapter 5. Process Optimization. (5.43) 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 opti- mized feedrates more closely, which results in increased productivity compared to filtering with- out 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 con- sidered. Bearing in mind that feed fluctuations along tool path will cause feed marks on the fin- ished 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. fmax i fs i( )2 D⋅ fe i( ) 2 A⋅ 2 A D si si 1––( )⋅ ⋅+ + A D+ -------------------------------------------------------------------------------------------------- TYPE 1 fs i( )2 2 A si si 1––( )⋅ ⋅+ 2 ---------------------------------------------------------- TYPE 2 fe i( )2 2 D si si 1––( )⋅ ⋅+ 2 ----------------------------------------------------------- TYPE 4⎝⎜ ⎜⎜ ⎜⎜ ⎜⎜ ⎜⎜ ⎜⎛ =179 Chapter 5. Process OptimizationOptimized feedrates given in Figure 5.36 are filtered out using the proposed algorithm. 1000 [mm/ ] ( 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. Optimized Optimized & Filtered 2 4 6 8 10 12 0 1 2 3 4 5 6 7 8 9 10 Path Length x 10 [mm] Fe ed ra te x 1 e3 [m m /m in ] B1 is realizable: (i) Move to next feed block, B2 B0 B1 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 maximum feed correction B3 Path Length Path Length Path Length Filtering is completed! STEP I STEP II STEP III STEP IV Objective: Filter feed blocks B0-B3 B0 B1B1 B2 B0 B1 B0 B2 Figure 5.43 : Feed filtering example with recursive updating & MFC. s2 ≅ 0 20 40 60 80 100 120 140 160 180 4 5 6 7 8 9 10 11 12 13 14 15 Fe ed ra te x 1 e3 [m m /m in ] Path Length [mm] Optimized Optimized & Filtered 0 20 40 60 80 100 120 140 160 180 4 5 6 7 8 9 10 11 12 13 14 15 Fe ed ra te x 1 e3 [m m /m in ] Path Length [mm] Optimized Optimized & Filtered (a) (b) Figure 5.44 : Filtered optimum-feed-commands for 1 [g] acc/dec; filtering (a) without MFC, (b) with MFC.180 Chapter 5. Process Optimization5.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 cir- cular. 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 engage- ment 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 previ- ous 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 con- tinuously 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, opti- mized 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 Optimizationpoints, 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/ 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 . . . 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 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 Optimization5.8. Verification for Post-Process Optimization The existing NC program for the die part described in the Chapter 3 was optimized using the pro- posed 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. (X,Y,Z)1 (X,Y,Z)2 (X,Y,Z)1.1 (X,Y,Z)1.2 (X,Y,Z)1.3 OFF1 OFF2 OFF3 OFF4 (X,Y,Z)1 (X,Y,Z)2(X,Y,Z)1.1 (X,Y,Z)1.2 (X,Y,Z)1.3 OFF2 OFF3 OFF4OFF1 (X,Y,Z)1 (X,Y,Z)2 OFF1 = min(OFF1,...,OFF4) G01: Linear Segment G02-G03: Circular Segment Feedrate can be varied along the path Uniform feedrate is required along the path : Tool coordinate from Original CL File : Intermediate point due to sampling for optimization Figure 5.46 : Feedrate update strategy for linear versus circular segments. Figure 5.47 : Estimated Torque and Power Curves.183 Chapter 5. Process OptimizationController parameters tabulated in Table 5.5 were also required for making reasonable estima- tions 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. 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, max- imum chip thickness was also specified for each cutter based on Sandvik Coromant’s Metal Cut- ting Technical Guide and Main Catalogue. Lastly, maximum spindle speed was set for some of Table 5.5 : Assumed Characteristics of the Heidenhain Controller. Table 5.6 : Other parameters of the machine tool used for verification.184 Chapter 5. Process Optimizationthe operations considering dynamic characteristics, speed limitations of the spindle, and recom- mended 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. 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 con- tained 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 diffi- culties 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 Table 5.7 : Optimization Criteria for Different Operations and Cutting Tools. Opt. Code Maximum Chip Thickness [mm] Maximum Spindle Speed [rpm] Comments 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, T4185 Chapter 5. Process Optimizationincluded in the comparison. Leaving operations 6-9 out of analysis do not affect results signifi- cantly 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). 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 maxi- mize 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 Table 5.8 : Cycle time results of operations. Cycle Times Cycle Time [hr:min:second] Real Predicted % Error Real Predicted %Error Difference Operation 1 1:58:00 1:57:00 -0.8% 1:28:00 1:17:24 -12.0% -25.4% Operation 2 0:03:10 0:03:13 0:02:29 -22.7% -21.5% Operation 3 2:58:33 2:05:00 2:28:00 18.4% -17.1% Operation 4 0:10:50 0:08:04 -25.0% 0:05:00 0:04:13 -15.6% -47.7% Operation 5 0:46:12 1:01:00 0:58:26 -4.2% 26.5% (Semi-Finishing) Operation 5 2:30:14 3:30:00 3:55:00 11.9% 37.7% (Finishing) Operation 6 0:11:00 Operation 7 0:19:27 Operation 8 0:14:06 Operation 9 0:18:50 Original VMS186 Chapter 5. Process OptimizationFigure 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 program- mer to foresee how the cutting conditions would change during machining, the feedrate was con- servatively 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 disap- peared due to a decrease in cutting force coefficients when chip thickness became larger. Compar- ison 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. Figure 5.48 : Chip formation with (a) original, (b) optimized NC programs. (a) (b)187 Chapter 5. Process OptimizationAn 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- O R IG IN A L O PT IM IZ E D Figure 5.49 : Pictures of T63 from front and bottom after execution of original and optimized CL files. Figure 5.50 : T20F - Part of the tool path from Operation 5, Semi-Finishing.188 Chapter 5. Process Optimizationness 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 fee- drate 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 con- straints. (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 Com- puter Aided Manufacturing (CAM) environment. The part-cutter engagement conditions are iden- tified 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, dimen- sional 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 automat- ically 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 contribu- tions of the thesis belong to the second component, which governs the physics of the machining 190 Chapter 6. Conclusionsprocesses. 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 signifi- cantly. 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 hav- ing arbitrary geometries. The model has two stages; pre-process planning and post-process simu- lation. 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. Conclusionscut. 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 pro- cess states, i.e. force, torque, power, tool deflection, and chip thickness, are used as process mea- sures, 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 geom- etries 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 con- straints 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 automati- cally adjusting feeds and speeds along the tool path. The novel algorithms are published in archi- val 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. Conclusionsence 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, ser- rated and variable pitch cutters. The algorithms of this work are finally integrated into a prototype industrial product: Virtual Mill- ing System. 6.2. Future Research Directions There are remaining key challenges before achieving a fully virtual machining of parts using dig- ital 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 computa- tional cost increases dramatically if the grid size is reduced. It is necessary to develop computa- tionally efficient yet accurate cutter-part engagement identification algorithms. The process mechanics and dynamics models must be able to cover 5-axis milling, turning, bor- ing, 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 prac- tical accuracy remain unsolved. The contact between the flank face of the tool and wavy finish 193 Chapter 6. Conclusionssurface 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. 1203- 1209. [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 sur- faces, 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 Engineer- ing 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 pulsa- tions 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 Engineer- ing Transactions, pp. 297. [26] Altintas, Y., and Spence, A., 1991, End Milling Force Algorithms for CAD Systems, Man- ufacturing 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 Applica- tions, 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 Engineer- ing, 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/prod- ucts/acis.html. [46] Saturley, P.V. and Spence, A.D., 2000, Integration of Milling Process Simulation with On- Line Monitoring and Control, International Journal of Advanced Manufacturing Technolo- gies, 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 Manufac- ture, 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 Tech- nology, 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, Interna- tional 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 sur- faces, 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 sur- faces, International Journal of Machine Tools and Manufacture, Vol. 46, n 7-8, pp. 747- 757. [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 Engineer- ing, 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.5- 14. [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 Engi- neering 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 mill- ing 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 Engineer- ing, 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, Interna- tional 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-dif- ferential 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 matri- ces, 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 gener- alized 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. 459- 469.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 sys- tems, Journal of Engineering Mechanics, Vol. 129-8, pp. 935-943. [86] Armarego, E.J.A., 1993, Material Removal Process: An Intermediate Course, The Univer- sity 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 Fee- drate Selection for 3-Axis Milling, Proc. IMECE Symp. on Dynamics and Control of Mate- rial 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 Applica- tions, Vol. 9-1, pp. 26-39. [92] W. H. Press et al., 1988, Numerical Recipes in C: The Art of Scientific Computing. Cam- bridge University Press. [93] Campomanes, M., 1998, Dynamics of Milling Flexible Structures, The University of Brit- ish 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 coeffi- cients, 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 Interna- tional Journal Series C: Mechanical Systems, Machine Elements and Manufacturing, Vol. 44-3-September, pp. 717-723.202
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Virtual three-axis milling process simulation and optimization
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Virtual three-axis milling process simulation and optimization Merdol, Doruk Sūkrū 2008
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Virtual three-axis milling process simulation and optimization |
Creator |
Merdol, Doruk Sūkrū |
Publisher | University of British Columbia |
Date Issued | 2008 |
Description | 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 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. |
Extent | 4559999 bytes |
Subject |
Virtual machining Milling optimization Virtual milling simulation |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-04-08 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0066338 |
URI | http://hdl.handle.net/2429/658 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2008-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
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}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0066338/manifest