A Transient Thermal Model for Machining by Coskun Islam B.Sc., Istanbul Technical University, 2008 M.Sc., Koc University, 2011 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) October 2018 © Coskun Islam, 2018 ii The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the dissertation entitled: A Transient Thermal Model for Machining submitted by Coskun Islam in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Mechanical Engineering Examining Committee: Dr. Yusuf Altintas Supervisor Dr. Hsi-Yung Feng Supervisory Committee Member Dr. Dana Grecov University Examiner Dr. Steve Cockcroft University Examiner iii Abstract Prediction of temperature in the tool, chip and workpiece surface layer is essential for tool design and the selection of most productive cutting conditions which yield the desired tool life and acceptable residual stresses left on the machined part. This thesis presents a comprehensive, finite difference method based numerical model on simulating the temperature distribution in the chip, tool, and finished workpiece surface layer as a function of material properties, cutting speed, feed rate and tool-workpiece engagement period. The heat is generated in the primary shear zone where the chip is sheared from the metal, in the secondary zone where the chip sticks and slides on the rake face, and in the tertiary zone where cutting-edge ploughs the workpiece surface. The chip, layer of the workpiece surface and the tool edge are meshed into discrete elements. The heat is transferred to the stationary tool, and dynamically moving chip and workpiece surface by conforming heat balance equations within each element. A finite difference technique with implicit time discretization is used to solve heat balance equations of the temperature fields on the tool, workpiece, and chip. Anisotropic material properties can be considered in the model which allows the inclusion of a coating layer on the tool. The proposed model allows two and three-dimensional heat transfer, hence it can be used to predict the temperature distribution in turning, drilling and milling operations. The continuous machining processes such as turning generate constant heat, so the temperature reaches a steady state after a transient period. The intermittent operations such as milling generate time-varying and periodic heat, hence the temperature variation is always in a transient state. The proposed model is experimentally validated with the data found in the literature and experiments conducted by the author at the industrial partner’s (Sandvik Coromant AB, iv Sweden) research facility. Experimental validations cover uncoated, single and multi-layer coated tools to simulate continuous turning, interrupted turning and milling operations. The proposed model is able to predict the temperature with less than 20% error in most of the validated cases. v Lay Summary The advancements in machining technology require increased productivity at a lower cost despite the challenges brought by new materials and machining methods. The generation of high heat and the resulting temperatures cause reduced tool life that leads to increased manufacturing costs. The cyclic thermal loading of the machined parts also reduces their fatigue life which is crucial especially in the aerospace industry. The thermal analysis of the cutting process is vital for tool design and planning of machining operations. This thesis proposes a novel and comprehensive simulation model to predict the temperature distribution in the tool, chip, and the machined surface layer as a function of tool and workpiece material properties, coating material and layer thickness, cutting speed, feed rate, and tool-workpiece engagement period. The developed method enables the tool manufacturers to design the most suitable cutting tool geometry and coating layer, and process planners to optimize cutting conditions in a virtual environment. vi Preface This Ph.D. dissertation proposes a transient thermal machining model that is capable of simulating continuous and interrupted machining processes. All the modeling work is conducted by the Ph.D. candidate in the Manufacturing Automation Laboratories at the University of British Columbia, under the supervision of Professor Yusuf Altintas. The research chapters of this dissertation are either already published or under preparation for submission for journal review. The contributions of the Ph.D. candidate for each chapter are explained in this section. A version of Chapters 3, 4 and 5.2 has been published in [1]. “Islam, C., Lazoglu, I. and Altintas, Y., 2016. A three-dimensional transient thermal model for machining. Journal of Manufacturing Science and Engineering, 138(2), p.021003.”. The manuscript was written by Ph.D. Candidate, and edited by his supervisor Dr. Altintas and visiting professor Dr. I. Lazoglu. Ph.D. Candidate was responsible for formulations of all the mathematical models and analysis of the simulation results. The mathematical model used in this paper does not include anisotropic material properties. A journal article is under preparation with the contents from Chapters 3, 4, 5.3, 5.4, 5.5, and 5.6. This article will include the temperature simulation with anisotropic material properties. Tool coating application, comparison of static contact and sliding contact, heat partition approximations will be included in the article. vii Table of Contents Abstract ........................................................................................................................................ iii Lay Summary ................................................................................................................................. v Preface ........................................................................................................................................... vi Table of Contents ......................................................................................................................... vii List of Tables .................................................................................................................................. x List of Figures .............................................................................................................................. xii List of Symbols ............................................................................................................................ xix List of Abbreviations ............................................................................................................. xxviii Acknowledgements ................................................................................................................... xxix Chapter 1: Introduction ................................................................................................................ 1 Chapter 2: Review of Temperature Prediction Literature in Machining ................................ 6 2.1 Overview ....................................................................................................................... 6 2.2 Analytical temperature prediction models ..................................................................... 6 2.3 Modeling of temperature in cutting with numerical models ......................................... 9 2.3.1 Finite element methods used in cutting temperature prediction ................................. 9 2.3.2 Finite difference methods used in cutting temperature prediction ........................... 13 2.4 Summary ...................................................................................................................... 15 Chapter 3: Heat Transfer Model of Metal Cutting .................................................................. 17 3.1 Overview ..................................................................................................................... 17 3.2 Modeling of heat sources in metal cutting ................................................................... 17 3.2.1 Heat Source in the primary deformation zone .......................................................... 18 viii 3.2.2 Heat source in the secondary deformation zone ....................................................... 21 3.2.3 Heat source in the tertiary deformation zone ........................................................... 26 3.3 Heat balance equations ................................................................................................ 26 3.3.1 Constant thermal property heat balance equations ................................................... 26 3.3.2 Variable thermal property heat balance equations ................................................... 29 3.3.3 Discretization of heat balance equations .................................................................. 30 Chapter 4: Application of Heat Transfer Model to Metal Cutting ......................................... 39 4.1 Overview ..................................................................................................................... 39 4.2 Computational domain transformation ........................................................................ 39 4.3 Discretization of cutting geometry .............................................................................. 46 4.4 Grid generation ............................................................................................................ 47 4.5 Chip temperature model .............................................................................................. 54 4.6 Tool temperature model ............................................................................................... 62 4.7 Workpiece temperature model ..................................................................................... 70 4.8 Simulation algorithm ................................................................................................... 76 Chapter 5: Simulations and Experimental Validations of Temperature Model ................... 81 5.1 Overview ..................................................................................................................... 81 5.2 Validation of constant thermal property approximation .............................................. 81 5.3 Validation of variable thermal property approximation .............................................. 97 5.3.1 Cutting Tests ............................................................................................................. 99 5.3.2 Temperature Simulations ........................................................................................ 106 5.4 Simulations and experimental results with tool coating ............................................ 121 5.4.1 Single-layer coated tool simulations ...................................................................... 123 ix 5.4.2 Multi-layer coated tool simulations ........................................................................ 139 5.5 Comparison between static contact and sliding heat partition assumptions .............. 151 5.6 Simulation input parameter sensitivity results ........................................................... 156 Chapter 6: Conclusion and Future Research Directions ....................................................... 164 6.1 Conclusions ............................................................................................................... 164 6.2 Future research directions .......................................................................................... 168 References ................................................................................................................................... 169 Appendix A : Oblique Cutting Mechanics ............................................................................ 175 x List of Tables Table 5-1-Thermal properties of the tool and workpiece materials for validations .................... 82 Table 5-2- Workpiece material orthogonal cutting constants ...................................................... 82 Table 5-3- Ti6Al4V orthogonal cutting constants ....................................................................... 95 Table 5-4- Insert grade information ........................................................................................... 100 Table 5-5-Tool and workpiece thermal properties for AISI316L simulations .......................... 106 Table 5-6- Workpiece material orthogonal cutting constants for AISI 316L ............................ 108 Table 5-7- Workpiece material orthogonal cutting constants for SS2541 extracted from [72] 124 Table 5-8-Tool substrate, tool coating and workpiece thermal properties for SS2541 simulations ................................................................................................................................................... 125 Table 5-9- Comparison of simulated and measured maximum rake face temperatures for SS2541 cutting with TiN coated tool ...................................................................................................... 125 Table 5-10- Comparison of simulated and measured maximum rake face temperatures for SS2541 cutting with the uncoated tool ...................................................................................... 126 Table 5-11- Example temperature measurements from [17] ..................................................... 141 Table 5-12- Workpiece material orthogonal cutting constants for AISI 1045 taken from CUTPROTM ............................................................................................................................... 141 Table 5-13-Tool substrate, tool coating and workpiece thermal properties for AISI 1045 simulations ................................................................................................................................. 142 Table 5-14-Variable tool substrate and workpiece thermal properties for AISI 1045 simulations ................................................................................................................................................... 143 Table 5-15- Comparison of simulated and measured mean rake face temperatures for AISI1045 cutting with TiC/Al2O3/TiN coated tool .................................................................................... 143 xi Table 5-16- Comparison of simulated and measured mean rake face temperatures for AISI1045 cutting with TiC/Al2O3/TiN coated tool (Temperature dependent material properties simulation result) ......................................................................................................................................... 144 Table 5-17- Constant thermal property static contact heat partition ratio approximation by using the equation (4.38) ..................................................................................................................... 152 Table 5-18- Static contact assumption heat partition ratios towards tool for AISI 316L simulation cases ......................................................................................................................... 153 Table 5-19- Static contact assumption heat partition ratios towards tool for TiN coated and uncoated tool simulation cases .................................................................................................. 154 Table 5-20- Static contact assumption heat partition ratios towards TiC/Al2O3/TiN coated tool for AISI1045 simulation cases................................................................................................... 155 Table 5-21- Ti6Al4V orthogonal cutting constants for sensitivity simulations ........................ 157 Table 5-22- Workpiece thermal properties for Ti6Al4V ........................................................... 157 xii List of Figures Figure 1-1- Deformation zones ...................................................................................................... 2 Figure 3-1 – Oblique cutting geometry ....................................................................................... 19 Figure 3-2 - Cutting force and velocity vectors in oblique cutting .............................................. 19 Figure 3-3- Assumed shear stress distribution in the secondary deformation zone .................... 23 Figure 3-4- Schematic of a control volume in cartesian coordinates .......................................... 27 Figure 3-5- Sample volume element in discretized chip domain ................................................ 31 Figure 3-6 Implicit solution representation in time ..................................................................... 34 Figure 3-7- Element to element varying conduction coefficient demonstration ......................... 36 Figure 4-1- Physical to computational space transformation ...................................................... 40 Figure 4-2 Discretized cutting geometry ..................................................................................... 47 Figure 4-3- Gridpoint mapping from physical to computational space ....................................... 48 Figure 4-4 Layout of chip substructures for a turning case ......................................................... 54 Figure 4-5 Discretization for the milling operation ..................................................................... 55 Figure 4-6 Chip substructure ....................................................................................................... 58 Figure 4-7 -Tool substructure layout ........................................................................................... 63 Figure 4-8 - Substructure shapes for tool ................................................................................... 66 Figure 4-9- Workpiece substructure layout ................................................................................. 70 Figure 4-10- Sample workpiece substructure .............................................................................. 72 Figure 4-11 – The flow chart of temperature prediction ............................................................. 79 Figure 5-1-Ti6Al6V-2Sn rake face simulation results compared to experimental results reported by Kitagawa et al. in [66] ............................................................................................................ 84 xiii Figure 5-2- Ti6Al6V-2Sn interrupted cutting rake face temperatures through one cutting and non-cutting period ........................................................................................................................ 85 Figure 5-3- Ti6Al6V-2Sn interrupted cutting tool and chip thermal gradients through one cutting and non-cutting period (time marker numbered with roman numbers are shown in Figure 5-2 . 86 Figure 5-4- Average heat partition to the tool during continuous turning of Ti6Al6V-2Sn ....... 88 Figure 5-5- Average heat partition to the tool during interrupted turning of Ti6Al6V-2Sn ....... 88 Figure 5-6- Turning Ti6Al6V-2Sn temperature distributions on the chip, tool, and workpiece at steady-state .................................................................................................................................. 89 Figure 5-7- The comparison of average rake face temperature prediction in the interrupted turning of gray cast iron with the experiments and computations reported in [6] ...................... 91 Figure 5-8- Average heat partition into the tool during interrupted turning of Gray Cast Iron .. 91 Figure 5-9- The comparison of average rake face temperature prediction in the interrupted turning of Al2024 with the experiments and computations reported in [6] for f=0.109 mm/rev 92 Figure 5-10- Average heat partition into the tool during interrupted turning of Al2024 for f=0.109 mm/rev ........................................................................................................................... 93 Figure 5-11- The comparison of average rake face temperature prediction in the interrupted turning of Al2024 with the experiments and computations reported in [6] for f=0.165 mm/rev 93 Figure 5-12 Average heat partition into the tool during interrupted turning of Al2024 for f=0.165 mm/rev ........................................................................................................................... 94 Figure 5-13- The comparison of simulated and measured forces for down milling of Ti6Al4V with V=3.67 m/s and f=0.105 mm/tooth, measurements from Sato et al. [19] .......................... 95 Figure 5-14- Simulated and measured resultant rake face temperature 0.1 mm below the tool-chip interface in down milling of Ti6Al4V, measurements from Sato et al. [19] ...................... 96 xiv Figure 5-15- Average heat partition into the tool during down milling of Ti6Al4V .................. 97 Figure 5-16-Orthogonal cutting specimen on the lathe ............................................................. 100 Figure 5-17-Drawing of insert geometry and thermocouple locations ...................................... 101 Figure 5-18- Example microscope measurement of an AISI 316L cutting insert (462 grade) . 102 Figure 5-19-Insert close-up picture with chip’s relative position .............................................. 102 Figure 5-20- Data acquisition system set-up for thermal measurements .................................. 103 Figure 5-21-Thermocouples on the tool holder without silica gel ............................................ 104 Figure 5-22-Simulation and thermocouple measurement comparison of AISI 316L turning at V=3 m/s and f=0.2 mm/rev ........................................................................................................ 110 Figure 5-23- Evolution of average rake face temperature with time for turning of AISI 316L at V=3 m/s and f=0.2 mm/rev ........................................................................................................ 110 Figure 5-24- Average heat partition to tool for turning of AISI 316L at V=3 m/s and f=0.2 mm/rev ....................................................................................................................................... 111 Figure 5-25-Simulated temperature distributions at steady-state for turning of AISI 316L at V=3 m/s and f=0.2 mm/rev ................................................................................................................ 112 Figure 5-26-Simulation and thermocouple measurement comparison of AISI 316L turning at V=1.67 m/s and f=0.2 mm/rev ................................................................................................... 113 Figure 5-27- Evolution of average rake face temperature with time for turning of AISI 316L at V=1.67 m/s and f=0.2 mm/rev ................................................................................................... 114 Figure 5-28- Average heat partition to tool for turning of AISI 316L at V=1.67 m/s and f=0.2 mm/rev ....................................................................................................................................... 114 Figure 5-29-Simulated temperature distributions at steady state for turning of AISI 316L at V=1.67 m/s and f=0.2 mm/rev ................................................................................................... 115 xv Figure 5-30-Simulation and thermocouple measurement comparison of AISI 316L turning at V=1.17 m/s and f=0.1 mm/rev ................................................................................................... 116 Figure 5-31- Evolution of average rake face temperature with time for turning of AISI 316L at V=1.17 m/s and f=0.1 mm/rev ................................................................................................... 117 Figure 5-32- Average heat partition to tool for turning of AISI 316L at V=1.17 m/s and f=0.1 mm/rev ....................................................................................................................................... 117 Figure 5-33-Simulated temperature distributions at steady state for turning of AISI 316L at V=1.17 m/s and f=0.1 mm/rev ................................................................................................... 118 Figure 5-34-Simulation and thermocouple measurement comparison of AISI 316L turning at V=1.17 m/s and f=0.2 mm/rev ................................................................................................... 119 Figure 5-35- Evolution of average rake face temperature with time for turning of AISI 316L at V=1.17 m/s and f=0.2 mm/rev ................................................................................................... 120 Figure 5-36- Average heat partition to tool for turning of AISI 316L at V=1.17 m/s and f=0.2 mm/rev ....................................................................................................................................... 120 Figure 5-37-Simulated temperature distributions at steady state for turning of AISI 316L at V=1.17 m/s and f=0.2 mm/rev ................................................................................................... 121 Figure 5-38- Simulation boundary conditions and coating layer application ............................ 122 Figure 5-39- Rake face temperature simulations compared to infrared measurements in [72] for turning of SS2541 at V=3.33 m/s, f=0.1 mm/rev with TiN coating after 100 ms cutting time . 127 Figure 5-40- Evolution of average rake face temperature with time for turning of SS2541 at V=3.33 m/s, f=0.1 mm/rev with TiN coating ............................................................................ 128 Figure 5-41- Average heat partition to tool for turning of SS2541 at V=3.33 m/s, f=0.1 mm/rev with TiN coating ........................................................................................................................ 128 xvi Figure 5-42-Simulated temperature distributions at steady state condition for turning of SS2541 at V 3.33 m/s, f=0.1 mm/rev with TiN coating after 100 ms cutting time ................................ 129 Figure 5-43- Rake face temperature simulation result for turning of SS2541 at V=3.33 m/s, f=0.1 mm/rev with the uncoated tool after 100 ms cutting time ............................................... 131 Figure 5-44- Evolution of average rake face temperature with time for turning of SS2541 at V=3.33 m/s, f=0.1 mm/rev with uncoated tool .......................................................................... 131 Figure 5-45- Average heat partition to tool for turning of SS2541 at V=3.33 m/s, f=0.1 mm/rev with uncoated tool ..................................................................................................................... 132 Figure 5-46-Simulated temperature distributions at steady state condition for turning of SS2541 at V 3.33 m/s, f=0.1 mm/rev with the uncoated tool after 100 ms cutting time ........................ 133 Figure 5-47- Evolution of average rake face temperature with time for turning of SS2541 at V=3.33 m/s, f=0.15 mm/rev with TiN coating .......................................................................... 134 Figure 5-48- Average heat partition to tool for turning of SS2541 at V=3.33 m/s, f=0.15 mm/rev with TiN coating ........................................................................................................................ 135 Figure 5-49-Simulated temperature distributions at steady state condition for turning of SS2541 at V 3.33 m/s, f=0.15 mm/rev with TiN coating after 100 ms cutting time .............................. 136 Figure 5-50- Evolution of average rake face temperature with time for turning of SS2541 at V=3.33 m/s, f=0.05 mm/rev with TiN coating .......................................................................... 137 Figure 5-51- Average heat partition to tool for turning of SS2541 at V=3.33 m/s, f=0.05 mm/rev with TiN coating ........................................................................................................................ 138 Figure 5-52-Simulated temperature distributions at steady state condition for turning of SS2541 at V 3.33 m/s, f=0.05 mm/rev with TiN coating after 100 ms cutting time .............................. 139 xvii Figure 5-53- Evolution of average rake face temperature with time for turning of AISI 1045 at V=2.4 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating .......................................................... 145 Figure 5-54- Average heat partition to tool for turning of AISI 1045 at V=2.4 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating......................................................................................... 145 Figure 5-55-Simulated temperature distributions at steady state condition for turning of AISI 1045 at V=2.4 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating after 100 ms cutting time .... 146 Figure 5-56- Evolution of average rake face temperature with time for turning of AISI 1045 at V=2.07 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating ........................................................ 147 Figure 5-57- Average heat partition to tool for turning of AISI 1045 at V=2.07 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating......................................................................................... 147 Figure 5-58-Simulated temperature distributions at steady state condition for turning of AISI 1045 at V=2.07 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating after 100 ms cutting time .. 148 Figure 5-59- Evolution of average rake face temperature with time for turning of AISI 1045 at V=1.48 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating ........................................................ 149 Figure 5-60- Average heat partition to tool for turning of AISI 1045 at V=1.48 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating......................................................................................... 149 Figure 5-61-Simulated temperature distributions at steady state condition for turning of AISI 1045 at V=1.48 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating after 100 ms cutting time .. 150 Figure 5-62-AISI316L heat partition change with cutting speed in sliding contact conditions 153 Figure 5-63- SS2541 heat partition change with feed rate in sliding contact conditions .......... 155 Figure 5-64- AISI1045 heat partition change with cutting speed in sliding contact conditions 156 Figure 5-65- Variable and constant thermal properties' effect on tool-chip interface temperature comparison for Ti6Al4V............................................................................................................ 158 xviii Figure 5-66- Sticking and sliding contact’s effect on tool-chip interface temperature comparison for Ti6Al4V ............................................................................................................................... 160 Figure 5-67- Sticking zone thickness’s effect on tool-chip interface temperature comparison for Ti6Al4V ..................................................................................................................................... 161 Figure 5-68- Convection coefficient sensitivity of the simulations for Ti6Al4V ..................... 162 Figure 5-69- Model size effect on rake face temperature comparison for AISI 316L .............. 163 Figure A- 1- Section view of tool and workpiece on the cutting edge normal plane ................ 175 Figure A- 2- Top view of the shear plane from Z4 direction .................................................... 175 xix List of Symbols A Square coefficient matrix for chip heat balance equations cA Tool-chip contact area (m2) sA Shear plane area (m2) itotal totaluHN X Ni xA The square coefficient matrix of steady-state terms with unit anisotropic mass transfer coefficients for chip heat balance equations , ,itotal totaluKN X Ni xx yy zzA The square coefficient matrix of steady-state terms with unit anisotropic heat conduction coefficients for chip heat balance equations total totalut t N X NA Transient part of the square coefficient matrix with unit material properties for the chip B Heat partition ratio for the secondary deformation zone 2B Heat partition ratio for the tertiary deformation zone C Heat generation array for chip heat balance equations itotal totaluHN X Ni xC The square coefficient matrix with unit anisotropic mass transfer coefficients to apply constant temperature boundary conditions for the chip , ,itotal totaluKN X Ni xx yy zzC The square coefficient matrix with unit anisotropic heat conduction coefficients to apply constant temperature boundary conditions for the chip xx 1totalut N XC Coefficient array of one time-step previous temperatures with unit thermal properties for the chip 1totalut t N XC Coefficient array of two time-steps previous temperatures with unit thermal properties for the chip pC Heat capacity (J/ Kg.K) D Square coefficient matrix for tool heat balance equations , ,itotal totaluKN X Ni xx yy zzD The square coefficient matrix of steady-state terms with unit anisotropic heat conduction coefficients for tool heat balance equations total totalut t N X ND Transient part of the square coefficient matrix with unit material properties for the tool E Heat generation array for tool heat balance equations , ,itotal totaluKN X Ni xx yy zzE The square coefficient matrix with unit anisotropic heat conduction coefficients to apply constant temperature boundary conditions for the tool 1totalut N XE Coefficient array of one time-step previous temperatures with unit thermal properties for the tool 1totalut t N XE Coefficient array of two time-steps previous temperatures with unit thermal properties for the tool inE Input energy to a control volume (W) generatedE Generated energy in a control volume (W) outE Output energy from a control volume (W) xxi storedE Stored energy in a control volume (W) nF Shear normal force (N) RF Resultant cutting force (N) sF Shear force (N) uF Friction force on the rake face (N) slidinguF Sliding contact portion of friction force on the rake face (N) stickinguF Sticiking contact portion of friction force on the rake face (N) vF Friction normal force on the rake face (N) teF Tangential edge force (N) consFV Flux volume constant , ,ii x y zH Orthogonal mass transfer term (W/m2.K) consIV Internal heat generation volume constant J Computational domain transformation Jacobian (1/m2 or 1/m3) teK Tangential edge force coefficient (N/mm) , ,iii x y zK Orthogonal heat conduction coefficients in Cartesian directions (W/m.K) , ,, ,iji x y zi jj x y zK Cross-heat conduction coefficients in Cartesian directions (W/m.K) M Iteration number xxii cN Chip substructure number tN Tool substructure number wN Workpiece substructure number totalN Total number of elements in a domain xN Number of elements in x-direction yN Number of elements in y-direction zN Number of elements in z-direction N Number of elements in curvilinear coordinate direction nN Number of elements in curvilinear coordinate n direction sP Shearing power (W) fP Friction power (W) stickingP Power generated in sticking contact zone (W) slidingfP Power generated in sliding contact zone (W) ploughingP Ploughing power (W) csliding frictionQ Sliding friction heat source array for the chip csticking frictionQ Sticking friction heat source array for the chip tsliding frictionQ Sliding friction heat source array for the tool tsticking frictionQ Sticking friction heat source array for the tool ttertiaryQ The tertiary zone heat source array for the tool xxiii wtertiaryQ The tertiary zone heat source array for the workpiece R Square coefficient matrix for workpiece heat balance equations itotal totaluHN X Ni x or zR The square coefficient matrix of steady-state terms with unit anisotropic mass transfer coefficients for workpiece heat balance equations , ,itotal totaluKN X Ni xx yy zzR The square coefficient matrix of steady-state terms with unit anisotropic heat conduction coefficients for workpiece heat balance equations total totalut t N X NR Transient part of the square coefficient matrix with unit material properties for the workpiece S Heat generation array for workpiece heat balance equations itotal totaluHN X Ni x or zS The square coefficient matrix with unit anisotropic mass transfer coefficients to apply constant temperature boundary conditions for the workpiece , ,itotal totaluKN X Ni xx yy zzS The square coefficient matrix with unit anisotropic heat conduction coefficients to apply constant temperature boundary conditions for the workpiece 1totalut N XS Coefficient array of one time-step previous temperatures with unit thermal properties for the workpiece 1totalut t N XS Coefficient array of two time-steps previous temperatures with unit thermal properties for the workpiece cT Chip temperature array (K) 1total Xboundaryc NT Chip constant temperature boundary condition array (K) xxiv toldcT Chip one time-step previous temperature array (K) t toldcT Chip two time-steps previous temperature array (K) tT Tool temperature array (K) 1total Xboundaryt NT Tool constant temperature boundary condition array (K) toldtT Tool one time-step previous temperature array (K) t toldtT Tool two time-steps previous temperature array (K) wT Workpiece temperature array (K) 1total Xboundaryw NT Workpiece constant temperature boundary condition array (K) toldwT Workpiece one time-step previous temperature array (K) t toldwT Workpiece two time-steps previous temperature array (K) shearT Transient average shear plane temperature (K) ssT Average shear plane steady state temperature (K) rT Room temperature (K) V Cutting velocity (m/s) cV Chip velocity (m/s) sV Shear velocity (m/s) xxv , ,iV i x y z Mass transfer velocity components in Cartesian directions (m/s) b Width of cut (m) dB Difference measure h Depth of cut (m) ch Chip thickness (m) k Thermal conductivity (W/m.K) cl Tool-chip contact length (m) pll Tool-chip plastic contact length (m) sl Shear plane length (m) m Element number cq Chip element number fq Frictional heat flux (W/m2) genq Internal heat generation rate for a control volume (W/m3) ploughingq Ploughing heat flux (W/m2) tq Tool element number wq Workpiece element number , ,iq i x y z Heat flux input to a control volume in Cartesian directions (W) , ,i diq i x y z Heat flux output from a control volume in Cartesian directions (W) xxvi er Tool edge radius (m) secondaryt Secondary deformation zone thickness (m) st Shear plane thickness (m) tn Time step number sQ Shear zone internal heat generation (W/ m3) stickingQ Sticking zone internal heat generation (W/ m3) Thermal diffusivity (m2/s) n Normal rake angle (rad) Taylor-Quinney coefficient a Friction angle (rad) Shear strain rate (1/s) K Maximum difference between heat conduction values from subsequent recursions (W/m.K) Convergence tolerance (oC) n Discrete element size in cuvilinear coordinate n direction Discrete element size in cuvilinear coordinate direction t Time step size (s) x Difference between successive iterations for grid point positions in x y Difference between successive iterations for grid point positions in y xxvii n Curvilinear coordinate direction c Chip flow angle (rad) i Resultant cutting force orientation angle with respect to cutting edge normal plane (rad) n Resultant cutting force orientation angle with respect to cut surface (rad) s Inclination angle (rad) Curvilinear coordinate direction Density (Kg/ m3) Time constant (s) av Average shear flow stress in the secondary deformation zone (Pa) sliding Average shear stress in the sliding contact region (Pa) s Shear stress (Pa) i Oblique shear angle (rad) n Normal shear angle (rad) w Curvilinear coordinate direction xxviii List of Abbreviations EMF Electro-motive force FDM Finite difference method FEM Finite element method IR Infra-red TC Thermocouple doc Depth of Cut xxix Acknowledgements First, I would like to express my deep gratitude to Dr. Altintas for his ongoing support, guidance, and patience throughout this whole work at every stage. He will always be an inspiration to me as an engineer, scientist, and most importantly as a person. Besides Dr. Altintas, I would like to thank Dr. Melkote as the external examiner and examination committee members, Dr. Feng, Dr. Grecov, and Dr. Cockcroft for their valuable questions and comments regarding the thesis. I also thank Dr. Lazoglu for introducing me to machining research during my MSc. studies and his valuable input related to edits of the first journal article during his sabbatical at UBC. My sincere thanks goes to Sandvik Coromant AB metal cutting research team for their support during the experimental part of this work. Especially, Dr. Jonas Östby for his valuable discussions and help. Mikael Lundblad for all the support and care during my time with his team. I was fortunate to work with brilliant people from all around the world during my time at Manufacturing Automation Laboratory. I would like to thank all the labmates for their friendship, support, and valuable input related to research. I especially thank to Oguzhan, Onur, Alptunc, Deniz, and Mohit. I also would like to thank Ozgur who has been a big brother to me. His sincere friendship and support have always been very important. Also, I want to thank all my friends and the Turkish community in Vancouver who made this journey enjoyable. Last but not least, I am grateful to my family: my mother, Hatice, father, Hasan, brother, Emre, sister in law, Zeynep, niece, Ayse, my aunt Vildan, uncle Grant, cousins, Deniz and Altan for their love, support and encouragement that they provided at every stage of this work. Finally, a special gratitude goes to all of my teachers till this point in my life. xxx to Alperen Kurtulus… 1 Chapter 1: Introduction Metal chip is sheared in primary plastic deformation zone and moves on the rake face by experiencing sticking and sliding friction in the secondary zone, and the finished surface contacts with the tool’s flank face which is the third deformation region (Figure 1-1). Cutting generates heat in all three zones which increases the temperature of the work material and cutting tool. The high temperature, in localized tool-workpiece contact zone, weakens the tool’s strength and causes accelerated tool wear and chipping of the cutting edge. High-temperature gradients lead to residual stresses on the finish surface of the product which shortens the fatigue life of the workpiece. Estimation of temperature distribution as a function of workpiece and tool material properties, cutting speed and uncut chip thickness (i.e. feed) has critical importance in developing coating and machining strategies especially in machining thermal resistant/high strength materials with low thermal conductivity like Titanium and Nickel alloys used in aerospace and biomedical industry. It is known that major tool wear mechanisms are based on the diffusion and the adhesion, which are related to temperature values at the tool-chip contact region. Steep thermal gradients also lead to thermal stresses which contribute to crack initiation along with high mechanical contact pressures on the tool surface. Surface integrity of the finished workpiece is also affected by thermal gradient resulted from machining. Induced residual stress field may contribute to the low fatigue life of parts and finished part distortion. In short, the thermal analysis is vital to understand the machining process physics which lead to improved tool design and higher productivity in industry. 2 Figure 1-1- Deformation zones The generation of heat and its dissipation to chip, tool and workpiece is a dynamic process which occurs in a small zone (i.e., 1 mm x 1 mm planar area). While the cutting process is continuous in single point machining operations such as turning and boring, the process is intermitted and periodic in multi-point cutting such as milling and drilling processes. The temperature rises initially in the transient period and reaches to steady state values in single point machining. The temperature fluctuates in intermitted operations depending on the tool-workpiece engagement periods hence cutting process is mostly in the transient phase. The transient, peak and periodic amplitudes of temperature need to be modeled to predict the cyclic thermal loading of the tool. The peak temperature must not exceed the diffusion limits of binding material to avoid premature chipping and accelerated wear of the tool. The cyclic thermal loading analysis is needed to assess the fatigue life of the tool and residual stress left on the cut surface. This thesis proposes two and three-dimensional thermal models of the machining operations for a given tool-workpiece material properties (thermal and mechanical), 3 tool geometry and cutting conditions such as speed and feed. The developed model has possible applications in different cutting processes by adopting geometry discretization specific to the process. It accepts variable and constant chip loads that may occur in single and multi-point cutting operations. The development of tool and coating materials require the information of peak and transient thermal loads. Diffusion between tool and workpiece at the tool-chip interface is one of the primary cause of wear mechanisms that is dependent on temperature. Cutting tool manufacturers use cobalt as a binder for tungsten-carbide tools, and it diffuses to workpiece material at elevated temperatures. Another form of common wear mechanism is thermal softening of the tool material; it increases susceptibility to deformation under machining loads. Cutting edge deformation, due to higher temperatures at the tool-chip contact region, leads to accelerated flank wear. In summary, the prediction of temperature in the tool and chip is required to plan the cutting conditions to avoid violating the thermal softening and diffusion limits of tool materials. The tool surface is coated with temperature resistant material to prolong tool life by thermally insulating the tool substrate. Coating layer reduces the heat transfer to the tool substrate and facilitates the removal of the heat with the chip from the cutting zone. Coatings can be applied as a single or multi-layer fashion depending on the workpiece material and process requirements. Multi-layer coatings can benefit from combination of different material properties. Some of the layers can keep their hardness at elevated temperatures; some can have increased thermal resistance during the cutting process. The prediction of temperature can benefit the tool designers to select the most suitable coating selection strategy for each machining application. This thesis provides a thermal model that is capable of modeling 4 anisotropic material properties within the modeled domain. The developed model is able to simulate multi-layer coated tools and varying material properties within tool, chip, and workpiece. Metal cutting is a high strain rate deformation process. Penetration of the tool into workpiece leads to material removal due to high amount of plastic deformation that takes place in the primary zone. Removed material travels over the rake face of the cutting tool and slightly deforms due to frictional effects in the secondary deformation zone. Rubbing between tool flank and machined surface creates some deformation and heat on the workpiece surface in the tertiary deformation zone. Mechanical and thermal loads are transferred to the newly generated surface via tertiary and primary deformation zones. This thesis models the heat generated in all three deformation zones. The temperature distribution is predicted by considering the heat transfer between the workpiece material, cutting tool and moving chip, and air. The thesis presents a hybrid model of the temperature distribution in the chip, tool, and workpiece. Principles of metal cutting mechanics are used to model the heat generated in the three deformation zones. A numerical solution is carried out to predict the thermal behavior of tool, chip, and workpiece by using finite difference method. Hybrid modeling allows solving temperature field without calculating material deformation at high strain rate, which has computational time advantage over integrated plastic deformation simulation approaches used in Finite Element method. The thesis is organized as follows: The existing thermal modeling approaches are reviewed in Chapter 2. The modeling of heat sources and transfer are presented in Chapter 3. Deformation zones, where the heat is generated, are modeled based on the mechanistic modeling of cutting process by considering the tool geometry, workpiece material, and cutting 5 conditions. Heat transfer equations are written based on general transient conduction with a possible mass transfer in three Cartesian directions. The heat balance equations are discretized using the finite difference method (FDM) with a second-order implicit time discretization. Chapter 4 presents the application of the heat transfer model to tool, chip and workpiece domains. Cartesian coordinate based heat transfer equations are converted to curvilinear coordinates to adapt the solution to arbitrary tool edge shapes. The derivation procedure of finite difference approximations is explained. Application of discrete heat balance equations to tool, chip and workpiece domains are presented. Chapter 5 consists of validation simulations of constant, temperature dependent thermal properties and tool coating application, and comparison against the experimental data found in the literature. In addition, experiments conducted with tool embedded thermocouples are presented and compared against the prediction results obtained from the proposed simulation model. The thesis is concluded in Chapter 6 by summarizing the contributions of the research and suggested future research directions. 6 Chapter 2: Review of Temperature Prediction Literature in Machining 2.1 Overview There are various thermal modeling approaches available in the machining literature that can be grouped under analytical and numerical approaches. While the analytical temperature models are computationally less intensive, they are limited to consider oblique cutting geometries and time-varying, interrupted milling operations. Although numerical methods can handle different cutting geometries and time variation of the process, they are computationally time-consuming. Furthermore, these methods do not allow for a clear understanding of the interrelationship between the cutting parameters and the resulting temperature changes in the tool and workpiece. This chapter gives an overview of the available analytical and numerical thermal models in the literature. The analytical, finite element and the finite difference methods based temperature prediction models are reviewed as follows. 2.2 Analytical temperature prediction models Analytical method’ s are mostly based on moving heat source solutions derived from Jaeger’ s [2] classical work. Loewen and Shaw published a work where they modeled orthogonal cutting [3]. They combined two moving band heat source solutions for primary deformation zone shearing as a uniform heat source moving with the shear velocity. The second source is approximated as a moving band heat source for the tool-chip interface. They assumed that chip boundaries are adiabatic. Therefore, the effect of tool-chip interface temperature on shear plane is not included as a heat source. Trigger and Chao later published [4] an analytical model for orthogonal cutting that considers both primary and secondary zone heat sources to calculate average tool-chip interface temperature. Venuvinod and Lau presented a model for oblique cutting which considers a heat source that moves with an orientation with respect to the 7 cutting edge [5]. Their work was one of the first that addressed the oblique chip movement. They determined rake face temperature distribution with their model. Venuvinod and Lau’ s model was the first one that gave two-dimensional temperature distributions on the rake face [5]. Stephenson and Ali gave the first example for transient temperature distribution on the rake face of the tool analytically [6]. They also integrated the instantaneous effect of point heat sources at given time steps to build the temperature field by using a steady-state solution [7]. They assumed a heat source distribution to reduce the computational load related to the convolution of the non-uniform heat source. The amount of frictional heat that goes into the tool is determined by using the approximation in [3]. Since their model worked for transient conditions, they were able to predict interrupted turning temperatures. Later Stephenson et al. extended this model for tool temperatures in contour turning operation [8]. They used a finite domain with a series solution as opposed to semi-infinite domain in the previous work [6]. The developed model did require experimentally determined contact area and frictional heat flux. Radulescu and Kapoor predicted tool temperature for continuous and interrupted cutting with an analytical model [9]. They used a three-dimensional transient heat conduction model. Quasi-steady state assumption was made to determine heat fluxes in the chip formation zone. Heat fluxes were calculated using a mechanistic approach. They forced heat balance around the chip formation zone through iterative solutions by considering boundary interaction between tool and workpiece. Komanduri and Hou proposed a comprehensive analytical model that predicts tool-chip and workpiece temperatures analytically for orthogonal cutting. They used non-uniform heat partition distribution between tool and chip and applied frictional heat generated based on 8 Jaeger’ s solution [2]. The shear plane heat source is modeled as an internal moving line heat source in a semi-infinite medium. Individual [10, 11] and combined [12] effects of these two heat sources are available in their work. Although semi-infinite medium assumption for tool and workpiece tend to result in low temperatures away from the heat sources, their work was significant regarding detailed temperature distribution it provides for the modeled domains. Other researchers used this reference work to build enhanced versions. Huang and Liang added a flank wear land heat source to Komanduri and Hou’ s work [13]. Li and Liang added tool flank convection effects to their model [14]. Karpat and Ozel applied modified Komanduri’s solution to include chamfered tool geometries [15]. Chou and Song extended it and introduced multiple rectangular heat sources to represent the complex geometry in turning with a wear-land heat source [16]. Tools with coating modeled analytically in literature as well by different researchers. Although analytical solutions are not valid for anisotropic mediums, material property changes were treated via the equivalent conductivity approach by Grzesik and Nieslony [17]. They determined the average and maximum rake face temperatures based on semi-analytical formulas of heat partition and temperature rise. Zhang and Li modeled monolayer coated tool temperatures by considering one-dimensional heat transfer [18]. Milling processes have been modeled to a limited degree with analytical methods. Generally, the transient behavior is solved via assuming quasi-steady-state conditions to determine heat partition and heat flux towards tool. Sato et al. [19] extended Stephenson and Ali’ s [6] approach to end milling temperature modeling. They considered the changing chip load and contact length via introducing time-dependent plane heat sources. However, heat flux and heat partition were approximated based on the assumptions in [6]. Richardson et al. 9 simulated the peripheral milling process by using a method of moving heat sources, and the heat flux value for the workpiece was determined via the temperature measurements from the workpiece [20]. Liu et al. modeled workpiece temperature in helical end milling by introducing heat sources along the flank and bottom cutting edges [21]. They used a point heat source Green function solution in their calculations. 2.3 Modeling of temperature in cutting with numerical models 2.3.1 Finite element methods used in cutting temperature prediction Finite Element Models have been used in temperature models since early 1970s. Initial finite element models were not entirely predictive. Tay et al. were one of the first who modeled the tool, chip, and workpiece steady-state temperatures in orthogonal cutting by the deformation field data gathered from quick stop experiments [22]. Later, Tay et al. updated this model with mechanistic model inputs to decrease the experimental input dependency [23]. Muraka et al. [24] added the tool flank heat source to Tay et al.’ s model [22]. They also included convection cooling effects in their analysis. The developed model required a substantial amount of experimental data to calculate deformation in primary and secondary zones. Dawson and Malkin modeled inclined moving heat source over workpiece [25]. Their approach included heat loss due to chip removal and convection cooling. The main aim was to determine workpiece temperatures for grinding and orthogonal cutting. It was a Eulerian solution that had fixed mesh and approximate shear plane geometry. Childs et al. modeled the two-dimensional tool, chip, and workpiece heat conduction for steady state conditions via finite element method [26]. They modeled deformation zone heat generation using mechanistic approximations. They investigated the effect of coolant on tool-chip interface temperatures within limits of numerical 10 experimentation. Later, Ueda et al. used similar finite element simulations in their experimental temperature measurements with diamond tools [27]. The predictive capacity of finite element models is increased by inclusion of thermomechanical formulations. Strenkowski and Moon developed a Eulerian finite element model for the orthogonal cutting process [28]. Although their model was a thermomechanical model that solves the deformation field, chip formation was not predicted due to Eulerian assumption. Chip geometry was approximated as an input to the model. Their model assumed sliding friction only, and the tertiary zone was neglected in the analysis. Lin and Lin proposed a model where temperature distributions were modeled by employing a finite difference approach, and the chip formation process was simulated by using a finite element method for orthogonal machining [29]. Shih developed an FEM model that uses elastic viscoplastic material behavior with sticking sliding friction on the rake face [30]. The model included a relaxation solution after the loading that gave residual stress distributions as well. Wu et al. proposed a thermomechanical model for orthogonal cutting [31]. They developed a new sticking sliding friction model that uses updated flow stresses along the tool-chip interface. Constant thermal properties were used in their temperature calculations. Marusich and Oritz modeled orthogonal cutting based on Lagrangian formulation [32]. Their model created the base of commercial machining finite element code AdvantEdge®. Movaheddy et al. [33] developed an updated Arbitrary Lagrangian-Eulerian finite element model for chamfered and blunt tools based on their previous work [34]. Their model solved the thermal problem on the Lagrangian formulation to include the mass transfer effects. Commercial finite element codes that can handle large deformations became widely used after late 1990s. They are used to simulate both continuous and interrupted cutting operations. 11 One of the main drawbacks of these combined solvers is their limitation to model tool-chip interface interaction. The results are sensitive to friction modeling. Therefore, the friction model and flow stress approximation for the secondary zone should be checked carefully. Using cutting tests for calibrating the flow stress prediction law is instrumental in reaching this goal. Ng et al. used FORGE2® finite element code to model orthogonal cutting [35]. A thermomechanical simulation approach was used with sticking-sliding contact at tool-chip interface. Two-dimensional heat transfer calculations were performed by determining heat partition via static contact assumption between tool and chip. They used thermal effusivities to calculate heat partition based on static contact assumption. Static contact assumption loses reliability when cutting speed increases due to increased mass transfer via chip removal. Ozel et al. used DEFORM2D® to simulate end milling [36]. They used Zorev’s friction model in the analysis. Moreover, temperature values for tool, chip, and workpiece were predicted. Soo et al. used Abaqus® Explicit to simulate ball end milling process [37]. They were able to simulate initial engagement only due to computational load limitations. Full engagement and chip load were not simulated in their work. Pittala and Monno modeled the temperatures in milling by employing two-dimensional analyses via DEFORM2D® and 3D [38]. First, they calculated the cutting temperatures using two-dimensional calculations. Then, they used that results as inputs to three-dimensional workpiece temperature simulations to decrease the computational burden. Tools with coating were also modeled in the past by different researchers. Yen et al. modeled orthogonal cutting of AISI 1045 with multi-layer coated tool via using DEFORM2D® [39]. They used shear strength based friction model with a shear friction factor. Friction factor was arbitrarily tuned by using experimental results. Arrazola et al. simulated orthogonal cutting of AISI 4140 with multi-layer coated tools by using Abaqus® Explicit [40]. They used infra-red 12 camera to measure cutting temperatures. Nemetz et al. followed a hybrid approach to simulate the 3D temperature field in milling [41]. They used two dimensional Abaqus® Explicit simulations (with deformation simulations) as an input to three-dimensional heat conduction only tool model. Thakare and Nordgren used cutting simulations of AdvantEdge® as an input for steady-state heat transfer analysis in MSC Marc® [42]. Two-dimensional models were run for AISI 4340 cutting with Ti (C, N) coated carbide tool. The reason behind separate analysis was the computational cost of deformation solutions in FEM. Wu and Zhang performed three dimensional Abaqus® simulations for end milling of Ti6Al4V with TiAlN coated carbide helical end-mill [43]. Four milliseconds long simulations took 180 hours according to reported results. Their force predictions were validated with experiments; however, temperature validation was not present in their work. Heat conduction analysis without solving the material deformation is another common finite element modeling approach. Model input can be gathered from experimental measurements or available deformation zone models in the literature. Anogonye and Stephenson modeled three-dimensional steady-state tool an tool holder temperatures by using Nastran® solver [44]. They derived correction factor model from FEM that was applied on Loewen and Shaw’ s [3] analytical solution to account for insert shape and nose radius. Dessoly et al. modeled tool temperatures using FEMLAB® (previous naming of COMSOL®) for rotary tool turning via the use of mechanistic analysis of deformation zones [45]. They used Shaw’s secondary zone heat partition calculation method in their analysis [46]. Mamedov and Lazoglu modeled micro-milling temperatures with COMSOL® [47]. The modeled deformation zone heat generation from the mechanistic analysis. Their model did not include calculation of heat partition towards tool. They did assign a fixed heat partition ratio before running their model. 13 2.3.2 Finite difference methods used in cutting temperature prediction Various researchers used the Finite difference method for solutions of partial differential equations in literature since early 1950s. The mathematical foundations of the method go back to 18th century. Use of finite difference method for metal cutting temperature distributions is less widely found compared to finite element method. Rapier proposed one of the first models that combine the finite difference method and analytical solution for orthogonal cutting temperature distributions [48]. Chip temperature field was modeled analytically, whereas, they determined tool and workpiece thermal fields with finite difference approximation. Rapier neglected heat conduction in chip flow direction by assuming that mass transfer effects are dominating the heat transfer. The Shear plane was assumed to have a uniform temperature distribution with an analytical solution of line heat source. The proposed model neglected the heat partition towards the tool. Dutt and Brewer proposed a finite difference method based thermal model for steady state orthogonal cutting [49]. They determined temperature distributions using finite difference method. Heat transfer between tool and chip was determined as opposed to Rapier [48]. Tool and chip thermal conduction was taken into account. Similarly, they considered the thermal interaction between chip and workpiece over the shear plane. Levy et al. proposed a two-dimensional transient thermal model that considers only tool and chip temperature distributions in orthogonal cutting [50]. Their model had chip velocity limit due to numerical solution stability. The primary deformation zone was modeled as internal volumetric heat generation, and secondary deformation zone heat flux was approximated via non-uniform distribution along the tool-chip contact length. 14 Usui et al. applied Bishop’s [51] “jerk” method to calculate temperature fields in turning [52] for chip and tool. In their assumption, the heat generation and transport occur instantaneously. Conduction occurs subsequently in a discrete time interval. They included the primary and secondary zone heat sources. Tool-chip interface heat balance was determined by assuming heat flux equivalence on coincident tool and chip surface locations. Smith and Armarego modeled steady-state three-dimensional heat transfer for chip and tool based on thin shear zone analysis in orthogonal cutting [53]. They assumed uniform temperature on shear plane which is identified from Chao and Trigger’s [4] average temperature rise relationship. The secondary zone heat source was assumed to be uniform in the tool-chip contact area. Lazoglu and Altintas developed a two-dimensional model of the tool and chip for orthogonal cutting and interrupted milling [54]. The primary zone heat source was included by introducing Oxley’s [55] energy partition function. Uniform strength heat source with a sliding contact assumption was employed to model the secondary zone heat generation. A quasi-steady state approximation was made to model the transient temperature of interrupted machining via considering the first order system dynamics. Later, Ulutan et al. [56] extended this model to 3D analysis with convection heat losses. Grzesik et al. developed a two-dimensional Cartesian coordinate model for coated and uncoated tools in orthogonal cutting [57]. They used an equivalent thermal conductivity formula to include multi-layer tool coating. Heat conduction along chip flow direction was neglected in their model due to chip flow velocity’s dominant effect over the heat conductance. Secondary zone thermal interaction was set by equating heat fluxes on tool and chip interface. 15 Lazoglu and Islam proposed a 3D temperature model for continuous oblique cutting operations [58]. This solution was for steady-state analysis without considering the transient effects. The developed model was capable of simulating tool, chip, and workpiece by considering primary and secondary zone heat sources. 2.4 Summary The heat transfer in metal cutting has been studied by the researchers extensively. The previous subsections presents the available thermal modeling approaches in machining. The limitations of prior work are stated here to clarify the benefits of the proposed work over the existing methods. One disadvantage related to the transient analytical solution is that heat flux towards tool assumed steady-state in them. This thesis provides an iterative heat partition determination through numerical solutions. Therefore, heat partition towards tool can be determined dynamically. Anisotropic medium modeling capability of analytical methods is limited due to constant material property assumption in available Green’s function or series expansion based solutions. Since a numerical model is used in this thesis, anisotropic mediums such as coated tools can be considered with the proposed approach. Introduction of multiple heat sources and boundary conditions is restricted due to the nature of analytical solutions. This thesis follows a Finite Difference Method based solution that does not have such limitation. This thesis proposes a finite difference method based heat transfer analysis model without solving the material deformation during metal cutting. The main advantage of the developed model over the existing finite element based deformation analysis is the computational time advantage. Since developed model uses thin shear zone analysis, secondary zone interaction defined based on experimental measurements which is an advantage over finite 16 element method based models. Another advantage of the developed model over heat transfer analysis with finite element method approaches is the heat partition calculation method. Heat partition is iteratively calculated based on the sliding contact assumption in the developed method. Most of the finite difference models (FDM) that are available in the literature are developed for a fixed coordinate system (i.e., Cartesian coordinates, polar coordinates) and orthogonal cutting. A three-dimensional transient temperature model for continuous and interrupted oblique cutting operations for coated and uncoated tools has not been developed using analytical or finite difference methods in the past. Since the developed model uses a curvilinear coordinate based formulation, different tool geometries can be represented. 17 Chapter 3: Heat Transfer Model of Metal Cutting 3.1 Overview Transient heat balance model for a general case of moving medium is introduced in this section. Heat balance equations are written based on both constant and varying thermal properties. Depending on the preference, the developed engine can switch between both options. First, the heat source types are presented. Heat generation in the primary, secondary, and tertiary zones are explained. Then, the heat balance equation and the discretization of governing equations with finite difference method in a Cartesian domain are shown to introduce the proposed numerical method applied on a sample case. 3.2 Modeling of heat sources in metal cutting The metal cutting process is shown in Figure 1-1, where the metal is plastically deformed with high strain and strain rates at the primary deformation zone. The deformed chip moves on the rake face under sticking and sliding friction conditions which constitute secondary deformation zone. The clearance face of the tool rubs against the finished surface in the tertiary deformation zone. Heat is generated and results in temperature gradients in all three zones. Heat generation in deformation zones depends on cutting tool geometry, workpiece material properties, and cutting conditions. In this thesis, the heat generation is estimated based on the orthogonal cutting theory with thin shear zone assumption. A finite volume method based approach is taken in modeling. Time-dependent behavior is treated using an implicit solution procedure due to mathematical restrictions. Model is developed to use multiple heat sources like internal heat generation, and surface heat fluxes. The heat sources in the primary shear and secondary deformation zones are considered in the thermal model. Since the tertiary zone contributes a negligible amount to heat generation, especially for 18 fresh tools, it is not taken into account in primary validated cases within this thesis although the developed model allows its inclusion. 3.2.1 Heat Source in the primary deformation zone The chip is formed in the primary shear zone as a result of high plastic shear deformation. Shearing power is modeled by assuming that strain and strain rate distributions are uniform along the shear plane which has a shear angle ( n ) with respect to cutting velocity direction (Figure 1-1). Shear angle is defined on the cutting edge normal plane ( nP ). Heat conduction in this deformation zone is neglected during the very first stages of chip formation due to high deformation rates. Hence, the heat generation is assumed to be uniform throughout the shear zone when chip starts to form. Shearing power is defined as follows; s s sP F V (3.1) where ,s sP F , and sV are shearing power, shearing force, and shearing velocity, respectively. Cutting velocity vector and tool cutting edge normal is oriented to each other with an inclination angle ( s ) in oblique cutting as shown in Figure 3-1. This generic case becomes orthogonal cutting when the inclination angle is zero. Cutting forces and velocity vectors are shown in Figure 3-2. Detailed information about oblique cutting mechanics is given in Appendix A . 19 Figure 3-1 – Oblique cutting geometry Figure 3-2 - Cutting force and velocity vectors in oblique cutting 20 Shearing force can be expressed by multiplying average shear flow stress ( s ) on the shear plane and shear plane area ( sA ) as, s s sF A (3.2) nnsin cossinsssb hhlA (3.3) where sl , h , b , n , and s are shear plane length, depth of cut, width of cut, normal shear angle, and inclination angle, respectively. Shearing velocity can be determined by using the vector relationship between the chip velocity and cutting velocity (see Appendix A ) as; nn n ccos cos sec( )(sin cos sec sin tan(η ))0s n ns s s nVV V (3.4) The magnitude of shear velocity vector is determined as; 2 2n n n ccos cos sec( ) (sin cos sec sin tan(η ))s s n n s s nV V V (3.5) Shearing power is introduced with two different options. It can be introduced as a uniform volumetric heat source by assigning a shear zone thickness as the first option. Average uniform temperature can be assigned at the shear zone as the second option. It is assumed that chip leaves the shear zone with average uniform temperature and then the secondary zone heat source further heats it as it moves over the rake face. Volumetric heat generation in the primary shear zone is expressed as; 21 tsss sPQA (3.6) where, sQ , sA , sP , and st are volumetric heat generation rate, shear plane area, shearing power and shear plane thickness respectively. Shearing power generates steady-state temperature that is approximated as; sss rp cPT TV C h (3.7) where, rT , , , pC , and ch are room temperature, Taylor-Quinney coefficient (ratio of shear power that is converted to heat), workpiece density, heat capacity, and cut chip thickness, respectively. Transient temperature is expressed as; 2 22 211tn tshear ssc sT t T eh t (3.8) where t is the time step size and is the time constant evaluated from the workpiece material’s thermal diffusivity constant ( ), the cut chip thickness ( ch ) and the shear zone thickness which is approximated as 0.15s ct h [59]. Since it is assumed that shearing happens in a thin zone, time constant of the shear zone is analytically determined similarly to a parallelogram with adiabatic boundaries. Note that the width of cut is neglected in time constant approximation. 3.2.2 Heat source in the secondary deformation zone Secondary deformation zone is modeled by using both sticking and sliding zone assumptions. Although the sliding zone assumption has been used for the most part in 22 validations, sticking zone assumption has been used within the context of checking model sensitivity to different input parameters. Uniform pressure distribution along the tool-chip contact length is assumed in sliding contact analysis. Force component along the chip flow direction ( uF ) is multiplied by the chip velocity ( cV ), and it is written as the power generated in the secondary deformation zone; nnc n.cos cos cos sin sin sin cosV cos sin Vcos η coss af u c ci i n n i i sscnb h sinP F V V (3.9) where, fP , uF , and cV are friction power, friction force along the chip flow direction and the chip velocity, respectively. Tool-chip contact length has been determined by using moment balance of frictional normal and shear normal forces with respect to cutting edge (see Appendix A ) as; ( ) cos (cos cos sin ( ))n ncc nn n n nh sinlsin sin (3.10) Uniform heat flux is applied to the tool-chip interface by dividing frictional power input to the contact area as; ,f cfc cc c c c cP NqA Nwhere A N b N l N (3.11) where, cA , cN are tool-chip contact area and discrete chip substructure number, respectively. Sticking zone assumption is made based on constant shear stress distribution in the sticking region. It is assumed that the shear stress has a constant value in sticking zone and then 23 decreases linearly to zero at the end of the tool-chip contact length (See Figure 3-3). This assumption changes the pressure distribution to non-uniform from a uniform distribution in sliding contact at the tool-chip interface. Hence, the energy dissipation in these regions are different. Figure 3-3- Assumed shear stress distribution in the secondary deformation zone Constant deformation zone thickness is assumed in sticking zone. Amount of deformation power is predicted by using shear strain rate and average shear flow stress as; secsticking av pl ondaryP b l t (3.12) where, stickingP , av , and are sticking power, average shear stress, and shear strain rate, respectively. Shear strain rate is determined by using average chip flow velocity and deformation zone thickness based on the constant shear strain rate assumption within the secondary deformation zone. secondarycVt (3.13) 24 Average shear flow stress ( av ) is determined by using the friction ( uF ) force balance distributed in sticking and sliding contact zones on the rake face as, 0bpl cpll lsticking slidingu u u av slidinglF F F dx b dx (3.14) where pll is sticking contact length, b is the width of cut, and sliding is the average shear stress in the sliding region (see Figure 3-3) of the contact. A ratio of plastic contact length ( pll ) to overall contact length ( cl ) is assumed (i.e., 0.8pl cl l , 0.5pl cl l etc.) in order to determine the average shear stress ( av ) in Figure 3-3 from the force balance given in equation (3.14). Friction force ( uF ) in equation (3.14) determined from the thin shear zone model by using the orthogonal cutting data base. Note that sticking contact analysis only used to check the model sensitivity. Contact length is determined in the same way as the sliding contact assumption via moment balance around the cutting edge. Equation (3.10) applies to fully sticking contact case without modification. Sticking-sliding contact assumption requires derivation of contact length due to change in moment arm length for friction normal force. Contact length is derived by using moment equilibrium between friction normal force ( vF ) and shear normal force ( nF ) with respect to cutting edge. Figure 3-2 shows the force components necessary for the derivation procedure. For example, if half of the total contact length is sticking contact. Moment arm length for friction normal force on rake face ( vF ) becomes 14.cos36ccl , and tool-chip contact length becomes; 25 14cos364 F2 2 223slidingc av cnv c nuu avccfrictil lF b bav l bstion shearhnormal force nockirmal forcesing slidingregion regionforcnF Fmmoment arm forfrictie forcon noremal forcel 18 sec ( ) 14 (cos cos sin ( )),cos ( )cos cos cos sin sin sin coscos cos cosc n ncn n n ns i n nni i n n i i n sns i nvoment arm forshear normal forceh sinlsin sinWhereb h sinFb hF sin sincos cos cos sin sin sin cosnn ni i n n i i n s Shear stress ( s ) in the friction force calculatios is calibrated at different cutting speeds for varying feed rates and cutting tool geometry. Therefore, orthogonal cutting coefficients indirectly includes temperature dependency of average shear stress values with a speed dependent calibration procedure followed. Power generated in the sliding zone is determined via multiplying the sliding contact portion of friction force with chip velocity as; cpllslidingsliding slidingf clu cbP V VdxF (3.15) 26 3.2.3 Heat source in the tertiary deformation zone Tertiary zone heat source occurs as a result of rubbing between the tool flank and newly generated workpiece surface (see Figure 1-1). It is applied within the context of sensitivity analysis in this thesis. Energy generation due to rubbing is directly predicted by using tangential edge force ( teF ) and cutting speed ( V ). It is applied over a distance equal to cutting edge radius ( er ). Total ploughing energy generation ploughingP and related heat flux ( ploughingq ) due to rubbing are approximated as; " " " "1000teploughing teploughingploughingeK b VP F V divisionis made for m to mm conversionPqr b (3.16) where teK is tangential edge force coefficient (see [59]). 3.3 Heat balance equations 3.3.1 Constant thermal property heat balance equations Heat balance equations are written based on the first law of thermodynamics at rate basis (equation (3.17)) for each control volume. Input heat rate values in x, y, and z orthogonal directions are represented by using 𝑞𝑥, 𝑞𝑦, 𝑞𝑧 terms. Output heat rate values at facing surfaces are represented by terms 𝑞𝑥+𝑑𝑥, 𝑞𝑦+𝑑𝑦, 𝑞𝑧+𝑑𝑧. Schematic of a control volume is shown in Figure 3-4. 27 Figure 3-4- Schematic of a control volume in cartesian coordinates The heat stored in a control volume is expressed as; in generated out storedE E E E (3.17) Taylor series expansion is used to define output heat flux terms without including second and higher order terms. xx dx xyy dy yzz dz zqq q dxxqq q dyyqq q dzz (3.18) Volumetric heat generation rate (?̇?𝑔𝑒𝑛𝑒𝑟𝑎𝑡𝑒𝑑) is expressed as; generated genE q dx dy dz (3.19) Stored energy (?̇?𝑠𝑡𝑜𝑟𝑒𝑑) term is expanded by assuming that there is no phase change in solid medium. After combining expanded expressions, the first law of thermodynamics leads to; 28 x y z x dx y dy z dz gen pyieldsyx zgen pTq q q q q q q dx dy dz C dx dy dztqq q Tdx dA dy dA dz dA q dx dy dz C dx dy dzx y z t (3.20) The equation (3.20) can be rewritten by introducing Fourier’s law for heat transfer rate expressions for an internal control volume. Note that heat conduction is assumed constant and equal in all directions. Heat conduction coefficient change with temperature is neglected. xyzTq kxTq kyTq kz (3.21) 2 2 22 2 2gen pq CT T T Tx y z k k t (3.22) A general form of the heat balance equation is written by assuming constant moving velocity as, 2 2 22 2 21 1 1(x, y, z, t), ,V V, ,Vgen p p p pxx y zx y zy zT T T x T y T zx x V t y y V t z z V tt t x t y t z tx yq C C C CT T T T T T Tx y z k k t k x k yzV V Vt tztk (3.23) Note that only one of these velocity components is activated for a given problem and others are kept zero. 29 3.3.2 Variable thermal property heat balance equations Use of varying thermal properties is important to increase thermal model accuracy through incorporating temperature dependence compared to constant thermal properties use. Another advantage is that it lets modeling composite mediums such as tools with the coating. The equation (3.20) is rewritten considering changing material properties. For a given anisotropic medium, the most generic type of heat conduction scenario will have a three by three heat conductance matrix. This matrix involves three Cartesian ( , ,xx yy and zz ) and cross-conduction terms , , , , ,xy xz yx yz zx and zy (equation (3.24)). Since a fully anisotropic material is not modeled within the scope of this thesis, cross-conduction terms are neglected. The whole conduction matrix reduced to a state similar to cubical structure (i.e., Orthorhombic crystal structure). Final heat balance equation (3.26) is written based on this assumption. For a generic anisotropic medium, Fourier conduction along three different Cartesian directions can be written as; , 00 00 00 0x xx xy xzy yx yy yzzx zy zzzxx xy xzyx yy yz xy yx xz zx zy yzzx zy zzxxyieldsyyzzTxq K K KTq K K KyK K KqTzK K KK K K K K K K K K KK K KKK KK (3.24) Which leads to the equation (3.25), 30 2 2 2 2 2 22 2 22 2 22 2 2xx yy zz xy yx xz zx yz zy gen preduces toxx yy zz gen pT T T T T T TK K K K K K K K K q Cx y z x y x z y z tT T T TK K K q Cx y z t (3.25) The equation (3.25) becomes final heat balance equation (3.26) by including moving medium terms, 2 2 22 2 2V V Vp x p y pxx yy zz gen p zT T TC C Cx yT T T TK K K q Cx y z t z (3.26) where all material properties are defined as a function of location and temperature, (x, y, z,T)K , (x, y, z,T)pC , and (x, y, z,T) . The equation (3.26) becomes the general heat balance equation which also encompasses constant thermal property case. The equation (3.23) can be derived by setting xx yy zzK K K k . 3.3.3 Discretization of heat balance equations Heat balance equations are written for the meshed tool, chip and workpiece geometries in a discrete form using finite volume method. Each geometry is discretized into finite number of non-overlapping volume elements around a center node , ,c c cx y z (Figure 3-5). Then, the resulting equation system is solved. Finite volume method is selected to discretize the problem due to its simplicity in discretization. Furthermore, the implementation of partial differential equations in thermal modeling is relatively easier compared to the finite element method. Since it is a control volume method, it has an integral formulation where overall energy balance in heat transfer problem is forced through volume integral. After the discretization and writing the volume integrals, the differential operators are approximated with piece-wise constant material property assumption. 31 Figure 3-5- Sample volume element in discretized chip domain Generated control volume size may change depending on the position of the nodes within the discretized domain (corner node, internal node, surface node, and borderline node). However, the interrelation between neighboring control volumes is forced through the resulting finite difference approximations. The equation (3.26) is written in an integral form to explain how the method is applied in a Cartesian domain in the equation (3.27). Integral limits are written around a central node , ,c c cx y z with halfway differential distances in x, y and z directions from neighboring nodes ,2 2 2x y zand as limits around the node (see Figure 3-5). Although curvilinear coordinate is used later in the thesis, Cartesian coordinates have been used to explain the application. 32 2 2 2 2 2 22 2 2 2 2 22 22 2 2c c c c c cc c c c c cc c cc c cdx dxz y x z y xt t t tt tz y x z y xdxz y xzdz dy dz dyxx yydz dy dx dz dy dxd dyzzdy xz dy dxT TK dxdydzdt K dxdydzdtx x y yTK dxdydzdtz z 2 2 2 22 2 22 2 2 2 2 22 2 2 2 2 2c c cc c cc c c c c cc c c c c cdxz y xt t t tt tz y xdxz dz dygendz dy dxdz dxz y x z y xt t t tdy dz dyxdzt tz y x z y xdy dx dz dy dxpyq dxdydzdtT Tdxdydzdt V dxdydzdtt xCTV 2 2 2 2 2 22 2 2 2 2 2c c c c c cc c c c c cdz dx dxz y x z y xt t t tt tz y x zdy dz dyzdz dy dx dz dy dxy xTdxdydzdt V dxdydzdty z (3.27) Volume integral formulation (equation (3.27)) is separately applied to the chip, tool, and workpiece domains. The interaction between them is defined as boundary conditions. Integral limits are modified depending on the location of the volume element. A sample chip geometry is used to explain the application of finite volume method in Figure 3-5. Volume integrals are written by considering the nodal location for tool-chip interface node , ,m m mx y z in Figure 3-5. Then, they are expanded into discrete version equation (3.30) by using finite difference approximations. Note that for illustration to be clear; node location is showed on the boundary of Figure 3-5. Although formulation is written for an internal node on the tool-chip interface (bottom face in Figure 3-5) The secondary zone heat source term is considered in the equation (3.30) as a heat flux. Since integral in z direction is limited at the tool-chip interface, derivative term is not expanded into finite difference approximation at lower level of integral. This term corresponds to heat flux for the node. 33 1, ,k , ,k2, , , ,22 2 222 2Δ Δ t Δ Δ t, Δ1 Δm m mmm mm m m m m mdz dyzz z zzzzdy dx i j i jHeat FluxTermxdxz y xt tt zy xy z z x y zzzfzT T TK dxdydzdt K x y K x yz z z zT TK i j x yzq B , ,kt ,1yields zzi jfTKzq B i j The time-dependent term in the heat balance equation is treated with a backward finite difference approximation due to stability requirements. Resulting equations generally show a sparse distribution of negative real Eigenvalues which dictates the use of an implicit method. Although there are semi-explicit methods available for the problems with negative real Eigenvalues, they are suitable for closely spaced Eigenvalues. The implicit methods require solution of a linear set of discrete governing equations for all nodes at each time step. Therefore, implicit methods are more time consuming compared to explicit methods (i.e., forward finite difference approximation of time term). It is assumed that , ,ti j kT rises to , ,t ti j kT at the beginning of the present time step and remains constant throughout the time step as follows; , , , , , ,1t t ti j k i j k i j kT wT w T where w is the time weighting factor and it is assumed equal to one. 34 Figure 3-6 Implicit solution representation in time Second order backward implicit finite difference approximation is used to treat time-dependent differential operator as; 2., , , , , ,3 42t t t ti j k i j k i j kT T TTt t (3.28) 35 First order backward implicit finite difference approximation is used for the first time step to start solutions. In the second time step and subsequent time steps, second order backward implicit approximation (equation (3.28)) is used. First order approximation is expressed as; , , , ,t ti j k i j kT TTt t (3.29) , , , , ,, , , ,, , ,2 22., , , , , ,1 ;23Δ;42Δ2.m m m m m mm m m m m m m m mm m m m m m m m m m m mstt tx y z x y zndt t t tx y z x y z xx x y z x y z x y zy zx x y zx xxx xxFor time stepT TT T T T zK K yxTIMEtAfter time stepT T TTIxMEt , , , , , , , ,2 2, , , ,2, , ,ΔΔ2Δ, Δ Δ Δ. 1 Δ2ΔΔ Δ2m m m m m m m m m m m mm m m m m mm m m mx y y z x y z x y z x y y zy yyy yyx y z z x y zz genzzx x y z x xfpxT T T T zK K xy yT T zK i j x y q x yq BT TTIMzzx yCE V ,, , , , , , , ,ΔΔ222ΔΔ Δ2Δ2m mm m m m m m m m m m m my zx y y z x y y z x yy zz z x y zzyzxT T T Tx yV V (3.30) Since Taylor series expansion is used in the derivation of finite difference approximations, spatial discretization error becomes first order (proportional to , ,x y z ). Although spatial discretization error is a second order for internal nodes, boundary node approximations have forward and backward finite difference approximations (which is the cause of first-order spatial error). Therefore, the element sizes are selected small enough to limit the discretization error’s 36 effect on the solutions. Since time discretization is second order, time discretization error is proportional to 2t . Further information can be found in [60]-[61]. Thermal properties in the equation (3.30) can be used as constant or varying depending on the preference. Conduction coefficient from one element to another is needed to be determined when using varying thermal properties in the simulations. Halfway terms that occur in equation (3.30) denote (i.e.2 2, ... .x yxx yyK K etc ) element to element heat conduction properties and they are determined via using heat flux equilibrium at two neighboring elements interface. Figure 3-7- Element to element varying conduction coefficient demonstration Assuming that element to element interface area does not change, heat conduction values from one element to another are written as (see Figure 3-7), 37 , , , , , , , ,'', , , , , ,2 2 2, ,2, ,, , , ,'' 2 2, , , ,, , , ,, ,2 2int ,20.5x y z x x y z x x y z x y zs x x xx y z x y z x y zxx y zx y z x x xx y z x y zs x y z x x x y z xx y z x y zx y zerface heat flux is written asT T d df K T A ddT T T Tf K T A K T Ad , ,2, ,, ,, ,2, , , ,, , , ,0.50.5 0.5xx y zx y zx x y zxx y zx y z x x y zx y z x x y zdleads todKd dK T K T (3.31) where, ,2xx y zA , , , ,, ,2, ,x y z x x x y zx y zd d d are element to element interface area and discretization sizes, respectively. The equation can be derived for other Cartesian directions in a similar manner. 38 , ,2, ,2 , , , ,, , , ,, ,2, ,2 , , , ,, , , ,, ,2, ,2 , , , ,, , , ,0.5 0.50.5 0.50.5 0.5xx y zxx y zx y z x x y zx y z x x y zyx y zyx y zx y z x y y zx y z x y y zyx y zyx y zx y z x y y zx y z x y y zdKd dK T K TdKd dK T K TdKd dK T K T , ,2, ,2 , , , ,, , , ,, ,2, ,2 , , , ,, , , ,0.5 0.50.5 0.5zx y zzx y zx y z x y z zx y z x y z zzx y zzx y zx y z x y z zx y z x y z zdKd dK T K TdKd dK T K T (3.32) Material density and heat capacity ( (x, y, z,T), (x, y, z,T)pC ) are taken as element centered average values during the solutions. 39 Chapter 4: Application of Heat Transfer Model to Metal Cutting 4.1 Overview This chapter explains the application of discretized heat balance equations in the tool, chip, and workpiece. First, curvilinear coordinate system transformation of the Cartesian coordinate frame is explained to create the computational domain. Then, the mesh generation method is introduced. Application of boundary conditions and setting the linear equation systems for each domain of interest is explained. Finally, the solution algorithm is given to represent the overall framework. The algorithm starts with the discretization of the tool, chip, and workpiece geometry to form substructures. Solutions are separately carried out for each substructure by considering boundary conditions between them. Then, the repetitive solutions are continued until a pre-set temperature convergence criterion is met for all discretized geometries. Detailed explanations are given in subsequent subsections. 4.2 Computational domain transformation Finite difference approximations require a rectangular computational domain. The derivation of finite volume formulations is presented in section 3.3.3. Cutting process geometries involve non-rectangular physical domains. The derivation of finite volume relations is cumbersome for non-rectangular domains due to unequally spaced grid distribution and non-uniform boundaries. Therefore, governing equation system in Cartesian coordinates is transformed into a rectangular computational domain by introducing curvilinear coordinates. A new coordinate frame with the following coordinate directions is defined. The generation of new grid points is explained in the grid generation section in detail. ( , , )( , , )( , , )x y zn n x y zw w x y z 40 Figure 4-1- Physical to computational space transformation Differential equation type remains unchanged after the coordinate system transformation. Thus, the second-order backward difference implicit solution scheme stays valid due to similar negative real Eigenvalues of the governing equations. Heat balance equations in Cartesian coordinates (equations (3.26) and (3.23)) are rewritten in computational space by using the following differential operators (equations from (4.1) to (4.7)). 2 2 22 2 22 2 22 2 22 2 22 2 2, , , , , , , , , , ,, , ,, , , x y z xx yy zzx y z xx yy zzx y z xx yy zzx y z x y zn n n n n nn n n n n nx y z x y zw w w w w ww w w w w wx y z x y z 2 2 2 2 2 22 22 2 22222 2 2.x x x x x x x xxx xx xx xT T T T T Tn n w w w nx n n w w w nT T T Tn wn w (4.1) 2 2 2 2 2 22 22 2 22222 2 2y y y y y y y yyy yy yy yT T T T T Tn n w w w ny n n w w w nT T T Twn w (4.2) 41 2 2 2 2 2 22 22 2 22222 2 2z z z z z z z zzz zz zz zT T T T T Tn n w w w nz n n w w w nT T T Tn wn w (4.3) x x xT T T Tn wx n w (4.4) y y yT T T Tn wy n w (4.5) z z zT T T Tn wz n w (4.6) t t tT T T T Tn wt n w (4.7) The equation (3.26) for variable thermal properties takes the following form after substituting differential operators, 42 2 2 2 2 22 22 22222 2222 2 22222 2 2.22 2x x x x x x x xxxxx xx xx xy y yyy y y y y y yy yyyyT T T T Tn n w w w nn n w w w nKT T T Tn wn wT Tn n wn n wT T T T TK w w n nw w n nT Tww 222 2 22 22 22 2 22222 2yz z z zzzz z z z zz zz zz zx x xgen p t t t pT T Tn n w wn n w wKT T T T T Tw n n ww n n wT T Tn wn wT T T Tq C n w Cn w y y yz z zxtT T T yn wn w tT T T zn wn w t (4.8) Constant thermal property heat balance equation can be written by substituting xx yy zzK K K k into the equation (4.8). The relationship between two coordinate frames is defined by relative differential mapping terms ( , ,... ,x xx z zzw w ) in the equation (4.8). These coefficients are identified within the section 4.4 later in detail. Finite volume formulation in the Cartesian coordinate frame can be applied similarly to curvilinear coordinates. Following relationships are derived for a reciprocating node , ,c c cn w in generalized coordinates of an example internal node , ,c c cx y z in the Cartesian coordinates (See Figure 3-5). Equations (4.9), (4.10), (4.11) show this derived relationships for an internal node. Since the integral limits change when deriving finite volume 43 formulas, similar derivations are made for 27 different total nodal locations on a generic domain (8 different corner node types, 12 different borderline node types, 6 different surface node type and an internal node type). Note that aforementioned discretization error sources (see Section 3.3.3) are reflected into computational domain transfer as well. 2 2 2 2 22 22 22222 2222 2 2222 2 222 2cx x x x x x x xt Vxx xx xx xy y yyy y y y y y yy yyxx cT T T T Tn n w w w nn n w w w nT T T Tn wK dV dtn wT Tn n wn n wT T T T TK w w n nw w n nTw 2222 2 22 22 22 2 22222 2cc ct Vyy yz z z zzz gent V t Vz z z z zz zz zz zp t tcc cTwT T Tn n w wn n w wK qT T T T T Tw n n ww n n wTdV dtdV dt dV dTCntTn 2 2 22 2 2,c c ccc c cc cx x xt p y y yt V t Vz z zd dndw dnc cw dw nt VwdnT T T xn wn w tT T T T yw C n ww n w tT T T zndV dt dV dtwwn w there 1t tctand dV d dndwJ (4.9) Assuming that domain grid point distribution does not change with time 0t t tn w , the discrete form can be written as; 44 2 2 2 2 2 21 1 1 1 1 12 2 2 2 2 22 2 2 2 2 21 1 1 1 1 12 2 2 2 2 2x y z x y zxx yy zz xx yy zzx y z x y zxx yy zz xx yy zzK K K w n K K K w nn K n K n K w n K n K n K wn n 2 2 2 2 2 21 1 1 1 1 12 2 2 2 2 22 21 12 2212, ,2 22c c cx y z x y zxx yy zz xx yy zzx xx y yyxx yyz zzzzT n ww K w K w K n w K w K w K nw wK w w w K w w wnK w w w H 1 1 12 2 22 21 12 221 1 1 12 2 2 2, , 12.2 2, , 122.28 w116c c cz y xz y xx xx y yyxx yyc c cz zz z y xzz z y xT n ww w H w w H w wwK w w w K w w wn T n wK w w w H w w H w w H w wwn 2 2 21 1 12 2 21 12 22121 1 12 2 22 21 4 318 w16x y zxx yy zzxx zzxx zzyyyyx y zx y zK K KK KK w wn H H H 2 2 21 1 12 2 21 12 22121 1 12 2 21, ,2 2 28 w1161 4 318 w16c c cx y zxx yy zzxx zzxx zzyyyyx y zx y zT n wK K KK KnK w wn H H H 1, ,c c cT n w (4.10) 45 2 21 12 221 1 1 12 2 2 22 21 12 221 1 12 2 22 2, 1,222 22x xx y yyxx yyc c cz zz z y xzz z y xx xx y yyxx yyz zz z yzz z y xK n n n K n n nw T n wK n n n H n n H n n H n nnK n n n K n n nwK n n n H n n H n n H 12, 1,211, 1,211, 1,212c c cxxxavg x x yy avg y y zz avg z z c c cxxavg x x yy avg y y zz avg z z c c cxxavg x x yy avg y y zz avg z z cT n wn nnw K n K n K n T n ww K n K n K n T n ww K n K n K n T 1, 1,11, 1,211, , 1211, , 1212c cxxavg x x yy avg y y zz avg z z c c cxxavg x x yy avg y y zz avg z z c c cxxavg x x yy avg y y zz avg z z c c cxxavg x x yy avgn ww K n K n K n T n wn K w K w K w T n wn K w K w K w T n wn K w K 1, , 111, , 121, 1, 121, 1,2y y zz avg z z c c cxxavg x x yy avg y y zz avg z z c c cxxavg x x yy avg y y zz avg z z c c cxxavg x x yy avg y y zz avg z z c c cw K w T n wn K w K w K w T n wK n w K n w K n w T n wK n w K n w K n w T n w 11, 1, 121,1 ;1, 1, ,22, ,xxavg x x yy avg y y zz avg z z c c cxxavg x x yy avg y y zz avg z z c c c gen pc cstt tndc c c cTIMEFor time stepT TTIMEtAfter time steK n w K n w K n w T n wK n w K n w K n w T n w q n w C n wn w n w 2.;3 , , , , , ,42c c c ct t t tc c c c cpT T TTn w n w n wIMEt (4.11) 46 where, , ,xxavg yy avg zz avgK K K are average thermal properties for discrete elements for an interior node , ,c c cx y z in Figure 3-5. ( 1, , ) ( , , ) ( 1, , )3( , 1, ) ( , , ) ( , 1, )3( , , 1) ( , , ) ( , , 1)3c c c c c c c c cxxavgc c c c c c c c cyy avgc c c c c c c c czz avgK x y z K x y z K x y zKK x y z K x y z K x y zKK x y z K x y z K x y zK (4.12) Directional thermal properties in equation (4.10) & (4.11) ( 1 1 1 12 2 2 2, ... ,xx xx zz zzK K K K ) are calculated by utilizing equations (3.31) and (3.32) (discrete element sizes are determined from the equation (4.16) when using them in generalized coordinates). All mass transfer terms in equations (4.10) and (4.11) are taken as individual element basis. 1 12 21 12 21 12 2( , , ) ( , , )( ,.., ) ( , , )( , , ) ( , , .)c c c c c cc c c cxavg p xx xy avg p yy yz avg p zz zc cc c c c c cx y z x yH H H zx y z x y zx y z x y zC VH H H C VH H H C V (4.13) Note that ( , , )c c cx y z corresponds to , ,c c cn w in the computational domain. 4.3 Discretization of cutting geometry The cutting geometry is modeled by discretizing the tool, workpiece, and chip into multiple oblique cutting substructures or a single orthogonal cutting depending on the case. They are generated by splitting the cutting edge into smaller linear segments that allow the application of cutting mechanics. The substructures of tool, workpiece, and chip have prismatic geometries as illustrated in Figure 4-2 for 3D simulations, and a sample geometry for a turning tool with a nose radius is given in Figure 4-4. After discretizing the cutting operation into 47 substructures, grid points are generated within each substructure. The generated grids allow to create discrete solution space to implement a numerical solution of equation (4.8). A temperature prediction model is constructed for each substructure by considering thermal interaction with its neighboring substructures. The temperature is predicted at each grid point of the finely meshed substructure via finite difference method. Figure 4-2 Discretized cutting geometry 4.4 Grid generation Grid generation procedure is applied to each domain of interest to build the relation between computational and physical spaces. The differential mapping terms ( , ,... ,x xx z zzw w ) 48 related to domain transformation in the equation (4.8) are determined via grid generation procedure. Elliptic grid generation method is utilized in this thesis. It falls under structured grid generation methods in which grid point locations are determined based on a general governing grid generation equation. This method is selected due to its suitability to domain shapes in this thesis. Elliptic grid generation methods are well established for 2D domains, and they result in one to one transformation between computational and physical spaces (see Figure 4-3). Although governing equations are written in 3D (equation (4.8)), 2D dimensional grid point generation and extruding the same cross-sectional area, to generate 3D shapes, found adequate within the application scope. Further information about grid point generation methods is available in the literature [62], [63]. Figure 4-3- Gridpoint mapping from physical to computational space 49 An elliptic differential equation system is solved to find the grid point locations in physical space by defining the boundary nodes on the periphery of tool, chip and workpiece domains. First boundary nodes are seeded in Cartesian coordinates. Then, the interior grid point distribution is determined based on the solution of the following differential equation [64]; 2 2 2 22 2 2 20 0n nandx y x y (4.14) where and n represent coordinate directions in the computational domain. The equation (4.14) is rewritten by using the following differential operators to get the solution. First, differentiation in Cartesian coordinates is written as; n nandx x x n y y y n (4.15) Then by using the relationship in the equation (4.15), differential distances in the computational and physical domain is linked via the following equation (4.16). 1x y nx y nx xd dx dx dandn n y ydn dy dy dnx xx y nn n y yx y n (4.16) Equation (4.16) leads to the following relation, 1, ,1y x n y n xJ J J anded Jx n y n x yx xnJx y y xy yn ntn (4.17) 50 where, J is Jacobian of transformation between computational and physical domains. It can be interpreted as the ratio of differential areas between two domains for 2D mesh and ratio of differential volumes for the 3D mesh. In a computational domain, , ,n w are dimensionless, and they vary between 0 and 1 as opposed to , ,x y z are in meters. Thus, Jacobian of the transformation interrelates the dimension in between two coordinates. It changes from element to element if element shapes are varying depending on grid point location. Governing grid generation equation (4.14) can be rewritten in the following form by using equations (4.15) to (4.17). 2 2 2 232 2 2 22 22 22 22 222 232 2 2 222 2 02 2y y y x x x x yJ a b c a b cn n n n nx y ny y y x x x x yJ a b c a b cnn ny n nx n 0 (4.18) where 2 2x yan n , x x y ybn n , and 2 2x yc . Solution of the equation (4.18) requires the following two conditions to be met. 2 2 22 22 2 22 22 02 0x x xa b cn ny y ya b cn n (4.19) The equation couple (Eq. (4.19)) can be solved numerically by using the finite difference method. Central difference approximation is used to solve governing grid generation equations after seeding grids on chip, tool or workpiece domain boundary. 51 221, 2 , ( 1, )1, 1 1, 1 1, 1 ( 1, 1)2, 1 2 , ( , 1)0y i j y i j y i jay i j y i j y i j y i jbny i j y i j y i jcn (4.20) 221, 2. , ( 1, )1, 1 1, 1 1, 1 ( 1, 1)2, 1 2 , ( , 1)0x i j x i j x i jax i j x i j x i j x i jbnx i j x i j x i jcn (4.21) where ,i j are used to define element location along and n directions, respectively. Coefficients , ,a b c in (4.20) and (4.21) are as the following; 2 2, 1 , 1 , 1 , 12 2x i j x i j y i j y i jan n (4.22) 1, 1, , 1 , 12 21, 1, , 1 , 12 2x i j x i j x i j x i jbny i j y i j y i j y i jn (4.23) 2 21, 1, 1, 1,2 2x i j x i j y i j y i jc (4.24) and n are discrete element sizes in the computational domain which are dimensionless. Grid size is determined based on numerical convergence test. Solutions converge to similar temperature distributions as the grid size becomes smaller. They are defined by using number discretization points along and n directions. 52 1 1nx n yand nN Nwhere N N and N N (4.25) The grid point locations are determined iteratively. First, boundary nodes are seeded. Then, an initial approximate distribution of grid points within the domain is set by connecting the seeded boundary nodes with straight lines to start the first iteration. Coefficients , ,a b c are approximated from initial grid distribution in the first iteration. These coefficients are recalculated as the iterations progress by utilizing the latest calculated grid point distribution to calculate the next iteration. Iterative solutions stop when the absolute summed difference between grid point locations in x and y coordinates decreases below a pre-set tolerance value as; 1, 111, 11, 111, 1[ , , ][ , , ]Mi Nx j NyM Mxi ji Nx j NyMyi jabs x i j x i jabs y i j y i j (4.26) The tolerance value is set to 1 m since it is cumulative sum of all errors. More detailed information on the elliptic grid generation procedure can be found in [64]. Once grid points are generated, , ,... ,x xx z zzw w differential mapping terms are determined by using equations (4.17) and (4.15). Second order derivative terms ( , , , , , , , ,xx yy zz xx yy zz xx yy zzn n n w w w ) are identified by using the following relationships. Since the same 2D grid is extruded to create 3D shapes in this thesis, z direction second order derivatives become zero. 53 2 2 22 22 2 22 2,,, , ,., , ., ,nn n nn nn nx y x xx y x y xnx y xxyy x y yn n n nyny 2 3 232 332 232 22 202n nn n nn n n n n n nxxn nnn nn n n n n nxxn nxxn nn n n n n nn n n nyyn nnyyx y y y y x y x y y x y y x x yy x x yx y y y x y y y x x y y x y x yny x x ywy x x x x x x x x y x x x y x yy x x yyn 3 2 232 200, 00, 00, 1n nn n n n n nn nyyzz zzz zzz zx y x x x x y x y x x y x x yy x x ywn nw w (4.27) where , , , , , , , , ,n nn n nn n nx x y y x x y y x y are approximated by using finite difference method numerically. 2 21, 1, 1, 1,,2 2, 1 , 1 , 1 , 1,2 21, 2 , ( 1, ) 1, 2 , ( 1, ),1, 1 1, 1 1, 1 ( 1, 1)21, 1 1,n nnnx i j x i j y i j y i jx yx i j x i j y i j y i jx yn nx i j x i j x i j y i j y i j y i jx yx i j x i j x i j x i jxny i j y i jy 1 1, 1 ( 1, 1)2y i j y i jn (4.28) 54 4.5 Chip temperature model Chip geometry is discretized into parallel substructures for 3D cases. Each substructure is modeled as an individual oblique or orthogonal cutting operation compared to global cutting geometry in either turning or milling operations. The global chip geometry is discretized into parallel substructures (Sc1 to Scn) as shown in Figure 4-4 and Figure 4-5. 2D simulations are sufficient for orthogonal cutting cases where z the direction in heat balance equation is omitted. Grid generation is applied to each substructure to write discretized heat balance equation with finite difference approximations. Figure 4-4 Layout of chip substructures for a turning case 55 Figure 4-5 Discretization for the milling operation Heat balance equations for the chip are written by using 3D heat and 1D mass transfer problem. The heat loss due to convection and radiation were neglected in the derived equations. Mass transfer rate is calculated based on the chip velocity cV . It is assumed that chip moves as a rigid body in the chip flow direction. The cartesian x direction is taken as chip flow direction. Hence, the velocity components in y and z direction are set to zero in equations (3.23) and (3.26). 0 V 0 30 2y zyy zV and for D Simulationt tyV for D Simulationt (4.29) Primary and secondary zone heat sources are used in modeling the chip temperature. It is assumed that chip leaves shear zone with an average temperature and then further heated by the 56 frictional heat source at the tool-chip interface. Shear zone temperature is determined by using the equation (3.8). In the sliding region, the uniform frictional heat flux is applied by introducing a local heat partition ratio between chip and tool within their contact zone. Sliding heat partition assumption is based on Blok’ s heat partition principle [65]. This principle dictates that average temperatures on reciprocal surfaces of one sliding and one stationary body should be equal. In the sticking region, frictional heat generation is applied as an internal volumetric heat generation within a predefined sticking region. Depending on the volume elements location, heat source formulations are updated. Related formulations are given as follows; 57 ( ) 1 {3 ,} , FV1 int1interface214csliding friction c z f t consconsconsconsQ q w nq B N i j in sliding regionFV for inernal nodes at tool chip erfaceFV for nodes onboundary lines at tool chipFFor DV for corner nodes at tooSimulations 1interface( ) 1 { }1int interface21interface418csticking friction c sticking t consconsconsconsl chipQ q Q B N i J n wIV in sticking regionIV for erior nodes at tool chipIV for nodes onboundary lines at tool chipIV for corner 2interface( ) 1 FV1 int1inter,face2(ycsliding friction c f consconsconscsticking friction cnodes at tool chipQ q n q B i in sliding regionFV for inernal nodes at tool chip erfaceFV forFor D Simucorner nodes at tlatiool chipo sQnq 1sec) 11interface line21interface4,sticking consconsconsstickingstickingpl ondaryQ B i J nIV in sticking regionIV for nodes ontool chipIV for corner nodes at tool chipPwhere Qbl t (4.30) where { } ,tB N i j is the proportion of frictional heat flux that flows into the tool. ,c tq N denote chip element number and tool substructure number respectively. ,i j or i is tool-chip interface element indices. Since finite volume relations are updated depending on element 58 location (i.e., corner node, borderline, interior node), the equation (4.30) has element volume modifiers. FVcons is flux volume constant; it modifies the flux application equation by considering related elements volume. consIV is the internal heat generation volume constant; it modifies the internal heat generation application equation by considering related elements volume. Irrelevant flux terms are not included in the heat source equations which arise when derivations are made for a general case. Figure 4-6 Chip substructure Chip geometry is discretized parallel to chip flow direction in a nested pattern where each substructure represents an oblique cutting operation in 3D modeling. A sample layout of substructures is shown in Figure 4-4. Thermal interaction between substructures is considered by applying temperature values at adjacent faces as a boundary condition. The length of consecutive substructure faces is kept equal to avoid the mesh skewness in the chip grid generation within each substructure. Hence, some chip substructure lengths are longer than local tool-chip contact length of the oblique cutting segment, and back end of chip substructures are connected with a conjoint face. The distance between P1 and P3 (Figure 4-6) is equal to 59 differential cutting edge length ( b in Figure 3-1) for chip substructures that are engaged with the tool. P1 to P2 distance (Figure 4-6) is equal to chip thickness. An average chip thickness value is used for each oblique cutting substructure. In 2D simulations, O1 to O2 distance is defined as tool-chip contact length, and O2 to O3 distance is taken equal to chip thickness. 3 , 1 3 , 1 22 , 1 2 , 2 3cc cIn D Simulations P P b P P hIn D Simulations O O l O O h (4.31) Grid points are generated after defining chip geometry boundaries. x direction is set parallel to chip flow direction, while z direction is selected perpendicular to the tool-chip interface for 3D substructure (Figure 4-4). In 2D cases, only x and y directions are modeled, and z direction is eliminated from the heat balance equation. Origin point is placed at P1 and O1 for 3D and 2D cases, respectively. Elements are numbered in a continuous manner by using cq index. Element index and the total number of elements are calculated as follows; 31 1 1 1 11 1 121 11 1, , , ,c x z xtotal x y zc xtotal x yD Simulationq N N j N k iN N N ND Simulationq N j iN N Nwherei j k are x y z indices respectively (4.32) Heat balance equations in discrete form (equations (4.10) and (4.11)) are assembled for all elements as the last step after grid generation with application of heat sources and defining boundary conditions. A linear set of equations is assembled for chip domain and solved at each 60 time step as follows; 11,1 ... 1, ... 1, 1,1 1,1,1 ... , ... , ,1 ,1,1 ... , ... , ,1total total totalcc total cc c c c total c c ctotal total c total total c total totaN X N N XA T CA A q A N T CA q A q q A q N T q C qA N A N q A N N T N C N 1,1N Xtotall (4.33) where A is the coefficient matrix of unknown nodal temperatures of the time step to be predicted. It includes all the coefficients that are coming from discretized heat balance equation (4.8).C is the heat source matrix which includes temperature history, known temperature values as boundary conditions, primary zone and secondary zone heat source effects. If simulations are to be made for constant thermal properties, coefficient matrix stays same at all time steps. If thermal properties change with temperature, the coefficient matrix is updated. In this case, both A and C are split into parts that correspond to conduction, mass transfer and time-dependent sub-matrices and sub-arrays. They are generated for unit material properties and saved. Later, they are multiplied by the corresponding material property (as it changes locally with temperature) at each time step to update them. They are represented as; 61 1 12 21 12 21 12 21 12 2xx avgxx xx total totaltotal total total totalyy avgyy yy total totaltotal total total totalu uuK K xxavg Kxx xx N X NN X N N X Nu uK K yy avg Kyy yy N X NN X N N X NA K A K A K AK A K A K A 1 12 21 12 21 12 21 12 2zz avgzz zz total totaltotal total total totalx avgx x total totaltotal total total totaluu uuK K zz avg Kzz zz N X NN X N N X Nu uuH H xavg H px x N X NN X N N X NK A K A K AH A H A H A C A 1 11 12 21 12 2( , , ) ( , , , ), ( , , , ), ( , , , )[ ]total totaltotal X total Xxx xxtotal total total totalut t N X Ni p i ii pu uboundary boundaryK c K cN Nxx xxN X N N X Nxwhere H C V i x y z K f x y z T f x y z T C f x y z TC K C T K C TK 11 121121212xx avgtotal X total Xyytotal totaltotal totalyy avgtotal Xyy total totaltotal totaluuboundary boundaryxavg K c K cN NyyN X NN X Nuuboundary bouK c yy avg K cNyy N X NN X NC T K C TK C T K C T 11 11 12 21121 12 212total Xtotal X total Xzz zztotal total total totalzz avgtotal X xtotal totalndaryNu uboundary boundaryK c K cN Nzz zzN X N N X Nuboundaryzzavg K c HN xN X NK C T K C TK C T H C 111 1212total Xtotal totalx avgtotal X total Xx total totaltotal totaltotal total total totaluboundaryc NN X Nuuboundary boundaryH c xavg H cN Nx N X NN X Ntu uold olp t c p t tN X N N X NTH C T H C TC C T C C T t td c cc sliding friction sticking frictionQ Q (4.34) 62 where total totalut t N X NA contains the unit coefficients of unknown temperature in the transient part of the heat balance equation (3.26). ,i itotal total total totalu uK HN X N N X NA A terms are square coefficient matrixes of steady-state terms with unit anisotropic heat conduction and mass transfer coefficients in the heat balance equation (3.26). ,i itotal total total totalu uK HN X N N X NC C terms are square coefficient matrixes with unit anisotropic heat conduction and mass transfer coefficients to apply constant temperature boundary conditions. 1totalut N XC and 1totalut t N XC are unit coefficients of previous time-step temperatures. ,c csliding friction sticking frictionQ Q are heat source arrays for the chip surface. 4.6 Tool temperature model Tool geometry is discretized into oblique cutting substructures as in chip discretization for 3D cases. The global tool geometry is discretized into parallel substructures (St1 to Stn) as shown in Figure 4-7 and Figure 4-5. Orthogonal cutting cases are modeled by using a 2D form of the heat balance equation. 63 Figure 4-7 -Tool substructure layout Heat balance equations are written based on the heat conduction only without mass transfer by considering adiabatic boundaries. The room temperature boundary condition is used at non-engaged ends of tool geometry. Although boundaries were modeled adiabatic, convection heat loss has been tested on tool geometry for sensitivity analysis. Since the mass transfer is not in effect for the tool, the velocity components in ,x y and z direction are set to zero in equations (3.23) and (3.26). 0, 0 V 0x y zx y zV V andt t t (4.35) Frictional heat sources on the rake and flank faces of the tool are considered. They are applied to elements depending on the contact pattern between the chip and workpiece. Tertiary zone (flank face) heat source is only modeled for sensitivity analysis purpose. Rake face heat source is modeled as either sliding, sticking or sticking-sliding assumptions. Sliding heat source 64 application differs from chip temperature modeling by the application of heat partition ratio. { } ,tB N i j percentage of frictional heat goes into the tool instead of 1 { } ,tB N i j in chip temperature model. The heat source in the sticking region is applied through the sticking heat flux since no internal tool deformation zone is assumed. Therefore, the calculated internal secondary zone volumetric heat generation rate stickingP is converted to surface flux stickingq for the tool. Heat source terms on the rake face are summarized as follows; ( ) { } , FV1 int1interface13 ,24tsliding friction t z f t consconsconsconsQ q w n q B N i j in sliding regionFV for inernal nodes at tool chip erfaceFV for nodes onboundary lines atFor D Simtool chipFV for corner nodesulatat tool cions interface( ) { } , FV,( ) F1,V2ytsticking friction t z sticking t consstickingstickingplasticcontacttsliding friction t f consconsFor D SimulhipQ q w n q B N i j in sticking regionPwhere qAQ q n q B i in sliding regionationFsV fo int1interface2( ) FV,yconststicking friction t sticking consstickingstickingplr inernal nodes at tool chip erfaceFV for corner nodes at tool chipQ q n q B i in sticking regionPwhere ql (4.36) where { } ,tB N i j is the proportion of frictional heat flux that flows into the tool. ,t tq N denote tool element number and tool substructure number, respectively. ,i j or i is tool-chip 65 interface element indices. Partition ratio array { } ,tB N i j is determined through recursive solutions at each time step. The detailed calculation procedure is presented in section 4.8. The tertiary zone heat source is modeled under sliding friction conditions for sensitivity analysis purpose only. It is applied to the newly generated surface by introducing ploughingq which is defined in section 3.2.3. Similar to the application of secondary zone heat source, following derivations are made. ( ) 2{ } FV1 int interface1interface2143 ,ttertiary t ploughing t consconsconsconsxQ q n wq mean B NFV for ernal nodes at tool workpieceFV for nodes onboundary lines at tool workpiForeceFV for corner nodes at tool wD Simulatoionsr interface( ) 2 FV1 int in2 ,terface1interface2xttertiary t ploughing consconsconskpieceQ q n q mean BFV for ernal nodes at tool workpieceFV for corner nodes at tool workpiecFor D Simulatioens (4.37) where 2{ }, 2tB N B are the proportions of ploughing heat flux that flows into the tool. ,t tq N denote tool element number and tool substructure number respectively. The partition value 2{ } ,tB N i j is approximated by using the ratio of thermal effusivity values of tool and workpiece material. This assumption is based on the static contact and steady-state conditions. Since the tertiary zone heat source is only considered for sensitivity analysis, it is not calculated dynamically under transient conditions by making a sliding contact assumption. 66 { } , { } , { } ,2{ } ,{ } , { } , { } , { } , { } , { },,,223tool tool toolt t p ttwpc wpc wpc tool tool toolt t p t t t p ttool tool toolpwpc wpc wpc toolpFor D SimulationsFor D SK N i j N i j C N i jB N i jK N i j N i j C N i j K Nii j N i j C N i jK i i C iB iK i i C i K imulations ::tool toolpPlease noti C iethatwpc workpiece (4.38) where ,i j or i is tool-workpiece contact interface element indices. Global tool geometry is discretized into multiple oblique cutting operations and assembled as adjacent substructures. Assembly layout is shown in Figure 4-7 and Figure 4-5. Thermal interaction between substructures is considered by applying temperature as a boundary condition at their common face from both adjacent faces. Example substructures for 3D and 2D cases are shown along with coordinate directions in Figure 4-8. Figure 4-8 - Substructure shapes for tool 67 Distance from P3 to P4 in Figure 4-8 is equal to differential cutting edge length ( b in Figure 3-1). Since room temperature boundary condition is applied on P1P2P5P6 (away from rake face), the distance between P5 and P7 (Figure 4-8) is generally taken multiples ( min ~ 10. cl) of contact length or equal to insert thickness. The distance between P5 and P7 distance is set based on the numerical experimentation; it is set long enough not to affect rake face temperature distribution. In 2D simulations, length of O2 to O3 is set similar to 3D case ( min ~ 10. cl ). O3 to O4 distance is set arbitrarily; generally, it is minimum four times the contact length. 3 , 3 4 , 5 7 min ~ 10.4 8 3 7 min ~ 4.2 , 3 4 min ~ 4. , 2 3 min ~ 10.ccc cIn D Simulations P P b P P insert thickness or lP P P P lIn D Simulations O O l O O insert thickness or l (4.39) Element locations are determined after tool geometry is defined. Coordinate direction for tool substructures is given in Figure 4-8. Origin point is placed to P1 for 3D cases, and O1 for 2D cases in Figure 4-8. tq index is defined to number tool elements. Total number of tool elements and element index is set by using the following; 31 1 1 1 11 1 121 11 1, , , ,t x z xtotal x y zt xtotal x yD Simulationq N N j N k iN N N ND Simulationq N j iN N Nwherei j k are x y z indices respectively (4.40) After grid generation, discrete heat balance equations (equation (3.26) in discrete form) are assembled for all elements. Assembled linear set of equations is solved for each tool 68 substructure at each time step to determine unknown temperature values (see equation (4.41)). 11,1 ... 1, ... 1, 1,1 1,1,1 ... , ... , ,1 ,1,1 ... , ... ,[ ],1total total totalt total tt t t t total t t ttotal total t total total t total toN X N N XtD D q D N T ED q D q q D q N T q E qD N D N q D N N T EDN NT E 1,1N Xtotaltal (4.41) where D is the coefficient matrix of unknown nodal temperatures to be predicted in current time step. E is the heat source matrix which includes temperature history, known temperature values as boundary conditions, secondary zone and tertiary zone heat source effects. Coefficient matrix stays invariant at all time steps for models with constant thermal properties. It is dynamically updated as the temperatures change when varying material properties are used. Similar to chip modeling, both &D E are split into varying time-dependent sub-matrices and sub-arrays. These are generated for unit material properties and saved. Later, they are multiplied by the corresponding material property (as it changes locally with temperature) at each time step to update them. The coefficient and heat source matrixes are split as; 69 1 12 21 12 21 12 21 12 2xx avgxx xx total totaltotal total total totalyy avgyy yy total totaltotal total total totalu uuK K xxavg Kxx xx N X NN X N N X Nu uK K yy avg Kyy yy N X NN X N N X ND K D K D K DK D K D K D 1 12 2121 12 212( , , , ) , ( , , , ), ( , , , )[ ]zz avgzz zz total totaltotal total total totaltotal totalxxtotaluu uuK K zz avg Kzz zz N X NN X N N X Nup t t N X Nii pKxxN XK D K D K DC Dwhere K f x y z T f x y z T C f x y z TE K E 11 1211 121212total X total Xxxtotal total totalxx avgtotal X total Xyytotal totaltotal totalu uboundary boundaryt K tN NxxN N X Nuuboundary boundaryxxavg K t K tN NyyN X NN X NyyT K E TK E T K E TK 11 121 112 2121 12 2yy avgtotal X total Xyy total totaltotal totaltotal Xzz zztotal total total totaluuboundary boundaryK t yy avg K tN NN X NN X Nu uboundary bouK t K tNzz zzN X N N X NE T K E TK E T K E T 11total Xzz avg total totaltotal Xtotal totaltotal totalndaryNu tuboundary oldzz avg K t p t tN X NNN X Nt tu old t t tp t t t sliding friction sticking friction tertiaryN X NK E T C E TC E T Q Q Q (4.42) where total totalut t N X ND contains the unit coefficients of unknown transient temperatures to be predicted in discrete form of heat balance equation (3.26). itotal totaluKN X ND terms are square coefficient matrix with unit heat conduction coefficients of steady-state terms in heat balance equation (3.26). itotal totaluKN X NE consists of coefficients to apply constant temperature boundary conditions with unit heat conduction terms. 1totalut N XE and 1totalut t N XE are coefficient arrays of 70 previous temperatures with unit material properties. , ,t t tsliding friction sticking friction tertiaryQ Q Q are heat sources on the tool surface. 4.7 Workpiece temperature model Global workpiece geometry is modeled by dividing the tool-workpiece engagement zone into parallel substructures (Sw1 to Swn) as shown in Figure 4-5 and Figure 4-9. Figure 4-9- Workpiece substructure layout Heat balance equations are written based on 3D heat conduction and 1D mass transfer similar to chip modeling without considering radiation and convection heat loss. Mass transfer rate is set based on the cutting velocity V . 2D form of heat balance equation is used for orthogonal cutting cases. All the governing heat balance equations are written based on the rigid body movement of the workpiece with cutting velocity. The cartesian direction z is taken as the cutting velocity direction for 3D modelling. The velocity components in x and y direction are 71 set to zero in equations (3.23) and (3.26). Whereas, the cutting velocity direction is selected to be the Cartesian x direction for 2D modeling. 0 V 0 320x yyD Simux yV and fort txVlationD Simulatiofortn (4.43) The primary zone heat source is considered in the analysis as the volumetric heat generation. Tertiary zone heat source between the tool flank and workpiece is also considered for sensitivity analysis purpose. Room temperature boundary condition is applied to below surface away from the heat sources for each workpiece substructure. It is assumed that neighboring chip elements at shear zone are interacting with corresponding workpiece elements. Therefore, the corresponding chip temperatures are applied as a boundary condition to the workpiece. The tertiary zone heat source is applied as; ( ) 1 2{3} FV1 int interface1interface2(,2 ,wtertiary w ploughing t consconysconswtertiaryQ q n wq mean B NFV for ernal nodes at tool workpieceFV for nodes onboundary lines at tool workpFor D SimulationsFor D SimulationsieceQ q ) 1 2 FV1 int interfacexw ploughing consconsn q mean BFV for ernal nodes at tool workpiece (4.44) where 2{ } ,tB N i j is the proportion of frictional heat flux that flows into the tool. ,w tq N denote workpiece element number and tool substructure number, respectively. 72 Sample workpiece substructures are presented in Figure 4-10 for 3D and 2D cases. Dimensions of these substructures are determined by considering the engagement profile. Figure 4-10- Sample workpiece substructure P8 to P12 length in Figure 4-10 is set to shear plane length or cut chip thickness. Since room temperature boundary condition is applied on P1P2P13P14, P15 to P13 distance (Figure 4-10) is generally taken multiples ( min ~ 10. cl ) of contact length. P15 to P13 distance is determined through numerical experimentation to find the distance that does not affect temperature values close to the surface. In 2D simulations, shear zone length is set equal to chip thickness. O1 to O4 distance is set to multiples of contact length ( min ~ 10. cl ) due to applied room temperature boundary condition between O1 & O2 in Figure 4-10. 3 , 8 12 , 13 15 min ~ 10.2 , 1 4 min ~ 10.c ccIn D Simulations P P h P P lIn D Simulations O O l (4.45) Grid generation procedure is applied after the workpiece substructure boundaries are defined. Origin point for the geometry is placed at P1 for 3D cases, and O1 for 2D cases in 73 Figure 4-10. wq index is used to number workpiece elements. Total number of elements and element index are set as; 31 1 1 1 11 1 121 11 1, , , ,w y z ztotal x y zt xtotal x yD Simulationq N N i N j kN N N ND Simulationq N j iN N Nwherei j k are x y z indices respectively (4.46) After defining heat generation terms and boundary conditions, equilibrium equation (equation (3.26) in discrete form) for all control volumes is assembled in a matrix form. Then, the linear equation system is solved for workpiece substructure at each time step. 11,1 ... 1, ... 1, 1,1 1,1,1 ... , ... , ,1 ,1,1 ... , ... ,[ ],1total total totalw total ww w w w total w w wtotal total w total total w total toN X N N XwR R q R N T SR q R q q R q N T q S qR N R N q R N N T SRN NT S 1,1N Xtotaltal (4.47) where R is the coefficient matrix of unknown nodal temperatures to be predicted in current time step. S is the heat source matrix which includes temperature history, known temperature values as boundary conditions, primary zone and tertiary zone heat source effects. In constant thermal properties model, the coefficient matrix is the same for all time steps. On the other hand, it is dynamically updated as the temperatures change when varying material properties are used. As 74 in chip and tool modeling, both &R S are split into varying time-dependent sub-matrices and sub-arrays when thermal properties change. The coefficient matrix and heat source matrix are generated for unit material properties and saved. Then, they are multiplied by the corresponding material property at each time step to update them with the new ones. Coefficient and heat source matrixes are split as; 75 1 12 21 12 21 12 21 12 2xx avgxx xx total totaltotal total total totalyy avgyy yy total totaltotal total total totalu uuK K xxavg Kxx xx N X NN X N N X Nu uK K yy avg Kyy yy N X NN X N N X NR K R K R K RK R K R K R 1 12 21 12 21 12 21 12 2zz avgzz zz total totaltotal total total totalz avgz z total totaltotal total total totaluu uuK K zz avg Kzz zz N X NN X N N X Nu uuH H z avg H pz z N X NN X N N X NK R K R K RH R H R H R C R 1 11 12 21 12 2( , , ) ( , , , ), ( , , , ), ( , , , )[ ]total totaltotal X total Xxx xxtotal total total totalut t N X Ni p i ii pu uboundary boundaryK w K wN Nxx xxN X N N X Nxwhere H C V i x y z K f x y z T f x y z T C f x y z TS K S T K S TK 11 121121212xx avgtotal X total Xyytotal totaltotal totalyy avgtotal Xyy total totaltotal totaluuboundary boundaryxavg K w K wN NyyN X NN X Nuuboundary bouK w yy avg K wNyy N X NN X NS T K S TK S T K S T 11 11 12 21121 12 212total Xtotal X total Xzz zztotal total total totalzz avgtotal X ztotal totalndaryNu uboundary boundaryK w K wN Nzz zzN X N N X Nuboundaryzzavg K w HN zN X NK S T K S TK S T H S 111 12121 11total Xtotal totalzavgtotal X total Xz total totaltotal totaltotal totaltotaluboundaryw NN X Nuuboundary boundaryH w z avg H wN Nz N X NN X Ntu uold oldp t w p t t wN X N XN XTH S T H S TC S T C S T 1totalt twtertiaryN XQ (4.48) where total totalut t N X NR includes the coefficients of discrete time derivative terms of heat balance equation (3.26). ,i itotal total total totalu uK HN X N N X NR R terms are unit material property based directional coefficient matrix components to be updated. Similarly, ,i itotal total total totalu uK HN X N N X NS S include 76 finite difference approximation coefficients to apply constant temperature boundary conditions. 1totalut N XS and 1totalut t N XS include the coefficients of the temperature history. wtertiaryQ is tertiary zone heat source array. 4.8 Simulation algorithm Simulation algorithm is categorized into two main sections. In the first one, input data to simulation is prepared. This initial step is carried out through the determination of the tool-chip contact area, power generation related to deformation zones, preparing chip, tool and workpiece geometries based on the cutting process. In the second section, iterative solution of discretized heat balance equations is carried out by applying boundary conditions. Grid points are generated first for each tool, chip and workpiece substructure based on the local cutting geometry by using the procedure explained in section 4.4. Then using the grid generation output, heat balance equations are written for each element position for unit material properties by using finite difference approximations (4.10) and (4.11) to be updated with material properties later. These equations are written depending on the element position (corner nodes, interior nodes, borderlines, and faces) and applied boundary condition (constant temperature, heat flux, convection). If variable material properties are used, coefficient matrixes are constructed from discrete heat balance equations by considering directional material property changes. They are generated for once for unit material properties and later used by updating them with corresponding material properties. The tool-chip contact region is determined in this step to tag heat flux application nodes at their interface. Chip, workpiece coefficient matrixes, and contact map are regenerated if chip thickness change with time such as in milling. Tool and chip contact map is calculated by using the relative distance between the tool and chip interface nodes for all substructures that may be in contact with each other. The 77 closest nodes from tool and chip sides at their interface are picked as reciprocal nodes to check convergence via temperature comparisons between them. Thermal properties at room temperature are assigned before the iterative solution phase to start the solutions. Later, thermal properties, coefficient and heat source matrixes, and arrays are updated as the thermal field is recalculated. If constant thermal properties are used, same assigned values are used throughout the recursions. Solution to the equation (4.8) is carried out recursively in the second stage of the algorithm. In this step, time-varying parts of the heat source and coefficient matrices are updated iteratively. The solution starts by assuming room temperature in all tool, chip and workpiece substructures. Main iteration parameter for each time step is the temperature difference ( dB ) between the reciprocating nodes of chip and tool at their interface. The temperature at the reciprocating nodes is iteratively converged via tuning frictional heat flux to chip and tool. Heat partition distribution is calculated as a result of the iterations. First, the initial values between 0 and 1 are assigned to all tool-chip interface nodes. After each iteration, heat partition values are updated by using dB value which is defined as; ( ) { }( )( ) { }( )2c c c t t tc c c t t tT N q T N qdBT N q T N q (4.49) where tq and cq are the tool and chip element numbers, respectively. Maximum absolute value of dB is checked in each iterated solution. If the difference is below a predefined set value, the solution is acceptable regarding the temperature difference between tool-chip interface reciprocal nodes. If there is a difference, the heat partition values are updated for each node at 78 the tool-chip interface to solve the heat balance equation again. Partitions are updated for contacting tool and chip nodes on the rake face by using: .t tt tB q B q Gai dBN nN (4.50) The gain value is used as a multiplier to adjust convergence speed. It can be set between 0 and 1 for convergence when the temperature difference between reciprocal nodes of the tool and chip that are in contact becomes high in some cases. Note that the heat partition ratio determination procedure is repeated at each time step. Each substructure uses the temperature values of adjacent faces as a boundary condition at each iteration. The boundary conditions are applied by imposing preceding recursion’s temperature in the solutions. Hence, the pre-set difference measure is defined to neglect the temperature distribution difference between successive iterations. This second convergence criterion within each time step is defined as; 1iteration iterationabs Norm T Norm T Tolerance (4.51) A third convergence criterion is required based on the thermal properties when temperature dependent thermal properties are used. The absolute difference of the maximum norm of thermal conduction values between successive recursions at each time step must be below pre-set tolerance to proceed to next time step. 111max , 2,iteration iterationchip chipiteration iterationtool tooliteration iterationwpc wpcNorm K Norm KK abs Norm K Norm K ToleranceNorm K Norm K (4.52) The flow chart of the algorithm to simulate temperature distribution is given in Figure 4-11. 79 Figure 4-11 – The flow chart of temperature prediction Generalized minimum residual method function, called “gmres” in MATLAB©, was used in solving of the resulting set of linear equations in sparse matrix form. After each solution recursion, an error refinement procedure is applied in which heat source matrices are reconstructed by direct multiplication of coefficient matrix and calculated temperatures for the chip as follows; 80 1error cerrorc errorupdatedc c cerrorwhile error pre set toleranceC C A Terror mean abs CT A CT T Tend (4.53) The same procedure is applied to workpiece and tool solutions as well. 81 Chapter 5: Simulations and Experimental Validations of Temperature Model 5.1 Overview The developed model has been used to simulate various conditions with constant and variable thermal properties. Simulation outputs have been compared against experimental data for both continuous and interrupted cutting operations. Tool coating application using anisotropic medium properties is explained, and benchmark cases against experimental data are presented. The primary functions of the developed model were tested with different workpiece material and tool material combinations. The simulations chapter is organized as follows; first, simulations with constant thermal property and their benchmarks against measurements from the literature are presented. Then, the variable thermal property simulations of cutting AISI 316L steel are compared against experimental measurements collected from the tool instrumented with thermocouples. Single layer and multi-layer coated tool simulations are compared against experimental measurements reported in the literature. A discussion about the comparison between static and sliding heat partition approximations is presented. The chapter is concluded with a simulation input sensitivity study. 5.2 Validation of constant thermal property approximation Ti6Al6V-2Sn, Al2025, Ti6Al4V, and Gray Cast Iron machining thermal measurements from the literature are used to validate 2D and 3D simulations with constant thermal property approximation. Simulations are carried out with 10 µm nominal element size and 50-microsecond discrete time steps ( t) . Temperature difference ( dB ) between tool and chip nodes at their interface was set to 1% and is set to 0.001 oC between successive recursive solutions. Note that all the secondary zone energy generation in this subsection is determined by 82 using the sliding contact assumption at the tool-chip interface. Some of the presented results are published in ASME Journal of Manufacturing Science [1]. Kitagawa et al. [66] embedded thin WC thermocouples (via sintering to tool matrix) with 25µm diameter to measure temperature on the rake face of the tool. They measured transient temperatures in both interrupted and continuous turning of Ti6Al6V-2Sn work material with K10 uncoated carbide tools. The rake angle of the tools was 6o degrees. The cutting speed was 100 m/min, and the feed rate was 0.1 mm/rev. The cutting and non-cutting lengths in the interrupted cutting experiments were set to 6mm and 30mm, respectively. Thermal properties and orthogonal cutting data are listed in Table 5-1 and Table 5-2. Table 5-1-Thermal properties of the tool and workpiece materials for validations Material Ti6Al6V-2Sn [66] K10 Carbide [66] Al2024 [67] Gray Cast Iron [67] C2 WC [67] Ti6Al4V [68] Thermal conductivity (W/m.K) 17.2 79.5 150.94 60.5 28.4 13.3 Heat capacity (J/Kg.K) 541 398 956.6 478.05 276 649 Density (Kg/m3) 4430 11800 2780 7150 11100 4346 Table 5-2- Workpiece material orthogonal cutting constants Material Ti6Al6V-2Sn [66] Al2024 [6] Gray Cast Iron [6] Shear Angle (o) 40 22.6 23 Friction Angle (o) 24 26.8 30.5 Shear Stress (MPa) 549 354.4 680 Rake Angle (o) 6 0 0 Cutting Speed (m/min.) 100 81.6 81.6 Feed Rate (mm/rev.) 0.1 0.109-0.165 0.109 Continuous and interrupted turning simulation results for Ti6Al6V-2Sn are shown in Figure 5-1. The simulated rake face temperatures are compared to experimental data given by 83 Kitagawa et al. [66]. Fresh chip geometry is regenerated when cutting starts in interrupted turning simulations. The chip and workpiece temperature fields are generated when the tool and workpiece are engaged. It is assumed that the newly generated chip geometry starts to form from room temperature (20oC) at each engagement cycle. After each cutting cycle, frictional heat flux on the tool is removed, and the tool is cooled down with the room temperature boundary conditions until the next engagement period. Interrupted cutting simulation results after five engagement periods are presented in Figure 5-1. Rake face temperature fluctuates between cutting (heating) and non-cutting (cooling) periods from 805 oC to 115 oC in interrupted cutting simulations. Steady-state tool temperature was experimentally measured 968 oC in continuous turning case. Predicted temperature value for the same case was 888 oC which corresponds to approximately -%8 difference. Cooling cycle in interrupted cutting tests results in lower temperature values at measurement location when compared to continuous cutting. Underprediction of average rake face temperature in continuous cutting can be related to primary shear zone modeling where Taylor-Quinney coefficient is taken around 0.9., which can be increased to allow more heat input to the tool. Another reason that may contribute to under prediction is fresh tool assumption; the worn edge penetrates and rubs against the newly generated surface, and transmits more heat into tool and workpiece. 84 Figure 5-1-Ti6Al6V-2Sn rake face simulation results compared to experimental results reported by Kitagawa et al. in [66] Heating and cooling cycles of the cutting tool for interrupted turning of Ti6Al6V-2Sn are shown in Figure 5-2 and Figure 5-3. Figure 5-2 shows average rake face temperature over first interrupted turning period in Figure 5-1. Tool heats up to 805 oC in 3.6 milliseconds of cutting, and then cools down to 115 oC in around 20 milliseconds. Figure 5-3 shows tool and chip thermal gradient at seven different moments within cutting and non-cutting periods. These seven moments are also marked with Roman numbers in Figure 5-2 for referencing purpose. Results show that there is a significant amount of temperature gradient change in the tool. 85 Figure 5-2- Ti6Al6V-2Sn interrupted cutting rake face temperatures through one cutting and non-cutting period 86 Figure 5-3- Ti6Al6V-2Sn interrupted cutting tool and chip thermal gradients through one cutting and non-cutting period (time marker numbered with roman numbers are shown in Figure 5-2 87 Heat flux that goes into tool fluctuates as a result of transient effects, and settles to a steady state value. Since tool body is colder in the beginning of cutting, heat flux towards to tool is initially higher compared to steady-state behavior (where it is almost constant). Evolution of heat partition ratio is shown in Figure 5-4 for the continuous turning of Ti6Al6V-2Sn. The heat partition ratio is around %90 in the beginning of cutting, and settles to %42 under steady-state conditions as the tool temperature increases. Similarly, the evolution of heat partition towards tool in interrupted cutting is displayed in Figure 5-5. Heat partition ratio varies between %76 and %51 for this case. Since partition change graph represents the results after five engagement periods (tool gets warmer as the tool engages and disengages), initial heat partition ratio is lower than the continuous case (%90). The initial high value of heat partition ratio is significant from tool wear perspective due to steep thermal gradient generation in regions close to the tool-chip interface. The low thermal conductivity of the workpiece material (i.e., titanium alloys) increases thermal load on the tool, and contributes to reduced tool-life due to increased heat transmission to tool [69]. The nature of thermal loading in interrupted cutting allows potentially higher cutting speeds in milling Titanium compared to continuous turning. Cycling between heating and cooling of the cutting tool may cause high thermal stresses that may lead to the initiation of thermal cracks, although there are lower temperatures in interrupted cutting. 88 Figure 5-4- Average heat partition to the tool during continuous turning of Ti6Al6V-2Sn Figure 5-5- Average heat partition to the tool during interrupted turning of Ti6Al6V-2Sn 89 Interrupted Ti6Al6V-2Sn simulations were carried out for 2D geometries by neglecting the third dimension in governing equations. 3D simulation capability is demonstrated only for continuous orthogonal turning case of Ti6Al6V-2Sn. Figure 5-6 shows 3D analysis results for the chip, tool, and workpiece. Oblique tool geometries can be modeled. However, current geometry discretization capacity of the simulation engine needs to be adapted to different insert and tool geometries in 3D. Figure 5-6- Turning Ti6Al6V-2Sn temperature distributions on the chip, tool, and workpiece at steady-state 90 Al2024-T351 and Gray cast iron cutting validations for interrupted cutting operations are carried out based on Stephenson and Ali’ s [6] experimental data. They used C2 tungsten carbide (WC) tools and measured temperature based on the tool-workpiece thermocouple technique. The measurement method is based on a change in EMF (Electro-Motive Force) at tool-chip contact depending on the temperature. Two types of slotted tubes were used with different cutting and non-cutting lengths in the experiments. First tube type had 51 mm cutting and 9 mm non-cutting length. The second had 25 mm cutting and 34 mm non-cutting lengths. Simulation parameters are listed in Table 5-1 and Table 5-2. Orthogonal cutting parameters are extracted from force measurements (see Table 5-2) reported in [6]. 10-micron nominal element size was used in the 2D simulations. Gray cast iron simulation results show that general prediction trends follow the experimental data with approximately 12% error. The differences between experimental and simulation results on the time scale (Figure 5-7) are due to the data digitization process from result plots in [6]. Temperatures fluctuate between 527 oC and 306 oC in Figure 5-7. Percentage of the heat that goes into the tool (the heat partition) is presented in Figure 5-8; it changes between %50 and %15. Since the thermal conductivity of grey cast iron is better than the titanium, the maximum amount of heat partition at the beginning of the cut is lesser than the one in cutting titanium. This behavior is particularly important because it gives an idea about the severity of thermal shock in the beginning of the cut. 91 Figure 5-7- The comparison of average rake face temperature prediction in the interrupted turning of gray cast iron with the experiments and computations reported in [6] Figure 5-8- Average heat partition into the tool during interrupted turning of Gray Cast Iron 92 Al2024 simulation results show less than %10 error for the two different feed rates at 0.109 mm/rev and 0.165 mm/rev. The comparison of average rake face temperatures is shown in Figure 5-9 and Figure 5-11. Heat partition ratios change between %50 and %10 (see Figure 5-10 and Figure 5-12). The reason behind the lower heat partition ratio compared to both Gray cast iron and Ti6Al6V-2Sn simulations are the higher thermal conductivity of Al2024. The lower heat partition ratio indicates more heat taken away by the chip. Heat partition ratios spike after the first cooling cycle to a higher value (Figure 5-10 and Figure 5-12) due to existing thermal gradient changes with time over different cutting cycles. As loading and unloading cycles progress, the average temperatures at the tool-chip interface tend to increase and heat builds up within the tool under simulated conditions. Figure 5-9- The comparison of average rake face temperature prediction in the interrupted turning of Al2024 with the experiments and computations reported in [6] for f=0.109 mm/rev 93 Figure 5-10- Average heat partition into the tool during interrupted turning of Al2024 for f=0.109 mm/rev Figure 5-11- The comparison of average rake face temperature prediction in the interrupted turning of Al2024 with the experiments and computations reported in [6] for f=0.165 mm/rev 94 Figure 5-12 Average heat partition into the tool during interrupted turning of Al2024 for f=0.165 mm/rev Sato et al. [19] made milling tests on Ti6Al4V with K10 uncoated carbide end mill. They measured the temperature with a custom build pyrometer made out of indium arsenide (InAs) and indium antimonide (InSb). Their developed pyrometer had 1 µs response time which allows high sampling rates. Down milling temperature measurements with a single flute end mill were used as a benchmark case. Thermal material properties and average orthogonal cutting coefficients used in the simulations are listed in Table 5-1 and Table 5-3, respectively. Calibration coefficients were determined from Altintas’ s book [59] for the workpiece material Ti6Al4V. Cutting forces were simulated to compare them against the experimental measurements reported in [19]. Force comparison is shown in Figure 5-13. It indicates that used orthogonal cutting coefficients properties are adequate to use in temperature simulations. 95 Table 5-3- Ti6Al4V orthogonal cutting constants Material Ti6Al4V Shear Angle (o) 34 Friction Angle (o) 29.52 Shear Stress (MPa) 613 Rake Angle (o) 0 Helix Angle (o) 0 Tool Diameter (mm) 52 Spindle Speed (rpm) 1300 Cutting Speed (m/min.) 214 Feed Rate (mm/tooth) 0.105 Radial Immersion (mm) 21 Figure 5-13- The comparison of simulated and measured forces for down milling of Ti6Al4V with V=3.67 m/s and f=0.105 mm/tooth, measurements from Sato et al. [19] Temperature simulations were performed by using 0.05 ms time steps and 10-micron nominal element size. Measurement point was reported to be located 0.1 mm below the rake face of the tool in [19]. The simulated and measured temperature results are given in Figure 5-14. The peak temperature prediction error is less than %10 at the measurement location. It is 96 observed that the maximum value of temperature increases in the beginning of the tool-workpiece engagement, and starts to reduce before the end of the engagement period as the chip thickness reduces. Heat partition ratio towards tool changes between %80 and zero heat flux through tool-workpiece engagement. Since the thermal conductivity of carbide is higher than Ti6Al4V, the heat generated at the tool-chip interface tends to flow into the tool. This effect is more apparent in the beginning of the tool-workpiece engagement. Figure 5-14- Simulated and measured resultant rake face temperature 0.1 mm below the tool-chip interface in down milling of Ti6Al4V, measurements from Sato et al. [19] 97 Figure 5-15- Average heat partition into the tool during down milling of Ti6Al4V 5.3 Validation of variable thermal property approximation Variable thermal property model is applied by considering temperature dependent material properties. An additional convergence criterion added to the simulation algorithm as defined by 2Tolerance in the equation (4.52). If 2Tolerance goes below a predefined value, the solution is accepted. This criterion checks the variation in thermal conduction coefficient between successive iterations. Thermal conduction is selected as a benchmark material property for the convergence criterion due to its primary effect on temperature distribution. Application of variable thermal properties to coefficient and heat source matrix assembly is explained in Section 3.3. Thermal properties are updated for each volumetric element by considering heat conduction towards neighboring element interfaces from the nodal location of interest. Equations (3.31) and (3.32) are used to calculate element to element heat 98 conduction coefficients. Heat capacity and material density values are based on the average value for all elements. One of the disadvantages of utilizing variable thermal properties in simulations is that it increases the cost of computation. Additional convergence requirement dictates an update of coefficient matrices with new thermal properties at each iteration. In addition, temperature-dependent material property use may require a denser meshing and smaller time steps depending on the workpiece material and tool combination. For example, if a low thermal conductivity workpiece material like titanium is simulated, thermal gradient close to heat sources will be steeper compared to an aluminum simulation case. This will require increased mesh density near heat sources. Smaller time steps are required when simulating high thermally conductive workpiece materials due to increased thermal diffusivity. Adjusting convergence rate is important when material properties are updated with changing temperatures. Heat partition calculation gain in the equation (4.50) is generally adjusted below one to approach the solution with smaller tuning steps to avoid overshoots. Overshooting the solution may trigger a convergence instability. Use of temperature dependent material properties increases prediction reliability compared to constant thermal property simulations. Typically, an average temperature (higher than room temperature) is selected for thermal properties when using constant material properties in the simulations. 99 5.3.1 Cutting Tests AISI 316L turning tests were performed in Sandvik Coromant AB facilities in Sandviken, Sweden. Temperature measurements were collected using embedded thermo-couples at two locations within the cutting insert. Tests were carried out on a George Fisher CNC Lathe, NDM-17/125 by using a Sandvik Coromant tool holder type STFCR2525M 16. Orthogonal cutting workpiece specimen was prepared by creating 3 mm thick circular disks with 25 mm depth on AISI316L steel bar (See Figure 5-16). Disk thickness variation was checked against any significant deviation by measuring with a caliper. It was observed that each disk had ± 0.05 mm deviation or less. Custom cutting inserts were prepared based on a standard commercial geometry. Two different uncoated substrate grades of cutting inserts were used in the tests. They are based on Sandvik’s TCMW 16T304 without additional chip breaker on them. The commercial version of these inserts is either coated via using chemical vapor deposition (CVD) or physical vapor deposition (PVD). Insert material composition, hardness and number of inserts instrumented are shown in Table 5-4. Eighty-nine experiments were carried out at different cutting speeds and feeds in total. Cutting speeds were varied from 70 m/min to 240 m/min, and the feed rates were varied between 0.1 mm to 0.3mm/rev. Only a limited portion of the results are shared in this text due to the confidentiality of the company data. 100 Figure 5-16-Orthogonal cutting specimen on the lathe Table 5-4- Insert grade information Carbide Grade 462 (PVD grade) 570 (CVD grade) Composition (% Volume) WC-85, Co-13.2 WC-83.6, Co -14.5 Hardness (Kg.mm2) 1814 1325 Number of Inserts 4 2 101 Inserts were instrumented with two thermocouples embedded under tool rake face. Insert geometry and thermocouple locations are shown in Figure 5-17. Thermocouple locations were verified by cutting inserts into the half on a diamond saw and taking pictures under the microscope. An example measurement is shown in Figure 5-18. It was found that measured thermocouple locations and hole dimensions are comparable to the planned ones. One drawback related to thermocouple hole dimensions is that it is hard to judge thermocouples reached till the end of the holes. Therefore, the average values of temperature around thermocouple tip locations are taken for benchmarks. Figure 5-17-Drawing of insert geometry and thermocouple locations 102 Figure 5-18- Example microscope measurement of an AISI 316L cutting insert (462 grade) Figure 5-19-Insert close-up picture with chip’s relative position Two unground Omega Brand K-type TJC36-CAXL-020U-SMPW-M thermocouples with 500 µm diameter and 150 mm length were fastened through a groove prepared on tool 103 holder with silica gel (See Figure 5-21). About 25 mm of the thermocouples were free to move at their tip, protruding through the insert seat to allow inserts to be attached and removed (See Figure 5-21). Thermocouples are flexible enough that they could reach to the bottom of each hole and then be pressed down through the insert seat to allow the insert to be tightened in place. A thick layer of HTSP50T silicone heat transfer paste by Electrolube was applied to thermocouples before they were mounted each time. An additional CO3-K thermocouple (foil type) was mounted on the underside of the tool holder as well. All three measurement channels were recorded at 100 Hz sampling frequency by using an HDM Spider 8 system (See Figure 5-20). During the cutting tests, disks on the workpiece are aligned about 1-4 mm from the corner of the tool (See Figure 5-19) to make thermocouples stay under the tool-chip contact area. Disks were aligned on approximately top of AA section line in Figure 5-17. Figure 5-20- Data acquisition system set-up for thermal measurements 104 One drawback related to selected sensors is that their time-dependent response is not fast enough to catch temperature rise time during machining because of having 500 µm thermocouple diameter. Another issue is that thermocouples are unground type which slow downs their time response compared to the ground thermocouples. Therefore, the measurements give an idea about temperature values rather than time-dependent response. According to the manufacturer’s data, thermocouple has around one second time constant in the open air. The time constant is expected to reduce when it is installed inside the cutting tool. However, it will be still high to catch the transient response. Selected thermocouples are reliable until 800 oC according to manufacturer’s specification sheet. Thermocouple information can be found in [70]. Figure 5-21-Thermocouples on the tool holder without silica gel 105 In parallel to the temperature measurements, cutting force data were collected during the tests to identify orthogonal cutting coefficients. A Kistler 9263A dynamometer was used to collect cutting forces with a sampling frequency of 1000 Hz. Since cutting force measurements have a dynamic component a 5th order Butterworth filter with 0.01 Hz cut-off frequency was employed in MATLAB environment to get the DC component of the force. Linear measurement drift was compensated after collecting each data set. Shear angle was determined from chip thickness measurements. Chip length and weight were used in the calculation of cut chip thickness. It was assumed that there is no side spread of chips (plain strain assumption). Therefore, chip thickness is determined by using the density of the workpiece material; chipcchip chipmhb l (5.1) where ch , chipm , chipl , b and chip are the thickness, mass, and the length of the collected chip, the width of cut, and density of the workpiece material, respectively. Tests were conducted until the end of tool life for each insert. Visual inspection of tool wear and increase in the cutting force were used to decide end of the tool life. Test duration for each cutting condition case was adjusted to take at least 5 seconds of cutting data to reach stable tool-chip interface temperature (although temperature values keep increasing at the measurement locations) and stable tool-chip contact conditions. In the experiments, depths of craters on rake face reached as high as 300 µm. The deepest craters were observed to be close to thermocouple locations. Although the crater depth was close to thermocouple tips, no damage on the sensors was observed after inspection at the end of tests. 106 Albeit it is not ideal, due to changing cutting edge geometry and friction conditions because of wear, same inserts used for successive tests due to a limited number of inserts. Since the number of inserts was limited, tool edges were worn out after few tests. This caused an increase in edge forces related to ploughing in the tertiary cutting zone. Therefore, the created orthogonal database is not based on the use of fresh tools. Validation cases were picked from the data where identified edge forces are relatively lower. 5.3.2 Temperature Simulations Four simulation cases with temperature dependent material properties are presented in this subsection. They are compared against experimental measurements at two different thermocouple locations. All simulations were carried out with 10 µm nominal element size and 0.05 milliseconds time steps. Note that all the secondary zone energy generation in this subsection is determined by using the sliding contact assumption at the tool-chip interface. The temperature dependent thermal properties are listed in Table 5-5 for AISI316L steel and 462 carbide tool material. Temperature-dependent thermal properties were identified for tool material within Sandvik Coromant AB facilities. Thermal properties for 462-grade is used to simulate both 570 and 462-grade carbide insert materials due to unavailability of data for 570-grade. AISI316L workpiece thermal material properties were taken from [71]. Table 5-5-Tool and workpiece thermal properties for AISI316L simulations Material AISI 316L [71] 462 Grade Uncoated Carbide Thermal conductivity (W/m.K) 214.307+ 0.0181 + 6×10-6TT (-0.0246) + 93.57T Heat capacity (J/Kg.K) 2 3440.79+ 0.5807 + 0.001 + 7×10-7TT T 0.0988 +192.4T Density (Kg/m3) 8000 1.478e+04- 0.318 T 107 Experimentally identified orthogonal cutting coefficients were used in simulations. Since the number of inserts was limited, each insert was utilized until the end of tool life as stated previously. Therefore, tool wear was an important factor that increases ploughing effects, which increased the edge forces. Orthogonal parameters are listed for four validation cases in Table 5-6. Average cutting coefficients were extracted from experimental data for 70, 100, and 180 m/min cutting speeds. Simulations were run at 0.1 and 0.2 mm/rev feed rates. Analysis showed that edge force coefficients are higher than usually expected amounts for a tool with an unworn edge. Although same cutting conditions repeated more than once for the given cases, a limited number of inserts prohibited to use fresh tools for each experiment case. It is seen that the tangential edge force coefficient is higher than the feed direction edge coefficients at 70 m/min cutting speed. Since normal force component on the workpiece surface is feed force, it is typically expected to be higher under sliding dominated contact conditions around the cutting edge. The reason behind this behavior possibly related to tool wear on the flank face of the tool. Besides, the use of only one insert for dry cutting experiments with 570 carbide grade did not allow the same testing conditions on a different sample. The second 570-grade insert was used for the tests at higher cutting speeds with coolant. 108 Table 5-6- Workpiece material orthogonal cutting constants for AISI 316L Case # 1 2 3 4 Shear Angle (o) 24.6 24.6 25.2 25.3 Friction Angle (o) 33.6 33.6 29.6 23.3 Shear Stress (MPa) 569 569 543 579 Rake Angle (o) 0 0 0 0 Tangential Edge Force Coefficient (N/mm) 52 52 69 65 Feed Edge Force Coefficient (N/mm) 48 48 87 103 Cutting Speed (m/min.) 70 70 100 180 Feed Rate (mm/rev.) 0.1 0.2 0.2 0.2 Cumulative Cutting Time for the insert (s) 33.4 13.8 35.8 4.5 Insert Grade 570 570 462 462 Results show that predicted temperatures at steady-state conditions are comparable to the measurements regarding the order of magnitude. However, time-dependent behavior of simulations is faster than the measurement time response. Since model size was limited for temperature calculations, the whole tool holder fixture assembly was not included in simulations. The faster response of simulations is partly due to the size effect. Approximated model size smaller compared to actual cutting tool and tool holder assembly. Although model size is adequate for predicting interface temperatures, current approximation with limited size along with the imposed boundary conditions causes a faster time response at points away from the heat sources in modeled domains. Another cause is that thermocouple diameter is large which hinders fast time response during measurements. Temperature values at thermocouple locations keep increasing within the measurement time in all four validation cases. This behavior can be attributed to ongoing tool wear on rake and flank faces during cutting. It causes an increase in contact surfaces on the tool; hence the frictional energy generation increases. 109 Since the number of inserts were limited in the tests (total nine inserts for eighty nine experiments), tests were run with worn tools. Figure 5-22 shows the average temperature at thermocouple locations compared to measurements at 180 m/min. (3 m/s) cutting speed and 0.2 mm/rev feed rate. This experiment was conducted with a 462-grade fresh insert. Average simulated temperatures along the diameter of the thermocouple are compared against the measurements. The average temperature is taken along the diameter for 150 µm of depth from thermocouple’s tip. Predicted steady state temperature at thermocouple location 1 is 421 oC and measured temperature is 402 oC after 4.5 seconds of cutting. Predicted temperature is 224 oC and thermocouple reads 330 oC at thermocouple location 2. Thermocouple readings keep increasing within the cutting period which suggests a higher actual temperature value compared to the simulation result. Evolution of the average tool-chip interface temperature with time is shown in Figure 5-23. It takes 40 milliseconds to reach %90 of the steady-state value of 707 oC. The rise time at the measurement locations has added delay due to their location. It takes around 350 milliseconds to reach %90 of steady-state value at the thermocouple location 1 in simulations. 560 milliseconds is spent to reach %90 of steady state value at thermocouple location 2 in simulations. The time-dependent response is sensitive to location as a result of the heat transmission related delay. Therefore, it adds up as an error source in the measurements. Thermal gradient plots for tool and chip at steady state condition are presented in Figure 5-25. 110 Figure 5-22-Simulation and thermocouple measurement comparison of AISI 316L turning at V=3 m/s and f=0.2 mm/rev Figure 5-23- Evolution of average rake face temperature with time for turning of AISI 316L at V=3 m/s and f=0.2 mm/rev 111 Figure 5-24- Average heat partition to tool for turning of AISI 316L at V=3 m/s and f=0.2 mm/rev 112 Figure 5-25-Simulated temperature distributions at steady-state for turning of AISI 316L at V=3 m/s and f=0.2 mm/rev Figure 5-26 shows the comparison of simulated and measured temperatures at 100 m/min (1.67 m/s) cutting speed and 0.2 mm/rev feed rate. The measurements after 9.6 seconds cutting indicated 468 oC at thermocouple location 1 whereas the simulated average temperature for the same location is 408 oC. Thermocouple 2 reads 340 oC and simulation result is 220 oC . Figure 5-27 displays the evolution of the average temperature at the tool-chip interface. It takes 113 48 milliseconds to reach %90 of the steady-state value (662 oC) which is higher than 180 m/min case (40 milliseconds Figure 5-24) due to reduced energy input and lower mass transfer rate. Heat partition ratio changes from %30 (180 m/min case) to %41 when cutting speed changes to 100 m/min (see Figure 5-24 and Figure 5-28). This behavior can be linked to the effect of mass transfer on heat conduction. Decreasing cutting speed allows more heat conduction towards the tool. Detailed steady-state temperature gradients of tool and chip are presented in Figure 5-29. Higher cutting speed increases the localization of thermal gradient both on tool and chip which causes a steeper thermal gradient. This behavior is particularly important from thermal stress evolution within the tool body. Steeper thermal gradient causes higher thermal stresses which contribute to crack initiation. Figure 5-26-Simulation and thermocouple measurement comparison of AISI 316L turning at V=1.67 m/s and f=0.2 mm/rev 114 Figure 5-27- Evolution of average rake face temperature with time for turning of AISI 316L at V=1.67 m/s and f=0.2 mm/rev Figure 5-28- Average heat partition to tool for turning of AISI 316L at V=1.67 m/s and f=0.2 mm/rev 115 Figure 5-29-Simulated temperature distributions at steady state for turning of AISI 316L at V=1.67 m/s and f=0.2 mm/rev Figure 5-30 shows the comparison between thermocouple measurements and simulation results at 70 m/min (1.17 m/s) cutting speed and 0.1 mm/rev feed rate. Simulations show 240 oC after 6.4 seconds at steady state for thermocouple location 1 and measurements reach 260.6 oC within the same amount of time. The measured temperature is 199.6 oC for thermocouple 116 location 2 where simulations show 123 oC average temperature for the same condition. Simulated average rake face temperature reaches to 495 oC under steady-state conditions. It takes 30.9 milliseconds to reach %90 of the steady-state rake face temperature. Lower chip thickness results in a faster time response compared to 100 m/min (48 milliseconds), and 180 m/min (40 milliseconds) cutting speed cases with 0.2 mm/rev feed rate. Figure 5-32 presents the evolution of heat partition ratio with time. Heat partition towards tool increases to %53 as a result of reduced cutting speed compared to previous two simulation cases (%30 for 180 m/min. and %43 for 100 m/min.). Tool and chip steady state thermal gradient plots are shown in Figure 5-33. Figure 5-30-Simulation and thermocouple measurement comparison of AISI 316L turning at V=1.17 m/s and f=0.1 mm/rev 117 Figure 5-31- Evolution of average rake face temperature with time for turning of AISI 316L at V=1.17 m/s and f=0.1 mm/rev Figure 5-32- Average heat partition to tool for turning of AISI 316L at V=1.17 m/s and f=0.1 mm/rev 118 Figure 5-33-Simulated temperature distributions at steady state for turning of AISI 316L at V=1.17 m/s and f=0.1 mm/rev The comparison between thermocouple measurement and simulation results at 70 m/min cutting speed and 0.2 mm/rev feed rate is shown in Figure 5-34. Thermocouple 1 reaches to 397 oC after 5.6 seconds where simulation results show 419 oC. Thermocouple 2 measurement is 276.4 oC and simulations gives an average of 234.1 oC after the same amount of cutting time. The simulated average rake face temperature is shown in Figure 5-35 with a steady-state value 119 of 654oC. It took 66.3 milliseconds to reach %90 of the steady-state value on rake face which is longer than the 180 m/min (40 milliseconds in Figure 5-24), and 100 m/min (48 milliseconds in Figure 5-28) cutting speed cases at the same feed rate. Compared to 0.1 mm/rev feed rate case at 70 m/min cutting speed, simulations resulted in slower response due to increased mass in the simulated chip domain. Heat partition ratio, on the other hand, is found to be lower %44 compared to 0.1 mm/rev feed rate case (%53) at the same cutting speed as the mass transfer increases with the feed rate. Simulated thermal gradient plots are shown in Figure 5-34. Figure 5-34-Simulation and thermocouple measurement comparison of AISI 316L turning at V=1.17 m/s and f=0.2 mm/rev 120 Figure 5-35- Evolution of average rake face temperature with time for turning of AISI 316L at V=1.17 m/s and f=0.2 mm/rev Figure 5-36- Average heat partition to tool for turning of AISI 316L at V=1.17 m/s and f=0.2 mm/rev 121 Figure 5-37-Simulated temperature distributions at steady state for turning of AISI 316L at V=1.17 m/s and f=0.2 mm/rev 5.4 Simulations and experimental results with tool coating Tool coating is applied by assigning anisotropic medium properties within the cutting tool. In the model, few layers of cutting tool grid points -that are close to coating thickness- is 122 assigned with coating material properties on top of substrate tool material as shown in Figure 5-38. After assigning the thermal properties, directional heat conduction coefficients are calculated by using the equations (3.31) and (3.32) at each volumetric element. Mesh size within the proximity of coating zone is adjusted in a denser form compared to rest of the tool geometry. Since thermal gradient is steep at the tool-chip contact zone, coarse mesh within the coating zone can cause numerical approximation errors. Discrepancies may lead to a degree where results do not show physical relevance regarding the resultant thermal gradient and heat partition. Figure 5-38- Simulation boundary conditions and coating layer application Tool coating model validations have been carried out via using the available data in the literature. Temperature simulations were run for transient conditions until steady-state is reached in simulated turning cases. Composite-layer and single-layer coating cases are benchmarked against the measurements By M’Saoubi & Chandrasekaran [72] in section 5.4.1. Experiments of Grzesik & Nieslony [17] were used for multi-layer coated tool validations in section 5.4.2. 123 5.4.1 Single-layer coated tool simulations M’Saoubi and Chandrasekaran [72] made infrared measurements on orthogonal turning of SS2541 steel, which has a similar composition to AISI4337, with P20 carbide cutting tools coated with TiN. Coating thickness is reported to be 5µm in their work. Tool geometry consists of 6o rake and 11o clearance angles. Cutting edge radius (sharp tool, denoted with “S” in [72]) was reported to be 2 ± 0.7 µm for uncoated and 5 ± 3 µm for coated tools (denoted with “SC” in [72]), respectively. Simulations were run by using the reported force and chip thickness measurements as inputs to determine orthogonal cutting coefficients (shear angle, friction angle, and shear stress). Then, tool-chip contact length, primary, and secondary zone heat generation rates were determined by using previously identified orthogonal cutting coefficients as described in Section 3.2. Since no force measurement data was reported at multiple feed rates other than 0.1 mm/rev., edge forces could not be identified by using the data provided. Therefore, edge forces were calculated by using slip-line field solution in CUTPROTM Advanced Machining Simulation software developed at UBC. Reported edge radius and workpiece material properties were used to determine the edge cutting forces. Details of the slip-line model used in CUTPROTM can be found in [73]. Infra-red (IR) temperature measurements are reported to be relatively reliable since tool material emissivity ( ) value does not fluctuate within the measurement range. M’Saoubi and Chandrasekaran [72] shared it to be 0.386 0.006748 between 500-1100 oC. Emissivity values change more below 500 oC which can cause a maximum measurement error of %8.5 in emissivity values between 0.3-0.5 for S6 uncoated insert [72]. The measured distribution could also be affected at a minor level by the IR measurement’s spatial resolution. Focal area is 124 reported to be 3.5x2.5 mm2 with 752x582 pixel resolution in [72] which corresponds to 4.7 µm/pixel to 4.3 µm/pixel spatial resolution. This spatial resolution level should be enough to catch gradient changes if there is no other source of error (i.e., measurement noise due to sensor heating, diffraction related to optical elements, sensor size effect). Simulations have been carried out for 3.33 m/s cutting speed, 0.05 mm/rev, 0.1 mm/rev and 0.15 mm/rev feed rates for the coated and sharp uncoated tool cases. Orthogonal cutting parameters of the workpiece-tool couple are listed in Table 5-7. Thermal properties of the cutting tool substrate, coating, and workpiece material are given in Table 5-8. Figure 5-38 shows the coating layer and simulation boundary conditions. 10 µm nominal element size was used in the simulations. 0.05 milliseconds long time step was set in all simulations. Element size perpendicular to rake face is set to 1 µm to represent coating with multiple sets of elements. 1 µm element size was kept within 10 µm thick coating layer in which 5 µm of this zone is assigned with coating material thermal properties while remaining 5 µm was assigned with carbide substrate material properties and used as a transition layer. Note that all the secondary zone energy generation in this subsection is determined by using the sliding contact assumption at the tool-chip interface. Table 5-7- Workpiece material orthogonal cutting constants for SS2541 extracted from [72] Material SS2541 TiN Coated SS2541 Uncoated Shear Angle (o) 27.7 32 Friction Angle (o) 26.4 33.1 Shear Stress (MPa) 648.8 560.6 Rake Angle (o) 6 6 Cutting Speed (m/s) 3.33 3.33 Feed Rate (mm/rev.) 0.05, 0.1, 0.15 0.1 125 Table 5-8-Tool substrate, tool coating and workpiece thermal properties for SS2541 simulations Material TiN Coating [74] P20 Carbide [75] SS2541 [76] Thermal conductivity (W/m.K) 22 67 34.3 Heat capacity (J/Kg.K) 783.4 398 494 Density (Kg/m3) 5420 11200 7850 The use of coarse element size with respect to coating thickness may create numerical convergence problems which was observed while numerically experimenting with the developed model. Since thermal gradient is steeper compared to rest of the tool geometry in the immediate contact region, it is generally a better practice to use a denser mesh in the proximity of tool-chip interface. Element size perpendicular to rake face selected in a way that coating layer is covered with an element stack. Room temperature boundary condition was applied to the far end of the insert geometry as shown in Figure 5-38. Simulation results are compared against experimentally measured maximum temperature values on the rake face. Table 5-9 shows the simulated results for the coated tool, and Table 5-10 lists the uncoated tool simulations. Table 5-9- Comparison of simulated and measured maximum rake face temperatures for SS2541 cutting with TiN coated tool Feed rate (mm/rev.) Simulation (oC) Measurement (oC) [72] Error (%) 0.05 673 790* -17 0.1 833 812 2.6 0.15 992 854 16 126 Table 5-10- Comparison of simulated and measured maximum rake face temperatures for SS2541 cutting with the uncoated tool Feed rate (mm/rev.) Simulation (oC) Measurement (oC) [72] Error (%) 0.1 997 841 18 Maximum rake face temperature for the coated tool predicted with %16 error at 0.15 mm/rev feed rate. Predictions were close to the measurements at the feed rate 0.1 mm/rev with a %2.6 error in maximum rake face temperatures. There was no direct comparison chance due to lack of measurement at 0.05 mm/rev feed rate case due to absence of coated tool measurement data available in [72]. Since coated and uncoated tool temperature measurements are close in terms of maximum rake face temperatures as shown in Table 5-9 and Table 5-10, a comparison is made by using uncoated cutting case measurement for the same condition. Measurement in uncoated 0.05 mm/rev feed rate case is reported to be 790 oC and simulation results show 673 oC (-%17 difference). Generally, it is expected to overpredict temperature compared to measurements due to adiabatic boundaries assumption. However, it is seen that the maximum temperature value was underpredicted for 0.05 mm/rev feed rate case. This can be related to different error source. Lack of chip thickness and force measurements at the specific feed rate to make orthogonal cutting coefficient identification is one of them. Also, the use of uncoated tool increases heat accumulation in tool-chip contact zone compared to the coated tool. This may have an additional contribution to the difference. Tool rake face temperature distribution is compared against the infrared measurement from [72] in Figure 5-39. Temperature distribution data points are extracted from the plot reported in [72]. It is observed that the rake face distribution pattern is different compared to 127 measurements. The possible reason behind the discrepancy can be related to uniform tool-chip contact pressure approximation in the model. In reality, there is a non-uniform pressure distribution which changes approximated contact length. Evolution of average rake face temperature at 0.1 mm/rev. cutting for the coated tool is shown in Figure 5-40. Average tool-chip contact temperature reaches 712oC after 100 milliseconds of cutting. It takes around 4.4 milliseconds to reach %90 of the steady-state temperatures. Heat partition also changes from %53 to % 20 after 100 milliseconds of simulation (see Figure 5-41). Thermal gradient plots are presented in Figure 5-42. Figure 5-39- Rake face temperature simulations compared to infrared measurements in [72] for turning of SS2541 at V=3.33 m/s, f=0.1 mm/rev with TiN coating after 100 ms cutting time 128 Figure 5-40- Evolution of average rake face temperature with time for turning of SS2541 at V=3.33 m/s, f=0.1 mm/rev with TiN coating Figure 5-41- Average heat partition to tool for turning of SS2541 at V=3.33 m/s, f=0.1 mm/rev with TiN coating 129 Figure 5-42-Simulated temperature distributions at steady state condition for turning of SS2541 at V 3.33 m/s, f=0.1 mm/rev with TiN coating after 100 ms cutting time 130 A case at 0.1 mm/rev feed rate is shown in Table 5-10 for benchmark purpose between uncoated and coated tools used under the same cutting condition. M’Saoubi provided additional chip thickness and force measurements for the uncoated carbide privately. Orthogonal coefficients are extracted based on the experimental data provided at different feed rates. Therefore, edge forces were identified based on the experimental data for the simulation of uncoated tool. Note that there is some amount of error related to data digitization from figures in [72] as a contributing factor in all validation cases related to data extraction. Maximum rake face temperature was reported to be 841oC from measurements in [72]. Simulations resulted in 997oC after 100 milliseconds. The difference between the measured and simulated maximum rake face temperature is %18. Figure 5-43 shows rake face temperature distribution comparison with respect to experimental profile. There are differences as in coated tool case. Evolution of the average rake face temperature is shared in Figure 5-44. Average temperature reaches 915oC after 100 milliseconds of cutting that is lower than the coated tool results for the same cutting condition. The main reason can be related to different friction conditions. Friction angle for the uncoated tool is identified as 33.1o, and it is 26.4o for the coated tool. The difference corresponds to %30 change in friction coefficient. Heat partition varies between %59 and %17. It settles slightly lower (about %3) than the coated tool case for the same condition. Usually, it is expected to observe a lower heat partition ratio when using coating. However, different orthogonal parameters are used in these two cases. Thus, the contact length and power inputs are different in two cases. This results in different temperatures and energy removal with mass transfer. Thermal gradient plots are shared in Figure 5-46 for the tool and chip domains after 100 milliseconds cutting time. 131 Figure 5-43- Rake face temperature simulation result for turning of SS2541 at V=3.33 m/s, f=0.1 mm/rev with the uncoated tool after 100 ms cutting time Figure 5-44- Evolution of average rake face temperature with time for turning of SS2541 at V=3.33 m/s, f=0.1 mm/rev with uncoated tool 132 Figure 5-45- Average heat partition to tool for turning of SS2541 at V=3.33 m/s, f=0.1 mm/rev with uncoated tool 133 Figure 5-46-Simulated temperature distributions at steady state condition for turning of SS2541 at V 3.33 m/s, f=0.1 mm/rev with the uncoated tool after 100 ms cutting time A TiN coated tool was simulated at 0.15 mm/rev feed rate with conditions given in Table 5-9. Average interface temperature (see Figure 5-47) increases slightly (from 712oC at 0.1 mm/rev to 845 oC at 0.15 mm/rev) with feed rate due to its effect on heat generation. Feed rate increase has a less profound effect on temperature values when compared to cutting velocity increase. Average heat partition towards tool decreases slightly compared to 0.1 mm/rev case 134 (from %20 see Figure 5-41) due to the increase in mass transfer rate with increased chip thickness. It starts at %56 and goes down to %17.5 (see Figure 5-48) as the simulation progresses. Thermal gradient plots for chip and tool after 100 milliseconds of cutting are shown in Figure 5-49. Figure 5-47- Evolution of average rake face temperature with time for turning of SS2541 at V=3.33 m/s, f=0.15 mm/rev with TiN coating 135 Figure 5-48- Average heat partition to tool for turning of SS2541 at V=3.33 m/s, f=0.15 mm/rev with TiN coating 136 Figure 5-49-Simulated temperature distributions at steady state condition for turning of SS2541 at V 3.33 m/s, f=0.15 mm/rev with TiN coating after 100 ms cutting time Tool rake face temperature for 0.05 mm/rev feed rate case is shown in Figure 5-50. Average rake face temperature within the contact zone goes down to 605oC from 712oC compared to 0.1 mm/rev feed rate coated tool case. Heat partition towards tool increases as a result of reduced mass transfer rate when compared to higher feed rate cases (%20 at 0.1 137 mm/rev see Figure 5-41, and %17.5 at 0.15 mm/rev see Figure 5-48). It varies from %52 to %30 (see Figure 5-51). Generally, it can be said that lower feed rates at the same cutting speed allow a higher percentage of heat to flow into the tool. Thermal gradient plots for tool and chip are shared in Figure 5-52. Figure 5-50- Evolution of average rake face temperature with time for turning of SS2541 at V=3.33 m/s, f=0.05 mm/rev with TiN coating 138 Figure 5-51- Average heat partition to tool for turning of SS2541 at V=3.33 m/s, f=0.05 mm/rev with TiN coating 139 Figure 5-52-Simulated temperature distributions at steady state condition for turning of SS2541 at V 3.33 m/s, f=0.05 mm/rev with TiN coating after 100 ms cutting time 5.4.2 Multi-layer coated tool simulations Grzesik and Nieslony modeled maximum and mean temperature rise at the tool-chip interface for multi-layer coated tools by using an equivalent conductivity approach [17]. Temperature measurements based on the tool-workpiece thermocouple method are presented for 140 the orthogonal turning of AISI 1045 steel in their work. The measurement method was based on changing EMF (Electro-Motive Force) values depending on tool-chip interface temperature. Reported validation tests of their model were for coated inserts with three layers (TiC/Al2O3/TiN), and four layers (TiC/ Ti(C,N) /Al2O3/TiN). Their measurement data is used as a benchmark reference. Grzesik and Nieslony’ s experimental conditions were simulated by using an orthogonal cutting database of UBC Manufacturing Automation Laboratory for AISI 1045 steel workpiece-uncoated carbide tool pair. Although using uncoated tool orthogonal database changes the friction properties between tool and chip, material deformation behavior will be similar in machining with both coated and uncoated tool. Generally, it is expected that uncoated tools experience higher temperature gradients under the same cutting conditions due to the less thermal barrier and changed friction at the tool-chip interface. Similar behavior is reported in experimental results of Grzesik and Nieslony [17]. Their results show that measured temperature values for uncoated P20 carbide tool experiments are higher than both three and four layer coated tool measurements. Some of their experimental results are shared in Table 5-11. It is seen that difference between uncoated and coated tool measured mean temperatures fluctuate roughly between -%6 and -%15 (coated tool measurements showed less mean temperature values) for the three-layer coated tool cases reported in [17]. 141 Table 5-11- Example temperature measurements from [17] Workpiece Material/Tool Combination AISI 1045/ P20 Carbide Uncoated AISI 1045/ TiC/Al2O3/TiN Coated P20 Carbide AISI 1045/ TiC/ Ti(C,N) /Al2O3/TiN Coated P20 Carbide Measured Mean Interface Temperature (oC) 709 616 568 Rake Angle (o) 0 0 0 Cutting Speed (m/s) 2.4 2.4 2.4 Feed Rate (mm/rev.) 0.16 0.16 0.16 Only three layer coated tool cases have been simulated at three different cutting speeds (2.4 m/s, 2.07 m/s, and 1.48 m/s) and 0.16 mm/rev feed rate. Orthogonal cutting coefficients and corresponding cutting conditions are shown in Table 5-12. Table 5-12- Workpiece material orthogonal cutting constants for AISI 1045 taken from CUTPROTM Workpiece Material/Tool Combination AISI 1045/ Uncoated carbide AISI 1045/ Uncoated carbide AISI 1045/ Uncoated carbide Shear Angle (o) 29.6 29.2 28.4 Friction Angle (o) 24.1 24.8 25.9 Shear Stress (MPa) 544.5 536.6 522.3 Rake Angle (o) 0 0 0 Cutting Speed (m/s) 2.4 2.07 1.48 Feed Rate (mm/rev.) 0.16 0.16 0.16 10 µm nominal element size and 0.05 milliseconds time steps were used in the simulations. 1 µm element size in perpendicular to rake face was used close to tool-chip interface. Coating thickness values for each layer was reported to be 6 µm, 3 µm, 1 µm for TiC, Al2O3, and TiN respectively with a total of 10 µm. Another five microns of material is meshed with 1 µm element size perpendicular to rake face within tool substrate to approximate material property changes gradually. Room temperature boundary condition was applied to the back of the insert 142 geometry. Simulations were run for 100 milliseconds long which was enough to reach steady state conditions at the tool-chip interface. Workpiece temperatures were not calculated. Note that all the secondary zone energy generation in this subsection is determined by using the sliding contact assumption at the tool-chip interface. Table 5-13-Tool substrate, tool coating and workpiece thermal properties for AISI 1045 simulations Material TiN Coating [74] Al2O3 Coating [75] TiC Coating [77] P20 Carbide [75] AISI 1045 [78] Thermal conductivity (W/m.K) 22 16 27 67 35.1 Heat capacity (J/Kg.K) 783.4 100 520 398 773 Density (Kg/m3) 5420 3980 4900 11200 7800 Predicted and measured mean tool-chip interface temperature values are compared in Table 5-15 after 100 milliseconds of cutting time. It is seen that the prediction error varies between -%11 to -%22 in benchmarked cases. Discrepancies between predictions and measurements can be related to various sources. Orthogonal cutting constants that are used in the simulations are for the uncoated carbide tool-AISI 1045 pair. Identified cutting constants can vary due to microstructural effects of workpiece material production steps (i.e., heat treatment, hardness variations). Use of constant thermal properties regardless of changing temperature is another source of error. AISI 1045 thermal conductivity decreases as the temperatures increase which forces more of heat flow towards the cutting tool. Tungsten carbide substrate, TiN coating layer, and TiC coating layer thermal conductivities increase with temperature which eases heat flow towards the tool as the temperatures increase. Al2O3 thermal conductivity 143 decreases as the temperature increases; then stays relatively stable after 750 oC. Model assumptions related to cutting mechanics can also contribute to the discrepancy. An example case result with variable thermal properties is shown in Table 5-16 where uncoated carbide and AISI 1045 were assigned with temperature dependent properties. Al2O3 and TiC were assigned with constant properties due to lack of temperature dependent heat capacity data, although their thermal conductivity information was found in [79]. TiN variable thermal property data was available from [74]. Simulated temperature results become closer to experimental measurements for the presented case when variable material properties were used. Table 5-14-Variable tool substrate and workpiece thermal properties for AISI 1045 simulations Material AISI 1045 [80] (up to 1000 oC, T is in oC) Uncoated Carbide [80] (up to 1000 oC, T is in oC) Thermal conductivity (W/m.K) 62.23- 0.0312 T (0.042) + 24.48T Heat capacity (J/Kg.K) 2282.33+ 0.504 + 1.04×10-18TT 2301.23+ 0.12 + 4.2×10-19TT Density (Kg/m3) 7844 11900 Table 5-15- Comparison of simulated and measured mean rake face temperatures for AISI1045 cutting with TiC/Al2O3/TiN coated tool Cutting speed (m/s) Simulation (oC) Measurement (oC) [17] Error (%) 2.4 548 616 -11 2.07 516 605 -14 1.48 450 582 -22 144 Table 5-16- Comparison of simulated and measured mean rake face temperatures for AISI1045 cutting with TiC/Al2O3/TiN coated tool (Temperature dependent material properties simulation result) Cutting speed (m/s) Simulation (oC) Measurement (oC) [17] Error (%) 1.48 562 582 -3 Average tool-chip interface temperatures increase with cutting speed gradually in simulations and experimental results [17] (see Table 5-15). Average temperature evolution with time is presented in Figure 5-53 for 2.4 m/s, Figure 5-56 for 2.07 m/s and Figure 5-59 for 1.48 m/s. Heat partition ratios change as the cutting speed changes due to change in mass transfer rate. Heat partition for 2.4 m/s cutting speed was simulated to be %14.6 after 100 milliseconds of cutting time (see Figure 5-54), and it slightly increases to %15.8 (see Figure 5-57) for 2.07 m/s. It reaches to %19.2 (see Figure 5-60) for 1.48 m/s case. Detailed thermal gradients are presented in Figure 5-55, Figure 5-58, and Figure 5-61. 145 Figure 5-53- Evolution of average rake face temperature with time for turning of AISI 1045 at V=2.4 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating Figure 5-54- Average heat partition to tool for turning of AISI 1045 at V=2.4 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating 146 Figure 5-55-Simulated temperature distributions at steady state condition for turning of AISI 1045 at V=2.4 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating after 100 ms cutting time 147 Figure 5-56- Evolution of average rake face temperature with time for turning of AISI 1045 at V=2.07 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating Figure 5-57- Average heat partition to tool for turning of AISI 1045 at V=2.07 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating 148 Figure 5-58-Simulated temperature distributions at steady state condition for turning of AISI 1045 at V=2.07 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating after 100 ms cutting time 149 Figure 5-59- Evolution of average rake face temperature with time for turning of AISI 1045 at V=1.48 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating Figure 5-60- Average heat partition to tool for turning of AISI 1045 at V=1.48 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating 150 Figure 5-61-Simulated temperature distributions at steady state condition for turning of AISI 1045 at V=1.48 m/s, f=0.16 mm/rev with TiC/Al2O3/TiN coating after 100 ms cutting time 151 5.5 Comparison between static contact and sliding heat partition assumptions Heat partition ratio in sliding systems is different from the static contact cases. Partition ratio is sometimes approximated by using a ratio of thermal effusivities, via static contact assumption, given in equation (4.38). This subsection adds a short discussion about the comparison between the static contact and Blok’s sliding contact assumptions used in the simulation cases presented in section 5.2, section 5.3, and section 5.4. There are experimental and theoretical evidences in the literature that indicate the importance of heat partition approximation. Heat partition ratio experimentally evaluated by Rech et al. [81] showed speed dependent behavior on pin sliding experiments on Ti6Al4V and couple of steel alloys with coated carbide pins. Heat partition ratios drop as the cutting speed increases in their results. Similar behavior is concluded by Abukhshim et al. [80] who did cutting tests on ASIS 4140 steel with uncoated carbide tools and measured temperature values by using a pyrometer. They did reverse analysis from experimental results and recreated the thermal loads by using finite element method. They observed speed dependent behavior as well. However, heat partition ratios dropped until the speed is increased to 600 m/min., and the tool-chip contact area increased significantly at the higher speeds which indicated the increased heat partition ratios. This behavior can be related to material softening at high cutting speeds and increased chip sticking to the tool which contributes to heat transmission. Table 5-17 shows the calculated values of heat partition to the tool by using the static contact assumption for constant thermal property validation cases (Section 3.3.1). These values are above the simulated results for a metal cutting case. As the cutting speed increases, the mass transfer effect on thermal behavior increases its effect. Mass transfer changes the energy balance between tool and chip via removing heat. Therefore, the use of static contact assumption 152 becomes less reliable. Contact body geometries, contact area, and thermal material properties also contribute to speed-dependent behavior. Heat flows more towards the direction, where there is lower thermal resistance. Table 5-17- Constant thermal property static contact heat partition ratio approximation by using the equation (4.38) Case Speed (m/min.) Feed Rate (mm/rev.) or (mm/tooth) Static Contact Heat Partition Ratio (%) Simulated Sliding Contact Heat Partition Ratio (%) Ti6Al6V-2Sn-Uncoated Carbide, Interrupted Turning 100 0.1 75 %51 Ti6Al6V-2Sn-Uncoated Carbide, Continous Turning 100 0.1 75 %42 Al2024-Uncoated Carbide, Interrupted Turning 81.6 0.109 32.8 %10 Al2024-Uncoated Carbide, Interrupted Turning 81.6 0.165 32.8 %10 Gray Cast Iron-Uncoated Carbide, Interrupted Turning 81.6 0.109 39 %15 Ti6Al4V-Uncoated Carbide, Down Milling 214 0.105 76 %80 to 0 The heat partition simulations for turning AISI316L with variable thermal properties (Section 5.3.2) are presented in Table 5-18. The equation (4.38) is used to determine the heat partition values with static contact assumption. They were determined by considering the local average temperature at the tool-chip interface. Although temperature dependent thermal properties were used in the calculations, the mass transfer effect of the chip forces a different thermal equilibrium compared to the static contact case. The estimated heat partition values are close to each other because the material properties dependent calculation in equation (4.38). Slight differences are due to temperature variation’s effect on material properties. Speed 153 dependence of sliding contact assumption can be seen in Figure 5-62. There is a linear relationship between heat partition change and cutting speed because of mass transfer rate’ s cutting speed dependence. Table 5-18- Static contact assumption heat partition ratios towards tool for AISI 316L simulation cases Speed (m/min.) Feed Rate (mm/rev.) Static Contact Heat Partition Ratio (%) Simulated Sliding Contact Heat Partition Ratio (%) 70 0.1 64.5 %53 70 0.2 63.4 %44 100 0.2 63.4 %41 180 0.2 62.9 %30 Figure 5-62-AISI316L heat partition change with cutting speed in sliding contact conditions Heat partition predictions in turning SS2541 with TiN coated tool are listed in Table 5-19 which are extracted from the simulation results in Section 5.4.1. Static contact heat partition values are calculated from equation (4.38). The static contact approximations are 154 higher than the sliding contact assumption results at 0.05 mm/rev feed rate. The TiN coated tool heat partition result (%30 towards tool) is closest to static approximation due to reduced mass transfer effect. Figure 5-63 shows simulation results of sliding heat partition ratio change with the feed rate. The increase of removed mass with higher feed rates results in a reduction of heat partition ratio towards the tool. Table 5-19- Static contact assumption heat partition ratios towards tool for TiN coated and uncoated tool simulation cases Tool Type Speed (m/min.) Feed Rate (mm/rev.) Static Contact Heat Partition Ratio (%) Simulated Sliding Contact Heat Partition Ratio (%) TiN Coated Tool 200 0.05 45.6 % 30 TiN Coated Tool 200 0.1 45.6 % 20 TiN Coated Tool 200 0.15 45.6 % 17.5 Uncoated Tool 200 0.1 60 % 17 155 Figure 5-63- SS2541 heat partition change with feed rate in sliding contact conditions The heat partition predictions in AISI 1045 with three-layer coated tool are shown in Table 5-20 which are extracted from the simulation results in Section 5.4.2. Static heat partition values were calculated from equation (4.38). Only the properties top coating layer (TiN) were considered in the calculation. The static contact heat partition ratio is %39.9 towards the tool (see Table 5-20) which is higher than the simulation results. Heat partition change with cutting speed is presented in Figure 5-64. Increased mass transfer’s effect on the heat partition ratio is observed in the simulated results. Table 5-20- Static contact assumption heat partition ratios towards TiC/Al2O3/TiN coated tool for AISI1045 simulation cases Speed (m/min.) Feed Rate (mm/rev.) Static Contact Heat Partition Ratio (%) Simulated Sliding Contact Heat Partition Ratio (%) 89 0.16 39.9 % 19.2 124 0.16 39.9 % 15.8 144 0.16 39.9 % 14.6 156 Figure 5-64- AISI1045 heat partition change with cutting speed in sliding contact conditions Generally, it can be concluded that the reduction in mass transfer rate makes sliding contact heat partition results to get closer to the static contact assumption. Mass transfer effect becomes more apparent when thermal diffusivity of the workpiece material is low such as Titanium and Inconel. 5.6 Simulation input parameter sensitivity results This subsection is intended to show specific input sensitivity behavior of the numerical model under steady-state conditions for two example cases. Ti6Al4V workpiece material is used in all simulations except one case with AISI316L. Variable workpiece and tool material properties were used in most of the simulations. 10 µm element size was used in all simulations. Orthogonal cutting coefficients for Ti6Al4V were taken from [59], and the material data is shown in Table 5-21. Shear angle can be determined from the given chip ratio by using equation 157 (5.2). Tool geometry was assumed sharp without any wear and rake angle was set to 0o in all simulations. Table 5-21- Ti6Al4V orthogonal cutting constants for sensitivity simulations Material Ti6Al4V [59] Chip Ratio 11.1.755 0.0280.337 0.0082Coo nnC hCC Friction Angle (o) 19.1 0.29 (deg)n Shear Stress (MPa) 613 Rake Angle (o) 0 1costan1 sin,:c nnc ncrrWherer chip ratio (5.2) Thermal properties of Ti6Al4V are listed in Table 5-22. Cutting tool material properties are the same with uncoated carbide material properties given in Table 5-5. Table 5-22- Workpiece thermal properties for Ti6Al4V Material Ti6Al4V [68] Thermal conductivity (W/m.K) 0.0157 5.52, 9900.0124 7.33, 990ooT forT CT for T C Heat capacity (J/Kg.K) 0.215 541.42, 9900.18 461.48, 990ooT forT CT for T C Density (Kg/m3) 4420 - T-25 0.154 The simulations with variable and constant thermal properties are compared in Figure 5-65. Constant thermal properties were used by assuming to have 400 oC average temperature in 158 tool and chip. Maximum and mean tool-chip interface temperatures were calculated for various cutting speeds from 20 m/min to 110 m/min. The difference between two approaches increased as the cutting speed and temperature values increase (i.e., material properties at 400 oC represents heat conduction less as the temperature increases). Difference reached to +%5.75 (Tmax difference) and +%11.6 (Tavg difference) at 110 m/min. Similarly, difference was calculated as -%3.8 (Tmax difference) and -%2.2 (Tavg difference) at 20 m/min. Note that calculated difference is referenced to simulation result for the same speed with variable thermal properties. Figure 5-65- Variable and constant thermal properties' effect on tool-chip interface temperature comparison for Ti6Al4V Sticking and sliding contact assumption at the tool chip-interface affects the predicted temperatures. A comparison case at two different cutting speeds for Ti6Al4V is presented in Figure 5-66. Fully sliding, fully sticking, and %80 sticking contact cases were simulated at 75 159 m/min, 90 m/min cutting speeds and 0.1 mm/rev feed rate. Contact length and average shear stress were determined as described in section 3.2.2. Secondary zone thickness was kept at %8 of the feed rate value which is 8 µm. Fully sticking contact resulted in lower temperature values (in terms maximum tool-chip interface temperatures) compared to fully sliding and %80 sticking contact cases. This is due to the redistribution of frictional power through differences in contact lengths and internal volumetric heat generation. Internal heat generation lets some of the generated heat carried via mass transfer as the chip is removed. Average tool-chip interface temperature results show that %80 sticking and fully sticking contact cases are close to each other (914 oC for %80 sticking and 934 oC for fully sticking contact at 90 m/min, 869 oC for %80 sticking and 889 oC for fully sticking contact at 75 m/min). Maximum tool-chip interface temperatures of %80 sticking and full sliding contact are close to each other (1183 oC for %80 sticking and 1147 oC for fully sliding contact at 90 m/min, 1098 oC for %80 sticking and 1079 oC for fully sliding contact at 75 m/min). These are only simulation results with simplified assumptions because actual interface pressure and shear stress distribution are not known. The intention is to give an example about the differences between sliding contact and sticking contact assumptions. Note that orthogonal cutting theory based on Merchant’s work [82] gives an average friction coefficient which is higher than the normal sliding conditions. Therefore, theoretical friction value is not based on regular sliding contact. 160 Figure 5-66- Sticking and sliding contact’s effect on tool-chip interface temperature comparison for Ti6Al4V Sticking zone thickness’s effect on tool-chip interface temperatures was simulated for Ti6Al4V at two different cutting speeds (45 m/min and 90 m/min) and 0.1 mm/rev feed rate. Sticking zone thickness values changing from %4 to %20 of the depth of cut are presented in Figure 5-67. It is observed that the sticking zone thickness has a minor impact on the predicted interface temperatures for the simulated cases with varying deformation zone thicknesses. Maximum temperature difference is observed to be close to %9.7 at 90 m/min at the simulated sticking zone thickness range. 161 Figure 5-67- Sticking zone thickness’s effect on tool-chip interface temperature comparison for Ti6Al4V Up to this section, convection cooling has been ignored within all presented validations. An example is presented to show the effect of ignored convection cooling. A Ti6Al4V case is simulated at 90 m/min cutting speed and 0.1 mm/rev feed rate. Convection coefficient is artificially changed between 1000 W/m2.K and 10000 W/m2.K. Since model size is limited, the convection area is small compared to the standard cutting tool and tool-holder assembly. Therefore, the convection coefficient was increased to outweigh model size effect. Change of heat that flows into the tool and average tool-chip interface temperature are shown in Figure 5-68 with respect to changing convection coefficients. The amount of heat goes into the tool increased as the convection coefficient increases in simulations. Ten times increase in convection coefficient resulted in %18 rise in transmitted heat (from 79.5 Watts to 94.2 Watts). Average tool-chip interface temperature dropped with an increase in convection coefficient from 924 oC to 887 oC (%4 reduction). As long as the tool-chip interface friction conditions are 162 unchanged, interface temperature at steady state is not affected significantly from convection in the simulations. Figure 5-68- Convection coefficient sensitivity of the simulations for Ti6Al4V Model size affects the tool-chip interface temperatures. One case was simulated with the same material properties in Table 5-5 for turning ASIS 316L at 100 m/min cutting speed and 0.1 mm/rev feed rate. Tool-chip interface temperature results are shown in Figure 5-69. An initial model with 2 mm by 2 mm size was selected for the tool. Then, the tool-chip interface temperatures were compared at different model sizes (between %25 and %250 of the initial model size.). Results show that tool-chip interface temperatures converge to one distribution as the model size is increased. This example is set to represent the effect of approximated model size on temperature gradients. Simulated domain sizes should always be selected carefully in order to get numerically converged solution at the deformation zones. 163 Figure 5-69- Model size effect on rake face temperature comparison for AISI 316L 164 Chapter 6: Conclusion and Future Research Directions 6.1 Conclusions Prediction of temperature distribution in the chip, machined surface, and the tool are important in designing cutting tool geometry, selecting cutting conditions that would not lead to accelerated tool wear and chipping, and predicting the residual stresses left on the finished surface. The cutting operations are classified either as continuous like in turning where the chip thickness remains unchanged within short time intervals or intermittent like in milling where the chip thickness periodically varies as the tool rotates. This thesis presents a comprehensive numerical model of temperature distribution for both continuous and intermittent metal cutting operations carried out with coated or uncoated tools. The temperature distribution in the chip, cutting tool, and the workpiece surface is predicted in two and three-dimensional space. The heat is generated in primary, secondary and tertiary deformation zones in cutting. The chip is sheared from the workpiece stock in the primary shear zone, and it moves on the rake face under sticking and sliding contact conditions in the secondary zone. The cutting-edge ploughs the workpiece surface in the tertiary zone. By assuming a thin shear zone model, the cutting forces and the resulting power that generates the heat are predicted based on the principles of metal cutting mechanics in the thesis. The three deformation zones meshed into finite volume elements. The heat is generated at the thin shear, tool-chip contact, and the tool edge-workpiece surface contact zones. The heat transfer from these zones to discrete elements in the cutting zones is solved numerically via the finite difference method. The time-dependent temperature terms are approximated with second-order backward implicit time discretization. The numerical stability is ensured by using the implicit solution of equations which contain scarcely distributed Eigenvalues of the tool, chip and workpiece coefficient matrices. 165 The model is able to handle both constant chip load (i.e. continuous machining) which generates time-invariant heat at the primary zone, as well as periodically fluctuating (i.e. intermittent) chip loads which lead to time-varying heat source. The anisotropic material properties are incorporated by varying them in three orthogonal directions at each nodal location. Therefore, temperature dependent material properties and coating layers on the cutting tool surface can be included in the simulations. The model allows both sticking and sliding friction conditions at the tool-chip contact zone. The predicted temperature distribution in the cutting zone is compared against the measurements found in the literature and experiments jointly conducted at the industrial research partner’s laboratory. In general, the predicted temperatures have about 20% discrepancy with the measurements which also have inherent errors due to the physical limitation of sensors which need to be placed away from the very small, actual cutting zone. The validations have been carried out for both continuous and intermittent cutting operations. While the temperature reaches at the steady state, constant value in continuous cutting, the temperature is always in a transient mode in interrupted cutting due to changes in the heat source amplitudes as a function of time. The steady-state maximum temperature values are important to select feed and speed that do not violate the tool’s thermal load limits. The transient changes in the temperature are essential to adjust the tool-workpiece contact periods, i.e. radial engagement and cutting speed, without overloading the tool thermally. Transient temperature predictions are also important in predicting cyclic thermal loading to predict the residual stresses left on the workpiece surface. The proposed model is incorporated to a MATLAB code which can predict the temperature in two and three dimensions. 166 The contributions of the thesis can be highlighted as follows : A novel, adaptable and comprehensive temperature prediction engine has been developed using the finite difference method specific to metal cutting. The model can handle both isotropic and temperature dependent anisotropic material properties. Coating material with varying properties at the tool-chip contact zone is also included. Tool coating formulation is embedded in the heat balance equations, which allows the analysis of different coating combinations in tool design. Time-varying heat sources as the chip thickness change with tool rotation are also incorporated into the model. In short, the proposed temperature prediction model is the most comprehensive finite difference based solution found in the literature, which has significantly shorter computation time than commonly used computationally intensive Finite Element Method based material deformation solvers. It is shown that the sliding and sticking contact zones affect the heat partition between the chip and tool significantly. The shorter sticking zone leads to the higher amount of heat transmitted into the tool, hence maximum and average rake face temperatures increase. Increased sticking zone length aids removal of heat with the disposed chip in the simulated cases. It is desired to keep the heat in the disposed chip as opposed to the cutting tool in order to prolong the tool life. The model predicts both the maximum and the transient temperature history. This allows the tool design engineer to test different substrate and coating material combinations cutting various materials at different speeds 167 and feeds. The model also allows the planner to select a feed and radial depth of cut which does not violate the thermal load limits of the tool and coating materials. Two and three-dimensional temperature prediction results are found to be very close to each other because the metal cutting process can be approximated as a plain strain deformation. Therefore, the computationally efficient two-dimensional temperature prediction models can be used in industrial applications. The transient thermal behavior of the machining operations is dependent on both thermal material properties and the cutting conditions (i.e., feeds and speeds). Reaching steady-state conditions takes longer when cutting speed is used due to the reduced mass transfer rate. Feed rate has an inversely proportional effect on rise-time to steady-state, lower feed rates lead to faster response because of reduced chip mass. Thermal diffusivity of the workpiece material has a direct effect on the time-dependent response. Lower thermal diffusivity (i.e. Titanium) results in longer rise-times compared to high thermally conductive materials (i.e. Aluminum) under the similar cutting conditions. Heat partition values are particularly affected by the mass transfer rate which is related to feed rate and cutting speed. It is demonstrated that there is an inversely proportional relation between heat partition and the cutting conditions (i.e., feeds and speeds). Increased feed rate and cutting speed cause a reduction in the percentage of heat that flows into the tool that conforms to the experimental observations in the literature. 168 6.2 Future research directions Prediction of temperature distribution in the tool and workpiece surface is the fundamental first step to conduct the following research in the future: Cyclic temperature variations on the workpiece surface can be integrated into cyclic elastic-plastic deformation model to predict the residual stresses left on the machined part. Residual stresses are correlated to fatigue failure of parts which are critical especially in the aerospace industry. Thin-walled aerospace parts are distorted severely after machining. Residual stresses can be used to predict and constrain the part distortions as a function of cutting speed, feed, temperature and tool geometry. Tool crater wear is dependent on the temperature variation on the rake face, while the flank wear is dependent on the temperature of tool’s flank face. Prediction of temperature is the core requirement to model the tool wear. Currently, torque, power, force, vibration and dimensional surface error limits are used in constraint-based process planning of machining operations. A database can be constructed to predict the maximum temperature as a function of speed, radial engagement, and feed for different tool and coating materials. Temperature constraint can be included in process planning of machining operations to prolong the tool life. 169 References [1] C. Islam, I. Lazoglu, and Y. Altintas, "A three-dimensional transient thermal model for machining," Journal of Manufacturing Science and Engineering, vol. 138, p. 021003, 2016. [2] J. C. JAEGER, "Moving sources of heat and the temperature of sliding contacts," in Proc. Roy. Soc., NSW, 1942, p. 203. [3] E. G. Loewen and M. C. Shaw, "On the Analysis of Cutting Tool Temperatures," Transactions of ASME, vol. 76, pp. 217-231, 1954. [4] B. Chao and K. Trigger, "Temperature distribution at the tool chip interface in metal cutting," Transactions of ASME, vol. 77, pp. 1107-1121, 1955. [5] P. K. Venuvinod and W. S. Lau, "Estimation of rake temperatures in free oblique cutting," International Journal of Machine Tool Design and Research, vol. 26, pp. 1-14, 1986. [6] D. A. Stephenson and A. Ali, "Tool temperatures in interrupted metal cutting," Journal of Engineering for Industry, vol. 114, pp. 127-136, 1992. [7] H. S. Carslaw and J. C. Jaeger, "Conduction of heat in solids," Oxford: Clarendon Press, 1959, 2nd ed., 1959. [8] D. A. Stephenson, T. C. Jen, and A. S. Lavine, "Cutting Tool Temperatures in Contour Turning: Transient Analysis and Experimental Verification," Journal of Manufacturing Science and Engineering, vol. 119, pp. 494-501, 1997. [9] R. Radulescu and S. G. Kapoor, "An Analytical Model for Prediction of Tool Temperature Fields during Continuous and Interrupted Cutting," Journal of Engineering for Industry, vol. 116, pp. 135-143, 1994. [10] R. Komanduri and Z. B. Hou, "Thermal modeling of the metal cutting process Part I : Temperature rise distribution due to shear plane heat source," International Journal of Mechanical Sciences, vol. 42, pp. 1715-1752, 2000. [11] R. Komanduri and Z. B. Hou, "Thermal modeling of the metal cutting process Part II: temperature rise distribution due to frictional heat source at the tool–chip interface," International Journal of Mechanical Sciences, vol. 43, pp. 57-88, 2001. [12] R. Komanduri and Z. B. Hou, "Thermal modeling of the metal cutting process Part III: temperature rise distribution due to the combined effects of shear plane heat source and the tool–chip interface frictional heat source," International Journal of Mechanical Sciences, vol. 43, pp. 89-107, 2001. [13] Y. Huang and S. Y. Liang, "Modelling of the cutting temperature distribution under the tool flank wear effect," Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 217, pp. 1195-1208, 2003. [14] K.-M. Li and S. Y. Liang, "Modeling of Cutting Temperature in Near Dry Machining," Journal of Manufacturing Science and Engineering, vol. 128, pp. 416-424, 2006. [15] Y. Karpat and T. Özel, "Analytical and thermal modeling of high-speed machining with chamfered tools," Journal of Manufacturing Science and Engineering, vol. 130, p. 011001, 2008. 170 [16] Y. K. Chou and H. Song, "Thermal modeling for white layer predictions in finish hard turning," International Journal of Machine Tools and Manufacture, vol. 45, pp. 481-495, 2005. [17] W. Grzesik and P. Nieslony, "Physics based modelling of interface temperatures in machining with multilayer coated tools at moderate cutting speeds," International Journal of Machine Tools and Manufacture, vol. 44, pp. 1451-1462, 2004. [18] S. Zhang and Z. Liu, "An analytical model for transient temperature distributions in coated carbide cutting tools," International Communications in Heat and Mass Transfer, vol. 35, pp. 1311-1315, 2008. [19] M. Sato, N. Tamura, and H. Tanaka, "Temperature Variation in the Cutting Tool in End Milling," Journal of Manufacturing Science and Engineering, vol. 133, pp. 021005-021005, 2011. [20] D. J. Richardson, M. A. Keavey, and F. Dailami, "Modelling of cutting induced workpiece temperatures for dry milling," International Journal of Machine Tools and Manufacture, vol. 46, pp. 1139-1145, 2006. [21] J. Liu, C. Ren, X. Qin, and H. Li, "Prediction of heat transfer process in helical milling," International Journal of Advanced Manufacturing Technology, vol. 72, pp. 693-705, 2014. [22] A. O. Tay, M. G. Stevenson, and G. D. V. Davis, "Using the finite element method to determine temperature distributions in orthogonal machining," ARCHIVE: Proceedings of the Institution of Mechanical Engineers 1847-1982 (vols 1-196), vol. 188, pp. 627-638, 1974. [23] A. Tay, M. Stevenson, G. de Vahl Davis, and P. Oxley, "A numerical method for calculating temperature distributions in machining, from force and shear angle measurements," International Journal of Machine Tool Design and Research, vol. 16, pp. 335-349, 1976. [24] P. D. Muraka, G. Barrow, and S. Hinduja, "Influence of the process variables on the temperature distribution in orthogonal machining using the finite element method," International Journal of Mechanical Sciences, vol. 21, pp. 445-456, 1979. [25] P. R. Dawson and S. Malkin, "Inclined Moving Heat Source Model for Calculating Metal Cutting Temperatures," Journal of Engineering for Industry, vol. 106, pp. 179-186, 1984. [26] T. H. C. Childs, K. Maekawa, and P. Maulik, "Effects of coolant on temperature distribution in metal machining," Materials Science and Technology, vol. 4, pp. 1006-1019, 1988. [27] T. Ueda, M. Sato, and K. Nakayama, "The Temperature of a Single Crystal Diamond Tool in Turning," CIRP Annals - Manufacturing Technology, vol. 47, pp. 41-44, 1998. [28] J. S. Strenkowski and K.-J. Moon, "Finite Element Prediction of Chip Geometry and Tool/Workpiece Temperature Distributions in Orthogonal Metal Cutting," Journal of Engineering for Industry, vol. 112, pp. 313-318, 1990. [29] Z. C. Lin and S. Y. Lin, "A Coupled Finite Element Model of Thermo-Elastic-Plastic Large Deformation for Orthogonal Cutting," Journal of Engineering Materials and Technology, vol. 114, pp. 218-226, 1992. [30] A. J. Shih, "Finite Element Simulation of Orthogonal Metal Cutting," Journal of Engineering for Industry, vol. 117, pp. 84-93, 1995. 171 [31] J.-S. Wu, O. W. Dillon, and W.-Y. Lu, "Thermo-Viscoplastic Modeling of Machining Process Using a Mixed Finite Element Method," Journal of Manufacturing Science and Engineering, vol. 118, pp. 470-482, 1996. [32] T. Marusich and M. Ortiz, "Modelling and simulation of high‐speed machining," International Journal for Numerical Methods in Engineering, vol. 38, pp. 3675-3694, 1995. [33] M. R. Movahhedy, Y. Altintas, and M. S. Gadala, "Numerical Analysis of Metal Cutting With Chamfered and Blunt Tools," Journal of Manufacturing Science and Engineering, vol. 124, pp. 178-178, 2002. [34] M. R. Movahhedy, M. S. Gadala, and Y. Altintas, "Simulation of chip formation in orthogonal metal cutting process: an ALE finite element approach," Machining Science and Technology, vol. 4, pp. 15-42, 2000. [35] E. Ng, D. K. Aspinwall, D. Brazil, and J. Monaghan, "Modelling of temperature and forces when orthogonally machining hardened steel," International Journal of Machine Tools and Manufacture, vol. 39, pp. 885-903, 1999. [36] T. Özel and T. Altan, "Process simulation using finite element method — prediction of cutting forces, tool stresses and temperatures in high-speed flat end milling," International Journal of Machine Tools and Manufacture, vol. 40, pp. 713-738, 2000. [37] S. L. Soo, D. K. Aspinwall, and R. C. Dewes, "Three-dimensional finite element modelling of high-speed milling of Inconel 718," Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture, vol. 218, pp. 1555-1561, 2004. [38] G. M. Pittalà and M. Monno, "A new approach to the prediction of temperature of the workpiece of face milling operations of Ti-6Al-4V," Applied Thermal Engineering, vol. 31, pp. 173-180, 2011. [39] Y. C. Yen, A. Jain, P. Chigurupati, W. T. Wu, and T. Altan, "Computer simulation of orthogonal cutting using a tool with multiple coatings," Machining Science and Technology, vol. 8, pp. 305-326, 2004. [40] P. J. Arrazola, I. Arriola, and M. A. Davies, "Analysis of the influence of tool type, coatings, and machinability on the thermal fields in orthogonal machining of AISI 4140 steels," CIRP Annals - Manufacturing Technology, vol. 58, pp. 85-88, 2009. [41] A. W. Nemetz, W. Daves, T. Klünsner, W. Ecker, T. Teppernegg, C. Czettl, et al., "FE temperature- and residual stress prediction in milling inserts and correlation with experimentally observed damage mechanisms," Journal of Materials Processing Technology, vol. 256, pp. 98-108, 2018. [42] A. Thakare and A. Nordgren, "Experimental Study and Modeling of Steady State Temperature Distributions in Coated Cemented Carbide Tools in Turning," Procedia CIRP, vol. 31, pp. 234-239, 2015. [43] H. B. Wu and S. J. Zhang, "3D FEM simulation of milling process for titanium alloy Ti6Al4V," The International Journal of Advanced Manufacturing Technology, vol. 71, pp. 1319-1326, 2014. [44] A. U. Anagonye and D. a. Stephenson, "Modeling Cutting Temperatures for Turning Inserts With Various Tool Geometries and Materials," Journal of Manufacturing Science and Engineering, vol. 124, pp. 544-544, 2002. 172 [45] V. Dessoly, S. N. Melkote, and C. Lescalier, "Modeling and verification of cutting tool temperatures in rotary tool turning of hardened steel," International Journal of Machine Tools and Manufacture, vol. 44, pp. 1463-1470, 2004. [46] M. C. Shaw and J. Cookson, Metal cutting principles: Clarendon press Oxford, 1984. [47] A. Mamedov and I. Lazoglu, "Thermal analysis of micro milling titanium alloy Ti-6Al-4V," Journal of Materials Processing Technology, vol. 229, pp. 659-667, 2016. [48] A. C. Rapier, "A theoretical investigation of the temperature distribution in the metal cutting process," British Journal of Applied Physics, vol. 5, pp. 400-405, 1954. [49] R. P. Dutt and R. C. Brewer, "On The Theoretical Determination of The Temperature Field in Orthogonal Machining," International Journal of Production Research, vol. 4, pp. 91 - 114, 1965. [50] E. K. Levy, C. L. Tsai, and M. P. Groover, "Analytical Investigation of the Effect of Tool Wear on the Temperature Variations in a Metal Cutting Tool," Journal of Engineering for Industry, vol. 98, pp. 251-257, 1976. [51] J. F. W. Bishop, "An Approximate Method for Determining The Temperatures Reached in Steady Motion Problems of Plane Plastic Strain," The Quarterly Journal of Mechanics and Applied Mathematics, vol. 9, pp. 236-246, January 1, 1956 1956. [52] E. Usui, T. Shirakashi, and T. Kitagawa, "Analytical prediction of three dimensional cutting process, part 3: cutting temperature and crater wear of carbide tool," Journal of Engineering for Industry, vol. 100, pp. 236-243, 1978. [53] A. Smith and E. Armarego, "Temperature Prediction in Orthogonal Cutting with a Finite Difference Approach," CIRP Annals - Manufacturing Technology, vol. 30, pp. 9-13, 1981. [54] I. Lazoglu and Y. Altintas, "Prediction of tool and chip temperature in continuous and interrupted machining," International Journal of Machine Tools and Manufacture, vol. 42, pp. 1011-1022, 2002. [55] P. L. B. Oxley, The mechanics of machining : an analytical approach to assessing machinability: E. Horwood ; Halsted Press, 1989. [56] D. Ulutan, I. Lazoglu, and C. Dinc, "Three-dimensional temperature predictions in machining processes using finite difference method," Journal of Materials Processing Technology, vol. 209, pp. 1111-1121, 2009. [57] W. Grzesik, M. Bartoszuk, and P. Nieslony, "Finite difference analysis of the thermal behaviour of coated tools in orthogonal cutting of steels," International Journal of Machine Tools and Manufacture, vol. 44, pp. 1451-1462, 2004. [58] I. Lazoglu and C. Islam, "Modeling of 3D temperature fields for oblique machining," CIRP Annals - Manufacturing Technology, vol. 61, pp. 127-130, 2012. [59] Y. Altintas, Manufacturing automation: metal cutting mechanics, machine tool vibrations, and CNC design: Cambridge university press, 2012. [60] J. Tannehill, D. Anderson, and R. Pletcher, Computational Fluid Mechanics and Heat Transfer: Taylor & Francis, 1997. [61] R. J. LeVeque, Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems vol. 98: Siam, 2007. [62] V. D. Liseikin, Grid Generation Methods, 2 ed.: Springer, 2010. [63] J. E. Thompson, Z. U. A. Warsi, and C. W. Mastin, Numerical Grid Generation Foundations and Applications: Elsevier Science Publishing 1985. 173 [64] K. A. Hoffmann and S. T. Chiang, Computational Fluid Dynamics vol. 1: Engineering Education System, 2000. [65] H. Blok, "Theoretical study of temperature rise at surfaces of actual contact under oiliness lubricating conditions," Proceedings of the General Discussion on Lubrication and Lubricants, vol. 2, pp. 222–235, 1937. [66] T. Kitagawa, A. Kubo, and K. Maekawa, "Temperature and wear of cutting tools in high-speed machining of Inconel 718 and Ti 6Al 6V 2Sn," Wear, pp. 142-148, 1997. [67] D. A. Stephenson, "Assessment of Steady-State Metal Cutting Temperature Models Based on Simultaneous Infrared and Thermocouple Data," Journal of Manufacturing Science and Engineering, vol. 113, pp. 121-128, 1991. [68] K. C. Mills, Recommended values of thermophysical properties for selected commercial alloys: Woodhead Publishing, 2002. [69] E. Ezugwu, "Key improvements in the machining of difficult-to-cut aerospace superalloys," International Journal of Machine Tools and Manufacture, vol. 45, pp. 1353-1367, 2005. [70] O. Engineering. (2018, 01-06). TJC36 Information Page [Website]. Available: https://www.omega.ca/pptst_eng/TJC36.html#description [71] M. Nasr, E. G. Ng, and M. A. Elbestawi, "A modified time-efficient FE approach for predicting machining-induced residual stresses," Finite Elements in Analysis and Design, vol. 44, pp. 149-161, 2008. [72] R. M'Saoubi and H. Chandrasekaran, "Investigation of the effects of tool micro-geometry and coating on tool temperature during orthogonal turning of quenched and tempered steel," International Journal of Machine Tools and Manufacture, vol. 44, pp. 213-224, 2004. [73] X. Jin and Y. Altintas, "Slip-line field model of micro-cutting process with round tool edge effect," Journal of Materials Processing Technology, vol. 211, pp. 339-355, 2011/03/01/ 2011. [74] F. Akbar, P. Mativenga, and M. Sheikh, "An evaluation of heat partition in the high-speed turning of AISI/SAE 4140 steel with uncoated and TiN-coated tools," Proceedings of the institution of mechanical engineers, Part B: Journal of Engineering Manufacture, vol. 222, pp. 759-771, 2008. [75] T. Obikawa, T. Matsumura, T. Shirakashi, and E. Usui, "Wear characteristic of alumina coated and alumina ceramic tools," Journal of materials processing technology, vol. 63, pp. 211-216, 1997. [76] O. Lesquois, J. Serra, P. Kapsa, S. Serror, and C. Boher, "Degradations in a high-speed sliding contact in transient regime," Wear, vol. 201, pp. 163-170, 1996. [77] A. Özel, V. Ucar, A. Mimaroglu, and I. Calli, "Comparison of the thermal stresses developed in diamond and advanced ceramic coating systems under thermal loading," Materials & Design, vol. 21, pp. 437-440, 2000. [78] J. Outeiro, A. Dias, and I. Jawahir, "On the effects of residual stresses induced by coated and uncoated cutting tools with finite edge radii in turning operations," CIRP Annals-Manufacturing Technology, vol. 55, pp. 111-116, 2006. [79] I. S. Jawahir and C. A. van Luttervelt, "Recent Developments in Chip Control Research and Applications," CIRP Annals - Manufacturing Technology, vol. 42, pp. 659-693, 1993. 174 [80] N. Abukhshim, P. Mativenga, and M. Sheikh, "Investigation of heat partition in high speed turning of high strength alloy steel," International Journal of Machine Tools and Manufacture, vol. 45, pp. 1687-1695, 2005. [81] J. Rech, P. J. Arrazola, C. Claudin, C. Courbon, F. Pusavec, and J. Kopac, "Characterisation of friction and heat partition coefficients at the tool-work material interface in cutting," CIRP Annals - Manufacturing Technology, vol. 62, pp. 79-82, 2013. [82] M. E. Merchant, "Mechanics of the metal cutting process. I. Orthogonal cutting and a type 2 chip," Journal of applied physics, vol. 16, pp. 267-275, 1945. [83] A. Moufki, A. Devillez, D. Dudzinski, and A. Molinari, "Thermomechanical modelling of oblique cutting and experimental validation," International Journal of Machine Tools and Manufacture, vol. 44, pp. 971-989, 2004. [84] E. Shamoto and Y. Altintas, "Prediction of Shear Angle in Oblique Cutting with Maximum Shear Stress and Minimum Energy Principles," Journal of Manufacturing Science and Engineering, vol. 121, pp. 399-407, 1999. 175 Appendix A : Oblique Cutting Mechanics The geometry of oblique cutting is demonstrated in Figure 3-1. Cutting velocity vector and tool cutting edge normal is oriented to each other with an inclination angle ( s ) in oblique cutting (see Figure 3-1). The orientation of cutting edge leads to the generation of cutting forces in three different directions (i.e., radial, tangential, feed forces). Figure A- 1- Section view of tool and workpiece on the cutting edge normal plane Figure A- 2- Top view of the shear plane from Z4 direction 176 The relationship between the cutting forces defined in Figure 3-2. The force balance on the primary and secondary zones can be analyzed via vector algebra. The shearing force ( sF ) and shear normal force ( nF ) exist on the shear plane due to plastic material deformation ( the primary deformation zone). The friction force on rake face ( uF ) and friction normal force ( vF ) are generated because of the contact between tool and chip. The resultant cutting force and the force balance on in the primary and secondary deformation zones are as follows; R s nR u vF F FF F F (A.1) The transformation frames are defined to make the derivations of cutting speed, force, and power relationships. Figure 3-1, Figure 3-2, Figure A- 1, and Figure A- 2 show the local coordinates used in the derivation procedure. 0 0 0, ,X Y Z is the main coordinate frame where 0X is in the cutting velocity direction, and 0Z is in the opposite of the feed direction. 1 1 1, ,X Y Z is defined on the cutting edge in which 1Y is along the cutting edge, and 1Z is on the cutting edge normal plane ( nP ). On the rake face, a local coordinate frame 2 2 2, ,X Y Z is defined by rotating the 1 1 1, ,X Y Z around 1Y counter-clockwise with the normal rake angle ( n ). A chip flow coordinate system 3 3 3, ,X Y Z is assigned in which 3Z is colinear with the chip flow direction. 4 4 4, ,X Y Z is defined on the shear plane via rotating 1 1 1, ,X Y Z . A shearing direction coordinate frame 5 5 5, ,X Y Z is created where 5X is aligned with the shear velocity direction. 177 6 6 6, ,X Y Z is the resultant cutting force coordinate frame where 6X is aligned with the resultant cutting force. The following rotation matrixes are defined between the coordinate frames; 01cos( ) sin( ) 0R sin( ) cos( ) 00 0 1s ss s (A.2) 12cos( ) 0 sin( )R 0 1 0sin( ) 0 cos( )n nn n (A.3) 23 c cc c1 0 0R 0 cos(η ) sin(η )0 sin(η ) cos(η ) (A.4) n n14n ncos( ) 0 sin( )R 0 1 0sin( ) 0 cos( ) (A.5) i i45 i icos( ) sin( ) 0R -sin( ) cos( ) 00 0 1 (A.6) n n i i16 i in ncos(θ ) 0 sin(θ ) cos(θ ) sin(θ ) 0R 0 1 0 -sin(θ ) cos(θ ) 0sin(θ ) 0 cos(θ ) 0 0 1 (A.7) where s , n , c , n , i are inclination angle, normal rake angle, chip flow angle, normal shear angle, and oblique shear angle, respectively. n and i are related to frictional contact and chip flow direction that is defined later in this section. 178 The primary deformation zone shear force ( sF ) is generated on the shear plane. It is calculated by multiplying the average shear flow stress ( s ) and shear plane area ( sA ). s s sF A (A.8) The shear plane area is approximated via conservation of material throughput in the primary zone volumetrically [83]. Material volume at the entry and exit of the primary shear zone is equated as follows; s nb h V A V (A.9) where b , h are the width of cut and depth of cut, respectively. The shear normal velocity component nV is written as (see Figure 3-2 , Figure A- 2) ; cos sin( )n s nV V (A.10) Combining (A.9) and (A.10) lead to, nsin cosssb hA (A.11) The shear force and the resultant cutting force are related by using vector transformations. The equation (A.12) shows the force component along the shear direction that is defined with the oblique shear angle i and the normal shear angle n . ncos cos cos sin sin ,cos cos cos sin sin sin coss R i i n n i isRi i n n i i sF F leads tob hF (A.12) 179 i is the angle between the resultant cutting force and cutting edge normal plane ( nP ). n is the angle between the newly generated surface and the resultant cutting force projected component on cutting edge normal plane ( nP ). The friction force and friction normal force on the rake face are written in terms of average friction angle ( a ) [84] (see Figure 3-2). 11sinsincos( )tantancossin sin sintan cos taniu R a Rcv R an nu v a vci a cn n c aleads toF F sin FF FF F F (A.13) The shear stress, friction angle and shear angle values are identified for each workpiece and tool material combination by conducting orthogonal cutting experiments. The further information about the mechanistic identification is found in [59]. The interrelation between cutting velocity, chip velocity, and shear velocity is determined from the following vector summation; s cV V V (A.14) where V , sV and cV are cutting, shearing and chip velocities, respectively. Vector summation in the equation (A.14) is written in shear plane coordinates via transforming chip and cutting velocity components to shear plane. 180 nT T14 01nV V cos cosR R 0 V sin0 V cos sinsss (A.15) c c n nT14 12 23 c cc c n nV cos η cos sin cos sin( )0R R 0 sin η VV V cos η cos cos sin sinn nc n nR (A.16) It is assumed that the velocity component perpendicular to the shear plane at the entrance and exit of the shear zone should be equal to each other [83]. The chip velocity magnitude is identified via this assumption as; ncc nV cos sin Vcos η cossn (A.17) nn n ccos cos sec( ) sin cos sec sin tan(η )0s n ns s s nVV V (A.18) The oblique shear angle on shear plane is determined by using equation (A.18) as the following; c n ntan sec tan η sin tan cos( )i n s n (A.19) The tool-chip contact area identification is vital to calculate the heat flux application area. It is evaluated by taking moments of normal shear force and friction normal force with respect to cutting edge by assuming uniform pressure distribution on the shear plane and tool-chip interface. cos cos cos cos sin sinv R a R i n nn nF F F (A.20) cos ( )n R i n nF F sin (A.21) 181 cos2 2sec ( ) cos cos sin ( )ncv c nc n ncn n n nnhsinlF Fh sinlsin sin (A.22) where ℎ, cl are uncut chip thickness and tool-chip contact length, respectively.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A transient thermal model for machining
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A transient thermal model for machining Islam, Coskun 2018
pdf
Page Metadata
Item Metadata
Title | A transient thermal model for machining |
Creator |
Islam, Coskun |
Publisher | University of British Columbia |
Date Issued | 2018 |
Description | Prediction of temperature in the tool, chip and workpiece surface layer is essential for tool design and the selection of most productive cutting conditions which yield the desired tool life and acceptable residual stresses left on the machined part. This thesis presents a comprehensive, finite difference method based numerical model on simulating the temperature distribution in the chip, tool, and finished workpiece surface layer as a function of material properties, cutting speed, feed rate and tool-workpiece engagement period. The heat is generated in the primary shear zone where the chip is sheared from the metal, in the secondary zone where the chip sticks and slides on the rake face, and in the tertiary zone where cutting-edge ploughs the workpiece surface. The chip, layer of the workpiece surface and the tool edge are meshed into discrete elements. The heat is transferred to the stationary tool, and dynamically moving chip and workpiece surface by conforming heat balance equations within each element. A finite difference technique with implicit time discretization is used to solve heat balance equations of the temperature fields on the tool, workpiece, and chip. Anisotropic material properties can be considered in the model which allows the inclusion of a coating layer on the tool. The proposed model allows two and three-dimensional heat transfer, hence it can be used to predict the temperature distribution in turning, drilling and milling operations. The continuous machining processes such as turning generate constant heat, so the temperature reaches a steady state after a transient period. The intermittent operations such as milling generate time-varying and periodic heat, hence the temperature variation is always in a transient state. The proposed model is experimentally validated with the data found in the literature and experiments conducted by the author at the industrial partner’s (Sandvik Coromant AB, Sweden) research facility. Experimental validations cover uncoated, single and multi-layer coated tools to simulate continuous turning, interrupted turning and milling operations. The proposed model is able to predict the temperature with less than 20% error in most of the validated cases. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2018-10-22 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0372960 |
URI | http://hdl.handle.net/2429/67650 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2018-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2018_november_islam_coskun.pdf [ 10.96MB ]
- Metadata
- JSON: 24-1.0372960.json
- JSON-LD: 24-1.0372960-ld.json
- RDF/XML (Pretty): 24-1.0372960-rdf.xml
- RDF/JSON: 24-1.0372960-rdf.json
- Turtle: 24-1.0372960-turtle.txt
- N-Triples: 24-1.0372960-rdf-ntriples.txt
- Original Record: 24-1.0372960-source.json
- Full Text
- 24-1.0372960-fulltext.txt
- Citation
- 24-1.0372960.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0372960/manifest