UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Modeling of machining thin-walled aerospace structures Tuysuz, Oguzhan 2018

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2019_february_tuysuz_oguzhan.pdf [ 12.41MB ]
Metadata
JSON: 24-1.0373471.json
JSON-LD: 24-1.0373471-ld.json
RDF/XML (Pretty): 24-1.0373471-rdf.xml
RDF/JSON: 24-1.0373471-rdf.json
Turtle: 24-1.0373471-turtle.txt
N-Triples: 24-1.0373471-rdf-ntriples.txt
Original Record: 24-1.0373471-source.json
Full Text
24-1.0373471-fulltext.txt
Citation
24-1.0373471.ris

Full Text

MODELING OF MACHINING THIN-WALLED AEROSPACE STRUCTURES by  Oguzhan Tuysuz  B.Sc., Mechanical Engineering, Istanbul Technical University, 2009 M.A.Sc., Mechanical Engineering, The University of British Columbia, 2012  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)  November 2018  © Oguzhan Tuysuz, 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: Modeling of Machining Thin-Walled Aerospace Structures  submitted by Oguzhan Tuysuz  in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Mechanical Engineering  Examining Committee: Dr. Yusuf Altintas Supervisor  Dr. Anasavarapu Srikantha Phani Supervisory Committee Member  Dr. Terje Haukaas University Examiner Dr. Gary S. Schajer University Examiner   iii  Abstract Thin-walled monolithic aerospace parts and rotor blades have high flexibilities causing severe static and dynamic deflections during machining. Poor dimensional accuracy and surface finish due to deflections scrap the costly parts. This thesis presents mathematical models to simulate thin-walled part machining in virtual environment considering the varying structural workpiece characteristics and process induced damping along complex toolpaths.  Traditionally having full order finite element (FE) models at several toolpath locations is prohibitive in the machining of curved blades. New, computationally efficient, reduced and full order workpiece dynamics update models are developed. Removed materials between discrete locations are defined as substructures of the initial workpiece. First, in-process workpiece frequency response functions (FRFs) are directly updated by coupling fictitiously negative dynamic stiffness of the removed materials. The model is improved by introducing substructure decoupling in time domain. The workpiece structure is modified by coupling fictitious substructures having the negative mass and stiffness of removed volumes. Mode shapes of the in-process workpiece are perturbed, and mode frequencies and workpiece FRFs are updated. The computed FRFs of the thin-walled parts are used to predict the chatter stability, static deflections, forced vibrations, and their effects on tolerance violations along the toolpath.    Unlike the conventional empirical process damping coefficients, a comprehensive analytical model to predict the machining process damping is proposed. The cutting edge is discretized in the chip width direction, and contact pressure between the edge element and workpiece surface is estimated using the tool geometry, vibration parameters, and work material properties. The specific process damping force of each element is evaluated by integrating the contact pressure. The damping force is linearized by representing it with equivalent viscous damper. iv  A generalized five-axis ball-end milling dynamics model is developed in frequency domain by incorporating the dynamics update and process damping models for flexible parts. Relative tool-workpiece vibrations are projected into the local chip thickness direction and the dynamic milling equation is derived. Milling stability is assessed at discrete locations using Nyquist criterion, and chatter regions and frequencies are predicted along cutting.  The proposed digital process models are experimentally verified and expected to guide engineers in process development.   v  Lay Summary Thin structural and functional components used in the aerospace industry are manufactured with milling operations by removing excess material from the initial geometry with a cutting tool. These highly flexible and costly parts severely vibrate during the operation. Unless avoided, these vibrations may become unstable, which in turn cause poor surface quality and scrapped parts. As opposed to the current trial and error based practice in the industry, this thesis develops the physics based mathematical methods to predict the vibrational characteristics of the part, and their interactions with the cutting tool and work material as cutting advances in the simulation environment. Although jet engine blades are used as the test bed, the introduced models are general, and they can identify the process related issues and guide the process planners to solve them before the actual part is made.   vi  Preface This Ph.D. study has been conducted at the Manufacturing Automation Laboratory under the supervision of Dr. Yusuf Altintas. The problem originated from the industrial collaborator, Pratt & Whitney Canada, and the candidate designed and his supervisor approved all the research steps including a six months long industrial internship with the collaborator. The candidate conducted all of the research and analyzed the collected data. The developed models and outcomes of the research have either been published, or submitted for review in journals with high standing in the field of manufacturing. Contributions of the candidate to each chapter are as follows: A concise version of Chapter 3 has been published [1] as, “Tuysuz, O., and Altintas, Y., 2017, Frequency domain updating of thin-walled workpiece dynamics using reduced order substructuring method in machining, Journal of Manufacturing Science and Engineering, 139(7), pp. 071013”. The manuscript was written by the first author and edited by the research supervisor. The first author developed the model, planned and conducted all the validation tests, and analyzed the collected data.    A concise version of Chapter 4 has been published [2] as, “Tuysuz, O., and Altintas, Y., 2017, Time domain modeling of varying dynamic characteristics in thin-wall machining using perturbation and reduced order substructuring, Journal of Manufacturing Science and Engineering, 140(1), pp. 011015”. The manuscript was written by the first author and edited by the research supervisor. The first author developed the model, planned and conducted all the validation tests, and analyzed the collected data.      Chapters 4 & 6 have also partially contributed to the study in [3] as, “Altintas, Y., Tuysuz, O., Habibi, M., and Li, Z. L., 2018, Virtual compensation of deflection errors in ball-end milling vii  of flexible blades, CIRP Annals, 67(1), pp. 365-368”. The candidate was responsible for updating the static stiffness, cutting mechanics and deflection calculations, and conducting the validation cutting tests. The candidate assisted the research supervisor by providing the procedures used in the computational models and experimental results. Dr. Habibi –a former post-doctoral fellow at Manufacturing Automation Laboratory– contributed to geometric calculations, error compensation solution, and conducted the machined surface measurements. Mr. Li –a former visiting Ph.D. student at Manufacturing Automation Laboratory– contributed to implementation of the deflection compensation algorithm.    A concise version of Chapters 5 & 6 has been submitted as a journal article as, “Tuysuz, O., and Altintas, Y., 2018, Analytical modeling of machining process damping”. The manuscript is written by the first author and edited by the research supervisor. The first author developed the model, used some existing experimental data at the Manufacturing Automation Laboratory, planned and conducted additional validation tests, and analyzed the collected data. viii  Table of Contents  Abstract ......................................................................................................................................... iii Lay Summary .................................................................................................................................v Preface ........................................................................................................................................... vi Table of Contents ....................................................................................................................... viii List of Tables ............................................................................................................................... xii List of Figures ............................................................................................................................. xiv List of Symbols ......................................................................................................................... xxiii List of Abbreviations ............................................................................................................. xxxiv Acknowledgments .................................................................................................................. xxxvi Dedication .............................................................................................................................. xxxvii Chapter 1: Introduction ................................................................................................................1 Chapter 2: Literature Review .......................................................................................................6 2.1 Overview ......................................................................................................................... 6 2.2 Varying Thin-Walled Workpiece Characteristics ........................................................... 6 2.3 Process Damping (PD) Phenomenon in Machining ..................................................... 11 2.4 Dynamics and Chatter Stability of Ball-End Milling Operations ................................. 14 2.5 Summary ....................................................................................................................... 16 Chapter 3: Modeling of Varying Workpiece Dynamics in Frequency Domain .....................18 3.1 Overview ....................................................................................................................... 18 3.2 Modeling of Position Dependent Dynamics of Thin-Walled Workpiece ..................... 19 ix  3.2.1 Updated Direct FRFs at the Interface DOF - iBccΗ  ................................................... 28 3.2.2 Updated Cross FRFs between the Interface and Internal DOF - iBcbH  ...................... 29 3.2.3 Updated Direct FRFs at the Internal DOF - iBbbH ...................................................... 30 3.3 Effect of the Truncated Dynamics ................................................................................ 31 3.4 Insertion of Damping .................................................................................................... 34 3.5 Reduced Order Dynamics Model of the Thin-Walled Workpiece ............................... 38 3.5.1 Selection of Master Degrees of Freedom.................................................................. 38 3.5.2 Model Order Reduction ............................................................................................ 40 3.6 Simulation Algorithm ................................................................................................... 43 3.7 Experimental Verification: Thin-Walled Rectangular Workpiece Cutting Test .......... 46 3.8 Summary ....................................................................................................................... 59 Chapter 4: Modeling of Varying Workpiece Dynamics in Time Domain ..............................60 4.1 Overview ....................................................................................................................... 60 4.2 Modeling of Position Dependent Mass and Stiffness of the Thin-Walled Workpiece . 61 4.3 Modeling of Varying Dynamic Characteristics of the In-Process Workpiece ............. 63 4.3.1 Evaluation of the Updated Modes and Mode Shapes ............................................... 64 4.3.2 Simulation Algorithm ............................................................................................... 72 4.4 Experimental Verification: Thin-Walled Twisted Blade Cutting Test ......................... 73 4.5 Summary ....................................................................................................................... 88 Chapter 5: Analytical Modeling of Process Damping in Machining ......................................90 5.1 Overview ....................................................................................................................... 90 5.2 Process Damping Phenomenon..................................................................................... 91 x  5.3 Indentation Geometry ................................................................................................... 93 5.4 Modeling of Contact Pressure under the Flank Face .................................................... 96 5.5 Evaluation of the Process Damping Force and Equivalent Viscous Damping ........... 103 5.6 Experimental Verification ........................................................................................... 106 5.6.1 Verification of Equivalent Viscous Damping (ceq) ................................................. 106 5.6.2 Parameter Uncertainty Analysis ............................................................................. 114 5.7 Summary ..................................................................................................................... 118 Chapter 6: Modeling of Chatter Stability in Frequency Domain..........................................120 6.1 Overview ..................................................................................................................... 120 6.2 Dynamics of Five-Axis Ball-End Milling................................................................... 121 6.2.1 Ball-End Mill Geometry ......................................................................................... 122 6.2.2 Static and Dynamic Chip Thickness ....................................................................... 125 6.2.3 Dynamic Cutting and Process Damping Forces ..................................................... 130 6.2.4 Relative Dynamic Displacements between the Tool and Thin-Walled Workpiece 136 6.2.5 Dynamic Milling Equation and Chatter Stability ................................................... 137 6.3 Experimental Verification ........................................................................................... 142 6.3.1 Stability Diagrams with Process Damping in Orthogonal Turning ........................ 143 6.3.2 Stability Diagrams with Process Damping in Flat-End Milling ............................. 146 6.3.3 Thin-Walled Rectangular Workpiece Chatter Test ................................................. 153 6.3.4 Impeller Blade Cutting Test .................................................................................... 156 6.3.5 Application: Prediction of Static Form Errors in Ball-End Milling of Blades ....... 165 6.3.6 Application: Prediction of Thin-Walled Rectangular Workpiece Vibrations......... 170 6.4 Summary ..................................................................................................................... 175 xi  Chapter 7: Conclusions and Future Research Directions ......................................................177 7.1 Summary and Contributions ....................................................................................... 177 7.2 Future Research Directions ......................................................................................... 180 References ...................................................................................................................................182 Appendices ..................................................................................................................................190 Appendix A Epsilon Algorithm .............................................................................................. 190 Appendix B Directional Coefficients for Dynamic Cutting and Process Damping Forces .... 192  xii  List of Tables Table 3.1: Cutting conditions in thin-walled rectangular workpiece verification experiments .... 47 Table 4.1: Cutting conditions in five-axis twisted blade experiments .......................................... 76 Table 4.2: Average prediction error and computation time comparisons of the developed thin-walled workpiece dynamics update models for the twisted blade case (EVP: Eigenvalue problem, TD: Time domain, FD: Frequency domain) ................................................................................. 86 Table 5.1: Cutting edge geometry, vibration and friction parameters, and the average separation angle for the orthogonal cutting tests in [48] .............................................................................. 107 Table 5.2: Cutting edge geometry, vibration and friction parameters, and the average separation angle used in the FE simulations in [56] ..................................................................................... 109 Table 5.3: Cutting edge geometry, vibration and friction parameters, and the average separation angle for the orthogonal cutting tests in [54] .............................................................................. 111 Table 6.1: Cutting edge geometry, vibration and friction parameters for the flat-end milling tests in [54] .......................................................................................................................................... 147 Table 6.2: Modal parameters of the flat-end cutter measured at the tool tip in feed and normal directions in [54], and the estimated simulation parameters ....................................................... 147 Table 6.3: Cutting edge geometry, vibration and friction parameters for the flat-end milling tests in [118] ........................................................................................................................................ 150 Table 6.4: Modal parameters of the flat-end cutter measured at the tool tip in feed and normal directions in [118], and the estimated simulation parameters..................................................... 150 Table 6.5: Cutting conditions in rectangular thin-walled workpiece chatter tests ...................... 153 Table 6.6: Cutting conditions in semi-finishing the blade top .................................................... 158 xiii  Table 6.7: Cutting edge geometry, vibration, friction, and simulation parameters for process damping calculations in semi-finishing the blade top ................................................................. 162 Table 6.8: Cutting conditions used in the stable cutting tests in [6] ........................................... 172  xiv  List of Figures Figure 1.1: Five-axis flank and point milling of blades .................................................................. 1 Figure 1.2 Outline of the thesis project ........................................................................................... 3 Figure 2.1: Effect of material removal on the in-process workpiece flexibility ............................. 7 Figure 2.2: Improved stable axial depth of cut limits in the presence of process damping .......... 11 Figure 2.3: Surface quality with and without chatter vibrations ................................................... 14 Figure 3.1: Time invariant tool dynamics and varying workpiece dynamics at different stages of machining ...................................................................................................................................... 20 Figure 3.2: (a) Geometrically removed material, (b) dynamically removed material equivalent to coupling of a fictitious substructure .............................................................................................. 20 Figure 3.3: (a) In-process workpiece i-1B  at ( )thi -1  machining step, (b) component level interaction between the substructure iA  and the rest of in-process workpiece i-1B  represented by rigid links, (c) component level interaction represented by interface coupling (reaction) forces, (d) opposite interface coupling forces generated by the fictitiously added iA , and interface force cancellation after addition of iA  ................................................................................................. 22 Figure 3.4: Three-step uniform material removal in flank milling of a rectangular thin-walled workpiece ...................................................................................................................................... 33 Figure 3.5: Comparison of the undamped FRFs with and without the truncated workpiece dynamics compensation in material removal steps ....................................................................... 34 Figure 3.6: Damping dominant regions on the real part of a lightly damped thin-walled workpiece FRF .............................................................................................................................. 36 xv  Figure 3.7: Nodal compatibility between the in-process workpiece and fictitious material, a) before, and b) after model order reduction ................................................................................... 38 Figure 3.8: Flow chart of the frequency domain in-process workpiece dynamics update model 44 Figure 3.9: Initial thin-walled rectangular workpiece (a) CAD geometry, (b) FE discretization . 47 Figure 3.10: (a) Experimental setup, (b) measurement locations and removed segments in roughing operation ........................................................................................................................ 48 Figure 3.11: Comparison of the measured and predicted direct FRFs along xW  at points #1, #2, #4, and #5 in roughing .................................................................................................................. 51 Figure 3.12: Comparison of the measured and predicted direct FRFs along xW  at points #7, #8, #10, and #11 in semi-finishing...................................................................................................... 52 Figure 3.13: Comparison of the measured and predicted direct FRFs along xW  at points #13, #14, #16, and #17 in finishing ...................................................................................................... 53 Figure 3.14: a) Variation of the first bending mode frequency at all FRF comparison points in the experiment and simulation, b) comparison of the prediction error in full and reduced order dynamics decoupling for the first bending mode .......................................................................... 54 Figure 3.15: a) Variation of the first torsional mode frequency at all FRF comparison points in the experiment and simulation, b) comparison of the prediction error in full and reduced order dynamics decoupling for the first torsional mode ......................................................................... 55 Figure 3.16: Schematic illustration of the master DOF selected due to, a) the first bending (shaded region), b) the first torsional (shaded region) modes of the initial workpiece ................ 56 Figure 3.17: Comparison of the experimentally measured and predicted static stiffness of the rectangular workpiece in xW  at all FRF measurement points ...................................................... 58 xvi  Figure 4.1: Varying mode frequencies and dynamic displacements in the cutting region ........... 65 Figure 4.2: Flow chart of the proposed in-process workpiece modal parameters update model.. 73 Figure 4.3: a) Initial workpiece geometry, machining sequence, and four removed segments in each operation, b) FE mesh of the initial thin-walled twisted blade ............................................. 75 Figure 4.4: a) Experimental setup, b) measurement (comparison) locations (points 1-9) in R1 .. 76 Figure 4.5: Comparison of the measured and predicted direct FRFs in xW  direction using the full order time domain (TD) model and repetitively solved eigenvalue problem (EVP) at points #1, #2, #6, and #7 in operation R1 ...................................................................................................... 77 Figure 4.6: Comparison of the measured and predicted direct FRFs in xW  direction using the full order time domain (TD) model and repetitively solved eigenvalue problem (EVP) at points #19, #20, #24, and #25 in operation SF1 .............................................................................................. 78 Figure 4.7: Comparison of the measured and predicted direct FRFs in xW  direction using the full order time domain (TD) model and repetitively solved eigenvalue problem (EVP) at points #37, #38, #42, and #43 in operation F1 ................................................................................................ 79 Figure 4.8: Comparison of the measured and predicted direct FRFs in xW  direction using the reduced order time domain (TD) and the reduced order frequency domain (FD) models at points #1, #2, #6, and #7 in operation R1 ................................................................................................ 81 Figure 4.9: Comparison of the measured and predicted direct FRFs in xW  direction using the reduced order time domain (TD) and the reduced order frequency domain (FD) models at points #19, #20, #24, and #25 in operation SF1 ...................................................................................... 82 xvii  Figure 4.10: Comparison of the measured and predicted direct FRFs in xW  direction using the reduced order time domain (TD) and the reduced order frequency domain (FD) models at points #37, #38, #42, and #43 in operation F1 ........................................................................................ 83 Figure 4.11: a) Variation of the predicted and measured first bending mode frequency with material removal, b) comparison of the absolute prediction error in the full and reduced order time domain (TD) and frequency domain (FD) models, and the full order repetitive eigenvalue (EVP) solution .............................................................................................................................. 85 Figure 4.12: a) Variation of the predicted and measured first torsional mode frequency with material removal, b) comparison of the absolute prediction error in the full and reduced order time domain (TD) and frequency domain (FD) models, and the full order repetitive eigenvalue (EVP) solution .............................................................................................................................. 87 Figure 4.13: Effect of number of the perturbed modes on the performance and prediction error in the reduced order time domain (TD) model ................................................................................. 87 Figure 5.1: (a) Two-dimensional cutting edge geometry, (b) static chip generation and material separation point (SP) ..................................................................................................................... 91 Figure 5.2: Machining process damping due to the cutting edge indentation and the associated forces ............................................................................................................................................. 92 Figure 5.3: Indentation geometry, surface wave form, and discrete process damping calculation points ............................................................................................................................................. 94 Figure 5.4: General contact pressure distribution at the interface between the cutting edge element and workpiece ................................................................................................................. 97 Figure 5.5: Flow chart for the analytical estimation of process damping ................................... 106 xviii  Figure 5.6: Comparison of the experimentally identified equivalent process damping against prediction of the proposed analytical model for AISI 1045 in (a) radial, and (b) tangential directions ..................................................................................................................................... 108 Figure 5.7: Comparison of the numerically (FE) identified equivalent process damping against prediction of the proposed analytical model for AISI 1045 in (a) radial, and (b) tangential directions ..................................................................................................................................... 109 Figure 5.8: a-b) Measured acceleration at the tool tip, and its frequency content using Fast Fourier Transform (FFT), c-d) corresponding displacement and its FFT for AISI 1018 at 2000 rpm .............................................................................................................................................. 111 Figure 5.9: a-b) Measured acceleration at the tool tip, and its frequency content using Fast Fourier Transform (FFT), c-d) corresponding displacement and its FFT for Ti6Al4V at 850 rpm..................................................................................................................................................... 112 Figure 5.10: Comparison of the experimentally identified equivalent process damping against prediction of the proposed analytical model for (a) AISI 1018, and (b) Ti6Al4V work materials..................................................................................................................................................... 113 Figure 5.11: Uncertainty analysis for the parameters: edge radius, material separation angle, and the plastic limit for the test case in [48] ...................................................................................... 115 Figure 5.12: Uncertainty analysis for the parameters: edge radius, material separation angle, and the plastic limit for the test case in [56] ...................................................................................... 116 Figure 5.13: Uncertainty analysis for the parameters: vibration amplitude, edge radius, material separation angle, and the plastic limit for the first test case in [54]............................................ 117 Figure 5.14: Uncertainty analysis for the parameters: vibration amplitude, edge radius, material separation angle, and the plastic limit for the second test case in [54] ....................................... 118 xix  Figure 6.1: A sample tool-workpiece engagement (TWE) in five-axis ball-end milling, its discretization, and variation of st  and ex  along z-axis ............................................................ 122 Figure 6.2: Ball-end mill envelope and tool coordinate system (TCS), a) front view, b) bottom view [112] ................................................................................................................................... 123 Figure 6.3: a) Representative variation in static chip thickness of the first flute with angular and axial positions, b) three dimensional dynamic chip formation in ball-end milling .................... 127 Figure 6.4: Chip regeneration and dynamic chip thickness on the cross section of the jth flute at elevation z ................................................................................................................................... 128 Figure 6.5: Flow chart to assess the chatter stability in five-axis ball-end milling of thin-walled parts ............................................................................................................................................. 142 Figure 6.6: Comparison of the stability lobes predicted using equivalent specific process damping (ceq) and the experimentally identified process damping force coefficient (Ci) for the parameters in Table 5.1. Work material: AISI 1045;Structural parameters: ɷn =450.7 Hz, mr =0.8081 kg,  kr =6.48x106 N/m,  cr =145 Ns/m; Cutting force coefficient: Kr =1384 MPa ........ 144 Figure 6.7: Comparison of the stability lobes predicted using equivalent specific process damping (ceq) and the experimentally identified process damping coefficient (cp,r) for the parameters in Table 5.3. Work material: AISI 1018; Structural parameters: ɷn =1403 Hz, kr =21x106 N/m, ξr=2.2%; Cutting force coefficient: Kr =1375 MPa ............................................. 145 Figure 6.8: a) Radial displacements (vibrations) of the three fluted cutter at marginally stable case (ns=1000 rpm, ap=3 mm) given in [54], and b) FFT of the vibrations ................................ 148 Figure 6.9: Simulated variation of the equivalent specific process damping (ceq) in flat-end milling of AISI 1018 for the conditions given in Table 6.1........................................................ 148 xx  Figure 6.10: Comparison of the stability lobes predicted using the equivalent specific process damping (ceq) and the experimentally identified process damping coefficient (cp,r) for the parameters in Table 6.1. Work material: AISI 1018; Cutter: three fluted flat-end tool with 19 mm diameter, 300 helix angle; Immersion: 100% radial; Cutting force coefficients:ctr K [1550 480] MPa ;Structural parameters: given in Table 6.2 ......................................................................... 149 Figure 6.11: a) Representative radial displacements (vibrations) of the cutter at two different marginally stable cases (ns=2071 rpm ap=0.75 mm, and ns=2519 rpm ap=0.35 mm) for Al 7075 given in [118], and b) FFT of the vibrations ............................................................................... 150 Figure 6.12: Simulated variation of the equivalent specific process damping (ceq) in flat-end milling of a) Al 7075 and b) AISI 1050 for the conditions in Table 6.3 .................................... 151 Figure 6.13 : Comparison of the predicted stability lobes using equivalent specific process damping (ceq)  at different vibration amplitudes for four fluted flat-end mill with 12 mm diameter, 300 helix angle with the parameters in Table 6.3 and Table 6.4, a) Work material: Al 7075; Cutting force coefficients: ctr K  [1160 440] MPa for 83% radial immersion down milling; b) Work material: AISI 1050; Cutting force coefficients: ctr K  [1800 350] MPa 50% radial immersion down milling ................................................................................................... 152 Figure 6.14: a) Initial dimensions of the workpiece, and WCS and TCS, b) experimental setup..................................................................................................................................................... 154 Figure 6.15: Predicted variation of the first two dominant mode frequencies of the workpiece 154 Figure 6.16: Predicted and experimental unstable (chatter) regions, and the experimental cutting sound ........................................................................................................................................... 155 Figure 6.17: a) Experimental, and b) predicted chatter frequencies along the toolpath ............. 156 xxi  Figure 6.18: CAD model of the impeller used in the verification experiment, and the nominal finished dimensions of the blade top and bottom ....................................................................... 157 Figure 6.19: a) Experimental setup, and b) the machined part ................................................... 157 Figure 6.20: Predicted variation of the first two dominant mode frequencies of the blade, and experimental modes in Strategy 2 ............................................................................................... 159 Figure 6.21: Predicted unstable (chatter) regions on a) front and back faces, and b) and experimentally observed chatter marks in semi-finishing the blade top (front face) at 8500 rev/min ........................................................................................................................................ 160 Figure 6.22: a) Experimentally identified chatter frequencies, b) predicted chatter frequencies, c) sound data of a sample toolpath segment and its FFT in semi-finishing the blade top at 8500 rev/min ........................................................................................................................................ 161 Figure 6.23: a) Predicted chatter regions and b) the corresponding chatter frequencies in semi-finishing the blade top with the process damping effect at 1500 rpm ........................................ 163 Figure 6.24: a) Experimentally identified chatter frequencies, b) sound data of a sample toolpath segment, and c) its FFT in semi-finishing the blade top at 1500 rpm ........................................ 164 Figure 6.25: Static deflections of the tool and flexible blade, and a sample cutter contact (CC) point ............................................................................................................................................ 166 Figure 6.26: Comparison of the measured and predicted dimensional form errors.................... 170 Figure 6.27: a) Workpiece geometry, b) comparison of the predicted and measured workpiece FRFs in normal (x )W direction. Measured modal parameters in [6]: natural frequency  , 731.2 1718.6n r  Hz, modal stiffness  62.04 8.67 10rk    N/m, modal damping ratio  0.01 0.0022r  , modal mass  0.097 0.074rm  kg ................................................................. 171 xxii  Figure 6.28: a) Simulated cutting force in Wx  and its FFT, b) simulated workpiece vibration in Wx and its FFT, c) experimental cutting sound data and its FFT around the point (A) analyzed at 10300 rpm ................................................................................................................................... 174 Figure 6.29: a) Simulated cutting force in Wx  and its FFT, b) simulated workpiece vibration in Wx  and its FFT, c) experimental cutting sound data and its FFT around the point (A) analyzed at 11000 rpm ................................................................................................................................... 175  xxiii  List of Symbols cdjA   directional dynamic cutting force coefficients matrix iA   removed material between  thi-1  and thi  machining steps iA   fictitiously added material between  thi-1  and thi  machining steps 0A   amplitude of vibration  uvA   mode residue at thu  and thv  DOF aDOF  internal DOF vector of iA  aDOF  internal DOF of iA   pa   axial depth of cut m mj ka  contribution coefficient of thmk  mode shape of the in-process workpiece to the first order mode shape perturbation terms of thmj  mode  pdjΒ   directional process damping force coefficients matrix iB   in-process workpiece structure at thi  machining step  bDOF  internal DOF vector of iB   effbDOF  effective bDOF vector within the cutting region along entire toolpath bDOF  internal DOF of iB  m mj kb  contribution coefficient of thmk  mode shape of the in-process workpiece to the second order mode shape perturbation terms of thmj  mode  pb   the width of cut in orthogonal cutting, radial depth of cut in milling xxiv  CD  computational DOF domain vector condC   condition number threshold iC  experimentally identified average process damping force coefficient in the radial direction pC   plastic contact pressure constant  cDOF  coupling/interface DOF vector cDOF  coupling/interface DOF eqc   equivalent specific viscous damping in the radial direction kc  contribution coefficient of the kth degree Chebyshev Polynomial of the first kind p,rc  experimentally/numerically identified average process damping coefficient in radial direction  rc   modal damping coefficient   cond   matrix condition number operator  D   dynamic stiffness matrix D   condensed dynamic stiffness matrix Td   tool deflection vector  Wd   workpiece deflection vector  , (t, z)csrta jdF  static cutting force acting on jth flute at axial position z  , (t, z)psrta jdF   static plowing force acting on jth flute at axial position z , (t, z)cdrta jdF   dynamic cutting force acting on jth flute at axial position z xxv  , (t, z)pdrta jdF  dynamic plowing/indentation force acting on jth flute at axial position z d  summation index dz   infinitesimal axial element height  eTd   dimensional form error due to tool deflection eWd   dimensional form error due to workpiece deflection db   chip width dS   length of the differential cutting element  diag  diagonal operator taking the main diagonal of the input matrix E  known coefficients matrix in least squares problem EI  effective independence distribution vector E  modulus of elasticity eqE   energy dissipated by the equivalent specific viscous damping PDE   energy dissipated by the indentation force e  row number of the epsilon table (t, z)AF  angular feed vector (t, z)LF  linear feed vector (t, z)TF  total feed vector PDF   process damping force ,d rF   specific indentation force  f   external force vector G  Guyan model order reduction xxvi  g   modal displacement vector H  FRF matrix uvH   FRF at thu  DOF due to excitation at thv  DOF (t, z)jh  total chip thickness of jth flute at axial position z , (t, z)s jh  static chip thickness of jth flute at axial position z d, (t, z)jh  dynamic chip thickness of jth flute at axial position z I  identity matrix i  machining step counter and imaginary number j  flute index ij    iteration index ,m mj k   mode indices K   stiffness matrix ,red typeK  reduced order stiffness matrix (type = SEREP, IRS, or G) K   uncoupled block stiffness matrix , (t, z)crta jK  vector of cutting force coefficients of jth flute at elevation z defined in , ,e e er t a   , (t, z)prta jK  vector of static plowing (edge) force coefficients jth flute at elevation z defined in , ,e e er t a    K   change in stiffness matrix of the in-process workpiece   spK   indentation force coefficient SK   static stiffness xxvii  k total dynamics update steps, degree of Chebyshev Polynomial of the first kind, axial disk number ck   clamping stiffness of the holder rk   modal stiffness  L  number of low frequency modes cL   gauge length of the tool measured from its tip to the collet/tool holder wL   equivalent wear length   engagement map index m  number of master (retained) DOF M   mass matrix ,red typeM  reduced order mass matrix (type = SEREP, IRS, or G) M   uncoupled block mass matrix M   change in mass matrix of the in-process workpiece   MAP   engagement map along the toolpath rm   modal mass CN   number of nodes of the Chebyshev Polynomial fN   number of flutes kN   total number of axial disk elements ON   number of the observable modes in the updated real FRF ( , )j t zn  unit surface normal direction vector in TCS xyn   unit surface normal direction vector in Tx - Ty  plane in TCS xxviii  n  total number of DOF  sn   spindle speed  P   position vector of the tool tip ccP   position vector of the cutter contact (CC) point  ctP   position vector of the tool center p  number of the considered significant modes ep   elastic contact pressure dp    overall contact pressure due to dynamic indentation pp   mean plastic contact pressure sp    overall contact pressure for static plowing Q   unitary matrix resulting from QR decomposition pq   perturbation order q   number of the circular frequencies away from the modes 1,Tjq  , ,Tjq  previous and current vibrations of the tool in radial direction 1,j Wq  , ,j Wq  previous and current vibrations of the workpiece in radial direction , (t,z)j T Wq   relative current surface vibration velocity in the radial direction R   upper triangular matrix resulting from QR decomposition Re()  real part operator (z)R   local radius of the ball part of the ball-end mill 0R   radius of the cylindrical part of the ball-end mill xxix  (z)effR  effective local radius of the ball part of the ball-end mill eR   cutting edge radius 21strstrr   interface coupling force vector felt by str2 due to connection to str1 r  mode index in FRF construction, harmonic number in Fourier series expansion  rank   matrix rank operator , ,e e er t a  local cutting edge coordinates in radial, tangential, and axial directions S  known right hand side matrix in least squares problem dS   designed surface  s   number of slave (eliminated) DOF and Laplace domain operator str  structure T  Boolean transformation matrix ,strred typeT  model order reduction basis matrix (type = SEREP, IRS, or G) rtaxyzT  transformation matrix between ( , , )T T Tx y z and ( , , )e e er t a coordinate frames TCSWCST  transformation matrix from TCS to WCS  MCSWCST  transformation matrix from MCS and WCS kT   kth degree Chebyshev Polynomial of the first kind  (z)jT   tooth passing period for the jth flute at elevation z t  time wt   thin-walled workpiece thickness t    sampling time interval xxx  kU   kth degree Chebyshev Polynomial of the second kind  u  displacement vector due to the external force vector , dV,V Vs  displaced total, static, and dynamic volumes in process damping cV   cutting speed x   displacement vector (total DOF vector) x   acceleration vector , ,T T Tx y z  unit vectors defining TCS in WCS , ,W W Wx y z  unit vectors defining WCS in MCS x j   roots of Chebyshev Polynomial of the first kind , ,T T Tx y z  tool’s Cartesian coordinates , ,W W Wx y z  workpiece’s Cartesian coordinates , ,M M Mx y z  machine’s Cartesian coordinates dx   discretization resolution of process damping calculation points along the wave sx   discretization resolution of static plowing calculation points  y  unknowns vector in least squares problem Wavey   equation of the surface undulation y   relative normal displacement (overlap amount) in the indentation region  ,k mZ   tool compliance (flexibility) at axial elevation kz  due to the cutting force at mz   z  column number of the epsilon table, axial position on the tool axis 1z   discrete time domain operator xxxi  n   normal clearance (flank) angle 0i    helix angle of the cylindrical part of the ball-end tool i   local helix angle of the ball part of ball-end tool  Γ   the matrix formed by the mode shapes and the associated perturbation terms Γ   partial column summation matrix  n   normal rake angle o   radial rake angle )( cylo   radial rake angle on the cylindrical part of the ball end mill , (t,z)j T WΔ  vector of relative vibrations between the workpiece and the tool from previous and current vibrations on  jth flute at elevation z , (t, z)j T WΔ     vector of current relative surface vibrations between the tool and workpiece on jth flute at elevation z , (t, z)j T WΔ  vector of current relative surface vibration velocities between the tool and workpiece on jth flute at elevation z ε   epsilon vector    perturbation parameter Tol   convergence tolerance Θ   general structural dynamics matrix (FRF, mass and stiffness matrices) ,A C    angular positions of the rotary drives around Mx  and Mz  axes    axial immersion angle  Λ   diagonal eigenvalue matrix xxxii     eigenvalue (square of the mode frequency) W   wavelength    eigenvalue perturbation term     coefficient of friction    Poisson’s ratio    total of modal process damping and modal structural damping ratios  ,p r   modal process damping ratio r   modal structural damping ratio    unit step function Y   yield strength (z)j  time delay between the current and previous vibrations of jth flute at axial position z Φ   mode shape matrix    mass normalized eigenvector (mode shape vector)    eigenvector perturbation term (t)   instantaneous rotation angle of the tool ex   cut exit angle st   cut start angle  , (z)p j   pitch angle of the jth flute at elevation z    material separation angle   ( )j z   lag angle of the jth flute at elevation z xxxiii  Ω   characteristic matrix equation    circular frequency c   chatter frequency n   natural frequency s   spindle frequency , (z)T j  tooth passing (forcing) frequency for jth flute at elevation z    frequency increment xxxiv  List of Abbreviations CAD  Computer Aided Design CAM  Computer Aided Manufacturing CC  Cutter Contact Point CL   Cutter Location Point CMM   Coordinate Measurement Machine CPU  Central Processing Unit  DOF   Degrees of Freedom EI  Effective Independence EMA   Experimental Modal Analysis EVP  Eigenvalue Problem FD  Frequency Domain FE  Finite Elements FFT   Fast Fourier Transform FRF  Frequency Response Function IBR  Integrally Bladed Rotor IRS  Improved Reduction System MCS   Machine Coordinate System        PD  Process Damping QR   Orthogonal-Triangular Decomposition SEREP System Equivalent Reduction and Expansion Process SP   Separation Point  SVD  Singular Value Decomposition xxxv  TCS   Tool Coordinate System TD  Time Domain TWE  Tool-Workpiece Engagement WCS   Workpiece Coordinate System                    xxxvi  Acknowledgments I am most indebted to my research supervisor Dr. Yusuf Altintas for his patience, mentorship, motivation, and sharing his research, professional, and life experience with me during the course of my M.A.Sc. and Ph.D. studies.  I also would like to extend my special thanks to my Ph.D. examining committee members Dr. E. Hall, Dr. A. S. Phani, Dr. T. Haukaas, Dr. G.S. Schajer, and Dr. M. A. Elbestawi for their constructive comments to reveal the value of my dissertation even more.       Many thanks are due to Mr. D. G. McIntosh, Mr. F. Perron, Dr. W. Ferry, Dr. S. Engin, Mr. F. Atabey, and Dr. J. Roukema in of Machining Technologies Development Department at Pratt & Whitney Canada for giving me the change of observing the real world problems in-situ during my 6-month long industrial internship.    I have always felt distinguished by being a member of Manufacturing Automation Laboratory (MAL) throughout my long journey in the graduate school. I want to thank all my friends at MAL for having time together at the laboratory and also outside. Among all, I would like to highlight the names of Coskun, Murat, Deniz, and Alptunc for their valuable discussions and friendships. It was a valuable experience to work in such a world class environment with brilliant colleagues.   My deepest appreciation extends to my beloved parents Hulya Tuysuz and Mehmet Tuysuz, and to my elder brothers Dr. M. A. Tuysuz and Dr. F. Tuysuz for their unconditional support and love from kilometers away, which is far beyond any words. I wish to express my special thanks to my dear wife, Tugce Tuysuz, who has been walking with and encouraging me at the most stressful moments throughout this journey. This dissertation would not have been accomplished without my family’s, Tugce’s and her family’s absolute confidence in me.    xxxvii  Dedication       I would like to dedicate this work to my mother Hulya Tuysuz …  1  Chapter 1: Introduction Thin-walled components having a thickness to length ratio of less than 1/10 are widely used in the aerospace industry due to increased weight efficiency requirements. The thin structural parts seen in the wing and fuselage have pockets that are machined with three-axis toolpaths. Highly flexible blades of the integrally bladed rotors (IBRs) in the jet engines have complex geometries with high curvatures, and they are machined in five-axis operations. In both cases, 90% of the material is removed from the initial forged blanks. Cycle times of the costly IBRs change from ~10 hours to ~50 hours depending on the size, material, and geometric complexity of the part. Modeling and simulation of these machining operations in the virtual environment to predict the cutting forces, machining dynamics and vibrations, and the tolerance violations before the actual part is made is essential to improve the efficiency of the process.     IBR blades composed of ruled profiles are machined with five-axis flank milling using the side (flank) of tapered ball-end cutters as illustrated in Figure 1.1(a), while more complex twisted blade profiles having free form surfaces are machined with five-axis point milling using the ball-end of the cutting tool (Figure 1.1(b)) for the surface integrity.   Figure 1.1: Five-axis flank and point milling of blades  2  The common practice in industry to determine the most productive cutting parameters to meet the tight tolerance requirements of the thin-walled structures is centered on the experience based knowledge and iteratively testing it at the process development stage. Current Computer Aided Design and Manufacturing (CAD/CAM) packages and the machining simulation systems used in industry are purely geometry based to visualize the geometric errors and tool collision locations along the toolpath. This thesis aims at minimizing the number of trials and errors in the process planning of the low rigidity aerospace components by developing the mathematical model of the thin-walled part machining. Although IBR blades are selected as the main test bed for the introduced models, the algorithms are applicable to machining of any thin-walled parts. As the material is removed, structural dynamics of the blade are calculated to model the dynamic interaction between the cutting mechanics and the thin-walled workpiece structure. The stability lobes (productivity charts) are accurately predicted considering the interference between the cutting edge and the newly generated workpiece surface for preprocess analysis. The chatter instability is detected along complex five-axis toolpaths for postprocess analysis. It is also shown how the developed fundamental models are used to predict the vibrations and the dimensional form errors left on the machined surface. The overall scope of the research conducted in this dissertation is schematically shown in Figure 1.2.   Process modeling starts with the cutting mechanics, i.e. prediction of the chip thickness and cutting forces. In this thesis, the cutting forces considering the multi-axis milling kinematics and its effect on the chip generation are predicted by employing the generalized five-axis point milling mechanics model previously developed by the author [4].  3   Figure 1.2 Outline of the thesis project This thesis develops the mathematical model of machining the thin-walled structures in digital environment to minimize and eventually avoid the costly physical trials. One of the main issues in thin-walled blade machining is its position dependent structural dynamics. The dominant vibration characteristics of the blade structure are different at distinct locations on the blade even before machining. The continuous removal of material causes the dominant natural frequencies and frequency response function (FRF) magnitudes of the flexible workpiece to vary so as the process stability along the toolpath. Therefore, dynamics of the thin-walled part have to 4  be modeled and included in the chatter stability prediction model. Contrary to the existing methods based on multiple finite element (FE) models of the in-process workpiece at discrete stations along the path, new semi-analytical frequency and time domain models are developed by discretizing the initial workpiece only, without updating its geometry and mesh as the material is removed. The proposed in-process workpiece dynamics update models are used to calculate the workpiece vibrations and to check the chatter stability of the process along the toolpath. Also, the developed models are used to evaluate the surface errors due to the combined workpiece and tool deflections. IBRs are usually made of thermally resistant alloys such as Ti6Al4V, and they are machined at low cutting speeds. Interaction between the cutting process and the structural dynamics of the thin-walled workpiece and slender cutting tool leads to the generation of undulated workpiece surface. Interference of the cutting edge with the newly generated wavy workpiece surface at low cutting speeds creates additional damping (process damping) in the flexible milling system and improves the chatter stability. Although all the existing studies are based on empirical calibration of a process damping coefficient, a novel analytical process damping force prediction model avoiding the experimental identification is developed considering the cutting edge geometry, work material properties, vibration characteristics, and the cutting conditions in this thesis. The proposed process damping model is used to accurately predict the stable cutting limits for preprocess cutting parameter selection, and to assess the process stability along the toolpath for postprocess analysis. Thin-walled, highly flexible parts both statically deflect and vibrate during machining. Dynamic interactions between the cutting tool and flexible workpiece may cause regeneration of the chip, which leads to unstable, self-excited (chatter) vibrations. Unless avoided, chatter leads 5  to exponentially growing chip loads, cutting forces and vibrations resulting in tool failure, overloading the machine, poor surface finish, and violation of the geometric tolerance and surface finish quality constraints. In this thesis, the generalized chatter stability model for IBRs having varying dynamics and varying process damping conditions along complex five-axis toolpaths is developed in the frequency domain. Nyquist stability criterion is used to detect the unstable regions and the associated chatter frequencies for the given cutting conditions, regular or variable pitch cutters, and complex and non-uniform tool-workpiece engagements (TWEs) along the path. The rest of the thesis is organized as follows: the relevant available literature is reviewed in Chapter 2. Two new, in-process thin-walled workpiece dynamics update models based on dynamic substructuring and matrix perturbation in frequency and time domains are presented along with their experimental validations in Chapters 3 and 4, respectively. A novel analytical process damping model is developed and experimentally verified in Chapter 5. Generalized modeling of five-axis point milling dynamics and stability of the thin-walled blades using the developed models in Chapters 3-5 are introduced along with the experimental verifications in Chapter 6. Workpiece vibrations and static deflection errors are also demonstrated with thin-walled part examples. Conclusions and future research direction are given in Chapter 7.  6  Chapter 2: Literature Review 2.1 Overview The objective of this thesis is to develop the mathematical models to simulate five-axis milling of thin-walled parts in the virtual environment. The physics based models can be used to avoid the scrapping of costly aerospace parts by selecting the cutting conditions that do not cause chatter vibrations and excessive deflection marks.  This chapter reviews the relevant literature in Section 2.2 for the varying thin-walled workpiece characteristics modeling in machining. The existing machining process damping models are reviewed in Section 2.3. Section 2.4 presents the past studies related to the dynamics and chatter stability of ball-end milling operations. The chapter is concluded with a summary of state-of-the-art literature in the field in Section 2.5.  2.2 Varying Thin-Walled Workpiece Characteristics Material removal causes change in the mass and stiffness of the in-process thin-walled workpiece, which becomes more flexible in time as illustrated in Figure 2.1. Variations in the static and dynamic structural properties of thin-walled parts due to the material removal have traditionally been obtained at discrete tool positions along the cut either using the experimental modal analysis (EMA), or the full order finite element (FE) models.  The EMA method was employed in the vibration suppression of five-axis milling of turbine blades at coarse discrete positions along the toolpath [5], and for the prediction of surface error and chatter stability in peripheral milling of a rectangular plate structure [6]. The time dependency, however, was ignored since EMA was not repeated as the material is removed in both studies [5, 6]. Measuring the frequency response functions (FRFs) of thin-webs as the 7  material is removed is cumbersome and infeasible for real industrial applications such as closely spaced integrally bladed rotor (IBR) blades.  Figure 2.1: Effect of material removal on the in-process workpiece flexibility Alternatively, change in the dynamic compliance of thin-walled tubular [7-10] and slender parts [11, 12] in the cylindrical (straight) turning was investigated by modeling the workpiece at different thickness, or cut length states using multiple FE models. Researchers [9], who assumed that the vibration mode frequencies do not change during cutting, obtained a large discrepancy between the experimentally identified and simulated chatter frequencies. Others [7, 11] reduced the order of the precalculated FE models in the modal domain and interpolated the dynamic parameters along the machined length. It was experimentally shown [7] how the part goes through different stability lobes during one pass due to the varying dynamics, and the spindle speed scheduling was proposed [9] to suppress the chatter vibrations in the turning of flexible parts. Precalculation of several FE models is time prohibitive and cannot be implemented in digital machining process models. Similarly, the cutting force induced static and dynamic deflections, and the form errors in peripheral milling of rectangular and tubular thin-walled structures were predicted [13-18] by 8  either updating the workpiece geometry and remeshing the associated FE models [13-17] along the cut, or by modifying the elements’ stiffness matrices proportional to their volume change [18] to consider the effect of material removal on the static and dynamic stiffness of the workpiece. Also, the authors [19] added the effect of the removed material from the final to initial workpiece geometry using the FE method for a simple cantilevered beam geometry and minimized its static deflections by orienting the cutting forces in the direction having the highest stiffness.  The effect of varying flexible workpiece dynamics on the surface quality, chatter vibrations, and the stability predictions have been widely studied by representing the different states of the part using separate full-order [20-26] or statically condensed [27, 28] FE models in single axis peripheral milling of the rectangular thin-walled parts with flat-end tools [20-28], and the thin-floor parts with bull-nose cutters [29, 30]. Instead of considering all the identified global modes, the ones effective in the cutting region were taken into account [21, 26, 27] for the stability predictions. Variation in the natural frequencies was considered by interpolating [20, 28] the dynamic characteristics of the intermediate workpiece states along the toolpath. In some studies, vibration characteristics of the part were assumed to be constant in the finishing operation due to slight material removal [23, 25, 27-29]. Either the cutting step number [21], or the tool position [26, 27, 29, 30] is added as the third axis to the classical stability lobe diagrams to emphasize the impact of varying workpiece dynamics on the chatter predictions. Spindle speed scheduling was proposed [24, 26, 27, 30] to reduce the vibrations and to have chatter free cutting along the toolpath due to the variation in dynamic parameters of the workpiece. Alternatively, the stability lobe diagrams were generated at discrete positions along the path and conservatively 9  superimposed [29] to find a single stability lobe diagram that is valid for the entire cutting operation, which is not usually possible in practice.  In another work [31], three different numerical methods were compared to update the varying thin-wall workpiece dynamics and to predict the workpiece deflections in five-axis milling of a flexible aerospace part. The FE and particle based models were found to be only applicable for the analysis of a small portion of the workpiece due to the complex mesh modifications required many times along the path. The oscillator based method was as time prohibitive as EMA to be used in the virtual machining systems due to the required modal parameters at different stages of cutting.  In place of repeating the computationally intense full order FE analyses and the restrictive EMA, the time-dependent modes and mode shapes of aerospace parts with thin-walled pockets were predicted using the finite strip method with re-meshing in [32]. The analytical multi-span beam method was also used by reducing the three-dimensional pocket to the equivalent two-dimensional geometry to consider the effect of simultaneous wall thickness reduction on both sides of rectangular pockets on part vibrations [33]. Both methods [32, 33] are only applicable to rectangular plate structures.  As an alternative, a structural dynamics modification method was combined with FE analysis for an aerodynamic workpiece surface with free form features in [34]. A bottom-up approach was followed and FRFs of the final workpiece structure were obtained from the FE analysis, and they were modified at sampled points along the toolpath by adding the removed material in the reversed order until the initial workpiece structure was recovered. The method [34] necessitates running the geometric simulation of the material removal process first to store the removed volume information at discrete points along the toolpath. In another study [35] , the Sherman-10  Morrison-Woodbury formula was used for inversion of the modified physical system matrices. Uniform modifications were applied to the in-process workpiece for each update step while the damping was kept constant throughout machining. This method needs the dynamic characteristics of both the initial and final workpiece states, and its accuracy depends on the number of predefined cutting steps. The dominant modes and mode shapes of rectangular and curved thin-walled parts were also updated at different cutting steps in the modal domain with inevitable modal truncations [36]. In order to minimize the error propagation, the authors updated the dynamics with respect to the initial workpiece (reference) state for each material removal by defining the removed volume as the one between the initial and the in-process workpiece states. This results in repetitive assembly of large mass and stiffness matrices along the toolpath. Although it gives accurate results for the first few slight material removals, results diverge as the cutting advances.     Previous studies either employed multiple full order FE models or experimental modal analyses to account for the change in the part dynamics. In these methods, workpiece models at discrete machining steps must be obtained, which is time consuming, and also impractical (due to the part setups) to be used in the production for process planning of the chatter free toolpaths with optimal material removal rates. It is preferred to predict the workpiece dynamics during the process planning and toolpath generation stages. This thesis develops two general and computationally efficient in-process workpiece dynamics update models [1, 2] using full and reduced order dynamic substructuring and perturbation in frequency and time domains. The models [1, 2] allow obtaining the updated workpiece FRFs anywhere within the cutting region.       11  2.3 Process Damping (PD) Phenomenon in Machining The cutting edge follows a wavy trajectory in the presence of structural vibrations, which are imprinted on the machined surface. Indentation of the cutting edge into the newly generated undulated workpiece surface has extensively been studied in machining literature. Change in the magnitude of the indentation force in the normal direction is shown [37] to be the dominant source of process damping, which improves the stable depths of cut at low cutting speeds as shown in Figure 2.2.   Figure 2.2: Improved stable axial depth of cut limits in the presence of process damping Due to the inconsistency in the experimental results reported for the same material by different researchers, modeling the process damping as a function of the tool geometry, vibration frequency, work material properties, and the cutting conditions is defined as one of the main challenges in machining [38]. Sisson and Kegg [39] recognized the increased process stability at low cutting speeds and experimentally studied the effect of the clearance angle, cutting edge (hone) radius, and vibration frequency on the process damping. It was concluded that the process damping effect is moved to the relatively higher speeds for high chatter frequencies. Numerous attempts have been reported in modeling the machining process damping using the experimentally identified indentation constants. Wu [40] presented a process damping (PD) model for orthogonal cutting considering the effect of the cutting edge radius only. It is assumed 12  that the indentation force ( )PDF  is proportional to the displaced volume ( )V  by a material specific indentation force coefficient ( )spK  as,  PD spF =K V   (2.1) spK  was evaluated for AISI 1018 using a constant indentation depth of 6.35 m along with the empirical relation between the indented volume and the indentation force originally derived for penetration of a punch into a flat surface [41]. Wu’s empirical model was modified for the tools with flank wear [42, 43] and sharp cutting edges [44] both in turning and end-milling by numerically calculating the indented volume, which was decomposed [45] into the static ( )sV  and dynamic ( )dV  components in the presence of the edge  radius. PDF  was found [42] to be more pronounced towards the tool tip of the ball-end cutter due to the decreasing cutting speed. Then, Chiou and Liang [46] simplified Wu’s model [40] for turning by reducing the numerical indentation volume calculation to evaluation of a triangular area for the sharp cutting edge with flank wear based on the small vibration amplitude assumption, which was later observed to be invalid for practical applications [47].     Alternatively, the average process damping force coefficient ( )iC , process damping coefficient ,(c )p r , and spK  were directly identified from the cutting and indentation tests. The indentation test was conducted [46] by axially pressing the turning insert into the stationary workpiece, or using a series of chatter free orthogonal plunge turning tests to identify the velocity dependent iC  for the vibration frequency range of 20-120 Hz with the constant vibration amplitude of 35 m  [48]. A similar set-up was also employed [49] to identify spK  directly by evaluating the interference volume from the wave-on-wave (stable) cutting conditions in peripheral milling by oscillating the cutter with different wavelengths and with the constant 13  vibration amplitude of 30 m  in the normal direction. p,rc  was also inversely identified from the experimentally found stability limits at various low and high cutting speeds both in orthogonal turning and milling [50, 51]. spK  was obtained [52] by equating the vibration energy dissipated by p,rc  to the one dissipated by spK  with numerically evaluated indentation volume over one vibration period. Recently [53], the process damping force coefficient ( )iC  was obtained by empirically fitting the stability prediction model to the experimentally identified stability boundaries, or the total modal damping ( )  from the stable cutting tests is identified for different cutting speeds both in orthogonal turning [54] and milling [55] using the operational modal analysis. Contribution of the process damping ,( )p r  was isolated [54, 55] by subtracting the structural damping ( )r  from the identified overall damping ( ) , and it was then converted to spK . Empirical values are not general in nature and found to be valid only for the identification test conditions [52]. Instead of physical cutting tests, iC  was identified from FE simulations both for macro and micro orthogonal turning at a single cutting speed due to its computational expense for the range of cutting speeds of interest [56], and spK  was then estimated from iC  using the small vibration amplitude assumption [46].       Existing literature represented indentation of the cutting edge into the wavy workpiece surface using the experimentally or numerically identified material dependent penetration coefficients. Researchers considered the effect of the edge radius, tool wear and clearance angle individually, or the combined effect of the edge radius and clearance angle only. Therefore, the additional damping forces resulting from the contact between the workpiece and the cutting 14  edge, and their effect on the process stability have not been fully modeled yet. In this thesis, the machining process damping has been analytically modeled by evaluating the overall contact pressure distribution and the associated indentation forces occurring at the interface of the cutting edge and the finished workpiece surface. 2.4 Dynamics and Chatter Stability of Ball-End Milling Operations Dynamic interactions between the cutting tool and the workpiece might cause the chip thickness to regenerate, which leads to unstable self-excited (chatter) vibrations causing exponentially growing chip loads, cutting forces, and poor surface quality as shown in Figure 2.3. Although there has been a great deal of studies predicting the chatter vibrations in machining, only a small fraction of them are for single and multi-axis ball-end milling operations.   Figure 2.3: Surface quality with and without chatter vibrations The analytical stability prediction model of flat-end milling [57] was adopted to single axis ball-end milling either by ignoring the axial forces [42], or by reducing the three-dimensional problem to an equivalent two-dimensional system [58]. Both studies [42, 58] linearized the local cutting edge geometry and the immersion dependent cutting force coefficients, and their non-15  linear dependence on the chip thickness. Then, the two-dimensional chatter stability theory was extended to three-dimensional problems [59] considering the chip regeneration affected by the flexibility in the tool axis (z) direction. In an attempt to migrate the chatter stability models to the ball-end milling operations having the feed direction not perpendicular to the tool axis, the tool was inclined, and its stability was investigated [60] in three-axis ball-end milling of a free-form surface by defining an empirical threshold for the ratio between the numerically calculated dynamic and static chip thicknesses. The method [60] requires numerous numerical analyses to generate the stability lobes. Later, the stability lobes were analytically evaluated for the inclined ball-end tool [61] by directly averaging the time varying directional force factors and calculating the dynamic forces over one tool revolution for unit dynamic displacements.  Stability of five-axis flank milling of an impeller blade with tapered ball-end mill was assessed in the frequency domain using Nyquist stability criterion in [62]. The high discrepancy between the predictions and experiments was mainly attributed to the varying blade dynamics and process damping, which were not taken into account. In another study [63], a computationally costly FE based time domain simulation concept for chip generation, cutting forces, and thin-walled workpiece vibrations in five-axis ball-end milling was presented. Instead of generating the stability lobes, the machining stability was qualitatively judged using the simulated surface quality. Next, the three-dimensional chatter stability model of ball-end milling [59] was extended [64] to five-axis ball-end milling by inclining the tool about two axes (lead and tilt angles), and the stability lobes were iteratively obtained both in frequency and time domains. The effect of multi-frequency solution on the lobes for a small radial depth was shown to be insignificant compared to an equivalent flat-end milling case due to the geometry of the 16  ball-end cutter [64]. Efficiency of the time domain solution was recently improved [65] by semi-analytically calculating the tool-workpiece engagements (TWEs).  The effect of tool orientation angles on the five-axis process stability was also studied [66, 67] due to their effect on TWE, feed direction, and the oriented dynamics of the cutting tool. Optimum chatter free ball-end mill orientations in five-axis machining were selected through checking the chatter stability of a predefined set of orientation angles [67]. The optimum stable tool orientation was set as the one causing minimum deviation from the original tool orientation at each cutter location (CL) point along a short tool path (~40 CL points).  Past studies regarding the ball-end milling dynamics focused on generating the stability lobes for rigid rectangular workpieces assuming sharp cutting tools without the effect of process damping. There is a lack of publications in the current literature on the chatter stability of five-axis ball-end milling of free form flexible surfaces with arbitrary tool-workpiece engagement (TWE) boundaries along the tool axis. These are all considered in this thesis, and a generalized five-axis ball-end milling dynamics model has been proposed for the regular and variable pitch cutters while the varying in-process thin-walled workpiece dynamics and process damping are simultaneously taken into account along the lengthy toolpath.  2.5 Summary Previous studies related to different aspects of this thesis are neither computationally efficient nor general in terms of the model capabilities and experimentally identified constants, and thus they are not suitable for machining process simulation and optimization for thin-walled aerospace structures in the digital environment. In this thesis, the generalized semi-analytical and analytical comprehensive models to predict the varying thin-walled workpiece dynamics, the process damping forces, the workpiece vibrations, the surface errors due to the static deflections, 17  and the chatter stability along complex and long toolpaths have been developed by considering the industrial applications. 18  Chapter 3: Modeling of Varying Workpiece Dynamics in Frequency Domain 3.1 Overview A large number of monolithic aerospace parts having a high strength to weight ratio are manufactured from blanks by milling more than 90% of the excess material. In spite of the advantages, thin-walled parts bring more challenges into the machining process due to lack of dynamic stiffness. The flexibility of these parts causes severe forced vibrations and chatter, which lead to poor surface finish and tolerance violations. Hence, the mathematical model of the process dynamics is needed to simulate the machining operations, and to predict the chatter free cutting conditions ahead of costly physical trials. However, prediction of the process dynamics in machining of the parts with low rigidity requires efficient algorithms that can handle the position dependent structural dynamics of the workpiece as the tool removes material along the toolpath. The position dependency results from excitation of different dominant vibration modes as the point of application of the cutting forces on the workpiece changes along the toolpath. Also, continuous material removal causes the dominant natural frequencies and corresponding frequency response function (FRF) magnitudes of the flexible workpiece to vary with time, which in turn affects the process stability and surface errors. In this chapter, a new frequency domain semi-analytical dynamic substructuring model is introduced to predict the varying dynamics of the thin-walled workpiece. It is a top-down hybrid model starting with evaluation of the undamped FRFs of the initial global flexible workpiece obtained from an FE system, and the FRFs are analytically updated along the toolpath as the material is removed (Section 3.2). Calculation of the workpiece FRFs is required only once for the unfinished part at the very beginning of cutting. Each removed volume is obtained by partitioning the initial global workpiece structure and defined as its substructure, which 19  eliminates the need for multiple FE models of the in-process workpiece to obtain its updated FRFs. Without updating the geometry and mesh of the part, the undamped dynamic stiffness of the removed substructure is analytically decoupled from the workpiece FRF at discrete toolpath locations (Sections 3.2.1-3.2.3) while compensating the effect of the truncated high frequency modes (Section 3.3). The small modal damping ratio is inserted into the updated, undamped workpiece FRFs (Section 3.4) afterwards. The model is further extended to the reduced order dynamics update to improve the computational efficiency (Section 3.5). The most compliant nodes are identified at the tool and part contact zone (Sections 3.5.1-3.5.2), and their FRFs in the directions affecting the dimensional accuracy and chatter stability the most are analytically predicted. The developed full and reduced order in-process workpiece dynamics update models are experimentally validated on a representative rectangular thin-walled workpiece (Section 3.7). The chapter is concluded with a summary in Section 3.8. 3.2 Modeling of Position Dependent Dynamics of Thin-Walled Workpiece Unlike the time invariant tool dynamics, the structural dynamics of the part change with the material removal as illustrated in Figure 3.1. The proposed model takes the FE discretization of the unfinished blank in the form of nodes and elements with the associated physical system (stiffness and mass) matrices. The eigenvalues (natural frequencies) of the initial undamped system are solved and a dominant range is retained. The part thickness dependent damping ratios are added to the system.  As shown in Figure 3.2(a), the geometric material removal between  1thi -  and thi  cutting steps ( 1... )i k  can be described by a decrease in the workpiece mass with an updated geometry, which is commonly used for validation of machining operations. In the classical part dynamics 20  update approach used in machining, mass and stiffness of the in-process workpiece at different cutting steps are used to obtain the associated dynamics by solving the eigenvalue problem of the full order FE models at multiple locations along the toolpath.  Figure 3.1: Time invariant tool dynamics and varying workpiece dynamics at different stages of machining  Figure 3.2: (a) Geometrically removed material, (b) dynamically removed material equivalent to coupling of a fictitious substructure 21  As opposed to the conventional method, the model presented in this study cancels the dynamics contribution of the removed material within the dynamics of the in-process workpiece without updating its geometry using the concept of negative structural modification [68, 69]. An equivalent fictitious substructure having the same geometry but the opposite dynamics of the removed material is added to the dynamics of the in-process workpiece. This transforms the material removal (dynamics decoupling) into a coupling problem as shown in Figure 3.2(b), which is dynamically equivalent to solve the updated dynamics of the remaining workpiece at any machining step, but with less effort. The model can be mathematically described in the time domain as,  1( ) ( ) ( )ii iBB At t t    Θ Θ Θ   (3.1) where t is time, 1Θ iB and Θ iB  represent the time dependent dynamics, or the physical system matrices of the machined workpiece at discrete cutting steps ( 1)i -  and i , respectively. Θ iA  is the negative dynamics of the removed material between the cutting steps ( 1)i -  and i . The removed volume ( )iA  between two consecutive dynamics update steps is defined as a substructure of the in-process workpiece 1( )i -B  as shown in Figure 3.3(a), where the nodes are partitioned as the coupling (i.e. interface), the internal nodes of iA , and the internal nodes of the remaining workpiece )( iB . Each node has three translational degrees of freedom (DOF), and the corresponding partitioned DOF sets are cDOF for coupling, aDOF and bDOF for the internal nodes of iA  and iB , respectively.  22   Figure 3.3: (a) In-process workpiece i-1B  at ( )thi -1  machining step, (b) component level interaction between the substructure iA  and the rest of in-process workpiece i-1B  represented by rigid links, (c) component level interaction represented by interface coupling (reaction) forces, (d) opposite interface coupling forces generated by the fictitiously added iA , and interface force cancellation after addition of iA  23  Interface interactions between the rest of workpiece i-1B  and substructure iA  within i-1B  can be represented by the rigid links connecting the interface DOF of both structures. This provides interface displacement compatibility and force equilibrium conditions as shown in Figure 3.3(b). The rigid link can be replaced with an interface coupling (reaction) force vector 21strstrr  at substructure level [70] denoting the forces felt by substructure 2 (str2) due to its connection to substructure 1 (str1) as illustrated in Figure 3.3(c) while maintaining the necessary coupling conditions. Although each node has three DOF, 21strstrr  is shown only in a single direction for brevity. Similarly, the dynamics decoupling is illustrated in Figure 3.3(d) with a fictitiously added substructure iA . Physically, coupling of the fictitious substructure with the current workpiece i-1B  provides the sufficient conditions for dynamics removal (decoupling) such that, the opposite dynamic stiffness eliminates the dynamics contribution coming from the removed volume iA  within the structure i-1B , and also the fictitious substructure generates negative (de)coupling forces at the interface to cancel the interaction between the removed volume iA  and the workpiece i-1B , which is detailed for an interface node in Figure 3.3(d) for all three directions.  Dynamics of the initial workpiece structure ( )0B  in Eq.(3.1) are expressed as receptance FRFs by solving the following eigenvalue problem,   0 0 00 0m m mB B BB Bj j jK M    (3.2) where the superscript denotes the structure, 0mBj  and 0mBj  are the thmj  eigenvalue (proportional to the square of the mode frequency, 0, mBn j ) and mass normalized eigenvector of 0B , in order. The 24  mass and stiffness matrices 0 0( , )B BM K  are obtained from the FE model of the initial workpiece. This step is required only once for the workpiece before material removal starts, and recursive implementation of Eq.(3.1) yields the updated receptance FRFs of the workpiece along cutting, which is detailed as follows. Henceforth, FRF refers to the receptance FRF throughout this thesis unless stated otherwise.   The undamped equations of motion of the in-process workpiece structure -1( )iB  can be written in the time domain as,  1 1 1 11 1( ) ( ) ( ) ( ) ( ) ( )i i i ii iB B B BB Bt t t t t t      M K f rx x   (3.3) where 1x iB  and 1iB x  are displacement (total DOF) and acceleration vectors, 1f iB is the external force vector, and 1r iB is the reaction (coupling) force vector acting only on the cDOF due to connection to the fictitious substructure. Eq.(3.3) can be rewritten in the frequency domain as,   1 1 1 11 1 1 1( ) ( ) ( ) ( )( ) ( ) ( ) ( )i i i ii i i iBB B BB B B B            x H f rx u H r  (3.4) where   is the circular frequency, 1 ( )H iB  is the FRF matrix, and the displacement vector due to the external force vector is,  1 1 1( ) ( ) ( )i i iBB B    u H f   (3.5) The undamped equations of motion for substructure iA  can be written using the physical system matrices of the substructure iA ,    ( ) ( ) ( ) ( ) ( ) ( )i i i ii iA A A AA At t t t t t       M K f rx x   (3.6) 25  where superscripts iA  and  iA  denote the removed and fictitious materials, respectively. MiA , K iA are also written as time dependent due to non-uniform material removal along the toolpath. iAf  is the external force vector acting on  iA , iAr  is the reaction (coupling) force vector acting only on the cDOF due to connection of  iA  to -1iB . Eq.(3.6) can be rewritten in the frequency domain as,  ( ) ( ) ( ) ( )i i i iA A A A       D x f r   (3.7) where 2( )i i iA A A   D M K  is the dynamic stiffness matrix. The last expression in Eq.(3.4) can be combined with Eq.(3.7) as,  1 1 11 ( )( )( ) ( ) ( )( ) ( ) ( )                                          r u xH 00 D x f ri i iii i i iB B BBA A A A  (3.8) Since the overall DOF 1( , )i iB A x x  and the interface force vectors 1( , )i iB A r r  appear on both sides of Eq.(3.8), the coupling with the interface compatibility and equilibrium cannot be carried out directly. Therefore, the structural dynamics of -1iB  and  iA  given in Eq.(3.8) are expanded to an uncoupled augmented form [71] with explicit DOF partitioning in the global set as, 26  11 111 11 1 1. . . .. . . .. . . . .. . . . . . . . .. . . . . . . . . .. . . . . . . . . .. . . . . . . . . .. . . . . . . . .. . . . . . .. . . . . .ii iii ii i ii ii iBB Baa ac aaabBB Bca cc cccbB B Bbbba bc bbA Aaa aa acA Acc ca ccI           H H H I 0 0H H H 0 I 0H H H I 000I D D0 D D1 11 11 11 11 11 1i ii ii ii ii ii ii ii ii ii iB Ba aB Bc cB Bb bB Ba aB Bc cB Bb bA Aa aA Ac cA Aa aA Ac c                                                             r ur ur ux 0x 0x 0r 0r 0x fx f (3.9) where I  is the identity matrix,   is dropped for clarity, and subindices of the matrices and the vectors show which DOF sets they are defined for. The global DOF set includes the interface displacements and the coupling forces ( , )c cx r  twice due to the augmented formulation given in Eq. (3.9). The duplicated interface information is eliminated from the global set by the following signed Boolean transformation [71] owing to partitioning a global workpiece structure, in turn the nodal compatibility between the in-process workpiece and the removed volumes. This transformation (T) also implicitly satisfies the interface compatibility and equilibrium required for the coupling. 27   111111.... . . . .. . . .. . . .. . . . . .. . . . . .. . . . .. . . . .. . . .. . . . . .. . . .iiiiiiiiiiBaBaacccBbbbBa aaBcccB bbbaaAaccAc aaAccaAc                    rI 0r0 I 0r0 I 0x I 0I 0x0 I 0x0 I 0rI 0 0 0r 0 I0 I 0 0xx11 11 11 11 ii ii ii ii ii iiB Ba ac cB Bb bB Ba ac cB Bb bA Aa aA Aa a                                                                             r rr rr rx xTx xx xr rx xT  (3.10) where 1i iB Ac c c   r r r  shows the distribution of coupling forces along the interface DOF and satisfies the interface equilibrium; and 1i iB Ac c c  x x x  is the interface displacement vector satisfying the interface compatibility. Insertion of the unique set of DOF defined by Eq.(3.10) into the augmented set of un-assembled equations (Eq.(3.9)) and the pre-multiplication with TTyields the coupled (i.e., assembled) set of equations as,  11111 111 11 1 1. . .. ... . . . . . .. 0 . .. . . . . . . .. . . . . . . .. . . . .iiiii iii ii i ii ii iBBB Baaa ac aaabBB B cca cc cccb BB B B bbbba bc bb BaA Acc cc caA Aac aa aa                     rH H H I 0rH H H 0 I 0rH H H 0 0 I 0x0I 0 D DD I D11111 1iiiiii ii ii iBaBcBbBaAc cB Bb bA Aa aA Aa a                                               000uuux fxrx f  (3.11) 28  Since the coupling is carried out through the interface DOF, the coupling forces act only on the cDOF, and thus the reaction forces on aDOF and bDOF in Eq.(3.11) can be replaced with zero vectors, i.e. 1i iaB Aa a  r r 0  and 1i bBb r 0 . Then, the final set of coupled equations becomes,  11111111.....iiiii iiiiiii ii iBBacac aaBBBcacc ccBBc bbbbcB AA Acbcc cc caAAA Aaaac aa                                                  urH I 0 0uxH 0 I 0x uH 0 I 0fxI D 0 Dfx0 D 0 D  (3.12) The linear set of matrix equations given in Eq.(3.12) is explicitly solved at the interface cDOF and the internal bDOF at each frequency to obtain the updated FRFs of the remaining workpiece structure ( )iB . By eliminating the interface forces ( )cr  using the Gaussian Elimination method combined with Singular Value Decomposition (SVD) as a stable matrix inversion technique, the numerical stability of the calculations are improved especially at the anti-resonances of the substructure iA  since the coupling is generally sensitive at these frequencies [68, 69]. Since the dynamic contribution of the removed substructure iA  is cancelled within the in-process workpiece structure 1iB  , and aDOF of the fictitious substructure are no longer within iB  after the decoupling, iAcf , iAaf , and 1iBaf  are set to zero vector in Eq.(3.12). Then, the isolated direct and cross dynamics of iB  are evaluated as follows. 3.2.1 Updated Direct FRFs at the Interface DOF - iBccΗ  Since the cutting forces act only at the interface between i -1B  and iA  where the material is removed, the updated FRFs at the interface cDOF are required for the prediction of chatter 29  stability and vibrations. The interface force vector can be evaluated from the fourth and fifth equations in Eq.(3.12) as,  1i i i i iA A A A Ac cc ca aa ac c c       r D D D D x D x   (3.13) where D  is the condensed dynamic stiffness matrix at cDOF. The updated direct FRFs at interface cDOF can be found by inserting Eq.(3.13) into the second equation in Eq.(3.12) and setting 1iBbb f 0  as,  1 1 11 111i i iii ii iB B BAcc cc cc c cB BB Acc cc cc cc         I H D H f xH I H D H  (3.14) 3.2.2 Updated Cross FRFs between the Interface and Internal DOF - iBcbH  The interface cDOF and the updated part’s bDOF in decoupling vary along the toolpath, and the current cDOF at machining step i might have common DOF with the subsequent dynamics update steps due to the continuous material removal. The updated cross FRFs between the interface cDOF and internal bDOF are evaluated by letting 1iBc c f 0 . Substitution of Eq.(3.13) into the second equation in Eq.(3.12) yields,   1 111111i ii iiii iB BB Acc cc ccb bBBB Acc cccb cb       I H D H f xH I H D H  (3.15) 30  3.2.3 Updated Direct FRFs at the Internal DOF - iBbbH  Since some bDOF at machining step i become cDOF after some material removal, the updated FRFs at the internal DOF of iB  are evaluated by keeping 1iBc c f 0 . Rearranging Eq.(3.13) and inserting it into the second equation in Eq.(3.12) gives the interface forces as,  1 11111ii ii iAc cB BB Ac cc cb b          x D rr H D H f  (3.16) Substituting Eq.(3.16) into the third equation in Eq.(3.12) yields,  1 1 1 1 111 1 111111i i i i ii ii i i ii iB B B B BB Accbb bc cb b bB B B BB Accbb bb bc cb                        H H H D H f xH H H H D H  (3.17) After some matrix manipulations, Eq.(3.17) can be written in a form similar to Eqs.(3.14)-(3.15) as,  1 1 111i i i iii iB B B BBA Acc ccbb bb bc cb      H H H D I H D H   (3.18)  The recursive implementation of Eqs.(3.14),(3.15), and (3.18) as the cutting advances leads to the updated FRFs of the machined workpiece at different stages of machining. Advantages of the proposed in-process workpiece dynamic update model can be listed as: (1) Unlike the existing methods, partitioning the global workpiece structure at the beginning and removing its own substructures as the cutting advances avoids multiple FE models, automatically results in nodal compatibility at the interfaces, and allows rigid coupling; (2) the proposed model 31  automatically provides mode discrimination in the cutting zone for vibrations and chatter stability predictions as the updated FRFs reveal the dominant dynamics.  Although the use of dynamic stiffness representation for the removed (or fictitious) substructures avoids dynamics truncation on them, evaluation of 0BH  ( 1)i   is still subject to modal truncations. Hence, 0BH  is compensated for the contribution of high frequency modes in Section 3.3. 3.3 Effect of the Truncated Dynamics Although the dynamics update model employs FE based FRFs of the initial workpiece structure ( )0B  and updates them as the material is removed, it uses the modal model to evaluate the FRFs of 0B  with a limited number of modes. Despite the computational advantage, the modal incompleteness adversely affects the accuracy of the coupling/decoupling [72], in turn the modes of the remaining structure [73]. Therefore, FRFs of the structure 0B  are compensated for the frequency dependent dynamics contribution of the higher order modes without explicitly solving them using the physical system matrices and the first L known low frequency modes and mode shapes of 0B  as follows.  FRF of the initial workpiece 0B  can be expressed as superposition of the low and high frequency undamped components as,  0 00 0 012B BB B BL H       H H H M K   (3.19) where 0BLH  and 0BHH  represent the contribution of the first L low frequency modes and the unknown truncated high frequency modes, respectively. 0BLH  can be expressed as,  32   0 0 0 012TB B B BL L L L       H Φ Λ I Φ   (3.20) where 0BLΦ  is the mode shape matrix, and 0BLΛ  is the corresponding diagonal matrix having the first L eigenvalues of 0B   0 0 02 2 2,1 ,2 ,( ) , ( ) , ..., ( )B B Bn n n L   . Inserting Eq.(3.20) into Eq.(3.19) gives,  0 0 0 00 0112 2TB B B BB BH L L L              H M K Φ Λ I Φ   (3.21) The matrix inversions in Eq.(3.21) can be replaced with their power series expansions, and the frequency dependent truncated contribution is obtained as [73, 74],  0 0 0 020dB B B BdH H0 H0d  H H M H   (3.22) where 0BH0H  is the static compensation term evaluated from Eq.(3.21) at 0   as,  0 00 0 01 1B BH0 LTB B BL L            H ΦK Λ Φ   (3.23) which cannot provide sufficient compensation for the entire frequency range of interest, and thus the higher order modes are considered using Eq.(3.22). The first two to five terms of the series are shown [74] to be sufficient when the following convergence condition [75] is satisfied,  ,upper n L    (3.24) where upper  is the upper boundary of the frequency range of interest for the workpiece FRF update. The number of low frequency modes (L) in the evaluation of 0BLH  can be determined using Eq.(3.24). The truncated FRF 0( )BHH  is calculated only for the DOF in the cutting region (Section 3.6), and it compensates for the higher order dynamics of 0BH  in Eqs.(3.14),(3.15), and (3.18). 33  A three-step uniform material removal in flank milling operation with a large axial depth of cut (150 mm) is shown in Figure 3.4, where the dynamics are updated after removal of a slice of material. The undamped in-process workpiece FRFs are obtained using a separate FE model for the workpiece state after each material removal, and also using the model introduced in Section 3.2 with and without compensation of the truncated workpiece dynamics. Comparison results are presented for the same position on the workpiece after each material removal step in Figure 3.5. While the model cannot accurately track the modes of the thin-walled workpiece particularly after the first material removal step due to the missing high frequency dynamics information of 0B , it accurately predicts the in-process workpiece FRFs when the truncated dynamics are included.   Figure 3.4: Three-step uniform material removal in flank milling of a rectangular thin-walled workpiece 34   Figure 3.5: Comparison of the undamped FRFs with and without the truncated workpiece dynamics compensation in material removal steps 3.4 Insertion of Damping Material removal does not only cause changes in the mass and stiffness, but also affects the damping. The researchers [28, 29, 33, 34] commonly assumed constant and uniform modal damping for all modes of a thin-walled part throughout machining. However, in this study, the cutting tests on the thin-walled workpiece showed that modal damping of the dominant modes change with the workpiece thickness ( )wt . FRF measurements before machining, after roughing, after semi-finishing, and after finishing workpiece states are used to identify the thickness dependent modal damping ratios ( ( ))r wt . The intermediate values of ( )r wt with respect to the material removal step are obtained with interpolation. 35   In order to insert the experimentally identified modal damping into the predicted (updated) undamped real FRFs given in Eqs. (3.14),(3.15), and (3.18), the mode residues at the DOF of interest are required. The updated damped FRFs can be expressed as,   , ,2 2 2 21 1, , , ,( )2 2O OuvuvN Nu r v r rr rn r r n r n r r n rAHi i                 (3.25) where ( )uvH   is the workpiece FRF at thu  DOF due to excitation at thv  DOF;   is the unknown mode shape vector; n  is the known undamped natural frequency from the updated FRFs; r  is the known modal damping ratio; and  uv rA  is the unknown residue of the thr  mode at thu  and thv  DOF. Eq.(3.25) can be expressed in real and imaginary form as,  22 2, ,2 2 2 2 2 2 21 1, , , ,( )( )( )( ) ( ) ( ) 2( ) (2 ) (2 )O OuvN Nuv r n r uv r r n rr rn r r n r n r r n rHA Ai                      (3.26) where i  is the imaginary number. The contribution of damping in FRF (Eq.(3.25)) decays rapidly at the frequencies away from the natural modes for the lightly damped thin-walled parts as schematically shown in Figure 3.6. Hence the predicted undamped workpiece FRFs are assumed to be equal to that of the real part of the damped workpiece FRFs (Eq.(3.26)) at these frequencies. Then, the real FRFs of the thin workpiece can be approximated at the frequencies away from the natural modes as,    2 2,2 2 2 2, ,2 2,1 1( )( ) (2 )( ) ( ) ( )O Ouvn rr rn r r n rN Nuv n r uvr rRe HA A             (3.27) 36   Figure 3.6: Damping dominant regions on the real part of a lightly damped thin-walled workpiece FRF The following empirical criterion is used for each mode in the updated undamped FRFs for the known r  to determine the frequencies where the approximation in Eq.(3.27) holds,    22 2,2,10 ,21.... On rr n rr N       (3.28) The number of frequencies away from the modes is usually much higher than the number of observable modes ( )ON  in the updated real FRF of the DOF of interest. An overdetermined set of equations is defined to evaluate the unknown residues  uv rA  in Eq.(3.27) as,          1 22 2 2 2 2 2,1 ,2 ,N2 2,1( ) .....OOOuv Nuv uv uvuvn n nn rNrrAA A ARe H               (3.29) which can be rewritten in matrix form as, 37  1 1 1 1 1211 12 13 11 2 3( ) ( ) ( ) ( ) ( )( ). . . ... ... . . ... ... . . ... ... . . ... .( )( ) ( ) ( ) ..... ( ).....OOOuvuvuv Nq q q qNq Nq q qE E E A RAAEE E E E                                          yE12( ( ))( ( )).,..( ( ))uvuvOuv qe HRe Hq NRe H           S  (3.30) where q  is the number of circular frequencies away from the modes obtained from the updated undamped FRFs using Eq.(3.28). Eq.(3.30) is solved by defining the following least squares problem,  22minyEy SEy S   (3.31) Once the solution for  uv rA  is obtained, the damped updated FRFs are evaluated using Eq.(3.25) at any stage of machining. Inserting the damping into the updated undamped FRFs as proposed in this section is advantageous: (1) the overall FRF update simulation becomes more computationally efficient as it is run only for the undamped conditions; (2) insertion of the small modal damping only for the FRF update DOF in the cutting region is allowed for any step of machining. Although the proposed flexible workpiece FRF update model analytically decouples the dynamics of the removed volumes along the toolpath, its efficiency is further improved by reducing the size of the initial workpiece FE model and its removed substructures as detailed in Section 3.5.       38  3.5 Reduced Order Dynamics Model of the Thin-Walled Workpiece The initial workpiece ( )0B  and the removed materials ( )iA  are reduced without modifying their original FE discretizations. The most flexible nodes and their DOF, which are mostly in the vicinity of the unsupported, free end of the thin-walled workpiece, are selected as the master (retained) degrees of freedom set (m) for 0B  and iA  using the following criterion.  3.5.1 Selection of Master Degrees of Freedom In the classical dynamic substructuring applications, the master DOF for each substructure is determined independently [76]. However, the removed volumes (fictitious substructures) in the developed top-down dynamics update model are obtained by partitioning the initial global workpiece structure. This makes the workpiece and the removed materials dependent structures with the nodal conformity (compatibility), which also needs to be preserved after the model order reduction as illustrated in Figure 3.7.  Figure 3.7: Nodal compatibility between the in-process workpiece and fictitious material, a) before, and b) after model order reduction 39  Due to preservation of the nodal compatibility, the master DOF for iA  are determined first, and they are kept in the master DOF set of 0B  using the iterative Effective Independence (EI) method [77], which is adapted by replacing the experimental mode shapes with the theoretical (FE based) ones. To eliminate the DOF having the least contribution to the modal kinetic energy and the spatial independence of the modes, the effective independence distribution vector ( )EI  [78] for the measurement noise free conditions is defined as,  1ii i i iT Tstr str str str strj np np np npj j j jdiag                      EI Φ Φ Φ Φ   (3.32) where ij  is the iteration number, the superscript str denotes the structure type  ,0 istr B A ;  diag  is the diagonal operator taking the main diagonal of the input matrix; the subindices n and p are structure’s total number of DOF and number of the significant modes considered, respectively. In order to simplify Eq.(3.32), the mode shape matrix in thij  iteration istrnpj     Φ  is decomposed [78] as,  i i istr str strnp np ppj j j          Φ Q R   (3.33) where strnpQ  and strppR  are the unitary and upper triangular matrices resulting from the orthogonal-triangular decomposition (QR) [75], which redistributes [78] the modal kinetic energy retained by the master DOF through reorthonormalizing the retained mode shape matrix in every iteration. Inserting Eq.(3.33) into Eq.(3.32) yields,  ii iTstr str strj np npj jdiag          EI Q Q   (3.34) 40  The selection scheme is started with treating all n DOF as potential masters, and strEI  is sorted in descending order after each iteration. The DOF constituting 90% of the overall strEI  are kept and carried to the next iteration as new potential masters. Unlike the original form of the EI method [77], the number of master DOF (m) is not known a priori for both 0B  and iA . Hence, the condition number and rank of the mode shape matrix istrnpj     Φ  are monitored after DOF elimination in each iteration and used as the stop criterion as,  ,i istr strnp cond npj jcond C rank p                Φ Φ   (3.35) where  cond  and  rank  are the operators evaluating the condition number and rank of the input matrix, and condC  is the predefined condition number threshold. Condition number and rank are measures of the linear independency between the columns of istrnpj  , and they physically qualify how well the remaining portions of the mode shapes in the iterative DOF elimination process are spatially independent. 3.5.2 Model Order Reduction The goal of model order reduction is to replace the full order displacement vectors strx  of the initial workpiece structure ( )0B  and the removed materials ( )iA  given in Eq.(3.3) and Eq.(3.6) with their subsets that are composed only of the selected master DOF (m) in Section 3.5.1. While retaining the significant dynamic characteristics of the original systems, the model order reduction is mathematically expressed as, 41   strmstr str strred mstrs     xx xxT   (3.36) where the subindex s is the number of the slave (eliminated) DOF, and strredT  is the reduction basis matrix of the size ( ) mm s  . The full physical system matrices ( , )M K  and the external force vector ( )f  are also partitioned to eliminate the slave DOF, and thus the undamped equations of motion become,   strmstrsstr str str str str strmm ms mm ms m mstr str str str str strsm ss sm ss s s                                        xxM M K K x fM M K K x f  (3.37)  The reduced equations of motion can be obtained by inserting Eq.(3.36) into Eq.(3.37) and premultiplying by  TstrredT , which projects the full DOF set (n) onto the subspace defined by the master DOF (m). The reduced physical system matrices of the size m m  become,  ,T Tred red red red red redmm ms mm mssm ss sm ss             M T T K T TM MM MK KK K  (3.38) The reduction method used in the dynamics update model is desired to be as accurate as possible while preserving the dynamic characteristics of the in-process workpiece at each update step. The projection basis for 0B   0BredT  is obtained by employing System Equivalent Reduction and Expansion Process (SEREP) [79] method. The original full DOF set of 0B  is rewritten in the modal coordinates as,  0000 0BBmpm BpB Bs sp           Φxgx Φ  (3.39) 42  where 0 0,B Bmp spΦ Φ  are modal matrices, and g  is the displacement vector in the modal coordinates. The modal displacements can be found in terms of the master DOF as,  0 0 0B B Bp mp m  g xΦ   (3.40) where mp  Φ  is the generalized inverse of the rectangular modal matrix, and it is evaluated as [80],  0 0 0 00 0 0 011,,T TB B B Bmp mp mp mpT TB B B Bmp mp mp mpm pm p                            Φ Φ Φ ΦΦ Φ Φ Φ  (3.41) Substitution of Eq.(3.40) into Eq.(3.39) gives 0SEREPBT  as,  00 00 00,red SEREPB mm mmm B BmB BBsp mp sp mps                           xx TΦ ΦxI IΦ Φ  (3.42) In order to overcome the dynamics truncation in Eq.(3.42), the iterated Improved Reduction System (IRS) [81] method with Guyan Reduction [82] is used for reduction of the removed substructures ( )iA  as,   11, ,1, ,,i i i i ii i ii ii i iA A A A AA A Ass sm ssIRS G G red G red GmmmmA Ared IRS red GA A Ass smIRS                   t t K M M t M KIIT Tt K K  (3.43) where ,iAred GT  and ,iAred IRST  are Guyan and IRS reduction bases; ,iAred GM  and ,iAred GK  are Guyan reduced system matrices of iA . Once IRS reduced matrices  , ,IRS,i iA Ared IRS redM K  are obtained 43  using Eq.(3.38) with ,iAred IRST  from Eq.(3.43), the reduction is iteratively improved by replacing ,iAred GM  and ,iAred GK  in Eq.(3.43) by ,iAred IRSM  and ,IRSiAredK  as [81],   11111,IRS ,IRS,i i i i ii i ij j j ji i i iiijjiA A A A AA A Ass sm ssGIRS IRS red redmmAAred IRSIRS             t t K M M t M KITt  (3.44) where ij  is the iteration number. Iterations continue until the convergence criterion is satisfied. The convergence of mode m of the reduced system is used as the criterion,  , 1 ,,( ) ( )( )i iin m j n m jToln m j     (3.45) where ,( ) in m j  is the mode m  of the reduced system in thij  iteration, and Tol  is the prespecified convergence tolerance.  The method leads to the reduced mass and stiffness matrices (Eq.(3.38)), which belong to the nodes contributing to the displacements most in the dominant mode shapes of the thin-walled workpiece at the tool and part contact zone. The dominant master (i.e. flexible) DOF are kept as they were in the original dense mesh, whereas the more rigid slave DOF are less densely populated for the milling process simulation.  3.6 Simulation Algorithm The overall solution procedure of the developed model has geometry, dynamics, and decoupling calculations as shown in Figure 3.8. Geometric calculations start with FE discretization of the initial workpiece ( )0B  in such a way that the different radial depths of cut in roughing, semi-finishing, and finishing operations are considered. The material between two 44  consecutive update locations is defined as a removed substructure, and its structural matrices used in Eqs. (3.14),(3.15), and (3.18) are obtained from the initial FE model of the blank.  Figure 3.8: Flow chart of the frequency domain in-process workpiece dynamics update model 45  The geometry module includes identification of DOF vectors of the removed substructures (and the corresponding fictitious substructures) as iAaDOF  and iAcDOF , as well as the machined workpiece as 1iB aDOF , 1iB cDOF , and 1iB bDOF  that vary along the toolpath. In order to increase the computational efficiency by localizing the dynamics decoupling to the cutting region, and to consider the interface DOF and the decrease in the number of total DOF with material removal along the entire toolpath, the varying computational DOF domain vector  CD  for each decoupling is defined as,      1 1 11 1 1...., 1...i i i i ki i iB B B B BeffB B Beffi k          bDOF bDOF cDOF cDOF cDOFCD cDOF bDOF  (3.46) where  1iBeffbDOF is the effective bDOF vector of the in-process workpiece within the cutting region along the whole path;  and   are union and intersection operators, respectively. CD  gives the minimum number of DOF whose FRFs need to be updated throughout machining. Since the interface DOF between the initial workpiece ( )0B  and the removed (fictitious) materials ( and )i iA A  have to be in the same order in both 0B  and iA  for Eqs.(3.9)-(3.18) to hold, the DOF of iA  and 0B  are rearranged as iAaDOF  and iAcDOF ; 0BcDOF  and  0BeffbDOF , respectively. Further rearrangements for the machined workpiece 1( )iB   are carried out at FRF level in decoupling calculations as cutting advances. In the dynamics module, the undamped FRFs of the initial workpiece ( )0B  are evaluated at each j  only once, and the frequency independent portion of the truncated dynamics 0( )BH0H  is evaluated for the full order simulation. For the reduced order dynamics decoupling, the first 0BN  46  and iAN  modes of 0B  and iA  are conventionally considered as significant modes in the iterative master DOF selection in Eq.(3.34). Model order reductions are performed as detailed in Section 3.5.2, and the required initial FRFs and, if needed, 0H0BH  are calculated. In the last phase, the frequency domain dynamics decoupling of all the removed substructures ( , 1... )iA i k  are carried out for one frequency ( )j  using Eqs. (3.14),(3.15), and (3.18), and the updated FRFs of the desired points are stored. The simulation then moves to the next frequency 1( )j   within the range of interest. Once the undamped updated FRFs are obtained the modal damping is inserted as detailed in Section 3.4. 3.7 Experimental Verification: Thin-Walled Rectangular Workpiece Cutting Test The proposed FRF update model has been verified with experimental measurements and simulations along the peripheral milling toolpath of a rectangular thin-walled workpiece made of Al7050. The thin-walled part thickness has been reduced from 8.25 mm initial uniform value to 4.75 mm after roughing, semi-finishing, and finishing operations with 1.75 mm, 1.1 mm, and 0.65 mm layers, respectively.  FE model of the initial workpiece is generated from its CAD geometry (Figure 3.9(a)) using 8-noded solid brick elements as illustrated in Figure 3.9(b). The assigned material properties for Al7050 are: modulus of elasticity of 71.7 GPa, density of 2700 kg/m3 , and Poisson’s ratio of 0.33. After checking FE mesh convergence, the initial full FE model is finalized with 46980 total DOF. 47   Figure 3.9: Initial thin-walled rectangular workpiece (a) CAD geometry, (b) FE discretization Table 3.1: Cutting conditions in thin-walled rectangular workpiece verification experiments  Roughing Semi-finishing Finishing Workpiece material Al7050 Spindle speed (rev/min) 10500 Axial depth of cut (mm) 20 Radial depth of cut (mm) 1.75 1.1 0.65 Feedrate (mm/tooth/rev) 0.1 Cutting Tool  19.05 mm, 3-fluted solid carbide flat-end mill with 85 mm flute length  Experiments have been conducted on a high speed machining center with the cutting conditions given in Table 3.1. The thin-walled workpiece is clamped on a vise that is rigidly mounted on the machine table (Figure 3.10(a)). Clamping boundary conditions at the base are reflected on the FE model as constrained x, y, and z motions on both sides of the initial workpiece as indicated with triangles in Figure 3.9 (b).  Although the proposed model can update the in-process workpiece dynamics at any time during cutting, overall machining is divided into 27 (k=27) material removal steps (9 segments 48  in each operation) with non-uniform cutting lengths. After a predetermined segment of the workpiece is removed, cutting is stopped for direct FRF measurements at the prespecified 6 points as exemplified for roughing operation in Figure 3.10(b). Semi-finishing and finishing operations also have the same material removal, FRF update, and measurement patterns. FRF measurements are carried out using Polytec CLV-2534 laser vibrometer with Kistler 9722A and PCB 086E80 impact hammers.  Figure 3.10: (a) Experimental setup, (b) measurement locations and removed segments in roughing operation The modal damping constants ( ( ))r wt  are experimentally identified for four different workpiece thickness states, i.e. (i) before machining, (ii) after roughing, (iii) after semi-finishing, and (iv) after finishing for the experimentally observable modes. In order to estimate the variation of modal damping ratios with the thin-walled workpiece thickness ( )wt , they are interpolated for the material removal steps. 49  In-process workpiece FRF update simulations are run for the full and reduced order models within the dominant frequency range of 50-2000 Hz.  For the full order simulations, defining the computational DOF domain for each decoupling 1( )iB CD  and rearranging the DOF sets accordingly decreased the full order model from 46980 DOF to 13821 computational DOF domain for the first decoupling ( )0BCD , and to smaller domains for the following decoupling stations as the material is removed. The initial workpiece FRFs are evaluated with the first 50 modes satisfying the convergence criterion in Eq.(3.24), and the truncated modes are compensated with the first two terms of the series expansion in Eq.(3.22). For the reduced order simulations, master DOF sets are selected based on the first 300 modes for 0B  and 50 modes foriA . The condition number threshold ( )condC  for master DOF selection in Eq.(3.35) is set to 210 . For SEREP reduction of 0B , the first 300 modes (p) are employed in Eq.(3.32), and all of these modes are used in the evaluation of initial FRFs. IRS reduction of iA  is based on 1% convergence in Eq.(3.45). Under these conditions, the initial workpiece having 46980 DOF is reduced to 6389 total DOF that has only 2396 DOF in the computational cutting zone for 0B , i.e., 23960B CD .  Once the undamped dynamics update simulations are completed, the damping is inserted into the model as explained in Section 3.4, and the updated FRFs are evaluated. Although the updated FRFs at 18 different machining stations were monitored in the simulation, comparison results are presented in the most flexible direction (x )W  for 12 points having different characteristics for clarity; Points #1, #2, #4, #5 in roughing, Points #7, #8, #10, #11 in semi-finishing, and Points #13, #14, #16, #17 in finishing operations in Figure 3.11-Figure 3.13. There are two dominant 50  modes (the first bending and the first torsional) within 50-2000 Hz and the dominancy changes between them depending on the tool position along the toolpath. Both frequencies and magnitudes of the modes show strong time dependency, and both modes become more flexible in time. 51   Figure 3.11: Comparison of the measured and predicted direct FRFs along xW  at points #1, #2, #4, and #5 in roughing 52   Figure 3.12: Comparison of the measured and predicted direct FRFs along xW  at points #7, #8, #10, and #11 in semi-finishing 53   Figure 3.13: Comparison of the measured and predicted direct FRFs along xW  at points #13, #14, #16, and #17 in finishing 54  Tracking the bending and torsional mode frequencies in the predicted FRFs is compared with the experimental measurements in Figure 3.14(a-b) and Figure 3.15(a-b). The full order dynamics update model predicts the variation of bending and torsional modes with maximum 2% overestimation at any stage of machining. On the other hand, the reduced order dynamics update model predicts both modes with less than 2% underestimation in roughing, with an average of ~3% overestimation in semi-finishing. In finishing, the average error becomes ~8% and ~6% overestimation for the bending and torsional modes, respectively. The increased deviation between the predicted and measured frequencies of modes in the reduced order model is due to the inevitable missing dynamics in the model order reduction, and also the dynamics modification cannot be exact with the missing dynamics [72]. The deviation propagates from roughing to finishing as the output of each dynamics update ( )iBH  is input for the next one.  Figure 3.14: a) Variation of the first bending mode frequency at all FRF comparison points in the experiment and simulation, b) comparison of the prediction error in full and reduced order dynamics decoupling for the first bending mode  55   Figure 3.15: a) Variation of the first torsional mode frequency at all FRF comparison points in the experiment and simulation, b) comparison of the prediction error in full and reduced order dynamics decoupling for the first torsional mode Moreover, all considered modes in the master DOF selection for 0B  and iA  are assumed to be equally significant; however, a weighting factor could improve the dynamics information carried to the reduced systems while avoiding overemphasis on some modes [83]. This explains why the prediction accuracy of the torsional mode is slightly better than the bending mode. In the master DOF selection of 0B , the first bending mode predominantly selects the DOF towards the tip of the initial workpiece structure due to high modal energy localization (Figure 3.16(a)), and the first torsional mode similarly selects the DOF towards the tip close to the side edges of the initial part (Figure 3.16(b)) since the nodal line of the torsional mode passes through the mid-section. Combining both masters while determining the overall DOF set of the workpiece even includes the DOF close to nodal line of the torsional mode due to the masters selected by the bending mode. This might have caused overemphasis of the torsional mode. 56   Figure 3.16: Schematic illustration of the master DOF selected due to, a) the first bending (shaded region), b) the first torsional (shaded region) modes of the initial workpiece Also, model order reduction methods, in general, do not reduce the interface DOF. In the proposed reduced order dynamics decoupling model, however, the interface DOF are inherently included in the reduction since iA  is a substructure of 0B , and most of the DOF in the cutting region are interface DOF between 0B  and iA . In experimental measurements, a negligibly small amplitude of the torsional mode appears at Points #2 and #8 in roughing and semi-finishing operations, but it is not excited at the corresponding measurement Point #14 in the finishing operation (Figure 3.15(a)) due to hitting on the nodal point of the mode. It is difficult to excite the same point precisely with an impact hammer. Some mismatch in mode magnitudes in Figure 3.11-Figure 3.13 can be attributed to the error in the experimental identification of damping at different thicknesses. In addition, the cut-off criterion used in Eq.(3.28) determines up to what level of the updated undamped dynamics can 57  be used in Eq.(3.30) to evaluate the unknown residues, which affect the mode magnitudes. The criterion in Eq.(3.28) needs to be different for each mode. Due to the thermally resistant nature of the materials used for the thin-walled monolithic aerospace parts in industry, they are cut at low cutting speeds, hence the cutting force excitations occur at low tooth passing frequencies. Therefore, static deflections of the thin-walled workpiece dominate the part accuracy. The full and reduced order dynamics update models are also used to estimate the static stiffness ( )SK  of the thin-walled part in the most flexible xW  direction at all 18 FRF measurement locations, and their comparisons with the experimentally extracted ones are presented in Figure 3.17. SK  decreases continuously at all measurement points, and the workpiece statically becomes more flexible with time as the mass is removed. The full order model can track the variation of SK  with ~20% average error at all stages of machining, while the average error is about 30% for the reduced order model that overshoots the last measurement location of each operation (Points #6, #12, #18). The increased deviation at Points #6, #12, #18 is due to the clustered master DOF close to the tip (Figure 3.16) of the part as the master DOF selection is based on modal energies (Section 3.5.1). However, the last measurement location of each operation is 8-10 times statically stiffer than the other locations and can be excluded in surface error calculations. 58   Figure 3.17: Comparison of the experimentally measured and predicted static stiffness of the rectangular workpiece in xW  at all FRF measurement points Increase in m for the model order reduction and number of the considered modes (p) for evaluation of the initial workpiece FRF have threshold effects on the prediction accuracy. For example, when a certain range is covered (300 modes and 6389 DOF), it retains sufficient dynamics information and modal energy of the dominant modes. Further increase does not improve the prediction accuracy significantly. On the other hand, the computational efficiency is adversely affected by the increase in p and m. The number of modes has a local effect on the computational expense as it is only required to calculate the FRF of the initial workpiece; however, m globally decreases the computational efficiency since the size of both the in-process workpiece and the removed (or fictitious) materials become larger. Although the overall mode prediction accuracy decreased slightly, the reduced dynamics update model has the same order of magnitude error as the full order model throughout machining of the thin-walled workpiece. It 59  uses ~5% of the total DOF of the full system in calculations and provides considerable computation time saving by taking ~8 minutes as compared to ~2.5 hours for the full order model on a workstation having 2.30 GHz central processing unit (CPU).   3.8 Summary A new in-process workpiece dynamics update model based on full and reduced order dynamic substructuring in the frequency domain is presented. Unlike the existing approaches using many FE models of the in-process workpiece, the removed volumes between the discrete dynamics update stations along the toolpath are defined as the substructures of the initial workpiece, and their dynamics are cancelled within the global workpiece FRFs. The effect of the high frequency dynamics on the updated FRFs is considered. Experimentally identified modal damping ratios are inserted into the updated FRFs whose magnitudes at the natural frequencies are highly dependent on the identified damping. The model has been validated with the measured FRFs and mode frequencies in peripheral milling of a thin-walled rectangular workpiece at different stages of cutting.  The proposed model can be used as a physics based tool to predict the varying dominant vibration modes and FRFs of the thin-walled parts in the virtual machining environment without going through the trial and error based cutting tests. The FRFs can be used in predicting the vibrations and chatter free cutting conditions along the toolpath as the material is removed, and the predicted vibration modes can be used for variable pitch cutter design.    60  Chapter 4: Modeling of Varying Workpiece Dynamics in Time Domain 4.1 Overview The prediction accuracy of the frequency domain model presented in the previous chapter depends on the frequency resolution ( )  in the simulation. Its computational efficiency is adversely affected by the number of frequencies in the analysis. These limitations are overcome by introducing the structural dynamics modification in the time domain.  Alternative to the frequency domain model proposed in Chapter 3, a new semi-analytical thin-walled workpiece dynamics update method is developed in the time domain. The model is based on the dynamic substructuring and dynamic characteristics perturbation to update the modal parameters of the workpiece as the material is removed. Similar to the frequency domain approach, the model takes the mass and stiffness matrices of the initial unmachined thin workpiece from its FE model and partitions the workpiece to obtain the removed substructures and their structural matrices. The equations of motion of the fictitious substructures having the negative mass and stiffness of the removed volumes are coupled with those of the in-process workpiece in the time domain to update its varying mass and stiffness matrices along the toolpath (Section 4.2). As opposed to conventionally solving the generalized eigenvalue problem (EVP) repetitively, the modal parameters of the machined workpiece are evaluated only once for the initial thin-walled workpiece. Then, the modal parameters are analytically updated by perturbing them from their initial state until the end of machining at discrete tool positions (Section 4.3.1). The flow chart of the model is detailed in Section 4.3.2. The proposed dynamic parameters update model is experimentally validated in five-axis machining of a twisted fan blade (Section 4.4). The chapter ends with a summary in Section 4.5. 61  4.2 Modeling of Position Dependent Mass and Stiffness of the Thin-Walled Workpiece Continuous material removal in milling of thin aerospace structures results in position dependent mass and stiffness of the machined workpiece, and makes the part more flexible as the cutting advances. Traditionally, the updated mass and stiffness of the newly generated workpiece at any time during machining are required to obtain the associated static and dynamic parameters used in the dimensional form error, vibration, and chatter stability predictions.  Time dependency of the workpiece mass and stiffness for the undamped conditions were expressed in Eq.(3.3) and Eq.(3.6) in Section 3.2. Mass and stiffness contributions of the substructure iA  ( , )i iA AM K  can be removed in the time domain by directly cancelling them within the in-process workpiece 1i-B . Writing Eqs.(3.3) and (3.6) in the block matrix form with DOF partitioning yields, 11 1 111 1 11 1 1 1ii i iii i ii i i iii iii iBB B Baa ac ab aBB B Bca cc cb cB B B Bba bc bb bAA A aaa acAA A cca cc                                           MM M M xM M M 0 xM M M xxM M0xM M11 1 111 1 11 1 1 1ii i iii i ii i i iii iii iBB B Baa ac ab aBB B Bca cc cb cB B B Bba bc bb bAA A aaa acAA A cca cc                                                KK K K xK K K 0 xK K K xxK K0xK K1 11 11 1i ii ii ii ii iB Ba aB Bc cB Bb bA Aa aA Ac c                                      f rf rf rf rf r  (4.1) where the time (t) dependency is dropped for clarity, and the uncoupled block mass and stiffness matrices are called M  and K , respectively. 62  Unlike the frequency domain model, the fictitious coupling is extended to involve both cDOF and aDOF of fictitiously added material iA  as interface DOF to fully decouple the mass and stiffness contribution of iA  from 1i-B , since cDOF and aDOF of iA  are all contained within the global workpiece structure ( )0B  prior to machining. Although aDOF degenerates and becomes a part of cDOF, the notation in Chapter 3 is preserved throughout this chapter for consistency. Structures iA  and 1i-B  can be directly coupled by selecting a unique set of DOF [71] satisfying the interface compatibility through a Boolean transformation ( )T  as,  1111 1iiii iiiBaaaBc a accBbb c cbB BAaaa b bA ccc                                                 xI 0 0x x x0 I 0x 0 0 I x T xI 0 0 x xx0 I 0x  (4.2) where ax  and cx  are the unique aDOF and cDOF vectors, respectively. The interface equilibrium is automatically satisfied by premultiplying the coupling force vector in Eq.(4.1) by TT  as,  1111iiiiiiBaBcaa aa aBcc cc cbBAbb a bAc                              rrI 0 0 I 0 00 I 0 0 I r 00 0 I 0 0 rrr  (4.3) where 1i iB Aa a  r r , 1i iB Ac c  r r  as the interface equilibrium is applied on both cDOF and aDOF, and 1iBb b r 0  since the coupling forces do not act on the internal bDOF of the 63  remaining workpiece structure ( )iB . Then, substituting Eq.(4.2) into Eq.(4.1) and premultiplying by TT  yield the coupled  (assembled) set of equations of motion for iB  as,  1111 1iiii i iiBaBca aBT T Tc c bB B Ab b aAc                                      ffx xT MT x T KT x T fx x ff  (4.4) where the updated mass and stiffness matrices of iB  are directly given as,  iiB TB TM T M TK T K T  (4.5) One can conventionally evaluate the updated modes and mode shapes of the in-process workpiece by recursively solving the generalized eigenvalue problem using the updated mass and stiffness matrices given in Eq.(4.5) for each material removal step along the toolpath. This is time consuming for large workpiece structures with long toolpaths and is prone to have ill-conditioned matrices due to the decoupling (Eq.(4.5)). Instead, the varying modal characteristics of the remaining workpiece ( )iB  are obtained by analytically perturbing the mode shapes of the in-process workpiece from its initial state 0( )B  to the final state ( )kB  as the material is removed, which is detailed in the next section. 4.3 Modeling of Varying Dynamic Characteristics of the In-Process Workpiece The dominant vibration modes and mode shapes in the cutting region vary with time both in longitudinal and transverse (tool axis) directions of the thin-walled parts due to the change in the mass and stiffness caused by the material removal. The proposed model accepts the removed substructures as perturbations to the dynamics of the machined part and analytically links the 64  eigensolution of the workpiece at cutting step i  ( )iB  to that of at cutting step 1i   1( )iB   to obtain the instantaneous modal parameters of the in-process workpiece.  4.3.1 Evaluation of the Updated Modes and Mode Shapes The cutting strategy determines how the dynamic characteristics change as illustrated for the first bending mode of the general case in Figure 4.1, and the variation should be known at the toolpath planning stage. The generalized eigenvalue problem for the in-process workpiece at the cutting step 1i   1( )iB   is,  1 1 11 1i i ii im m mB B BB Bj j j   K M    (4.6) where 1imBj   and 1imBj  are the thmj  eigenvalue (square of the mode frequency) and the mass normalized eigenvector of 1iB  . Eq.(4.6) is subject to the following modal orthogonality conditions,    1 1 111 11i i iim mm m mi iim mm mTB B BBj kj k jTB BBj kj k    KM    (4.7) where m mj k  is Kronecker Delta, and it is given as,  1 ,0 , otherwisem mm mj kj k      (4.8) 65   Figure 4.1: Varying mode frequencies and dynamic displacements in the cutting region The updated eigenvalue problem for the machined workpiece at cutting step i  ( )iB  becomes,  i i ii im m mB B BB Bj j jK M    (4.9) which is also subjected to Eq.(4.7). The change in mass and stiffness of 1iB   due to the material removal (Eq.(3.1)) can be mathematically rewritten using the mass and stiffness of the removed material as perturbation terms,  11i i ii i iB B AB B A    K Κ KM M M  (4.10) where the updated iBM  and iBK  are same as the ones in Eq.(4.5),   is the perturbation parameter showing the order of the terms, iAK  and iAM  are the change in stiffness and mass matrices of the in-process workpiece due to the non-uniform material removal, and they are equal to the structural matrices of the removed material ( )iA . Instead of solving Eq.(4.9) repeatedly for each material removal step, the eigenvalues and eigenvectors of the in-process workpiece structure ( )iB  can be expressed as a series expansion [84, 85], 66   112 31, 2, 3, q ,2 31, 2, 3, q ,......pi i i i i im m m m m p mpi i i i i im m m m m p mqB B A A A Aj j j j j jqB B A A A Aj j j j j j                                     (4.11) where q ,ip mAj  and q ,ip mAj  ( 1,2,3...)pq   are the thpq  order perturbation terms for the thmj  eigenvalue and eigenvector of the in-process workpiece 1iB   due to removal of the substructure iA  (or addition of iA ). Truncating Eq.(4.11) after the second order perturbation term [86], and inserting it along with Eq.(4.10) into Eq.(4.9) give the updated (modified) eigenvalue problem as,     1 11 112 21, 2, 1, 2,21, 2,i i i i i ii i i im m m m m mi i im m mB A A B A AB A B Aj j j j j jB A Aj j j                                 Κ K M M     (4.12) Expanding Eq.(4.12) by neglecting the higher order terms yields, 1 11 11 1 11 11 11 121, 2,2 21, 2, 1,2 21, 1, 1, 1,i i i ii im m m mi i i i i ii i i im m m m m mi i i i i ii i i im m m m m mB A B AB Bj j j jA A B B B AB B A Aj j j j j jB A A B A AA A B Bj j j j j j                                         M MΚ Κ M MΚ Κ M M       1 112 21, 2,i i i ii im m m mA B A BA Bj j j j                M M   (4.13) Grouping the same order of   on both sides of Eq.(4.13) leads to,    1 1 1 1 11 1 11, 1, 1,i i i i i i i ii i i i im m m m m m m mA B B A B B A BB A B A Bj j j j j j j j                Κ Κ M M M       (4.14) and,   1 1111 1112, 1,2 22, 1, 1, 1, 1,2,i i i ii im m m mi i i i i ii i i im m m m m mi iim mB A B AB Aj j j jA A A A A BB A B Aj j j j j jA BBj j                         M MΚ Κ M MM      (4.15) 67  The unperturbed eigenvectors of 1iB   are used as basis for the first and second order mode shape perturbation terms since they define an orthogonal basis [87],  1 11, 2,1 1, , (1,....., )O Oi i i im m m mm m m mm mN NA B A Bj k j k m Oj k j kk ka b j N             (4.16) where m mj ka  and m mj kb  are the contribution coefficients of thmk  mode shape of the workpiece structure 1iB   to the first and second order mode shape perturbation terms of its thmj  mode. Substitution of 1,imAj  in Eq.(4.16) into Eq.(4.14), and premultiplying, respectively, by  1imTBj  and    1imTBk   m mk j  yield the following coefficient equations,     1 111 1 1 1 1 1111111,Oi iim mm mmOi i i i i ii i im mm m m m m mmi iim mNB BBj kj kkNT TB B B B B BB A Aj kj k j j j jkA BBj jaa                               MΚ Κ MM (4.17) where    1 1 1 1 1 1 11 11 1,O Oi i i i i i ii im m m m m m m mm m m m m m mm mN NT TB B B B B B BB Bj k j j j k j jj k j j j k jk ka a a a            Κ M     due to the orthogonality conditions (Eq.(4.7)) and,     1 111 1 1 1 1 1111111,Oi iim mm mmOi i i i i ii i im mm m m m m mmi iim mNB BBj kj kkNT TB B B B B BB A Aj kk k j k j jkA BBj jaa                           MΚ Κ MM     (4.18) 68  where  1 1 111Oi i iim m m mm m mmNTB B BBj k j kk k kka a   Κ  ,   1 1 1 111Oi i i iim m m mm m m mmNTB B B BBj k j kk j k jka a    M   due to Eq.(4.7). The first order perturbation term of the thmj  eigenvalue is evaluated from Eq.(4.17) as,     1 1 11, i i i ii im m m mTA B B BA Aj j j j       Κ M    (4.19) Similarly, the coefficient m mj ka  in Eq.(4.16) is calculated from Eq.(4.18) as,     1 1 11 1,i i ii im m mm m i im mTB B BA Ak j jj k m mB Bj ka k j       Κ M   (4.20) Analogously, inserting 2,imAj  given in Eq.(4.16) into Eq.(4.15) and premultiplying, respectively, by  1imTBj  and  1imTBk   m mk j  yield the following set of coefficient equations,    1 111 1 1 11 11 1112 21, 1, 1, 1,11, 2,Oi iim mm mmOi i i i i i i ii i i im mm m m m m m m mmi i i ii im m m mNB BBj kj kkNT TB B A B B A A AB A A Bj kj k j j j j j jkA B A BA Bj j j jbb                                     MΚ Κ M MM M        (4.21) and, 69     1 111 1 1 11 11 1112 21, 1, 1, 1,11, 2,Oi iim mm mmOi i i i i i i ii i i im mm m m m m m m mmi i i ii im m m mNB BBj kj kkNT TB B A B B A A AB A A Bj kk k j k j j j jkA B A BA Bj j j jbb                                     MΚ Κ M MM M             (4.22) Due to the modal orthogonality (Eq.(4.7)), the second order perturbation term of the thmj  eigenvalue is evaluated from Eq.(4.21) as,         1 1 1 112, 1, 1, 1,i i i i i i i ii i i im m m m m m m mT TA B B A A B A BA A B Aj j j j j j j j               Κ M M M     (4.23) Similarly, the contribution coefficient m mj kb  in Eq.(4.16) is calculated from Eq.(4.22) as,       1 1 1 111 11, 1, 1,,i i i i i i ii i i im m m m m m mm m i im mT TB B A A B A BA A B Ak j j j k j jj k m mB Bj kb k j                  Κ M M M  (4.24) In order to obtain the contribution coefficients of thmj  mode shape to its own first and second order perturbation terms (m mj ja and m mj jb ), the mass orthogonality condition in Eq.(4.7) is rewritten for the in-process workpiece structure at time step i  ( )iB  as,    1i iim mTB BBj jM    (4.25) Inserting the second order perturbed mode shape imBj  from Eq.(4.11) and the perturbed mass from Eq.(4.10) into Eq.(4.25) gives,     1 112 21, 2, 1, 2, 1i i i i i ii im m m m m mTB A A B A AB Aj j j j j j                 M M        (4.26) Expanding Eq.(4.26) to its terms and neglecting the higher order terms lead to, 70            1 1 11 11 1 111 11 11, 2, 1,21, 1, 1,1, 2,i i i i i ii i im m m m m mi i i i i ii i im m m m m mi i i ii im m m mT T TB A B A B AB B Aj j j j j jT T TB B A A A BA B Aj j j j j jT TA B A BB Bj j j j                                        M M MM M MM M             0   (4.27) Substituting the mode shape perturbation terms 1, 2,( , )i im mA Aj j    given in Eq.(4.16) into Eq.(4.27) and setting the coefficients of   and 2  to zero yield,      1 1 1 1 1 11 11 11 1012O Oi i i i i ii i im m m mm m m m m mm mi iim m m mTN NT TB B B B B BB A Bj k j kj k j j k jk kTB BAj j j ja aa                         M M MM (4.28) and,           1 1 11 11 1 111 11, 1, 1,11,11, 1,012Oi i i i i ii i im mm m m m m mmOi i i ii im mm m m mmi i ii im m m m mNT T TB B B A A AB A Bj kj k j j j jkTNTA B B BA Bj kj j k jkT TB A AA Bj j j j jbbb                                 M M MM MM M             11, 1,i i iim m mTA A BAj j j      M    (4.29) Then, the first 1,( )imAj  and second 2,( )imAj  order mode shape perturbation terms for the thmj   mode shape in Eq.(4.11) become,      1 1 11 1 11 11,112i i ii iOm m mi i i iim m m mi im m mm mTB B BA ANT k j jA B B BAj j j kB Bk j kk j            Κ MM       (4.30) and, 71            1 111 1 1 1111 11, 1, 1, 1,2,1, 1, 1,112i i i i i ii i im m m m m mii i i i i i ii i i iOmm m m m m m m imi im m mm mT T TB A A A A BA B Aj j j j j jA T TB B A A B A BA A B ANjk j j j k j j BkB Bk j kk j                             M M MΚ M M M           (4.31) Inserting Eqs.(4.30)-(4.31) into Eq.(4.11) evaluates thmj  mode shape of the in-process workpiece structure at machining step i  ( )iB  as a function of the modes and mode shapes from the previous cutting step, the mass and stiffness of the removed material in the current cutting step, and the updated mass of the in-process workpiece in the previous cutting step. The recursive implementation of the Eqs.(4.1)-(4.5),(4.11),(4.30)-(4.31) from the initial workpiece structure 0( )B  until the finished workpiece ( )kB  leads to the updated varying mode shapes of the in-process workpiece without solving the eigenvalue problem repetitively as the cutting advances.  The updated mode frequencies are evaluated using Rayleigh Quotient [88] with the updated mass and stiffness matrices obtained using Eq.(4.5) and the perturbed mode shapes ( )imBj  obtained from Eqs.(4.11), (4.30), (4.31) as,  2,( )( )( )i iim mi imm i iim mB BBTj jB Bn jj B BBTj j  KM    (4.32) where ,( )imBn j  is the updated thmj  mode frequency of the remaining workpiece structure iB . Eqs.(4.5), (4.11), (4.30)-(4.32) are applied as the tool removes material along the toolpath to update the dominant modes and mode shapes of the in-process workpiece without explicitly solving the generalized eigenvalue problem after each material removal step. Eq.(4.1) is the 72  general form to couple the in-process workpiece 1( )iB   with the fictitious substructures ( )iA , and both structures’ mass and stiffness matrices can be approximated with their dynamically equivalent reduced order counterparts 1iBredM  and 1iBredK ; and iAredM  and iAredK . The developed time domain in-process workpiece dynamic characteristics update model is extended to the dynamically reduced order model on a twisted blade example in Section 4.4 to have a fair comparison with the frequency domain model developed in Chapter 3.  4.3.2 Simulation Algorithm The overall calculation scheme of the proposed model is shown in Figure 4.2, which starts with discretizing the initial global workpiece 0( )B  to obtain its structural matrices. After planning the update locations along the toolpath ( 1... )i k , the structural matrices of the fictitiously added substructures (removed material) are obtained from 0B ’s mesh. Mass and stiffness contributions of the removed material are directly decoupled (Eq.(4.5)) from those of the in-process workpiece 1( )iB  . After solving the first few dominant modes and mode shapes (Eq.(4.6)) only for 0B  in the decoupling module, perturbations of the mode shapes starting from 0B  for all the dynamics update stations along the toolpath are carried out using Eqs. (4.11), (4.30)-(4.32). Without loss of generality, the model can be implemented with dynamically reduced order structures by selecting the master DOF using Eqs.(3.32)-(3.35) as outlined in Section 3.5.1, and the model order reduction can be carried out for iA  and 0B  using Eqs.(3.38), (3.42), and (3.44) as presented in Section 3.5.2.  73  The obtained modal parameters  ,, ( )i immB Bn jj  are stored and used in frequency response function (FRF) evaluation, which can be employed in chatter stability and workpiece vibration calculations.  Figure 4.2: Flow chart of the proposed in-process workpiece modal parameters update model 4.4 Experimental Verification: Thin-Walled Twisted Blade Cutting Test The developed in-process workpiece modal parameters update model has been verified in five-axis point (ball-end) milling of a thin-walled twisted fan blade made of Al7050. The 74  experimentally measured FRFs are obtained through impact hammer tests by interrupting the cutting at different tool positions along the toolpath, and they are compared to the FRFs evaluated using the updated mode frequencies and mode shapes of the in-process blade in Eq.(3.25). The full order time domain model predictions have been compared against those of the repeatedly solved eigenvalue problem (EVP) along the path. Also, the time domain model is extended to the reduced order dynamics perturbation using the methods detailed in Section 3.5, and its results are compared with the reduced order frequency domain dynamics update model of Chapter 3. An initial rectangular block is first machined to a rough thin blade geometry having a uniform thickness of 10 mm, followed by the machining sequence applied to both sides; Roughing-1 (R1), Roughing-2 (R2), Semi-Finishing-1 (SF1), Semi-Finishing-2 (SF2), Finishing-1 (F1), and Finishing-2 (F2) as illustrated in Figure 4.3(a). Each operation is divided into four segments (k=24 in total), and the last segment includes milling of the blade’s root. Thickness of the finished blade is non-uniform, and it varies between 4.5 and 4.8 mm. The maximum reduction in the thickness in R1 and R2 is 1.5 mm, while it is 0.75 mm in SF1 and SF2, and 0.5 mm in F1 and F2. FE model of the initial blade geometry is generated with 10-node tetrahedral structural solid elements as shown in Figure 4.3(b) by assigning the same material properties for Al7050 as Section 3.7. The mesh convergence for 0B  is obtained with 118,002 total DOF for the full order FE model considering different radial depths of cut in all operations as given in the close-up view of Figure 4.3(b).  75   Figure 4.3: a) Initial workpiece geometry, machining sequence, and four removed segments in each operation, b) FE mesh of the initial thin-walled twisted blade Five-axis ball-end milling tests have been conducted on a machining center with the cutting conditions given in Table 4.1. The workpiece is clamped on a vise that is bolted on the machine table as shown in Figure 4.4(a). Unlike the peripheral cutting, a small amount of material is removed in each ball-end milling pass. Hence, the mode frequencies and the mode shapes are updated after the last pass of each segment, and the updated FRFs are evaluated at 3 different tool positions in each segment using Eq.(3.25) for the frequency range of 100-2500 Hz. Due to the limited access, direct FRFs of the workpiece in the xW  direction are measured at the first 9 tool positions in each operation (54 update positions in total) as shown for the operation R1 in 76  Figure 4.4(b). The remaining operations (R2, SF1, SF2, F1, F2) have the same measurement pattern.  Table 4.1: Cutting conditions in five-axis twisted blade experiments  Pre-machining 1 Pre-machining 2 R1 R2 SF1 SF2 F1 F2 Spindle speed (rev/min) 4000 10000 5000 5000 1000 1000 1000 1000 Maximum radial depth of cut (mm) max. 50% of diameter max. 50% of diameter 1.5  1.5  0.75  0.75  0.5  0.5  Feedrate (mm/tooth/rev) 0.2 0.1 Cutting tool  50 mm, indexable face mill  12.7 mm, 4-fluted tapered ball-end mill with 3 0  taper angle and 62 mm flute length Workpiece material Al7050   Figure 4.4: a) Experimental setup, b) measurement (comparison) locations (points 1-9) in R1 Although variations of the FRFs at all 54 points are tracked in the simulation, their comparisons against the experimentally measured FRFs are presented only for the Points #1, #2, 77  #6, #7 in the first operation (R1); the Points #19, #20, #24, #25 in the third operation (SF1); and the Points #37, #38, #42, #43 in the fifth operation (F1) in Figure 4.5-Figure 4.7 for clarity.  Figure 4.5: Comparison of the measured and predicted direct FRFs in xW  direction using the full order time domain (TD) model and repetitively solved eigenvalue problem (EVP) at points #1, #2, #6, and #7 in operation R1 78   Figure 4.6: Comparison of the measured and predicted direct FRFs in xW  direction using the full order time domain (TD) model and repetitively solved eigenvalue problem (EVP) at points #19, #20, #24, and #25 in operation SF1 79   Figure 4.7: Comparison of the measured and predicted direct FRFs in xW  direction using the full order time domain (TD) model and repetitively solved eigenvalue problem (EVP) at points #37, #38, #42, and #43 in operation F1 80  The first bending and the first torsional modes dominate the frequency range of interest, and their frequencies and magnitudes show time dependency. Both the time domain (TD) dynamics update model and repetitive EVP solution can predict the variation of the FRFs along the toolpath well with less than 20% prediction error in the dynamic stiffness of the modes. The discrepancy between the measured and predicted FRF magnitudes at some points along the path is mainly caused by the experimentally identified modal damping ratios ( )r  of the observable (first three) modes, which are identified from the measured workpiece FRFs for different thin-walled blade states (before machining, after operation R2, after operation SF2, and after operation F2) as explained in the previous chapter. They are used both in the perturbation and EVP based calculations by interpolating them for the intermediate FRF comparison points. Also, the blade has a root blend section where thickness of the removed material rapidly decreases from the radial depth of cut (point A in Figure 4.4(b)) to zero at the hub (point B in Figure 4.4(b)). Since it is hard to accurately represent this region with the FE mesh, and it increases the model size significantly, the material removal in this region is not considered in the simulations. Though the clamping has a finite stiffness, the DOF on the clamping interface are fixed in the initial FE model, which contributes to the prediction error as well.  Variations of the first five mode frequencies and the mode shapes of the workpiece are also predicted using the reduced order time domain (TD) and the reduced order frequency domain (FD) models having similar prediction errors as the full order models. The same reduction settings of Section 3.7 are used. The initial full order workpiece is reduced to 33,962 DOF, which has a total of 9657 interface DOF (~8% of the full system) in the cutting zone along the toolpath. Frequency domain FRF update simulations are run with the frequency increment ( )  81  of 2.5 Hz for the same frequency range. Comparisons of the reduced order simulation results with the experimentally measured FRFs are given in Figure 4.8-Figure 4.10.  Figure 4.8: Comparison of the measured and predicted direct FRFs in xW  direction using the reduced order time domain (TD) and the reduced order frequency domain (FD) models at points #1, #2, #6, and #7 in operation R1 82   Figure 4.9: Comparison of the measured and predicted direct FRFs in xW  direction using the reduced order time domain (TD) and the reduced order frequency domain (FD) models at points #19, #20, #24, and #25 in operation SF1 83   Figure 4.10: Comparison of the measured and predicted direct FRFs in xW  direction using the reduced order time domain (TD) and the reduced order frequency domain (FD) models at points #37, #38, #42, and #43 in operation F1 84  In addition to the error introduced by the identified modal damping, model order reduction of the structures, and the idealized clamping conditions, the magnitude discrepancy in Figure 4.8-Figure 4.10 can be partially attributed to the 2.5  Hz used in the frequency domain model as it limits how well the sharply erupting and lightly damped modes are captured. Smaller   provides more accurate predictions with longer computation times. This explains why the reduced order time domain model provides more accurate mode amplitudes than its frequency domain counterpart.  The reduced order time domain dynamics update model predicts the dynamic stiffness of both modes with ~26% average error throughout machining. On the other hand, the frequency domain model predicts the dynamic stiffness of both modes with ~32% average error excluding the tool positions #37-#38, where the predicted FRF magnitude of the bending mode has the highest mismatch with the experimentally measured one due to improper   selection for that mode frequency, which is not known a priori. This issue can be circumvented by the simulation with a coarse   to bracket the mode frequencies first, and then running the simulation with a finer   around the identified modes from the initial simulation. The better prediction of the time domain model could be attributed to allowance of fine frequency increment (0.1 Hz) in Eq.(3.25) as it does not affect the computational efficiency of the model.   The perturbed first bending and torsional mode frequencies predicted using the full and reduced order time and frequency domain models are compared with the measurements in each material removal segment as shown in Figure 4.11 and Figure 4.12, respectively.   85   Figure 4.11: a) Variation of the predicted and measured first bending mode frequency with material removal, b) comparison of the absolute prediction error in the full and reduced order time domain (TD) and frequency domain (FD) models, and the full order repetitive eigenvalue (EVP) solution  All the developed models keep the same order of magnitude error for the dominant bending and torsional modes of the thin-walled twisted blade throughout machining. The average prediction accuracy and the computational performance of the proposed models are compared with the repetitive eigenvalue solution (EVP) in Table 4.2. The full order time domain perturbation model utilizes the sparsity of the physical system matrices, and its considerable computational advantage over the full order EVP model is due to the ill-conditioned matrices after each material removal step in Eq.(4.5). The model order reduction apparently does not provide any computational savings in the time domain model due to the dense matrices. In order 86  to keep the same accuracy level as the other developed models, the reduced order frequency domain model uses ~8% of the full DOF set and takes the largest amount of time. Its computational burden increases with the number of frequencies in the analysis. Also, behaviors of the reduced order models fluctuate depending on the position along the toolpath since the DOF located towards the free end of the blade are predominantly selected as master DOF.  Table 4.2: Average prediction error and computation time comparisons of the developed thin-walled workpiece dynamics update models for the twisted blade case (EVP: Eigenvalue problem, TD: Time domain, FD: Frequency domain)  Average Absolute Prediction Error for Mode 1 Average Absolute Prediction Error for Mode 2 Computation Time Full Order EVP 3.7% 2.7% ~30 minutes Full Order TD 4.9% 3.1% ~2 minutes Reduced Order TD 5% 3.6% ~30 minutes Reduced Order FD 4% 4.9% ~120 minutes  Different error characteristics in different operations might have resulted from using the same perturbation order for all operations. This could be circumvented either by using adaptive (different) perturbation orders for different operations, or by using a separate initial FE model at the beginning of each operation type. Also, the number of tracked modes in the time domain dynamics update model is constant and cannot be changed in the perturbation calculations as the material is removed. Perturbing more mode shapes could possibly improve the time domain predictions as the number of perturbed modes affects the accuracy of the perturbation terms in Eq.(4.11). Sensitivity analysis on the effect of number of perturbed modes (Figure 4.13) showed that further increase in the number of modes marginally improves the mode prediction accuracy of the first two dominant modes after the perturbation converges.   87   Figure 4.12: a) Variation of the predicted and measured first torsional mode frequency with material removal, b) comparison of the absolute prediction error in the full and reduced order time domain (TD) and frequency domain (FD) models, and the full order repetitive eigenvalue (EVP) solution  Figure 4.13: Effect of number of the perturbed modes on the performance and prediction error in the reduced order time domain (TD) model 88  The proposed perturbation based time domain dynamic characteristics update model eliminates the need for repeatedly solving the computationally intense generalized eigenvalue problem for the full and the reduced order systems along the toolpath. Its accuracy slightly improves with the increase in the number of perturbed mode shapes at the expense of increased computation time. When the thin-walled workpiece has closely spaced mode frequencies, the presented method might become inaccurate due to the perturbation terms given in Eqs.(4.30)-(4.31). Thus, it might require more perturbation terms in Eq.(4.11), or simulation of the material decoupling with smaller size iA . In order to improve the convergence of the truncated expansion in Eq.(4.11), the algorithm given in [85, 89] is used as shown in Appendix A. Similarly, when the stiffness and mass contributions of the removed material become significant compared to the overall mass and stiffness of the in-process workpiece, the proposed method might need to use higher order perturbation terms. Since the blade geometries generally have distinct mode frequencies and relatively small material removals in ball-end milling passes, using the second order perturbation method for the complex blade geometries is found to be sufficient in this study. 4.5 Summary A new and efficient instantaneous modal parameter update model is presented using the substructure decoupling and mode shape perturbation for the in-process workpiece in place of currently used full order FE models for each cutting step. The updated mass and stiffness of the machined part are obtained by structurally modifying the workpiece with the fictitious negative of the removed material at different tool positions along the toolpath. The change in the distinct modes and mode shapes induced by the material removal is evaluated by perturbing the mode shapes without repeating the eigensolution for the newly generated workpiece at discrete stations 89  along the toolpath. The proposed model has been verified with the cutting tests and measured FRFs in all operations. Evidently, the proposed time domain model is computationally more efficient than both the previously developed frequency domain model (Chapter 3) and the repeated eigensolution method. The proposed model can be used for virtual simulation and optimization of five-axis milling of fan, gas turbine, and compressor blades by predicting their time varying dynamic parameters at the toolpath planning phase.              90  Chapter 5: Analytical Modeling of Process Damping in Machining 5.1 Overview High strength and difficult to cut aerospace materials are machined at low cutting speeds due to tool wear concerns. Although the machining process induced damping due to indentation of the cutting edge into the undulated cut surface has a strong effect on the process stability, it has not been fully modeled yet. As opposed to the empirical methods reviewed in Section 2.3, an analytical model is derived to estimate the process damping without the cumbersome experimental calibrations [46, 48-55] and the time prohibitive numerical identification methods [56].  In this chapter, a comprehensive model of process damping has been developed using two-dimensional (2D) representation of the cutting edge element. The indentation of the rigid cutting edge into the work material and the associated elastic-plastic contact pressure are analytically modeled by assuming the workpiece as a half-space with plane strain deformation. The indentation geometry and boundaries are evaluated using the edge geometry and the surface wave form in Section 5.3. In Section 5.4, the elastic contact pressure is obtained for the rounded and straight sections of the cutting edge at discrete positions along the wavelength of the surface undulations. The overall contact pressure is obtained as a function of the cutting edge geometry (edge radius, clearance angle, equivalent wear length), vibration frequency and amplitude, and the material properties (modulus of elasticity, Poisson’s ratio, yield strength) of the work material using Tabor’s equation as the plastic limit. The resulting specific indentation force is evaluated by integrating the overall contact pressure along the contact length in Section 5.5. Then, the non-linearity of the penetration force is represented by an equivalent specific viscous damping, which is used in the chatter stability calculations in the next chapter. The newly 91  proposed analytical process damping model is compared against the experimentally and numerically identified process damping force coefficients for different materials in Section 5.6. The chapter is completed with a summary in Section 5.7.       5.2 Process Damping Phenomenon The 2D cutting edge element is obtained by intersecting the flute with a plane in its normal direction. The cutting edge geometry is defined by the normal rake angle ( )n , normal clearance (flank) angle ( )n , and the edge (hone) radius ( )eR  as shown in Figure 5.1(a). In the absence of the structural vibrations, a constant chip thickness is generated with material shearing as illustrated in Figure 5.1(b). Due to the curved (honed) tip of the cutting edge element, some material is sheared away from the workpiece and forms the chip, while the rest moves beneath the flank face and is plowed by the cutting edge (Figure 5.1(b)). The angle between the radial direction and the line connecting the material separation point (SP) and the hone center is named as the separation angle ( ) . Since determination of   is still a fundamental research question and is not within the main scope of the current work, it is directly obtained from the available literature in this thesis.           Figure 5.1: (a) Two-dimensional cutting edge geometry, (b) static chip generation and material separation point (SP)  92  In the presence of relative structural vibrations due to the flexible workpiece and cutting tool- tool holder-spindle assembly, the cutting edge follows a wavy trajectory as illustrated in Figure 5.2. While the tangential (horizontal) motion of the cutting edge generates the chip, its radial (vertical) motion causes interference between the freshly generated workpiece surface and the cutting edge element. When the tool moves downward the rigid cutting edge extrudes the work material, and the resulting plowing forces are shown both in experimental [37, 43, 44, 46, 54] and numerical (FE analysis) [56] studies to be responsible for the machining induced damping called process damping.        Figure 5.2: Machining process damping due to the cutting edge indentation and the associated forces The cutting edge element penetration changes the interference between the edge and the work material depending on the position along the wave. Throughout the thesis, the plowing in Figure 5.1(b) and the change in the interference (Figure 5.2) will be named as static and dynamic indentations, respectively. The dynamic indentation causes process damping, which improves the chatter stability of the machining operations especially at low cutting speeds. An analytical process damping force prediction model eliminating the need for cumbersome empirical identifications is presented in the following sections.   93  5.3 Indentation Geometry Since the dominant source of  process damping is shown to be the movement of the cutting edge into the workpiece [43, 44, 46, 47, 50-52, 56], half the wavelength of the undulated surface is discretized, and the material separation point (SP) is moved along the wave due to the tool trajectory as shown in Figure 5.3. The overall contact pressure ( ( ))dp x  is evaluated at each discrete point for the cutting edge element and the associated indentation geometry. Indentation boundaries are governed by the wavelength of the undulated workpiece surface ( )W , the normal clearance angle ( )n , the material separation angle ( ) , the edge radius ( )eR , and the positions of the points EO , O, F, J, SP when the effect of material pile-up or sink-in is not considered. As the objective is to predict the chatter stability of the machining operations, W  can be expressed as [47],  W ccV   (5.1) where cV  is the cutting speed (mm/s) and c  is the chatter frequency (Hz), which is close to the dominant natural frequency of the flexible milling system.   94   Figure 5.3: Indentation geometry, surface wave form, and discrete process damping calculation points Profiles of the wavy workpiece surface  Wavey , the rounded  J SPy   and the flank  EO Jy   sections of the cutting edge element are expressed as,  22( ) sin( )( ) , ,2( ) tan( ) , ,EWWave0 WJ SP xx xeO J En x xxxy x Ax Jy x x J SPRy x x n x O J              (5.2) where 0A  is the amplitude of the undulated workpiece surface; Wxn  shows the x-coordinate of the intersection of the flank face with Wx  axis; xJ , xSP , and ExO  are x-coordinates of the associated points in Figure 5.3. The equivalent wear length ( )wL  in Figure 5.3 can be calculated using the cutting edge geometry and the material separation angle as [54],   sin sin cot (cos cos )w e n n nL R           (5.3) 95  The positions of points F and the cutting edge roundness center (O) are analytically expressed for the known discrete position of the separation point (SP) as,  02 3, sin , ,4 4W Wxx x w y xWSPF SP L F A SP              (5.4)  02sin , sin cosxx x e y eWSPO SP R O A R         (5.5) where the subindices x and y denote x and y coordinates of the points. Position of the point J, where the clearance face meets the honed section of the cutting edge element, is analytically calculated as,  sin , cosx x e n y y e nJ O R J O R       (5.6) Equation of the line connecting the points EO  and F can be expressed as,    ( ) tanEO Fy n xg x F x F     (5.7) The difference between the Eq.(5.7) and the wave equation in Eq.(5.2) is,   02( ) ( ) ( ) sin tanEWave O Fy n xWxf x y x g x A F x F          (5.8) Then, the position of EO  is evaluated using Newton’s Method [90] at the cutting edge position nx x  as,  '1 0'( ) 2( ) 2( ) , ( ) cos tan( )nn nn n n nW Wx xnf x xdf xf x x f x Adxf x         (5.9) where 1nx   iteratively approaches to ExO  when the predefined solution tolerance is satisfied, 1n n Tolx x    , and the corresponding EyO  is obtained from the wave form given in Eq.(5.2). 96  Eqs.(5.1)-(5.9) fully define the indentation geometry that is used to calculate the contact pressure distribution at the interface of the cutting edge element and wavy workpiece surface as detailed in next section. 5.4 Modeling of Contact Pressure under the Flank Face  The process induced interaction between the rigid cutting edge and the elastic-plastic workpiece is modeled as quasi-static indentation assuming plane strain conditions with line loading. Depending on the position of the cutting edge element on the surface undulation, the clearance (flank) face and the rounded section indent the workpiece as illustrated in Figure 5.3.  The relative normal displacement in the indentation region  (x)y  is written as,  ( ) ( ) , ,( )( ) ( ) ( ) , ,EEWave O J Ex xWave O J J SPx xy x y x x O Jy xy x y x y x x J SP               (5.10) Substitution of the profiles given by Eq.(5.2) into Eq.(5.10) yields,    22sin tan , ,( )( )2sin tan , ,2WWE0 n x xW xx0 n x xW xexA x n x O Jy xx JxA x n x J SPR                         (5.11)  The normal surface displacement gradient is obtained from Eq.(5.11) as,     2 2cos tan , ,( )( )2 2cos tan , ,E0 n x xW Wx0 n x xW WexA x O Jy xx x JxA x J SPR                             (5.12) Contact pressure at the interface between the cutting edge element and the work material (Figure 5.4) needs to be evaluated to calculate the process damping forces ( )PDF  shown in Figure 5.2. Calculation of the contact pressure is simplified by neglecting the effect of friction on 97  the interface pressure. This effect is shown to be ~5% for the typical indenter geometries [91, 92], and also the contact area is practically found to be unchanged due to the friction [92].   Figure 5.4: General contact pressure distribution at the interface between the cutting edge element and workpiece   Hence the relation between the normal surface displacement gradient (Eq.(5.12)) and the elastic contact pressure ( )ep  at the interface is written as [93],  2( )d( ( ))4(1 )xExSPeOpE y xx x       (5.13) where E and   are modulus of elasticity and Poisson’s ratio of the work material, respectively. The elastic contact pressure is obtained by inverting Eq.(5.13) [4, 94],         2(x) d( )4 1 ( )xExSPEe x xEO x xyEp x x O SP xx O SP x             (5.14) where [ , ]Ex xx O SP . Inserting Eq.(5.12) into Eq.(5.14) gives the elastic contact pressure distribution as, 98              2dtan( )d( )( )1 d4 12cos d2( )xExxxxxxExSPnEO x xSPxEe J x xESPe x xEe J x xSP W0W EO x xO SP xx JR O SP xEp x x O SP xR O SPAO SP x                                           (5.15) The first three integrals on the right hand side of Eq.(5.15) are the same as the governing equations for the indentation of an inclined punch with a rounded tip into a flat surface, and their solutions are given in [94]. The last integral equation in Eq.(5.15) represents the contribution of the wavy workpiece surface to the elastic contact pressure. Then, Eq.(5.15) can be rewritten as,  ( ) ( ) ( )P We e ep x p x p x    (5.16) where Pep  and Wep  stand for the elastic contact pressure contribution of the punch and wave, respectively. ( )Pep x  can be evaluated using the contact geometry in Figure 5.3 as [4, 94],          1 12 1 112 2 2111 111( ) 1 tan ln4 14 1 1 12 21 1( ) tan sin , tan sin2 2Ex xPeeE Ex x x x xE Ex x x xE O SPp xR yx O SP J O SPxSP O SP O                                                           (5.17) where ( )x  and 1  are temporary variables. The elastic contact pressure due to the wavy workpiece surface ( )Wep  can be explicitly written from Eq.(5.15) as, 99      22cos d( ) (x)( x)2 1xExSP WW 0e WOEAp x ww         (5.18) where ( ) ( )( )Ex xw x x O SP x   . The integral in Eq.(5.18) could be regular or principal value depending on the x position along the undulated workpiece surface. Since the singular integral results are usually given for the interval of [-1,1], the integral in Eq.(5.18) is normalized by introducing the following change of variables,     ,2 2E Ex x x xO SP O SPx C C          (5.19) Then the new interval becomes,  2121EE xx Ex xxx Ex xOx O CSP OSPx SP CSP O            (5.20) where C can be evaluated as,  Ex xEx xO SPCSP O  (5.21) Introducing Eq.(5.21) into Eq.(5.19), and then inserting Eq.(5.19) into (x)w  yields,  2( ( )) 12Ex xSP Ow x      (5.22) Similarly,  w   becomes,  2( ( )) 12Ex xSP Ow x    (5.23) 100  Introducing Eqs.(5.19),(5.22), and (5.23) into Eq.(5.18) along with the derivative of   i.e.2Ex xO SPd d      gives,    122 21cos (SP O ) O SP( ( )) 12 1 1 ( )E Ex x x xWW 0e WEAp x d               (5.24) Numerator of the integrand in Eq.(5.24) can be approximated using the Chebyshev expansion [95],    '0( ) = cos (SP O ) O SP ( )CNE Ex x x x k kWkf c T          (5.25) where kc  is the contribution coefficient of the kth degree Chebyshev Polynomial of the first kind ( )kT , CN  is the number of polynomial nodes considered, the prime sign indicates that the first term in the series expansion is halved. ( )kT  is defined as,  1( ) cos( cos ( ))kT k   (5.26)  The contribution coefficients in Eq.(5.25) are evaluated as [96],  02(x ) (x )1CNk j k jC jc f TN    (5.27) where x j  are the roots of the polynomial and given as,   0.5cos1jCjxN      (5.28) Inserting the approximation in Eq.(5.25) into Eq.(5.24) yields the elastic pressure contribution of the surface wave as, 101    12 '2 20 1( )( ( )) 12 1 1 ( )CNW 0 k ke WkEA c Tp x d           (5.29) The integral in Eq.(5.29) can be analytically evaluated using Chebyshev Polynomial of the second kind ( )kU  [95],   2 '120( ( )) 1 ( )2 1CNW 0e k kWkEAp x c U              (5.30) where 1( )kU    is,    11 1sin cos ( )( )sin cos ( )kkU    (5.31) Contribution of the wavy surface to the elastic contact pressure is calculated by setting the number of nodes of the Chebyshev Polynomial ( )CN  in Eq.(5.30). The overall elastic contact pressure is obtained by superposing Eq.(5.17) and Eq.(5.30),             1 12 1 112 2 2112 '1201( ) 1 tan ln4 14 1 1 112 1CEx xeeN0k kWkE O SPp xR yEAc U                                          (5.32) where   is obtained by combining Eq.(5.21) with the first expression in Eq.(5.19),  2 Ex xEx xx O SPSP O   (5.33) Indentation of the cutting edge into the newly generated, undulated workpiece surface might result in both elastic and plastic deformations. The contact pressure contributed by the material yield should also be included in the analytical process damping model. As explained by Johnson 102  [97], when the elastic contact pressure exceeds the yield point of the workpiece material, initially elastic contact of the cutting edge element might become elastic-plastic which is contained by the surrounding elastic material. If the indentation becomes more severe, then the deformation becomes fully plastic. Since the boundary between the fully elastic, elastic-plastic, and the fully plastic regions are not analytically known [97], the mean plastic contact pressure ( )pp  sets the upper boundary for the elastic contact pressure [98], and it is given by Tabor [97, 99] in the following general form,    , 1 3p Y ppp C C     (5.34) where pC  is the geometry dependent plastic contact pressure constant and Y  is the yield strength of the work material. pC  quantitatively characterizes the plastic deformation, and it is either obtained from the indentation tests [97, 99, 100], or from extensive FE simulations [4, 101]. In this study, pC  obtained from slip-line field solution of the wedge indentation [98] is adapted to the cutting edge penetration due to their geometric similarity as,  0 02.57 0.023 , 0 30p n nC        (5.35) which covers the practical range for the normal clearance angle (i.e. 0 05 12n  ). The total elastic contact pressure distribution ( (x))ep  in Eq.(5.32) is combined with the plastic limit ( )pp  defined by Eq.(5.34) to obtain the overall contact pressure distribution ( ( ))dp x  and, in turn, to calculate the process damping forces as explained in the following section. 103  5.5 Evaluation of the Process Damping Force and Equivalent Viscous Damping The specific indentation force ,( )d rF  acting in the radial direction of the cutting edge element due to its penetration into the undulated workpiece surface is evaluated by integrating the overall contact pressure distribution ( ( ))dp x  obtained in Section 5.4 over the contact length,  , ( )xExSPd r d dOF p x dx   (5.36) where ,d rF  is in N/m, and ddx  is the discretization resolution for pressure calculation. The cutting edge statically plows some material in the presence of the edge radius ( )eR  even at the chatter free cutting conditions [37, 45, 54, 55]. The static plowing does not change along the vibration cycle, hence its effect should be eliminated by removing the static indentation force from Eq.(5.36). The elastic contact pressure and its plastic limit in the case of static indentation of the cutting edge element can be respectively obtained by replacing ExO  in Eq.(5.17) by xF , and directly from Eqs.(5.34)-(5.35). Then, the specific process damping force equation is written as,    ( ) ( )x xExxSP SPPD d d s srFOF p x dx p x dx     (5.37) where ( )sp x  and sdx  are the overall contact pressure distribution and discretization resolution in static plowing, accordingly. The process damping force acting in the tangential direction of the cutting edge element can be traditionally written as [40, 42-47, 49-51, 54, 55, 102],     PD PDt rF F   (5.38) where   is the Coulomb friction coefficient.  104  Eqs.(5.37)-(5.38) lead to non-linear process damping forces as they only act along half the wavelength of the vibration cycle. Although this non-linearity can be considered in computationally prohibitive time domain chatter stability analysis, it can be circumvented by replacing the damping force by an equivalent linear viscous damper acting in the radial direction [50, 52] for the frequency domain analysis. The approach given in [52] is adapted by replacing the experimentally identified specific indentation force coefficient ( )spK  and the displaced volume ( )V  based process damping force calculation (Eq.(2.1)) with the analytically evaluated indentation force in the radial direction (Eq.(5.37)) as follows.  The energy dissipated by the indentation force (Eq.(5.37)) over half the wavelength can be written as,   344WWWavePD PD rdyE F dxdx    (5.39) and inserting the derivative of the wave form (Eq.(5.2)) into Eq.(5.39) gives,   34042 2cosWWPD PDW WrA xE F dx         (5.40) The energy dissipated by the equivalent specific viscous damping ( )eqc  over the entire wavelength in the radial direction is expressed as,  20 0W WWave Waveeq eq eq cdy dyE c dy c V dxdt dx            (5.41) Similarly, substituting the derivative of the wave form (Eq.(5.2)) into Eq.(5.41) yields, 105   22002 2cosWeq eq c W WA xE c V dx             (5.42) The integral in Eq.(5.42) is analytically defined, and eqE  becomes,  2 202eq eq c WAE c V   (5.43) Equating the energy dissipated by the equivalent damping (Eq.(5.43)) and the indentation force (Eq.(5.40)) gives eqc  (Ns/m2) as,   34402cosWWPD r WeqcxF dxcV A     (5.44) The overall calculation scheme for the analytical estimation of equivalent specific process damping ( )eqc  is summarized in Figure 5.5. The model accepts the cutting edge geometry, vibration parameters, cutting conditions, and the work material properties as inputs. Half wavelength of the surface undulation is defined and discretized as the process damping region, and the corresponding contact pressure distribution (Eqs.(5.32),(5.34)) is evaluated. The overall contact pressure is integrated along the contact length, and the contribution of static plowing is removed (Eqs.(5.37)). The process damping force is represented with the equivalent viscous damping (Eq.(5.44)) to be integrated into the chatter stability analysis that is presented in Chapter 6, while the direct experimental validation results for estimation of the equivalent specific viscous damping ( )eqc  are presented in the next section. 106   Figure 5.5: Flow chart for the analytical estimation of process damping 5.6 Experimental Verification 5.6.1 Verification of Equivalent Viscous Damping (ceq) The proposed analytical process damping model has been verified by comparing the predicted equivalent specific viscous damping ( )eqc  against the experimentally and numerically identified process damping force coefficients for different work materials in the previously conducted cutting tests at UBC-MAL.   Altintas et al. [48] experimentally identified the cutting speed ( )cV  dependent average process damping force coefficient in the radial direction ( )iC  for AISI 1045 work material using the chatter free orthogonal plunge turning tests. A tube with a 35 mm diameter is machined with 0.05 mm/rev feedrate using a piezo actuator providing vibratory motion to the cutting tool. iC  was identified as 0.611610  (N/m), which was related to its tangential component via the friction coefficient ( )  as iC . The cutting edge geometry, vibration and friction parameters 107  used in the tests are given in Table 5.1. The edge radius of the fresh tool is assumed to be 20 m , and the average material separation angle ( ) , which is required to calculate the elastic contact pressure (Eq.(5.32)) and the contact geometry (Eqs.(5.2)-(5.9)), is obtained from the machining literature as 500 [103-107].      Table 5.1: Cutting edge geometry, vibration and friction parameters, and the average separation angle for the orthogonal cutting tests in [48] Material Edge Radius ( )eR m  Rake Angle 0( )n  Clearance Angle 0( )n  Spindle Speed sn  (rpm)  Vibration Amplitude 0A  ( )m  Vibration Frequency c  (Hz) Friction Coefficient,   Average Separation Angle, 0( )  AISI 1045 20 0 7 200-4000  35  20-120 0.3 50  The proposed analytical process damping model employs the temperature dependent properties ( , )YE   of the work material in the contact pressure calculations. Prediction of the workpiece temperature during machining is an active research topic and is not within the scope of the current study. Hence, the workpiece temperature rise beneath the flank face of the cutting edge element due to shearing is estimated using a previously developed [108] finite difference based workpiece temperature prediction code at UBC-MAL to consider its effect on E and Y . The average temperature at 20eR m  below the workpiece surface is used as the workpiece temperature, and it is estimated as 2200C for the given spindle (cutting) speed range in Table 5.1.   The experimentally identified process damping force coefficient is compared against the analytically estimated equivalent specific viscous damping ( )eqc  given in Eq.(5.44) and its corresponding tangential component ( )eqc . Comparison results are presented in Figure 5.6 for 108  lower and upper boundaries of the vibration frequency range (Table 5.1) using the following relation [48],  ieqcCcV   (5.45)  Figure 5.6: Comparison of the experimentally identified equivalent process damping against prediction of the proposed analytical model for AISI 1045 in (a) radial, and (b) tangential directions There is a good agreement between the measured and predicted values. Higher vibration frequency decreases the wavelength of the surface undulations, which leads to more rubbing of the tool against the machined surface, and in turn more process damping particularly towards the low cutting speeds in Figure 5.6.  Jin and Altintas [56] modeled the chatter free orthogonal macro-cutting of AISI 1045 workpiece using the FE method. The average process damping force coefficient ( )iC  and its tangential component are respectively identified as 1.54 510  (N/m) and 0.462 510  (N/m) for 109  the cutting edge geometry, vibration and friction parameters shown in Table 5.2. The corresponding edge radius ( )eR  for the given wear length in Table 5.2 is obtained using Eq.(5.3) as 21.4 m . Numerically identified [56] process damping force coefficients are compared against the prediction of the proposed analytical model in Figure 5.7 for the given cutting speed range (Table 5.2) using Eq.(5.45).  Table 5.2: Cutting edge geometry, vibration and friction parameters, and the average separation angle used in the FE simulations in [56] Material Wear Length ( )wL m  Rake Angle 0( )n  Clearance Angle 0( )n  Cutting Speed cV  (m/min) Vibration Amplitude ( )0A m  Natural Frequency (Hz)n  Friction Coefficient  Separation Angle, 0( )  AISI 1045 80  0 7 50-400 ~3 1500 0.3 50   Figure 5.7: Comparison of the numerically (FE) identified equivalent process damping against prediction of the proposed analytical model for AISI 1045 in (a) radial, and (b) tangential directions 110  The average workpiece temperature at 21.4 m  below the surface for the cutting speed range in Table 5.2 is estimated using the method in [108] as 2020C, and E and Y  are selected accordingly. On top of the input parameter uncertainty in the proposed model, the increasing error towards the low cutting speed region can be partially attributed to the way the temperature effect is considered in the analytical model as compared to the coupled thermo-mechanical FE analysis used in [56]. Also, the wavelength of the surface undulations decreases with the decreasing cutting speed, which increases the local slope of the wave causing more dynamic indentation. This might have amplified the prediction error, which is more observable at low cutting speeds since the predicted process damping variation does not reach to its vertical asymptote in Figure 5.7.  Overall, the analytically approximated equivalent specific viscous damping ( )eqc  agrees well with the experimentally and numerically identified values in radial and tangential directions in Figure 5.6 and Figure 5.7. This corroborates that the proposed process damping model can be reasonably employed in place of impractical empirical and computationally intense FE based identifications considering that each point is a separate experiment and FE analysis in Figure 5.6 and Figure 5.7, respectively. Ahmadi and Altintas [54] experimentally identified the additional modal damping due to interference of the cutting edge with the finished workpiece surface from the chatter free orthogonal cutting tests using the operational modal analysis method. AISI 1018 and Ti6Al4V were used as work materials. The cutting edge geometry, friction and vibration parameters used in [54], and the average material separation angles obtained for AISI 1018 and Ti6Al4V from literature [105, 109, 110] are given in Table 5.3. 111  Table 5.3: Cutting edge geometry, vibration and friction parameters, and the average separation angle for the orthogonal cutting tests in [54] Material Edge Radius ( )e mR   Rake Angle 0( )n  Clearance Angle 0( )n  Cutting Speed cV  (m/min) Vibration Amplitude ( )0A m  Natural Frequency (Hz)n  Friction Coefficient  Separation Angle, 0( )  AISI 1018 60 0 7 55-300 3-3.6 1403 0.72 48 Ti6Al4V 30 0 6 85-210 2.5-4.5 1397 - 40  The tool tip vibrations were measured with the attached accelerometer at different cutting speeds in [54]. In this study, the tool tip vibrations are calculated from the measured accelerations [54], and the representative results for AISI 1018 and Ti6Al4V are given in Figure 5.8 and Figure 5.9, respectively. The average of the vibration amplitudes given in Table 5.3 are employed in the simulations.    Figure 5.8: a-b) Measured acceleration at the tool tip, and its frequency content using Fast Fourier Transform (FFT), c-d) corresponding displacement and its FFT for AISI 1018 at 2000 rpm  112   Figure 5.9: a-b) Measured acceleration at the tool tip, and its frequency content using Fast Fourier Transform (FFT), c-d) corresponding displacement and its FFT for Ti6Al4V at 850 rpm The process damping ratio ,( )p r  is obtained [54] by subtracting the measured structural damping ( )r  from the identified total damping ( ) , and it is converted to the process damping coefficient ,( )p rc  using the measured modal parameters. ,p rc  is related to the equivalent specific process damping ( )eqc  through the width of cut ( )pb  as,  , ,2( )p r r n r r p eqc m b c       (5.46) where  rm  is the modal mass. The identified ,p rc  is converted to the same form as eqc  (Eq.(5.46)) and compared against the analytically obtained eqc  in Figure 5.10. In order to consider the effect of workpiece temperature beneath the flank face of the cutting edge, the average workpiece temperature for the 113  cutting speed range in Table 5.3 is estimated in dry cutting conditions [108] as 2900C and 4600C for AISI 1018 and Ti6Al4V, respectively.   Figure 5.10: Comparison of the experimentally identified equivalent process damping against prediction of the proposed analytical model for (a) AISI 1018, and (b) Ti6Al4V work materials While the proposed model accurately predicts the process damping at different cutting speeds for AISI 1018, it is overpredicted for Ti6Al4V although the general trend is captured. The uncertainty in the geometric input parameters becomes more important for Ti6Al4V since the cutting edge profile (especially eR ) changes due to the fast tool wear during dry machining of Ti6Al4V. In [54], however, the measurement condition (before or after cutting) of the cutting edge radius is not given, which might have affected the analytical prediction in Figure 5.10(b). Also, the thermally resistant Ti6Al4V makes it more sensitive to the coupled thermo-mechanical effect, which is only considered through the approximated workpiece temperature in the current work. 114  Although satisfactory results are obtained for the experimental and numerical comparison cases presented in this section, the effect of uncertainty in the input parameters (Figure 5.5) needs to be studied. Sensitivity of the predicted eqc  to the individual uncertainty of the effective parameters is investigated by defining the simulation results of this section as the nominal predictions as follows.      5.6.2 Parameter Uncertainty Analysis The proposed analytical model is advantageous for a quick analysis of the model sensitivity to the individual input parameter uncertainties. While the basic cutting edge geometry ( , )n n   is provided by the cutting tool manufacturers with high accuracy, the cutting edge radius ( )eR  is usually measured. The material separation angle ( )  obtained from the literature for the validation cases in the previous section has a range of values, and an average value is employed. Therefore, eR  and   are selected as the geometric parameters having the most uncertainty. Also, the vibration amplitude 0( )A  is obtained from the measured acceleration data in some of the experiments [54]. Hence, uncertainty of 0A  and its effect on the specific equivalent process damping ( )eqc  is also investigated. Finally, the mean plastic contact pressure ( )pp  is obtained using an empirical constant ( )pC . In order to account for the error induced into the estimation of pp , its uncertainty is also studied in this section.     Due to the piezo actuator used for vibratory motion of the tool in [48], the amplitude of vibration 0( )A  is not considered as an uncertain variable. 20% , 20% , and 30%  variations are respectively introduced into the previously used values of cutting edge radius ( )eR , separation angle ( ) , and the mean plastic contact pressure ( )pp . The reason for the higher 115  uncertainty in pp  is due to the higher number of parameters affecting it. Comparison of the nominally predicted equivalent specific process damping ( )eqc  with the simulations having different uncertainties is presented for c =120 Hz in Figure 5.11. As evident, eqc  is more sensitive to the unit change in   than pp , and larger eR  and   increase the indentation force, which in turn results in higher eqc  as expected.                 Figure 5.11: Uncertainty analysis for the parameters: edge radius, material separation angle, and the plastic limit for the test case in [48] A similar analysis is conducted for the numerical case given in [56] with the same uncertainties for the same input parameters. Comparison of the nominal prediction and the simulations with different uncertainties is shown in Figure 5.12. While eqc  is more sensitive to the unit change in the separation angle ( )  than the mean plastic contact pressure ( )pp  for the cutting speeds cV   60 m/min, it becomes equally sensitive to the unit change in   and pp  for  116  cV  60 m/min. This could be attributed to the fact that the increase in the indentation area and the associated penetration force with increasing   becomes less pronounced at low cutting speeds due to the shorter wavelengths. Therefore, the dominancy of   on eqc  diminishes with the decreasing cutting speed, and the effect of pp  on eqc  becomes equally significant.          Figure 5.12: Uncertainty analysis for the parameters: edge radius, material separation angle, and the plastic limit for the test case in [56] Besides the uncertain inputs in the previous cases, 20%  uncertainty is also added to the vibration amplitude ( 0A ) obtained from the experimentally measured acceleration of the tool tip to simulate the process induced additional damping for the cases given in [54]. Comparison of the nominal prediction and the simulations with different uncertainties for AISI 1018 and Ti6Al4V work materials are shown in Figure 5.13 and Figure 5.14, respectively. For both materials, equivalent specific process damping ( )eqc  is almost insensitive to the change in 0A  due to the small nominal vibration amplitudes (3.3 m  and 3.5 m ). In Figure 5.13, unit change 117  in the material separation angle ( )  is the main source of the variation in eqc  for the cutting speeds cV   100 m/min, while the mean plastic contact pressure ( )pp  becomes responsible for the largest change in eqc  for the low cutting speeds ( cV  100 m/min) due to the same reasoning as the previous case. On the other hand, eqc  is most sensitive to the unit change in   for the entire cutting speed range for Ti6Al4V in Figure 5.14.  Figure 5.13: Uncertainty analysis for the parameters: vibration amplitude, edge radius, material separation angle, and the plastic limit for the first test case in [54] Overall, the material separation angle ( )  and the plastic contact pressure limit ( )pp  are found to be the most significant geometric and deformation parameters affecting the proposed equivalent viscous damping ( )eqc  for the materials and the cutting conditions used in the validation tests. 118   Figure 5.14: Uncertainty analysis for the parameters: vibration amplitude, edge radius, material separation angle, and the plastic limit for the second test case in [54]  5.7 Summary This chapter presents a novel process damping model that considers the cutting edge geometry, cutting conditions, material properties of the workpiece, and the properties of the wave imprinted on the part surface. The elastic-plastic contact pressure at the interface between the cutting edge and workpiece surface is analytically evaluated at the discrete locations along half of the wavelength. The non-linear specific contact force is obtained by integrating the overall contact pressure distribution over the contact length, and it is linearized for the chatter stability analysis by converting it to a representative viscous damping. It is shown that the proposed analytical model can predict the additional damping caused by interference of the cutting edge with the undulations left on the newly generated surface, and it can be used in place of the present cumbersome empirical and computationally intense numerical models. The 119  developed varying workpiece dynamics update and process damping models are integrated into a new five-axis chatter stability model in the following chapter.    120  Chapter 6: Modeling of Chatter Stability in Frequency Domain  6.1 Overview The productivity of multi-axis thin-walled part machining is usually limited by the unstable chatter (self-excited) vibrations due to the interaction between the cutting forces and the flexibilities of the workpiece and the tool. Chatter vibrations cause the chip thickness to modulate (regeneration) between the consecutive flutes of the tool. Unless avoided, they lead to exponentially growing cutting forces and vibrations, which result in a damaged tool and part, overloading the machine, poor part surface quality, and violation of the part dimensions. Therefore, chatter needs to be predicted and avoided by selecting stable cutting conditions before the thin-walled part is machined.  This chapter develops a general dynamics and chatter stability prediction model for five-axis point (ball end) milling of thin-walled parts with free from surfaces using the varying thin-walled workpiece dynamics (Chapters 3 & 4) and process damping (Chapter 5) along the complex toolpaths. The ball-end mill is axially discretized to consider the local geometry and tool-workpiece engagement (TWE) variations along the cutter axis (Section 6.2.1). The static chip thickness is calculated by considering the five-axis milling kinematics while the dynamic chip thickness is modeled using the cutting tool and flexible workpiece vibrations (Section 6.2.2). The static and dynamic cutting forces are evaluated using the modeled chip thickness (Section 6.2.3). The process damping forces are estimated through the introduced equivalent viscous damping. The characteristic equation of the flexible milling system is obtained by summing the contribution of each axial disk to the total dynamic milling forces, which takes the relative dynamic displacements between the workpiece and tool into account (Section 6.2.4). Chatter stability of the ball-end milling is checked (stable or unstable cut), and the chatter frequency is 121  identified for the unstable locations using the Nyquist stability criterion along the toolpath (Section 6.2.5). The proposed analytical process damping model is further validated by comparing the predicted stability lobe diagrams with the experimental data for flat-end milling and orthogonal turning operations (Section 6.3.1-6.3.2). The developed chatter stability model is experimentally verified on a thin-walled plate and an impeller blade with and without process damping (Sections 6.3.3-6.3.4). Also, the practical applications of the developed models are demonstrated by predicting the dimensional form errors and vibrations of thin-walled parts (Sections 6.3.5-6.3.6). The chapter concludes with a summary in Section 6.4.       6.2 Dynamics of Five-Axis Ball-End Milling Five-axis ball-end milling of flexible parts have the following specific challenges: (i) dynamics of the machined part change as the material is removed along the toolpath, and also the cutter axis, (ii) varying local tool geometry (helix, normal rake and clearance angles) changes the cutting speed, and in turn the process damping conditions along the tool axis (z), (iii) cutting force coefficients are both chip thickness and cutting speed dependent, and they change both in tool axis (z) and tangential (te) directions, (iv) the engagement boundaries vary along the tool axis due to the five-axis machining and machined surface features, (v) the updated workpiece dynamics is defined in the workpiece coordinate system (WCS), and it needs to be transformed into the tool coordinate system (TCS), (vi) the tool dynamics is usually measured in the fixed machine coordinate system (MCS), and it needs to be transformed into TCS at discrete points along the toolpath.  In order to accommodate the variations in the local tool geometry and the tool-workpiece engagement (TWE) along the tool axis, the cutting tool is divided into infinitesimal axial elements (disks) having the height of dz in TCS, i.e., kth disk has the axial position of kz kdz . 122  Complex TWEs throughout the five-axis toolpath are obtained from MACHproTM-Virtual Machining System developed in the Manufacturing Automation Laboratory at UBC. The three dimensional engagement region on the cutting tool is projected into TCS at discrete stations (maps) along the toolpath, and TWE is assumed unchanged between two successive map points MAP  and 1MAP  . TWE is represented in TCS with the cut start ( )st  and cut exit ( )ex  angles of the cutting flutes at different elevations along the tool axis as exemplified in Figure 6.1.  Figure 6.1: A sample tool-workpiece engagement (TWE) in five-axis ball-end milling, its discretization, and variation of st  and ex  along z-axis  6.2.1 Ball-End Mill Geometry The geometric model of the general end mill presented in [111] is adapted in this thesis. The required geometric features of the ball-end mill is briefly summarized from [4] here for the sake of notation compatibility and completeness.  The helical ball-end mill is composed of a spherical bottom and a cylindrical shank as shown in Figure 6.2(a). Tool’s Cartesian coordinate system ( , , )T T Tx y z  is located at the tool tip, and the rotating local radial, tangential, and axial coordinate frame ( , , )e e er t a  of the infinitesimal cutting edge element with the height of dz has a different orientation at each z=zk  along the tool axis. 123   Figure 6.2: Ball-end mill envelope and tool coordinate system (TCS), a) front view, b) bottom view [112]  The transformation between these two coordinate frames is carried out using the following transformation matrix ( )rtaxyzT ,  sin (t, )sin ( ) cos (t, ) sin (t, )cos ( )(t, ) cos (t, )sin ( ) sin (t, ) cos (t, )cos ( )cos ( ) 0 sin ( )j j j j jrtaxyz j j j j jj jz z z z zz z z z z zz z                  T   (6.1)  (t, ) , (t, )j j j jrta xyzj xyz j j rta jj j j jx r r xy z t t z yz a a z                                    T T   (6.2) where rtaxyzT  is an orthogonal matrix 1[ ] [ ]rta rta T xyzxyz xyz rta  T T T . The local cutter radius (R(z)) and the local helix angle ( i z ( ) ) are expressed as,  20 000 0,( ),1 1>zR zzzR RRR R         (6.3) 124   10 000 0tan tan>iiiR zz RRzz R        ( ),( ),  (6.4) where 0R  and 0i  are the radius and the helix angle of the cylindrical part of the tool. The local radius dependent cutting speed (mm/s) is,  2 (z)(z)60scR nV   (6.5) where sn  is the spindle speed (rev/min). Normal orientation of the two dimensional cutting edge element (Chapter 5) is defined in a plane normal to the cutting edge. Then, the local normal rake angle ( ))( n z  is defined for the ball-end tool as,  ( ) ( )tan tan cos ( )n o iz z z     (6.6) where ( )o z  is the local radial rake angle and defined as,  00))( sin ( ) ,( )( ,cylcyloooz z Rzz R   (6.7) where )( cylo  is the radial rake angle on the cylindrical part of the tool. The axial orientation of the cutting edge element is defined by the axial immersion angle ( ( ))z  of the jth flute at elevation z using the local radius as,  10( )( ) sinR zzR       (6.8) In order to determine the immersion conditions for each flute, the instantaneous angular position (immersion angle) of a point on the jth flute at elevation z is evaluated in the Tx - Ty  plane (Figure 6.2(b)) as, 125   ,1(t, ) (t) (z) ( )j p j jjz z            (6.9) where (t)  is the instantaneous rotation angle of the tool measured from Ty -axis to Tx -axis, , (z)p j  is the angular spacing between the consecutive flutes (pitch angle), ( )j z  is the lag angle due to the helix angle. (t)  is evaluated as,  2(t)60sn t    (6.10) The pitch angle for the uniformly and non-uniformly spaced flutes is defined as,  ,,uniform pitchnon-uniform pitch(z)2,,fp jp jN   (6.11) where fN  is the total number of flutes. Points along the jth flute get into cut with a delay due to the lag angle, and it is defined for constant lead tools as,   tan (z)(z)ijzzR    (6.12) Having mathematically defined the ball-end mill geometry, the dynamic chip thickness is derived in the following section. 6.2.2 Static and Dynamic Chip Thickness The total chip thickness of the infinitesimal cutting edge element (disk) on the jth flute ( (t, z))jh  is consisted of a static component ,( (t, z))s jh  generated by the relative rigid body (feed) motion between the tool and workpiece, and also a dynamic component d,( (t,z))jh  due to the structural vibrations resulting from the excited modes of the thin-walled workpiece and the slender ball-end mill. Then, (t, z)jh  can be expressed as, 126   , d,(t,z) ( (t,z) (t,z)) ( (t,z))j s j j j jh h h      (6.13) where ( (t,z))j j   is the unit step function determining if the jth flute at elevation z is in or out of the cut,  1 , (t, z)( (t, z))0 ,st j exj jelse       (6.14) Since the cutting forces and the chatter stability are evaluated in TCS, the chip thickness is measured in the surface normal (radial) direction ( ( , ))j t zn  of the jth flute at elevation z, and it is defined in TCS as,    00sin (t, z)sin (z) cos (t, z)sin (z) cos (z) ,( , )sin (t, z) cos (t, z) 0 , ( 2)j j j j jjj j jz Rt zz R           n   (6.15) The static chip thickness ,( (t, z))s jh  is calculated by projecting the total feed vector ( (z))TF  onto ( , )j t zn  in TCS as [4],   ,,(z)(t, z) ( , ) (t, z)2p js j j Tsh t zn n F   (6.16) where (t, z)TF  is superposition of the angular feed ( (t, z))AF  resulting from the rotational motion between two consecutive discrete locations (maps) along the five-axis toolpath, and the linear feed ( (t, z))LF  due to the relative rigid body motion between the tool and the workpiece. Detailed calculation of (t, z)TF  at discrete stations along complex five-axis toolpaths can be found in [4].  The static chip thickness varies (Figure 6.3(a)) both in axial and tangential directions of the ball-end tool due to its geometry. Since the cutting forces are equally felt by the thin-walled 127  workpiece and the slender ball-end mill in opposite directions, the forces may excite their structural vibration modes in all three Cartesian coordinates. Then, the chip thickness dynamically varies in the presence of the relative vibrations between the flexible workpiece and tool as depicted in Figure 6.3(b).   Figure 6.3: a) Representative variation in static chip thickness of the first flute with angular and axial positions, b) three dimensional dynamic chip formation in ball-end milling The dynamic chip thickness d,( (t,z))jh  is composed of the present relative vibrations and the vibrations imprinted on the workpiece surface one tooth period before as illustrated in Figure 6.4 for a cutting edge element at elevation z in Tx - Ty  plane. Then, the original expression [57] for d, (t, z)jh  can be extended by considering the relative vibrations between the tool and workpiece as,     d, 1,T 1, ,T ,(t, (z), z) (t (z), z) (t (z), z) (t, z) (t, z)j j j j j W j j j Wh q q q q           (6.17) 128  where 1,Tjq  , ,Tjq  and 1,j Wq  , ,j Wq  are respectively the previous and current vibrations of the jth flute and the flexible workpiece in the radial direction. (z)j  (sec) is the time delay between the current and previous vibration marks left on the surface, and it is equal to,  ,60 (z)(z)2p jjsn   (6.18)  Figure 6.4: Chip regeneration and dynamic chip thickness on the cross section of the jth flute at elevation z Current and previous dynamic displacements of the slender ball-end mill are respectively obtained by projecting its vibrations into the chip thickness (radial) direction using the transformation in Eq.(6.1) as,  1,T 1,T1,T1,T(t (z), z) (t (z), z) sin (t, z) sin (z)(t (z), z) cos (t, z) sin (z)(t (z), z) cos (z)j j j j j jj j j jj j jq xyz             (6.19) where 1(t,z) (t (z),z)j j j    , and 129   ,T ,T ,T,T(t, z) (t, z) sin (t, z) sin (z) (t, z) cos (t, z) sin (z)(t, z) cos (z)j j j j j j jj jq x yz      (6.20) Similarly, 1,j Wq   and ,j Wq  are,  1, 1,1,1,(t (z), z) (t (z), z) sin (t, z) sin (z)(t (z), z) cos (t, z) sin (z)(t (z), z) cos (z)j W j j W j j jj W j j jj W j jq xyz             (6.21) and  , , ,,(t, z) (t, z) sin (t, z) sin (z) (t, z) cos (t, z) sin (z)(t, z) cos (z)j W j W j j j W j jj W jq x yz      (6.22) Inserting Eqs.(6.19)-(6.22) into Eq.(6.17) and rearranging it yields,           d, , 1, , 1,, 1, , 1,, 1, ,(t, , z) (t, z) (t (z), z) (t, z) (t (z), z) sin (t, z) sin (z)(t, z) (t (z), z) (t, z) (t (z), z) cos (t, z) sin (z)(t, z) (t (z), z) (t, z)j j j T j T j j W j W j j jj T j T j j W j W j j jj T j T j j Wh x x x xy y y yz z z z                         1, (t (z), z) cos (z)j W j j    (6.23) which can be rewritten in a more compact form using Eq.(6.15),  , ,h (t, (z),z) ( , ) (t, (z),z)d j j j j T W jt z  n Δ   (6.24) where,  ,, ,,(t, (z), z)(t, (z), z) (t, (z), z)(t, (z), z)j T W jj T W j j T W jj T W jxyz          Δ   (6.25) where ,j T Wx  , ,j T Wy  , and ,j T Wz   are the relative vibrations between the thin-walled workpiece and the tool from previous and current vibrations in TCS. 130  6.2.3 Dynamic Cutting and Process Damping Forces The cutting force model based on the orthogonal to oblique transformation presented in [4] is employed for five-axis ball-end milling in this study. Cutting forces are linearly dependent on the chip thickness and the cutting force coefficients. The differential forces ,( (t, z))rta jdF  acting on the cutting element at axial position z in , ,e e er t a  frame are composed of the static cutting ,( (t, z))csrta jdF  and static plowing ,( (t, z))psrta jdF  forces, dynamic cutting ,( (t, z))cdrta jdF  force due to the dynamic chip thickness, and dynamic plowing/indentation ,( (t, z))pdrta jdF  forces that are caused by the interference between the cutting edge element’s flank face and the wavy workpiece surface. All force components are superposed as,  , , , ,1(t, z) ( (t, z)) ( (t, z) (t, z) (t, z) (t, z))fNcs cd ps pdrta j rta j rta j rta j rta jjd d d d d    F F F F F   (6.26) Since the static forces do not take any role in chatter stability of the milling system [57, 113], Eq.(6.26) is reduced to the following equation for the dynamic milling case,  , ,1(t, z) (t, z) (t, z)fNd cd pdrta rta j rta jjd d d F F F   (6.27) The dynamic cutting forces are expressed as,  , , ,(t, z) (t, z) h (t, (z), z) (z)cd crta j rta j d j jd dbF K   (6.28) where  , (t, z) (t, z) (t, z) (t, z)Tcrta j rc tc acK K KK  is the vector consisting of the cutting force coefficients that are obtained from orthogonal to oblique transformation [114] in the radial, tangential, and axial directions, (z)db  is the chip width of the oblique flute segment and evaluated as [4], 131   (z)sin (z)jdzdb   (6.29) Substitution of Eqs.(6.24) and (6.29) into Eq.(6.28) gives,   ,, , ,,(t, (z), z)(t, z) (t, z) sin (t, z) cos (t, z) cot (z) (t, (z), z)(t, (z), z)j T W jcd crta j rta j j j j j T W jj T W jd dz            xF K yz (6.30) or, more compactly,  , , ,(t,z) (t,z) (t,z) (t, (z),z)cd cjrta j rta j j T W jd dzF K n Δ   (6.31) The second component in Eq.(6.27) is the non-cutting (damping) force, which is modeled in Chapter 5 and  integrated into Eq.(6.27) as follows. The process damping force acting on the cutting edge element at elevation z is written by replacing the indentation volume based process damping force calculation with the analytically estimated equivalent viscous damping ( )eqc  as,  r, , ,t, r,(t, z) (z) (t, z) (z)(t, z) (t, z)pdj eq j j T Wpd pdj jdF c q dSdF dF  (6.32) where , (z)eq jc is the position dependent equivalent process damping coefficient due to varying cutting speed along the ball-end tool axis (Eq.(6.5)), , (t,z)j T Wq   is the current relative surface vibration velocity in the radial direction of the cutting edge, and (z)dS  is the differential length of the cutting edge element. Although the exact expression of (z)dS  is given in [4], it approaches to (z)db  for infinitesimal dz ( (z) (z))dS db  [62, 115, 116]. Similar to the relative vibrations in the radial direction (Eq.(6.23)), , (t,z)j T Wq   is written as, 132   , , ,, ,, ,(t, z) ( (t, z) (t, z))sin (t, z)sin (z)( (t, z) (t, z))cos (t, z)sin (z)( (t, z) (t, z))cos (z)j T W j T j W j jj T j W j jj T j W jq x xy yz z        (6.33) Introducing Eqs.(6.29) and (6.33) into the first expression in Eq.(6.32) gives,  , ,r, , , ,, ,( (t, z) (t, z))sin (t, z)(t, z) (z) ( (t, z) (t, z))cos (t, z)( (t, z) (t, z))cot (z)j T j W jpdj eq j j T j W jj T j W jx xdF c y y dzz z          (6.34) Then the overall process damping forces acting on the cutting edge element at elevation z can be written using Eq.(6.32) as,  , ,rta, , , ,, ,( (t, z) (t, z))sin (t, z)1(t, z) (z) ( (t, z) (t, z))cos (t, z)0 ( (t, z) (t, z))cot (z)j T j W jpdj eq j j T j W jj T j W jx xdF c y y dzz z                (6.35) which can be rewritten in a more compact form similar to Eq.(6.31) as,   rta, , ,(t,z) (z) 1 0 (t,z) (t,z)Tpdj eq j j j T WdF c dz  n Δ   (6.36) where , (t, z)j T WΔ  is the current relative surface vibration velocity vector,  , ,, , ,, ,(t, z) (t, z)(t, z) (t, z) (t, z)(t, z) (t, z)j T j Wj T W j T j Wj T j Wx xy yz z      Δ   (6.37) The total dynamic force exerted on the cutting edge element of the jth flute at elevation z is obtained by inserting Eqs.(6.31) and (6.36) into Eq.(6.27) as,  , , , , ,(t, z) (t, z) (t, z) (t, (z), z) (z) 1 0 (t, z) (t, z)Td crta j rta j j j T W j eq j j j T Wd c dz      F K n Δ n Δ  (6.38) 133  The cutting forces in the rotating frame are transformed into TCS by premultiplying Eq.(6.38) by the transformation matrix (t, )rtaxyz zT  given in Eq.(6.1) as,  xyz, , , ,(t, z) (t, z) (t, (z), z) (z) (t, z) (t, z)d cd pdj j j T W j eq j j j T Wd c dz     F A Δ B Δ   (6.39) where (t, z)cdjA  and (t, z)pdjB  are time and z position dependent periodic coefficient matrices for the dynamic cutting and process damping forces for the jth cutting flute at elevation z,  , , ,, , ,, , ,, , ,, , ,(t, z) (t, z) (t, z)(t, z) ( (t, z)) (t, z) (t, z) (t, z)(t, z) (t, z) (t, z)(t, z) (t, z) (t, z)(t, z) ( (t, z)) (t, z) (t, z) (t,xx j xy j xz jcdj j j yx j yy j yz jzx j zy j zz jxx j xy j xz jpdj j j yx j yy j yz ja a aa a aa a ab b bb b b        AB, , ,z)(t, z) (t, z) (t, z)zx j zy j zz jb b b       (6.40) After the necessary matrix manipulations, the individual entries of (t, z)cdjA  can be found as,       ,,,1(t, z) (t, z)sin 2 (t, z) (1 cos 2 (t, z)) (t, z)sin (z) (t, z)cos (z)21(t, z) (t, z) 1 cos 2 (t, z) sin 2 (t, z) sin (z) cos (z)2(t, z) (t, z)sin (t, z)cos (z) sin (t, z)xx j ct j j cr j ca jxy j ct j j cr j ca jxz j cr j j ja K K Ka K K Ka K                        ,,cot (z) (t, z)cot (t, z) (t, z)cos (z)1(t, z) (t, z) cos 2 (t, z) 1 sin 2 (t, z) (t, z)sin (z) (t, z)cos (z)21(t, z) 1 cos 2 (t, z) (t, z)sin (z) (t, z)cos (z) (t, z)sin 2 (t2j ct j ca jyx j ct j j cr j ca jyy j j cr j ca j ct jK Ka K K Ka K K K                  ,,,, z)(t, z) (t, z)cos (t, z)cos (z) sin (t, z)cot (z) (t, z) (t, z)cot (t, z)cos (z)(t, z) sin (t, z) (t, z)cos (z) (t, z)sin (t, z)(t, z) cos (t, z) (t, z)cos (z) (t, z)sinyz j cr j j j j ct ca j jzx j j cr j ca jzy j j cr j caa K K Ka K Ka K K               ,(z)(t, z) cos (z) (t, z) (t, z)cot (z)jzz j j ca cr ja K K    (6.41) 134  Similarly, (t, z)pdjB  is composed of the following equations,        ,,,yx,yy,1(t, z) 1 cos 2 (t, z) sin (z) sin 2 (t, z)21(t, z) sin 2 (t, z)sin (z) 1 cos 2 (t, z)2(t, z) sin (t, z)cos (z) cos (t, z)cot (z)1(t, z) sin 2 (t, z)sin (z) 1 cos 2 (t, z)2(txx j j j jxy j j j jxz j j j j jj j j jjbbbbb                         ,zx,zy,zz,1, z) 1 cos 2 (t, z) sin (z) sin 2 (t, z)2(t, z) cos (t, z)cos (z) sin (t, z)cot (z)(t, z) sin (t, z)cos (z)(t, z) cos (t, z)cos (z)1 cos 2 (z)1(t, z)2 sin (z)j j jyz j j j j jj j jj j jjjjbbbb                     (6.42) Although the axial position dependency of (t, z)cdjA  and (t, z)pdjB  is taken care of by axially discretizing the cutting tool, the time dependency cannot be included in the frequency domain analysis. Since (t, z)cdjA  (t, z)pdjB  include the switching function ( ( (t,z)))j  , they are both periodic functions and can be represented with Fourier series expansion at the uniform tooth passing frequency ( )T  for the regular pitch cutters, or at a different tooth passing frequency for each flute ,( (z))T j  for the variable pitch cutters. For the most general case of variable pitch cutters, the Fourier Series expansion yields,  , ,(z) (z), ,(t, z) , (t, z)T j T jir t ir tcd cd pd pdj r j j r jk k k kr re e                   A A B B   (6.43) and the corresponding Fourier coefficients are, 135  , ,(z) (z)(z) (z), ,0 01 1(t, z) , (t, z)(z) (z)j jT j T jT Tir t ir tcd cd pd pdr j j r j jk k k kj je dt e dtT T                  A A B B  (6.44) where (z)jT  is the tooth passing period for the jth flute at elevation z. The time variable in Eq.(6.44) can be replaced with the rotation angle of the jth flute ( )j  as shown in Appendix B, then Eq.(6.44) becomes,    ,,,,2,02,01(t, z)(z)1(t, z)(z)p jp j jp jp j jircd cdr j j jk ks jirpd pdr j j jk ks je dTe dT                A AB B  (6.45) Since the directional coefficients are non-zero only if the flute is in cut, Eq.(6.45) can be written as,    ,,, , ,2, , , ,,, , ,,2,,1 1(t, z)(z)1 1(t, z)(z)exp j jstexp j jstr r rxx j xy j xz jircd cd r r rr j j j yx j yy j yz jk ks j p jr r rzx j zy j zz jkrxx j xyirpd pdr j j jk ks j p ja a ae d a a aTa a ab be dT                            A AB B, ,, , ,, , ,r rj xz jr r ryx j yy j yz jr r rzx j zy j zz jkbb b bb b b        (6.46) where elements of the matrices (t, z)cdjk  A  and (t, z)pdjk  B  are given in Eqs.(6.41)-(6.42). Solutions of the general integral equations in Eq.(6.46) are given in Appendix B. The results in Appendix B do not consider the change in the cutting force coefficients ,( (h (t, z)))crta j sK  due to the chip thickness variation in the tangential direction (Figure 6.3(a)). Since the amplitude of the dynamic chip thickness is not known in frequency domain analysis [64], the static chip 136  thickness is employed to average out , (h (t, z))crta j sK  on the jth flute at axial position z by adapting the method given in [62] as,  ,st,(z), ,, st, (z)1(z) (h (t, z))(z) (z)ex jjc crta j rta j sex j jd  K K   (6.47) It is shown in [57, 64, 113] and also in Section 6.3 of this study that taking only 0r    , ,,cd pd0 j 0 jk k      A B  is sufficient to detect the chatter occurrence and to predict the fundamental chatter frequency ( )c  for a given set of cutting conditions, tool and wokpiece pair, and the boundary conditions st( , )ex   at any location along the five-axis ball-end milling toolpath. 6.2.4 Relative Dynamic Displacements between the Tool and Thin-Walled Workpiece The relative dynamic displacements (vibrations) between the cutting tool and low rigidity workpiece are explicitly given in Eqs.(6.23) and (6.25). The relative vibrations with the time delay ( (z))j  between the current and previous vibrations on the machined surface can be represented in the Laplace domain,          ,, ,, , ,, ,( ) (s)(s) 1 ( ) (s)( ) (s)j kj T j Wj T W j T j Wkj T j Wk ksx s xe y s yz s z                              Δ   (6.48) where subindex k is used to show the axial position dependency. The tool and the workpiece vibrations in Eq.(6.48) can be expressed using their frequency response functions (FRFs),  x, , , , , , ,, , , , , , ,, , , , , ,,(s)(s) (s) (s) (s)(s) (s) (s) (s) (s)(s) (s) (s)(s) (s)dj T xx j T xy j T xz j Tdj T yx j T yy j T yz j T ydzx j T zy j T zz j Tj T zkk kdFx H H Hy H H H dFH H Hz dF                              (6.49) 137  and,  x, , ,W , ,W , ,W, , ,W , ,W , ,W, ,W , ,W , ,W,(s)(s) (s) (s) (s)(s) (s) (s) (s) (s)(s) (s) (s)(s) (s)dj W xx j xy j xz jdj W yx j yy j yz j ydzx j zy j zz jj W zkk kdFx H H Hy H H H dFH H Hz dF                               (6.50) where a negative sign shows the direction of the cutting forces on the workpiece. Inserting Eqs.(6.49)-(6.50) into Eq.(6.48) yields,       , ,x,, , ,(s)( ) 1 ( ) ( ) (s) 1 ( )(s)j k j kds dj T Wj T W j W j T y kk k kdz ksdFs e s s dF e sdF                    Δ H H Δ  (6.51) where  , ( )j T WksΔ  is the current relative vibrations vector, , ( )j W ks  H  and , ( )j T ks  H  are 3-by-3 FRF matrices of the thin-walled workpiece and the ball-end cutter for the kth axial disk of the jth cutting flute. Although , ( )j T ks  H  is time independent, , ( )j W ks  H  varies with position as the material is removed along the toolpath as modeled in detail in Chapters 3 & 4. Therefore, the updated workpiece FRFs at discrete locations along the toolpath are employed, and they are assumed to be unchanged between two successive FRF update points.    6.2.5 Dynamic Milling Equation and Chatter Stability This section evaluates the cutting and process damping forces acting on the kth axial disk of the jth flute and sums the contribution of all the engaged flutes to calculate the total dynamic forces exerted on the cutting tool and flexible workpiece. The governing equation is derived by writing Eq.(6.39) with the time invariant, axial position and angular immersion dependent directional factors ,cd0 jk  A  and ,pd0 jk  B  as, 138        xyz, , , , , ,,(t) (t, ) ( ) (t)d cd pdj 0 j j T W j k eq j k j T W0 jkk kkkd c dz           F A Δ B Δ   (6.52) Eq.(6.52) can be rewritten in the Laplace domain as,       xyz, , , , ,,(s) (s) ( ) (s)d cd pdj 0 j j T W eq j k j T W0 jkk kkkd c dz           F A Δ B Δ   (6.53) where  , (s)j T W kΔ  is defined in Eq.(6.51), and the current relative surface vibrations velocity vector can be represented in the Laplace domain by,     ,, (s) (s)j T Wj T Wkks  Δ Δ   (6.54) Substituting Eqs.(6.51) and (6.54) into Eq.(6.53) leads to,            , , ,xyz, , , ,(s) (1 ) ( ) ( ) (s)j kd cd pdj T W j T Wj 0 j eq j k 0 jk kk k ksd e s c s dz         F A Δ B Δ   (6.55) Inserting  , (s)j T WkΔ  from Eq.(6.51) into the dynamic milling equation in Eq.(6.55), and rearranging it yields,       ,xyz, , , , , xyz,,(s) 1 ( ) ( ) ( ) (s)j kd cd pd dj 0 j eq j k j W j T j0 j k kk kk ksd e c s s s d dz                    F A B H H F       (6.56) Then the dynamic milling equation becomes,    , , , , xyz,,(z)1 ( ) ( ) ( ) (s)jcd pd d0 j eq j k j W j T j0 j k kk k kkse c s s s d dz                         I A B H H F 0    (6.57) where I  is 3-by-3 identity matrix. Summing the contributions of all kN  axial disks and fN  cutting flutes gives, 139    ,, , , , , xyz1 11 ( ) ( ) ( )fkj kNNcd pd d0 j eq j k 0 j j W j Tk kk kk jse c s s s dz                          I A B H H F 0  (6.58) Eq.(6.58) has nontrivial solution only if,    ,, , , , ,1 1det 1 ( ) ( ) ( ) 0fkj kNNcd pd0 j eq j k 0 j j W j Tk kk kk jse c s s s dz                          I A B H H  (6.59) Eq.(6.59) yields the following characteristic equation for the milling system when s is replaced with i  at the border of stability,   ,, , , , ,1 1( ) det 1 ( ) ( ) ( )fkj kNNcd pd0 j eq j k 0 j j W j Tk kk kk jii e c i i i dz                             Ω I A B H H  (6.60) where i is the unit imaginary number.  Eq.(6.60) considers variations of the cutting force coefficients, engagement conditions, cutting edge geometry, and also includes the effects of varying workpiece FRFs and process damping conditions both in tool axis direction and along the five-axis toolpath. Using the Nyquist criterion [117], chatter stability is assessed (stable/unstable) by writing Eq.(6.60) at discrete locations along the toolpath and scanning the chatter frequencies ( )c  around the dominant natural frequencies of the tool and workpiece. Eq.(6.60) can also be used to obtain the stability lobes (productivity charts) in ball-end milling, flat-end milling, and orthogonal turning operations by checking the stability of the system at spindle speed and axial depth of cut pairs.    Although Eq.(6.60) is defined in TCS, , ( )j W ki  H  is defined in WCS in Chapters 3 & 4,  and , ( )j T ki  H  is measured in the fixed machine coordinate system (MCS). Hence, both FRF 140  components need to be transformed into TCS while assessing the chatter stability at any discrete point along the toolpath.  As a convention, MACHpro defines [4] ( )Tx  axis in the plane formed by ( )Tz  axis and the linear feed vector ( )LF  at th  map point ( )MAP  in WCS while projecting TWEs on TCS, and the same convention is followed in this thesis. The unit vectors defining TCS can be obtained in WCS as [4],  1 , ( ) ( ) , ( ) ( ) ( )LL T T T T TL     FF P P y z x y zF  (6.61) where ( )Tz  is the unit orientation vector of the tool axis, ( )  is cross product of vectors, P  and 1P  are the position vectors of the tool tip at MAP  and 1MAP  , respectively. ( )Tz , P , and 1P  are available from the five-axis toolpath. Then, the coordinate transformation matrix from TCS to WCS at MAP  is,  , ,x ,, , ,, , ,( ) ( ) ( )( ) ( ) ( )( ) ( ) ( )T x T T xTCSWCS T y T y T yT z T z T z      x y zT x y zx y z  (6.62)       A similar transformation can be carried out between MCS and WCS depending on the kinematic configuration of the five-axis machine. The machine used in the verification experiments in Section 6.3 has two rotary motions ( , )A C   around Mx  and Mz  axes on the table side. For this configuration, MCSWCST  can be obtained as [67],  cos sin 0 1 0 0sin cos 0 0 cos sin0 sin cos0 0 1C CMCSWCS C C A AA A                    T   (6.63) 141  where A  and C  are either directly obtained from the post processed toolpath file, or through the inverse kinematics model of the machine configuration.   The updated FRFs of the thin-walled workpiece obtained in WCS using the methods developed in Chapters 3 & 4 are transformed into TCS at MAP  as,   , ,( ) ( )TTCS WCS WCS WCSj W TCS j W TCSk ki i       H T H T   (6.64) where 1WCS TCSTCS WCST T . Similarly, the measured tool FRFs in MCS can be transformed into TCS at MAP ,  ,T ,T( ) ( )TTCS MCS WCS MCS MCS WCSj WCS TCS j WCS TCSk ki i               H T T H T T   (6.65) The overall flow chart for assessing the chatter stability in five-axis ball-end milling of thin-walled parts is summarized in Figure 6.5. The stability is checked at discrete stations ( )MAP  along the toolpath for the given tool geometry, cutting conditions, tool tip position and tool orientation obtained from the NC toolpath. The cutting tool is axially discretized (dz) to consider the variation in the local ball-end tool geometry and its effect on the chip thickness (Eqs.(6.16) -(6.25)), directional coefficients (Eqs.(6.40)-(6.46)), and in turn the dynamic milling forces (Eq.(6.52)). The varying thin-walled workpiece dynamics update models developed in Chapters 3 & 4 are employed to estimate the workpiece FRFs, and the proposed analytical model in Chapter 5 is used to approximate the process induced damping forces. Both models are integrated into the characteristic equation (Eqs.(6.52)-(6.65)) calculations. The chatter stability is checked through tracking the characteristic matrix equation’s (Eq.(6.60)) encirclements the origin of the complex plane.          142   Figure 6.5: Flow chart to assess the chatter stability in five-axis ball-end milling of thin-walled parts 6.3 Experimental Verification Experimental verifications presented in this section have two purposes: (i) further validation of the proposed analytical process damping model in orthogonal turning and flat-end milling of rigid workpieces by comparing the predicted chatter stability diagrams (lobes) with the experimental results, (ii) verification of the proposed chatter stability prediction model in ball-end milling of a thin-walled rectangular workpiece and a five-axis ball-end milling of a complex impeller blade with and without the process damping effect. Also, the developed models’ ability to evaluate the static deflection errors and part vibrations are demonstrated on two thin-walled workpiece machining cases.  143  6.3.1 Stability Diagrams with Process Damping in Orthogonal Turning The governing equation of chatter stability with process damping in orthogonal cutting of a rigid workpiece can be deduced from the general characteristic equation given in Eq.(6.60). The stable width of cut ( )pb  limit is sought while the chip regeneration takes place in feed (radial) direction only (one dimensional cutting). Also, cutting forces do not change their directions in orthogonal cutting, and hence the directional coefficients for the dynamic cutting and process damping forces vanish. Then, Eq.(6.60) reduces to,    1 (1 ) ( ) 0cr eq c T c piK e c i H i b         (6.66) In Figure 6.6, stability lobes prediction using the analytically calculated equivalent specific process damping ( )eqc  in Section 5.6.1 is compared with the experimentally identified process damping force coefficient ( )iC  based prediction of [48]. The experimentally stable and unstable points are indicated with circle and cross markers, respectively. Stability lobes predicted with the experimental iC  misses the asymptotic stability line at low spindle speeds while the upper and lower stability lobes predicted with eqc  for the vibration frequencies of c =120 Hz and c =20 Hz lead to better estimations. There is still some discrepancy between the predicted lobes with eqc  and the experimental stable and unstable points. In addition to the reasons discussed in Section 5.6, it could be due to the vibration frequency ( [20,120]Hz)c   used in the simulations since the tool vibrates at the chatter frequency that is close to the dominant mode frequency ( 450.7 Hz)n   of the structure at the border of stability. As the vibration amplitude at the stability boundary is not known in [48], eqc  is not calculated at the chatter frequency.  144   Figure 6.6: Comparison of the stability lobes predicted using equivalent specific process damping (ceq) and the experimentally identified process damping force coefficient (Ci) for the parameters in Table 5.1. Work material: AISI 1045;Structural parameters: ɷn =450.7 Hz, mr =0.8081 kg,  kr =6.48x106 N/m,  cr =145 Ns/m; Cutting force coefficient: Kr =1384 MPa Similarly, stability lobes are predicted using the equivalent specific process damping ( )eqc  from Section 5.6.1, and they are compared with the experimentally identified process damping coefficient ,( )p rc  based predictions of [54] in Figure 6.7. There is some discrepancy between their stability lobe predictions in Figure 6.7. The vibration amplitude at the border of stability is different at different cutting/spindle speeds due to the non-linear relation between the process damping and stable cutting width/depth limit [47, 50]. However, the vibration data at the stability boundary is only available at few experimental points in Figure 6.7, and thus the analytical eqc  is evaluated for the average vibration amplitude of 0 3.3A m  as given in Table 5.3.  145   Figure 6.7: Comparison of the stability lobes predicted using equivalent specific process damping (ceq) and the experimentally identified process damping coefficient (cp,r) for the parameters in Table 5.3. Work material: AISI 1018; Structural parameters: ɷn =1403 Hz, kr =21x106 N/m, ξr=2.2%; Cutting force coefficient: Kr =1375 MPa Furthermore, the inner and outer diameters of the tubular workpiece used in [54] are given but the effective diameter of the workpiece at different speeds is not known. Hence, the same average tube diameter is used for the entire spindle speed range both in eqc  and stability lobe simulations. Also, the stability simulation in [54] is based on the correlation between the indented volume ( )dV  and the indentation force coefficient ( )spK  as given in Eq.(2.1), which is not the case in the proposed chatter stability model (Section 6.2). Overall, it can be said that the proposed analytical process damping model yields the chatter stability chart with an acceptable accuracy when all these differences and the involved uncertainties (Section 5.6.2) are considered.    146  6.3.2 Stability Diagrams with Process Damping in Flat-End Milling Eq.(6.60) can be simplified to govern the stability of two dimensional flat-end milling operations having constant tool-workpiece engagement (TWE) boundaries, cutter geometry, cutting and process damping conditions along the tool axis. Further simplifications can be achieved for the flat-end milling of a rigid workpiece with the regular pitch ( 2 )p fN   tool. When the tool, tool holder, and spindle assembly is rigid in the tool axis (z) direction, the stable axial depth of cut ( )pa  limit is sought due to the chip regeneration contributed by the vibrations of the tool in Tx - Ty  plane. Therefore, the characteristic equation of the milling system (Eq.(6.60)) can be reduced to the following form,       det (1 ) ( ) 0ccd pd0 eq c T c p0ie c i i a      I A B H   (6.67) where I , ,cd0 jA , ,pd0 jB , TH  become 2-by-2 matrices. Experimental stability lobe results given for the flat-end milling of a rigid workpiece made of AISI 1018 in [54] are used to validate the lobe predictions with the proposed analytical process damping model. Unlike the orthogonal cutting, vibrations at the chatter frequency on the stability border are either directly measured using non-contact vibrometers [118], or are estimated [54] from the measured cutting forces at the marginal stability in milling. The measured/estimated resultant vibrations in the Tx - Ty  plane need to be projected on the radial direction of the engaged cutting flutes. Since the immersion of each flute is not known when the cutting starts, flute #1 is assigned as the reference flute, and its angular position in TCS is set to 00 at the beginning of the cut. The angular positions of the other flutes are calculated using the pitch angle (Eq.(6.11)) for the full immersion conditions given in [54]. The cutting edge geometry, vibration 147  and friction parameters given in [54] are presented in Table 6.1. The cutter has a 19 mm diameter and 300 helix angle, and its measured modal parameters along with the estimated simulation parameters are given in Table 6.2. Tool vibrations (displacements) are projected on the radial direction of the engaged flutes and filtered around the natural frequencies to estimate the vibration amplitude at the border of stability as shown in Figure 6.8, which is used in the equivalent specific process damping ( )eqc  calculations. Simulated variation of eqc  for the given spindle speed range and the average vibration amplitude is shown in Figure 6.9.   Table 6.1: Cutting edge geometry, vibration and friction parameters for the flat-end milling tests in [54] Material Edge Radius, ( )eR m  Rake Angle, 0( )n  Clearance Angle, 0( )n  Spindle Speed,  sn  (rpm)  Vibration Amplitude, 0A  ( )m  Vibration Frequency,c  (Hz) Friction Coefficient,   AISI 1018 15 5 8 500-2750  ~9  ~765 0.72  Table 6.2: Modal parameters of the flat-end cutter measured at the tool tip in feed and normal directions in [54], and the estimated simulation parameters Natural Frequency,  n  (Hz) Modal Stiffness,  rk   (N/m) Modal Damping,    (%) Average Separation Angle, 0( )  Estimated Workpiece Temperature (0C) Feed (x): 748 Normal (y): 729 Feed (x): 98106 Normal (y): 124106 Feed (x): 1.9 Normal (y): 1.9 48 234  148   Figure 6.8: a) Radial displacements (vibrations) of the three fluted cutter at marginally stable case (ns=1000 rpm, ap=3 mm) given in [54], and b) FFT of the vibrations  Figure 6.9: Simulated variation of the equivalent specific process damping (ceq) in flat-end milling of AISI 1018 for the conditions given in Table 6.1 The stability diagram for full immersion flat-end milling of AISI 1018 workpiece is predicted using the analytically estimated eqc , and it is compared with the experimental calibration based prediction of [54] in Figure 6.10. The simulation result without process damping is not given for clarity. Prediction with the proposed eqc  has the same level of accuracy as the empirical 149  identification based prediction. Since there is only one identified marginally stable experimental point, its vibration amplitude is used for the entire spindle speed range in the simulations.   Figure 6.10: Comparison of the stability lobes predicted using the equivalent specific process damping (ceq) and the experimentally identified process damping coefficient (cp,r) for the parameters in Table 6.1. Work material: AISI 1018; Cutter: three fluted flat-end tool with 19 mm diameter, 300 helix angle; Immersion: 100% radial; Cutting force coefficients:ctr K [1550 480] MPa ;Structural parameters: given in Table 6.2 The stability lobe predictions with the proposed process damping model is also validated with two representative experimental flat-end milling cases for rigid Al 7075 and AISI 1050 workpieces given in [118]. Edge geometry, vibration and friction parameters are given in Table 6.3, and the measured modal parameters of the cutting tool and the estimated simulation variables are given in Table 6.4. Unlike the previous example, different chatter vibration amplitudes are obtained from multiple marginally stable points identified experimentally. This yields a range of vibration amplitudes, which is exemplified for Al 7075 in Figure 6.11(a-b).      150  Table 6.3: Cutting edge geometry, vibration and friction parameters for the flat-end milling tests in [118] Material Edge Radius, ( )eR m  Rake Angle, 0( )n  Clearance Angle, 0( )n  Spindle Speed,  sn  (rpm)  Vibration Amplitude, 0A  ( )m  Vibration Frequency,c  (Hz) Friction Coefficient,   Al7075 15 5 10 1750-4500  3.6-7.5  ~3200 0.7 AISI 1050 1500-3250  3.6-9   Table 6.4: Modal parameters of the flat-end cutter measured at the tool tip in feed and normal directions in [118], and the estimated simulation parameters Natural Frequency,  n  (Hz) Modal Stiffness,  rk   (N/m) Modal Damping,    (%) Average Separation Angle [105, 119], 0( )  Estimated Workpiece Temperature (0C) Feed (x): 3110 Normal (y): 3106 Feed (x): 9.73106 Normal (y): 10.35106 Feed (x): 1.43 Normal (y): 1.68 45  174 Feed (x): 3210 Normal (y): 3217 Feed (x): 9.94106 Normal (y): 9.06106 Feed (x): 1.15 Normal (y): 1.45 50  191   Figure 6.11: a) Representative radial displacements (vibrations) of the cutter at two different marginally stable cases (ns=2071 rpm ap=0.75 mm, and ns=2519 rpm ap=0.35 mm) for Al 7075 given in [118], and b) FFT of the vibrations 151   Figure 6.12: Simulated variation of the equivalent specific process damping (ceq) in flat-end milling of a) Al 7075 and b) AISI 1050 for the conditions in Table 6.3 Higher vibration amplitudes are obtained towards lower spindle speeds, while the lower amplitudes are observed at relatively higher spindle speeds in the range of interest. Therefore, the equivalent specific process damping ( )eqc  is calculated for the maximum, minimum, and average amplitudes as shown in Figure 6.12(a-b). The corresponding predicted stability lobes along with the experimentally observed points are presented in Figure 6.13(a-b). Lobes are better represented with higher vibration amplitudes at relatively low spindle speeds, while the lower amplitudes give better predictions at the higher spindle speeds. In order to estimate the stability lobes more accurately for the entire spindle speed range, vibration amplitudes at the stability border need to be either experimentally known, or calculated with the iterative time domain simulation of tool vibrations, which is beyond the scope of the current study and left as a future research topic.  152   Figure 6.13 : Comparison of the predicted stability lobes using equivalent specific process damping (ceq)  at different vibration amplitudes for four fluted flat-end mill with 12 mm diameter, 300 helix angle with the parameters in Table 6.3 and Table 6.4, a) Work material: Al 7075; Cutting force coefficients: ctr K  [1160 440] MPa for 83% radial immersion down milling; b) Work material: AISI 1050; Cutting force coefficients: ctr K  [1800 350] MPa 50% radial immersion down milling Considering all the uncertainties in the model and use of the average vibration amplitude for the entire spindle speed range of interest, it can be said that the proposed analytical process 153  damping model can avoid the need for costly experimental identifications and guide the process planners to select cutting conditions and tool geometry for the improved stability. 6.3.3 Thin-Walled Rectangular Workpiece Chatter Test The proposed chatter stability prediction model for ball-end milling of flexible structures has been first validated in single axis machining of a thin-walled rectangular workpiece similar to the one used in Section 3.7. For the model verification purpose, the toolpath is composed only of two passes and the cutting conditions are given in Table 6.5. Machining is conducted at a relatively high spindle speed to avoid process damping. The workpiece is made of Al7050 and dimensions of the unmachined part are shown in Figure 6.14(a). Experiments have been conducted on a high speed machining center by clamping the thin-walled workpiece on a vise that is rigidly mounted on the machine table (Figure 6.14(b)). The stability of the cut is assessed from the recorded cutting sound with a microphone and by inspecting the machined surface quality. Although the path is short, variation in the dominant mode frequencies is predicted using the model developed in Chapter 4 by discretizing the initial workpiece geometry as detailed in Section 3.7, and the results are presented in Figure 6.15.  Table 6.5: Cutting conditions in rectangular thin-walled workpiece chatter tests  Spindle Speed (rpm) Feed (mm/rev/tooth) Radial Depth of cut (mm) Cutting Tool 1st Pass 5250 0.1 2 2 fluted solid ball-end mill with 16 mm shank diameter 2nd Pass  154   Figure 6.14: a) Initial dimensions of the workpiece, and WCS and TCS, b) experimental setup  Figure 6.15: Predicted variation of the first two dominant mode frequencies of the workpiece  Chatter stability of the cut is evaluated by scanning the frequencies around the predicted dominant natural frequencies (Figure 6.15) of the part at discrete tool positions along the toolpath using Eq.(6.60) without the process damping terms. The model damping ratios identified for a very similar geometry in Section 3.7 are used with the predicted modes to calculate the varying workpiece FRF in the most flexible direction of the workpiece ( Wy ). The predicted and experimentally observed unstable (chatter) regions are compared in Figure 6.16 for 155  both passes. The cutting is unstable along the entire toolpath both in the simulation and experiment.      Figure 6.16: Predicted and experimental unstable (chatter) regions, and the experimental cutting sound In order to compare the predicted and measured chatter frequencies, a time-frequency analysis is carried out on the collected cutting sound data by taking its FFT every 0.1 seconds along cutting. The comparison results are given in Figure 6.17, where the dominant frequencies in the cutting sound are shown in Figure 6.17(a) while the predicted chatter frequencies are given in Figure 6.17(b). As seen, there is a good match between the simulated and measured chatter frequencies as well. Since the second dominant mode (first torsional mode) of the workpiece diminishes around half way through each pass, it does not contribute to chatter in the middle of both passes, which is also captured in Figure 6.17(b) as the predicted workpiece FRF automatically includes this dynamics information.  156   Figure 6.17: a) Experimental, and b) predicted chatter frequencies along the toolpath 6.3.4 Impeller Blade Cutting Test The proposed thin-walled workpiece dynamics update and chatter stability models with and without the process damping effect are experimentally verified in five-axis ball-end milling of an impeller having six evenly spaced identical twisted blades. The workpiece is made of Al 7050, and its height, bottom and top diameters are 115 mm, 145 mm, and 44 mm, respectively. The nominal dimensions of a single blade are shown in Figure 6.18. Five-axis ball-end milling tests have been conducted on a machining center by clamping the workpiece on a three-jaw fixture that is fastened on the machine table as shown in Figure 6.19(a). 157   Figure 6.18: CAD model of the impeller used in the verification experiment, and the nominal finished dimensions of the blade top and bottom  Figure 6.19: a) Experimental setup, and b) the machined part The nominal blade thickness after the opening pass, which includes the blade roughing, is 2.2 mm. It then goes through the semi-finishing and finishing operations with 0.24 mm and 0.11 mm stock thicknesses, respectively. Blade geometry is split into “top” and “bottom” portions at 50% length of the leading edge (excluding the root blend) as illustrated in Figure 6.18 to show the effect of the machining sequence on the blade dynamics by following two different machining strategies: (i) semi-finishing the blade top and bottom, and then finishing the blade top and bottom; (ii) semi-finishing and finishing the blade top, and then semi-finishing and finishing the 158  blade bottom. Also, in order to reveal the effect of process damping on the five-axis ball-end milling dynamics, semi-finishing the blade top in these sequences is conducted at 8500 rev/min and 1500 rev/min, respectively. The test at 8500 rpm is used to validate the chatter stability prediction model for complex five-axis toolpaths with varying dynamics, while the test with 1500 rpm is used to further validate the chatter stability prediction with process damping. The finished impeller is shown in Figure 6.19(b), and the cutting conditions for semi-finishing the blade top are given in Table 6.6. Table 6.6: Cutting conditions in semi-finishing the blade top  Spindle Speed (rpm) Feed (mm/rev/tooth) Stock Thickness (mm) Cutting Tool Strategy 1 8500 0.055-0.16 0.24 2 fluted solid ball-end mill with 16 mm shank diameter Strategy 2 1500  Variation of the first two (first bending and first torsional) dominant mode frequencies in two strategies are presented in Figure 6.20, where experimental measurements are done only for Strategy 2. As seen, the proposed in-process workpiece dynamics update model accurately predicts the mode variations. Prediction errors in the mode frequencies are respectively ~6.8% and 3.5% for the first two dominant modes, which result from the FE model of the initial blade, and they are preserved throughout machining the blade top. The first operation (semi-finishing the blade top) is common in both sequences, hence both strategies yield the same mode variations. However, next operations in the above strategies, i.e. semi-finishing the blade bottom and finishing the blade top, result in completely different behaviors, which clearly show the significance of the machining sequence for in-process blade dynamics.  159   Figure 6.20: Predicted variation of the first two dominant mode frequencies of the blade, and experimental modes in Strategy 2 Chatter stability of semi-finishing the blade top with 8500 rpm is simulated to predict the unstable toolpath segments and the associated chatter frequencies. Figure 6.21(a-b) compares the predicted chatter (unstable) regions on the impeller blade with the experimentally observed surface quality. As seen, the proposed five-axis chatter stability prediction model accurately determines the unstable toolpath segments, and chatter predominantly occurs towards the free end of the leading edge as it is the most flexible section of the blade.   In order to compare the predicted and measured chatter frequencies, time-frequency analysis is carried out on the collected cutting sound data similar to the rectangular thin-walled part machining case. Comparison results are given in Figure 6.22, where the measured dominant frequencies in the cutting sound are shown in Figure 6.22(a) while the predicted chatter frequencies are given at discrete toolpath points in Figure 6.22(b). The details of the cutting 160  sound for a sample toolpath segment are also shown in Figure 6.22(c) with the FFT spectrum. This segment partially covers the front face, leading edge and the back face machining. The first mode of the blade is responsible for the dominant chatter frequency 1( )c  along the entire semi-finishing operation.    Figure 6.21: Predicted unstable (chatter) regions on a) front and back faces, and b) and experimentally observed chatter marks in semi-finishing the blade top (front face) at 8500 rev/min 161  Although there is a good match between the variation trends of the simulated and measured chatter frequencies in the whole operation, there are some discrepancies in the frequency value, which is due to the error (~6.8%) in the FE model of the initial blade. Chatter frequency prediction can be improved by refining the initial FE mesh.  Figure 6.22: a) Experimentally identified chatter frequencies, b) predicted chatter frequencies, c) sound data of a sample toolpath segment and its FFT in semi-finishing the blade top at 8500 rev/min Lastly, semi-finishing the blade top test is repeated at 1500 rev/min to further validate the introduced process damping model on the impeller blade with the parameters given in Table 6.7.   162  Table 6.7: Cutting edge geometry, vibration, friction, and simulation parameters for process damping calculations in semi-finishing the blade top Material Edge Radius, eR( )m  Normal Rake Angle, n  0( )  Normal Clearance Angle, n  0( )  Spindle Speed,  sn   (rpm)  Vibration Amp.,  0A   ( )m  Vibration Freq., c  (Hz) Average Sep. Angle, 0( )  Friction Coefficient,   Al7050 49 Varying 5 1500  Varying Varying 45 0.3  In Table 6.7, edge radii of the flutes are measured at three different points and the average value is used for all flutes in the simulations. The normal rake angle of the ball-end cutter is calculated from the radial rake angle at the cylindrical part (013.5 ) of the tool using Eq.(6.6). The vibration (chatter) frequency is already known from Figure 6.22 as ~5kHz, the same average separation angle of 045   given for Al7075 in Table 6.4 is used for Al7050, and a general average friction coefficient [48] is employed. Since the vibration amplitude is not known a priori and cannot be calculated in the frequency domain analysis, the critical vibration amplitude given in [120] before the tool jumps out of the cut due to vibrations is used to obtain a representative amplitude,  0, ,(z) (z) 2sA h   (6.68) where ,sh  is the  average static chip thickness at th  discrete point along the toolpath. Since each point has a varying static chip thickness ( )sh  both in tangential and axial directions of the tool (Section 6.2.2), the vibration amplitude used at each discretized cutting edge element is different. The workpiece temperature, whose prediction along complex five-axis toolpaths is still an active research area, is not considered. This is a valid assumption since aluminum has a high thermal diffusivity, both sides of the blade are open, and the removed material thickness in semi-163  finishing the blade top is small. Chatter stability prediction of the operation is repeated by taking the varying process damping conditions both in tool axis and along the toolpath into account with the varying blade dynamics along cutting. The predicted chatter regions and the corresponding chatter frequencies are validated with the collected sound data and its FFT analysis as follows. Figure 6.23 shows the predicted unstable (chatter) regions on the blade and the corresponding predicted chatter frequencies when the effect of process damping is taken into account. In the presence of process damping, extend of the chatter region on the front face is remarkably decreased as compared to the one in Figure 6.21.  Figure 6.23: a) Predicted chatter regions and b) the corresponding chatter frequencies in semi-finishing the blade top with the process damping effect at 1500 rpm The experimentally identified chatter frequencies along the entire semi-finishing of the blade top is obtained from FFT analysis of the recorded cutting sound and presented in Figure 6.24(a). Chatter is still present at certain segments of the toolpath, but with diminished amplitude. The same toolpath segment and its cutting sound data as Figure 6.22(a-b) are analyzed in detail in 164  Figure 6.24(b). The chatter frequency 1( )c  and the tooth passing frequency ( )T  equally dominate (marginally stable) the power spectrum of the cutting sound while it was completely dominated by 1c  in Figure 6.22(c). However, chatter is still present around the leading edge and on the back face of the blade with the diminished effect. These changes conclude that the predicted stable zones on the front face in Figure 6.23(a) are in a good agreement with the experimental results given in Figure 6.24(c). Since the spindle speed is the only parameter changed, both predicted (Figure 6.23(b)) and experimentally identified (Figure 6.24(a)) chatter frequencies stay similar to those in Figure 6.22(a-b) both in trend and magnitude.         Figure 6.24: a) Experimentally identified chatter frequencies, b) sound data of a sample toolpath segment, and c) its FFT in semi-finishing the blade top at 1500 rpm  165  6.3.5 Application: Prediction of Static Form Errors in Ball-End Milling of Blades In the machining of thermally resistant metals at low cutting speeds, the accuracy of thin-walled parts is not dominated by the chatter but by the combined static deflections of the flexible workpiece and slender ball-end mill [3, 13, 18, 62, 121] as illustrated in Figure 6.25. Instead of obtaining the static stiffness of the thin-walled workpiece from the predicted FRFs as shown in Chapter 3, Eq.(4.5) derived in Chapter 4 can be directly used to estimate the static stiffness of the thin part in the cutting zone as the material is removed.  In ball-end milling, the designed part surface ( )dS  is generated at the cutter-workpiece contact (CC) point, and the dimensional form error is evaluated in the surface normal direction while the flute generates that point. The position of CC at a discrete point along the path ( )ccP  is calculated in the workpiece coordinate system (WCS) as the closest point on the designed surface to the ball center,  0 minWCScc T dR  P P z S   (6.69) where P  is the position vector of the tool tip, Tz  is the unit tool orientation vector in WCS, 0R  is the ball-end mill radius at its cylinderical section. The surface normal is evaluated from WCSccP  to the ball center 0( )TRP z .   166   Figure 6.25: Static deflections of the tool and flexible blade, and a sample cutter contact (CC) point The static cutting forces in the radial, tangential, and axial directions of the local cutting edge are predicted using the author’s previously developed five-axis mechanics model in his masters thesis [4]. The engagement region (Figure 6.25) is divided into discrete elements along the tool axis in the tool coordinate system (TCS) as,  , s, ,1(t, z) ( (t, z)) ( (t, z) (t, z) (z) (t, z) (z))fNs c prta j rta j j rta jjd h db dS  F K K   (6.70) where fN  is the number of flutes, srtadF  is the differential static cutting forces; ,crta jK  and ,prta jK  are the cutting and edge (static plowing) force coefficients, respectively; s, jh  is the local static chip thickness. The radial compliance (flexibility) of the tool at axial elevation kz  due to the cutting force at mz  ,( )k mZ  is calculated using the cantilevered beam formulation given in [3, 13, 113, 121]. ,k mZ  is used to predict the tool deflections in the direction perpendicular to the tool axis and written as, 167   2, 2( ) (2 3 ) 1,6( ) (2 3 ) 1,6c k c m km kck mc m c k mk mcL z L z zz zEI kZL z L z zz zEI k          (6.71) where cL  is the gauge length of the tool measured from its tip to the collet or tool holder, ck  is the measured clamping stiffness of the holder, E is the modulus of elasticity of the tool material, I is the area moment of inertia for the effective local tool radius ( (z))effR , which is approximated as (z) 0.8 (z)effR R  due to the cutting flutes [13, 113]. Elemental forces are projected into the xyn  (Figure 6.25) direction at each disk element and multiplied with the corresponding compliance to evaluate the components of the tool deflection vector ( )Td  as,  , , ,1(( ) )kNrta sT k k m xyz rta m xymZ dd T F n   (6.72) where kN  is the total number of disk elements, rtaxyzT  is the transformation matrix from the local cutting edge coordinate system to tool’s Cartesian coordinate system (see Section 6.2.1). The portion of Td  belonging to the axial disk containing the CC point ( Td ) is projected into the surface normal direction to obtain the deflection error due to the tool ( )eTd  as,   0, 0 0TCS TCSTe TCS TCS WCS WCSct ccT T ct cc TCS ccTCS TCSct ccd R  P Pd P P T PP P  (6.73) where TCSctP  is the position vector of the ball-end tool center in TCS, WCSTCST  is the transformation matrix from the workpiece coordinate frame to the tool frame located at the tool tip (see Section 6.2.5). 168  In order to calculate the contribution of the workpiece deflections on the surface error, the predicted elemental forces are transformed from TCS to WCS as,  ( )s WCS TCS sxyz WCS xyzd d F T F   (6.74) where the negative sign shows the direction of the cutting forces on the workpiece. Mesh elements containing the cutting forces (Eq.(6.74)) are identified, and the elemental forces are distributed over their nodes using their Barycentric coordinates [3, 121] as weighting. The nodal displacements vector for the workpiece ( )Wd  is evaluated using the updated stiffness matrix ( )iBK  from Eq.(4.5) and the cutting forces (Eq.(6.74)) as,  ( )iB s WCSW xyzK d F   (6.75) where sxyzF  is the total forcing vector for the workpiece. Deflection of the closest node ( )Wd  to the CC point is projected into the surface normal direction to evaluate the surface error due to the workpiece,  00( )( )WCSe T ccW W WCST ccRdR  P z PdP z P  (6.76) Then the static form error left on the surface at th  discrete station along the toolpath is obtained by superposing the contributions of the workpiece and tool,  e e eT Wd d d    (6.77) TWE is updated at each axial disk considering the combined deflection of the workpiece and tool, and the surface error is recalculated until the error between two consecutive iterative calculations converges. The engagement boundaries are updated in down milling as [3],    169   W* 1W* 1WW* 1(z) (z) (z)(1 cos (z));(z) cos ( cos ( (z) (z)) (z))(z) cos ( cos ( (z) (z)) (z))(z)(1 cos (z)) (z) (z) (z)(1 cos (z));(z) cos ( cos ( (zT exst st Tex ex Tst T exst st Tif RRRif R R                          d dd dd dd dd W*) (z)) (z))(z)exR   d  (6.78) where superscript * denotes the updated engagement boundaries. The static form error prediction is implemented in three-axis ball-end milling of a blade made of Al7050. The workpiece is cut at 1000 rpm with 0.1 mm chip load under stable cutting conditions. A four-fluted ball-end mill with 12 mm diameter is used in the tests. The blade height and width are 50 mm and 54 mm, respectively. The finished blade thickness varies between 4 mm and 6 mm along its longitudinal direction, and the removed stock thickness is uniform at 1 mm. The form errors left on the machined surface are measured under the coordinate measurement machine (CMM), and they are compared with the predicted values for the first three passes of the toolpath in Figure 6.26. As seen, the form errors in the most flexible sections of the toolpath are predicted with 15% and 22% average prediction errors on the front and back sides of the blade. Prediction error could be attributed to error in the estimated part stiffness, the cutting force coefficients obtained from orthogonal to oblique transformations, and the error in TWE calculations. Different errors on the front and back faces result from the spiral toolpath. The predicted values can be used to compensate for the errors in the toolpath generation stage.  170   Figure 6.26: Comparison of the measured and predicted dimensional form errors 6.3.6 Application: Prediction of Thin-Walled Rectangular Workpiece Vibrations The predicted in-process workpiece dynamics ( ( ))uvH   in Chapter 4 can be used to predict the forced vibrations of thin-walled parts at stable spindle speeds. Prediction of the forced vibrations is demonstrated on the machining of a rectangular thin-walled workpiece presented by Eksioglu [6]. The cutting forces are predicted using the mechanics model given in [4], and the forced vibrations of the part are evaluated using the predicted in-process workpiece FRF as follows. The cutting positions from the free end of the rectangular part (Figure 6.27(a)), stable cutting conditions, and the cutter geometry are given in Table 6.8. Both tests are single pass toolpaths, 171  and the measured workpiece FRF before machining is taken from [6], while the workpiece FRF is predicted as the material is removed using the proposed model in Chapter 4 with the given modal damping ratios in [6]. Due to small material removal, the part dynamics are updated only for Test 1 and used in simulations of both tests. Comparison of the measured and predicted direct FRFs of the workpiece in its normal direction ( )Wx  is presented in Figure 6.27(b), where two dominant vibration frequencies are the first bending and torsional modes, respectively. Since the workpiece is found to be much more flexible than the tool [6], only the part dynamics is considered.      Figure 6.27: a) Workpiece geometry, b) comparison of the predicted and measured workpiece FRFs in normal (x )W  direction. Measured modal parameters in [6]: natural frequency  , 731.2 1718.6n r  Hz, modal stiffness   62.04 8.67 10rk    N/m, modal damping ratio  0.01 0.0022r  , modal mass  0.097 0.074rm  kg     172   Table 6.8: Cutting conditions used in the stable cutting tests in [6]  Tool Tip Position (mm) Axial Depth,  pa  (mm) Spindle Speed, n (rev/min) Feed (mm/tooth) Radial  Depth,  pb (mm)  Tool Test 1 0.5 mm from the free end 0.5 10300 0.1 2 4-fluted flat-end mill with  19.05 mm diameter,  350 helix angle Test 2 2 mm from the free end 2 11000  In order to be consistent with [6], and due to the low axial depth of cut, the predicted cutting forces are lumped at the tool tip in TCS (( ) )s TCSxyzF  and transformed into WCS using TCSWCST  (Eq.(6.62)). The workpiece dynamics are represented at the corresponding workpiece node at point A (Figure 6.27(a)), which has Wy =21.4 mm. Dynamic displacements of the workpiece are expressed in the Laplace domain as,  , , , , ,y, ,, , ,,z,, , ,(s) (s) (s) (s) (s)(s) (s) (s) (s) (s)(s)( )(s) (s) (s)i i ii i ii i iB B Bxx W xy W xz W x W d WB B BW d Wyx W yy W yz WB B Bd WW Azx W zy W zz W AAH H H F xH H H F yzF sH H H                           (6.79) where (s)iBWH  is the predicted workpiece FRF and can be rewritten from Eq.(3.25) as,   uv, 2u, v,21 , ,(s) , u, v = x , , z2OiBW W WWNr rr r n r n rH yss i       (6.80) Discrete time displacements of the part in its normal direction ( )Wx  are evaluated by transforming the predicted workpiece FRFs from the Laplace domain to the discrete time domain using the trapezoidal approximation [113] as, 173   2 1 11 2 1 0uv, 2 1 12 1 012(1 )( ) ,(1 )OiBWNrz z zH z sa z a z a t zb b b            (6.81) where t  is the sampling time interval. Eq.(6.79) can be rewritten by sampling the simulated cutting forces at the same interval,  1 1 1, , , , ,1 1 1y, ,, , ,1 1 1,z,, , ,( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( )( )( )( ) ( ) ( )i i ii i ii i iB B Bxx W xy W xz W x W d WB B BW d Wyx W yy W yz WB B Bd WW Azx W zy W zz W AAH z H z H z F k t x k tH z H z H z F k t y k tz k tF k tH z H z H z                                    (6.82) where k is the sampling interval counter. Dynamic displacements of the workpiece are evaluated using Eq.(6.82) for the cases given in Table 6.8. Simulated cutting forces and the dynamic displacements of the flexible part for both test cases are shown in Figure 6.28(a-b) and Figure 6.29(a-b) along with their frequency spectrums.  Although the predicted first bending mode frequency 1n   744 Hz is excited in Test 1, amplitude of the forced vibration at Point A is ~50 m  at the forcing frequency of T  688 Hz, which is well correlated with the dominant frequency in the measured cutting sound (Figure 6.28(c)). In Test 2, however, T  733 Hz almost matches with the predicted first bending mode frequency 1n   744 Hz and causes the mode to resonate under heavy vibrations reaching to ~1 mm of vibration amplitude. The frequency spectrum of the collected sound data also confirms its resonance in Figure 6.29(c).  174   Figure 6.28: a) Simulated cutting force in Wx  and its FFT, b) simulated workpiece vibration in Wx and its FFT, c) experimental cutting sound data and its FFT around the point (A) analyzed at 10300 rpm The introduced models can be used to calculate the vibrations in the flexible milling system before any actual cutting is conducted. Vibration calculations demonstrated in this example can be taken one step further by simulating them using the previously developed semi-discrete time domain solution by Kilic [122] at UBC-MAL to evaluate both the forced and self-excited vibrations, and the surface quality of the part ahead of physical cutting trials.    175   Figure 6.29: a) Simulated cutting force in Wx  and its FFT, b) simulated workpiece vibration in Wx  and its FFT, c) experimental cutting sound data and its FFT around the point (A) analyzed at 11000 rpm 6.4 Summary This chapter presents a new, general dynamics and chatter stability model for the five-axis point milling of thin-walled structures considering the position dependent workpiece dynamics (Chapters 3 & 4) and varying process damping forces (Chapter 5) along the complex toolpaths. The cutting tool is axially discretized to accommodate the change both in the local cutting edge geometry and tool-workpiece engagements (TWEs). The dynamic cutting and process damping forces are evaluated for each axial cutting edge element considering the relative vibrations between the slender tool and flexible part. The overall dynamic milling equation of a discrete point along the toolpath is obtained by integrating individual contributions of all engaged cutting edge elements at that point. The chatter stability of the governing equation at these discrete 176  stations is checked using the Nyquist stability criterion. The proposed model is general and covers regular and variable pitch cutters, which are commonly used in industry.  Additional experimental validations for the introduced process damping model to accurately predict the stability charts (lobes) in orthogonal turning and flat-end milling cases are presented by deducing their dynamic equations from the general five-axis ball-end milling case. The presented general dynamics model is also experimentally verified in single-axis ball-end milling of a rectangular thin-walled workpiece and in five-axis ball-end milling of an impeller blade. It is shown that the proposed model detects the chatter regions and the corresponding chatter frequencies without going through any costly physical testing prior to the actual part machining. It is further shown how the developed models form a foundation for the digital milling of thin-walled aerospace structures by predicting the surface errors and workpiece vibrations as well. 177  Chapter 7: Conclusions and Future Research Directions 7.1 Summary and Contributions This thesis contributes to virtual, five-axis machining of thin-walled, highly flexible parts. The cutting forces cause static deflections and forced vibrations leading to tolerance violations; and chatter vibrations creating poor surface finish and tool breakage, which scrap the costly aerospace components.     Although the five-axis ball-end milling mechanics to predict cutting forces, workpiece and tool deflections have been studied in literature, they lack the consideration of varying workpiece characteristics due to the material removal along the toolpath. The position dependent part characteristics need to be included in predicting the chatter stability, static and dynamic deflections in a computationally efficient manner. Furthermore, process damping resulting from the contact between the cutting edge and wavy surface finish needs to be taken into account in low speed machining. The process damping forces greatly affect the chatter stability of the machining process and should be incorporated into the five-axis ball-end milling dynamics.  This thesis has developed the physics based and efficient process simulation models for five-axis ball-end milling of low rigidity parts to minimize the number of iterations in the process development stage. The process related issues can be identified ahead of time using the introduced predictive models having the following scholarly contributions:  A new, frequency domain dynamics substructuring based workpiece dynamics update model avoiding multiple FE models and remeshing of the in-process workpiece structure has been developed. The position dependent workpiece FRFs are analytically updated along the toolpath. FRFs of the part are constructed only once for the initial thin-walled part before machining, which is systematically discretized considering the radial depth of 178  cut in the operations. Since the model partitions the initial global structure, the removed materials are its own substructures with automatic nodal compatibility. Thus the model allows rigidly decoupling the dynamics of removed substructures from the workpiece FRFs at discrete toolpath locations. Efficiency of the model is improved by updating the in-process workpiece dynamics only at the interface degrees of freedom (DOF) of all material removals and also reducing the system matrices by selecting the most compliant DOF.     Alternative to the proposed reduced order frequency domain method, a new time domain dynamics update model for the in-process thin-walled workpiece has been introduced. The model is based on the dynamic substructuring and matrix perturbation to track the variations in modal frequencies and mode shapes of the thin-walled workpiece as it is machined. The model solves the eigenvalue problem (EVP) only once for the initial global workpiece and eliminates the need for repetitively solving EVPs for the reduced order and dense, or full order and ill-conditioned updated physical system matrices along the path. Unlike its frequency domain counterpart, the model potentially eliminates the need for model order reduction for better performance. The model can also be used to update the dynamics of flexible parts in turning and turn-milling operations.   The experimentally validated proposed dynamics update models in the frequency and time domains form the foundation for the analytical process damping and chatter stability models in thin-wall milling. A novel model has been proposed to analytically estimate the process damping forces in machining. As opposed to the existing experimentally calibrated empirical models, the proposed analytical model is general and considers the effect of the cutting edge geometry, cutting conditions, vibration parameters, and 179  workpiece material properties on the process damping. The non-linear specific contact force is represented with an equivalent damping for frequency domain analysis, while the calculated non-linear process damping forces can be directly used in time domain simulations. The model can be used to evaluate the machining dynamics behavior of the newly designed tools by generating the tool and workpiece specific stability charts without doing the costly machining tests.   Sensitivity of the proposed analytical process damping model to the uncertainties in the material separation angle, cutting edge radius, vibration amplitude, and the mean plastic yielding pressure threshold has been investigated for the given test cases. The model is found to be sensitive to the material separation angle and plastic yielding threshold for the analyzed cases.     Generalized dynamics and chatter stability prediction models for multi-axis ball-end milling of thin and rigid parts have been developed. The model considers both varying thin part dynamics and varying process damping forces at discrete stations along the five-axis toolpath. The model also considers different cutting, tool-workpiece engagement (TWE), and dynamics conditions at each discrete point. The chatter stability of any point on the path is checked and the chatter frequencies are identified for the unstable points. The model can simulate the process stability for regular and variable pitch cutters with arbitrary TWE boundaries along the tool axis either by lumping the blade and tool dynamics at a single point, or distributing them over the engagement region along the tool axis. The model can be used for variable pitch cutter design and spindle speed scheduling to avoid the chatter vibrations without resorting to costly physical machining trials. 180  The proposed models have been experimentally verified and prepared to be integrated into the digital machining of thin-walled aerospace parts and blades in Computer Aided Manufacturing (CAM) packages. The system will be capable of predicting the chatter stability, process forces, and the resulting static deflection and forced vibration errors left on the part surface to guide the process planners in achieving the most productive cutting conditions while avoiding part and tool damages. 7.2 Future Research Directions The fundamental mathematical models are presented in this thesis, and the following studies can improve the accuracy and efficiency of the proposed methods.     The calculation efficiency of the introduced frequency domain FRF update model is adversely affected by the number of decoupling frequencies in the simulation. Its computation time should be improved with efficient parallel programming schemes.   Though it has not been encountered in this study, the perturbation based dynamics parameters update model becomes inaccurate for closely spaced modes. Therefore, the model needs to be extended to close mode cases.   The machining sequence has a significant effect on how the dominant mode frequencies of the thin workpieces evolve as shown in Figure 6.20. Optimization of the material removal sequence and planning can be studied by considering the changes in the structural dynamics of the thin walled parts.  The proposed analytical process damping model yields acceptable results for the given test cases. However, its accuracy is directly related to the accuracy of the material separation angle, vibration amplitude, and the geometry dependent plastic contact pressure constant ( )pC . The separation angle should be identified for different cutting 181  edge geometries and workpiece materials from finite element (FE) simulations of the cutting with different vibration amplitudes. Finally, a comprehensive thermal model to predict the workpiece temperature should be integrated into the process damping model to obtain accurate material properties under the clearance face of the cutting edge.     Stiffness gain from the relatively rigid cutter at the most flexible sections of the thin part should be investigated. Although it is not observed in point milling tests in this study, it could cause shifts in the in-process dominant vibration mode frequencies of the parts with low rigidity in the peripheral/flank milling applications due to the long contact lengths between the tool and workpiece.  182  References [1] Tuysuz, O., and Altintas, Y., 2017, "Frequency Domain Updating of Thin-Walled Workpiece Dynamics Using Reduced Order Substructuring Method in Machining," Journal of Manufacturing Science and Engineering, 139(7), p. 071013. [2] Tuysuz, O., and Altintas, Y., 2017, "Time-Domain Modeling of Varying Dynamic Characteristics in Thin-Wall Machining Using Perturbation and Reduced-Order Substructuring Methods," Journal of Manufacturing Science and Engineering, 140(1), p. 011015. [3] Altintas, Y., Tuysuz, O., Habibi, M., and Li, Z. L., 2018, "Virtual compensation of deflection errors in ball end milling of flexible blades," CIRP Annals, 67(1), pp. 365-368. [4] Tuysuz, O., Altintas, Y., and Feng, H.-Y., 2013, "Prediction of cutting forces in three and five-axis ball-end milling with tool indentation effect," International Journal of Machine Tools and Manufacture, 66, pp. 66-81. [5] Ismail, F., and Ziaei, R., 2002, "Chatter suppression in five-axis machining of flexible parts," International Journal of Machine Tools and Manufacture, 42(1), pp. 115-122. [6] Eksioglu, C., 2011, "Mechanics and Dynamics of Thin-Wall Machining," MASc. Thesis, The University of British Columbia, Vancouver, BC, Canada. [7] Fischer, A., Eberhard, P., and Ambrósio, J., 2013, "Parametric Flexible Multibody Model for Material Removal During Turning," Journal of Computational and Nonlinear Dynamics, 9(1), p. 011007. [8] Gerasimenko, A., Guskov, M., Duchemin, J., Lorong, P., and Gouskov, A., 2015, "Variable Compliance-related Aspects of Chatter in Turning Thin-walled Tubular Parts," Procedia CIRP, 31, pp. 58-63. [9] Lorong, P. L., Arnaud; Duarte, Alexis Perez, 2011, "Dynamic Study of Thin Wall Part Turning," 13th CIRP International Conference on Modeling of Machining Operations, Advances Materials Research, Sintra, Portugal, pp. 591-599. [10] Mehdi, K., Rigal, J. F., and Play, D., 2002, "Dynamic Behavior of a Thin-Walled Cylindrical Workpiece During the Turning Process, Part 1: Cutting Process Simulation," Journal of Manufacturing Science and Engineering, 124(3), pp. 562-568. [11] Stepan, G., Kiss, A. K., Ghalamchi, B., Sopanen, J., and Bachrathy, D., 2017, "Chatter avoidance in cutting highly flexible workpieces," CIRP Annals, 66(1), pp. 377-380. [12] Urbikain, G., López de Lacalle, L. N., Campa, F. J., Fernández, A., and Elías, A., 2012, "Stability prediction in straight turning of a flexible workpiece by collocation method," International Journal of Machine Tools and Manufacture, 54–55, pp. 73-81. [13] Budak, E., and Altintas, Y., 1995, "Modeling and avoidance of static form errors in peripheral milling of plates," International Journal of Machine Tools and Manufacture, 35(3), pp. 459-476. [14] Elbestawi, M. A., and Sagherian, R., 1991, "Dynamic modeling for the prediction of surface errors in the milling of thin-walled sections," Journal of Materials Processing Technology, 25(2), pp. 215-228. [15] Kline, W. A., DeVor, R. E., and Shareef, I. A., 1982, "The Prediction of Surface Accuracy in End Milling," Journal of Engineering for Industry, 104(3), pp. 272-278. [16] Tsai, J.-S., and Liao, C.-L., 1999, "Finite element modeling of static surface errors in the peripheral milling of thin-walled workpieces," Journal of Materials Processing Technology, 94(2), pp. 235-246. 183  [17] Tsai, M.-P., Tsai, N.-C., and Yeh, C.-W., 2016, "On milling of thin-wall conical and tubular workpieces," Mechanical Systems and Signal Processing, 72–73, pp. 395-408. [18] Wan, M., Zhang, W., Qiu, K., Gao, T., and Yang, Y., 2005, "Numerical Prediction of Static Form Errors in Peripheral Milling of Thin-Walled Workpieces With Irregular Meshes," Journal of Manufacturing Science and Engineering, 127(1), pp. 13-22. [19] Koike, Y., Matsubara, A., and Yamaji, I., 2013, "Design method of material removal process for minimizing workpiece displacement at cutting point," CIRP Annals - Manufacturing Technology, 62(1), pp. 419-422. [20] Arnaud, L., Gonzalo, O., Seguy, S., Jauregi, H., and Peigné, G., 2011, "Simulation of low rigidity part machining applied to thin-walled structures," The International Journal of Advanced Manufacturing Technology, 54(5), pp. 479-488. [21] Bravo, U., Altuzarra, O., López de Lacalle, L. N., Sánchez, J. A., and Campa, F. J., 2005, "Stability limits of milling considering the flexibility of the workpiece and the machine," International Journal of Machine Tools and Manufacture, 45(15), pp. 1669-1680. [22] Ding, Y., and Zhu, L., 2016, "Investigation on chatter stability of thin-walled parts considering its flexibility based on finite element analysis," The International Journal of Advanced Manufacturing Technology, pp. 1-15. [23] Liu, Y., Wu, B., Ma, J., and Zhang, D., 2016, "Chatter identification of the milling process considering dynamics of the thin-walled workpiece," The International Journal of Advanced Manufacturing Technology, 89(5-8), pp. 1765-1773. [24] Mundim, R. B., and Borille, A. V., 2017, "An approach for reducing undesired vibrations in milling of low rigidity structures," The International Journal of Advanced Manufacturing Technology, 88(1), pp. 971-983. [25] Seguy, S., Dessein, G., and Arnaud, L., 2008, "Surface roughness variation of thin wall milling, related to modal interactions," International Journal of Machine Tools and Manufacture, 48(3–4), pp. 261-274. [26] Thevenot, V., Arnaud, L., Dessein, G., and Cazenave–Larroche, G., 2006, "Influence of material removal on the dynamic behavior of thin-walled structures in peripheral milling," Machining Science and Technology, 10(3), pp. 275-287. [27] Mane, I., Gagnol, V., Bouzgarrou, B. C., and Ray, P., 2008, "Stability-based spindle speed control during flexible workpiece high-speed milling," International Journal of Machine Tools and Manufacture, 48(2), pp. 184-194. [28] Song, Q., Ai, X., and Tang, W., 2011, "Prediction of simultaneous dynamic stability limit of time–variable parameters system in thin-walled workpiece high-speed milling processes," The International Journal of Advanced Manufacturing Technology, 55(9), pp. 883-889. [29] Adetoro, O. B., Sim, W. M., and Wen, P. H., 2010, "An improved prediction of stability lobes using nonlinear thin wall dynamics," Journal of Materials Processing Technology, 210(6), pp. 969-979. [30] Campa, F. J., Lopez de Lacalle, L. N., and Celaya, A., 2011, "Chatter avoidance in the milling of thin floors with bull-nose end mills: Model and stability diagrams," International Journal of Machine Tools and Manufacture, 51(1), pp. 43-53. [31] Kersting, P., and Biermann, D., 2014, "Modeling techniques for simulating workpiece deflections in NC milling," CIRP Journal of Manufacturing Science and Technology, 7(1), pp. 48-54. 184  [32] Ahmadi, K., 2017, "Finite strip modeling of the varying dynamics of thin-walled pocket structures during machining," The International Journal of Advanced Manufacturing Technology, 89(9-12), pp. 2691-2699. [33] Meshreki, M., Attia, H., and Kövecses, J., 2011, "Development of a New Model for the Varying Dynamics of Flexible Pocket-Structures During Machining," Journal of Manufacturing Science and Engineering, 133(4), p. 041002. [34] Budak, E., Tunç, L. T., Alan, S., and Özgüven, H. N., 2012, "Prediction of workpiece dynamics and its effects on chatter stability in milling," CIRP Annals, 61(1), pp. 339-342. [35] Song, Q., Liu, Z., Wan, Y., Ju, G., and Shi, J., 2015, "Application of Sherman–Morrison–Woodbury formulas in instantaneous dynamic of peripheral milling for thin-walled component," International Journal of Mechanical Sciences, 96–97, pp. 79-90. [36] Yang, Y., Zhang, W.-H., Ma, Y.-C., and Wan, M., 2016, "Chatter prediction for the peripheral milling of thin-walled workpieces with curved surfaces," International Journal of Machine Tools and Manufacture, 109, pp. 36-48. [37] Huang, C. Y., and Wang, J. J. J., 2006, "Mechanistic Modeling of Process Damping in Peripheral Milling," Journal of Manufacturing Science and Engineering, 129(1), pp. 12-20. [38] Altintas, Y., and Weck, M., 2004, "Chatter Stability of Metal Cutting and Grinding," CIRP Annals, 53(2), pp. 619-642. [39] Sisson, T. R., and Kegg, R. L., 1969, "An Explanation of Low-Speed Chatter Effects," Journal of Engineering for Industry, 91(4), pp. 951-958. [40] Wu, D. W., 1988, "Application of a comprehensive dynamic cutting force model to orthogonal wave-generating processes," International Journal of Mechanical Sciences, 30(8), pp. 581-600. [41] Shaw, M. C., and Desalvo, G. J., 1970, "A new approach to plasticity and its application to blunt two dimensional indenters," Journal of Engineering for Industry, 92(2), pp. 469-479. [42] Abrari, F., Elbestawi, M. A., and Spence, A. D., 1998, "On the dynamics of ball end milling: modeling of cutting forces and stability analysis," International Journal of Machine Tools and Manufacture, 38(3), pp. 215-237. [43] Elbestawi, M. A., Ismail, F., and Ullagaddi, B. C., 1994, "Modelling Machining Dynamics Including Damping in the Tool-Workpiece Interface," Journal of Engineering for Industry, 116(4), pp. 435-439. [44] Lee, B. Y., Tarng, Y. S.,Ma, S.C., 1995, "Modeling of the process damping force in chatter vibration," International Journal of Machine Tools and Manufacture, 35(7), pp. 951-962. [45] Shawky, A. M., and Elbestawi, M. A., 1997, "An Enhanced Dynamic Model in Turning Including the Effect of Ploughing Forces," Journal of Manufacturing Science and Engineering, 119(1), pp. 10-20. [46] Chiou, R. Y., and Liang, S. Y., 1998, "Chatter stability of a slender cutting tool in turning with tool wear effect," International Journal of Machine Tools and Manufacture, 38(4), pp. 315-327. [47] Ahmadi, K., and Ismail, F., 2010, "Experimental investigation of process damping nonlinearity in machining chatter," International Journal of Machine Tools and Manufacture, 50(11), pp. 1006-1014. [48] Altintas, Y., Eynian, M., and Onozuka, H., 2008, "Identification of dynamic cutting force coefficients and chatter stability with process damping," CIRP Annals - Manufacturing Technology, 57(1), pp. 371-374. 185  [49] Sellmeier, V., and Denkena, B., 2012, "High speed process damping in milling," CIRP Journal of Manufacturing Science and Technology, 5(1), pp. 8-19. [50] Budak, E., and Tunc, L. T., 2009, "A New Method for Identification and Modeling of Process Damping in Machining," Journal of Manufacturing Science and Engineering, 131(5), p. 051019. [51] Budak, E., and Tunc, L. T., 2010, "Identification and modeling of process damping in turning and milling using a new approach," CIRP Annals - Manufacturing Technology, 59(1), pp. 403-408. [52] Ahmadi, K., and Ismail, F., 2011, "Analytical stability lobes including nonlinear process damping effect on machining chatter," International Journal of Machine Tools and Manufacture, 51(4), pp. 296-308. [53] Tyler, C. T., Troutman, J. R., and Schmitz, T. L., 2016, "A coupled dynamics, multiple degree of freedom process damping model, Part 2: Milling," Precision Engineering, 46, pp. 73-80. [54] Ahmadi, K., and Altintas, Y., 2014, "Identification of Machining Process Damping Using Output-Only Modal Analysis," Journal of Manufacturing Science and Engineering, 136(5), p. 051017. [55] Wan, M., Feng, J., Ma, Y.-C., and Zhang, W.-H., 2017, "Identification of milling process damping using operational modal analysis," International Journal of Machine Tools and Manufacture, 122, pp. 120-131. [56] Jin, X., and Altintas, Y., 2013, "Chatter Stability Model of Micro-Milling With Process Damping," Journal of Manufacturing Science and Engineering, 135(3), p. 031011. [57] Altintas, Y., and Budak, E., 1995, "Analytical Prediction of Stability Lobes in Milling," CIRP Annals, 44(1), pp. 357-362. [58] Altintas, Y., Shamoto, E., Lee, P., and Budak, E., 1999, "Analytical Prediction of Stability Lobes in Ball End Milling," Journal of Manufacturing Science and Engineering, 121(4), pp. 586-592. [59] Altintas, Y., 2001, "Analytical Prediction of Three Dimensional Chatter Stability in Milling," JSME International Journal Series C Mechanical Systems, Machine Elements and Manufacturing, 44(3), pp. 717-723. [60] Liu, X., Li, R., Wu, S., Yang, L., and Yue, C., 2017, "A prediction method of milling chatter stability for complex surface mold," The International Journal of Advanced Manufacturing Technology, 89(9), pp. 2637-2648. [61] Shamoto, E., and Akazawa, K., 2009, "Analytical prediction of chatter stability in ball end milling with tool inclination," CIRP Annals, 58(1), pp. 351-354. [62] Ferry, W. B., 2008, "Virtual five-axis flank mlling of jet engine impellers," PhD Thesis, The University of British Columbia, Vancouver, BC. [63] Kersting, P., and Biermann, D., 2009, "Simulation concept for predicting workpiece vibrations in five-axis milling," Machining Science and Technology, 13(2), pp. 196-209. [64] Ozturk, E., and Budak, E., 2010, "Dynamics and Stability of Five-Axis Ball-End Milling," Journal of Manufacturing Science and Engineering, 132(2), p. 021003. [65] Dai, Y., Li, H., Wei, Z., and Zhang, H., 2018, "Chatter stability prediction for five-axis ball end milling with precise integration method," Journal of Manufacturing Processes, 32, pp. 20-31. 186  [66] Ozturk, E., Tunc, L. T., and Budak, E., 2009, "Investigation of lead and tilt angle effects in 5-axis ball-end milling processes," International Journal of Machine Tools and Manufacture, 49(14), pp. 1053-1062. [67] Sun, C., and Altintas, Y., 2016, "Chatter free tool orientations in 5-axis ball-end milling," International Journal of Machine Tools and Manufacture, 106, pp. 89-97. [68] D’Ambrogio, W., and Fregolent, A., 2010, "The role of interface DoFs in decoupling of substructures based on the dual domain decomposition," Mechanical Systems and Signal Processing, 24(7), pp. 2035-2048. [69] Voormeeren, S. N., and Rixen, D. J., 2012, "A family of substructure decoupling techniques based on a dual assembly approach," Mechanical Systems and Signal Processing, 27, pp. 379-396. [70] Craig, R. R., Jr., and Kurdila, A. J., 2006, Fundamentals of Structural Dynamics, Wiley, Hoboken, NJ. [71] Valk, P. L. C. V. d., 2010, "Model Reduction and Interface Modeling in Dynamic Substructuring: Application to a Multi-megawatt Wind Turbine," M.Sc. Thesis, Delft University of Technology, Delft, The Netherlands. [72] Braun, S. G., and Ram, Y. M., 2001, "Modal modification of vibrating systems: some problems and their solutions," Mechanical Systems and Signal Processing, 15(1), pp. 101-119. [73] Duarte, M. L. M., 1996, "Experimentally-derived structural models for use in further dynamic analysis," Ph.D. Thesis, Imperial College of Science, Technology, and Medicine, London, UK. [74] Zhao, Y.-q., Chen, S.-h., Chai, S., and Qu, Q.-w., 2002, "An improved modal truncation method for responses to harmonic excitation," Computers & Structures, 80(1), pp. 99-103. [75] Wilkinson, J. H., 1965, The algebraic eigenvalue problem, Oxford University Press,, London, UK. [76] Law, M., Phani, A. S., and Altintas, Y., 2013, "Position-Dependent Multibody Dynamic Modeling of Machine Tools Based on Improved Reduced Order Models," Journal of Manufacturing Science and Engineering, 135(2), p. 021008. [77] Kammer, D. C., 1991, "Sensor placement for on-orbit modal identification and correlation of large space structures," Journal of Guidance, Control, and Dynamics, 14(2), pp. 251-259. [78] Li, D. S., Li, H. N., and Fritzen, C. P., 2007, "The connection between effective independence and modal kinetic energy methods for sensor placement," Journal of Sound and Vibration, 305(4), pp. 945-955. [79] John O'Callahan, Peter Avitabile, and Riemer, R., 1989, "System Equivalent Reduction Expansion Process (SEREP)," Proc. 7th International Modal Analysis Conference (IMAC), Society of Experimental Mechanics, pp. 29-37. [80] Sastry, C. V. S., Roy Mahapatra, D., Gopalakrishnan, S., and Ramamurthy, T. S., 2003, "An iterative system equivalent reduction expansion process for extraction of high frequency response from reduced order finite element model," Computer Methods in Applied Mechanics and Engineering, 192(15), pp. 1821-1840. [81] Friswell, M. I., Garvey, S. D., and Penny, J. E. T., 1995, "Model reduction using dynamic and iterated IRS techniques," Journal of Sound and Vibration, 186(2), pp. 311-323. [82] Guyan, R. J., 1965, "Reduction of stiffness and mass matrices," AIAA Journal, 3(2), pp. 380-380. 187  [83] Cho, M., and Kim, H., 2004, "Element-Based Node Selection Method for Reduction of Eigenvalue Problems," AIAA Journal, 42(8), pp. 1677-1684. [84] Chen, J. C., and Wada, B. K., 1977, "Matrix Perturbation for Structural Dynamic Analysis," AIAA Journal, 15(8), pp. 1095-1100. [85] Chen, S. H., Wu, X. M., and Yang, Z. J., 2006, "Eigensolution reanalysis of modified structures using epsilon-algorithm," International Journal for Numerical Methods in Engineering, 66(13), pp. 2115-2130. [86] Jezequel, L., 1990, "Procedure to reduce the effects of modal truncation in eigensolution reanalysis," AIAA Journal, 28(5), pp. 896-902. [87] Brandon, J. A., 1984, "Derivation and significance of second-order modal design sensitivities," AIAA Journal, 22(5), pp. 723-724. [88] Parlett, B., 1998, The Symmetric Eigenvalue Problem, Society for Industrial and Applied Mathematics. [89] Chen, S. H., Ma, L., and Meng, G., 2008, "Dynamic response reanalysis for modified structures under arbitrary excitation using epsilon-algorithm," Computers & Structures, 86(23–24), pp. 2095-2101. [90] Press, H. W., Teukolsky, A. S., Vetterling, T. W., and Flannery, P. B., 1992, Numerical Recipes in C-The Art of Scientific Computing, Cambridge University Press, USA. [91] Adams, G. G., and Nosonovsky, M., 2000, "Contact modeling — forces," Tribology International, 33(5), pp. 431-442. [92] Carlsson, S., Biwa, S., and Larsson, P. L., 2000, "On frictional effects at inelastic contact between spherical bodies," International Journal of Mechanical Sciences, 42(1), pp. 107-128. [93] Hills, D. A., Nowell, D., and Sackfield, A., 1993, Mechanics of Elastic Contacts, Butterworth-Heinemann, Great Britain. [94] Goryacheva, I. G., Murthy, H., and Farris, T. N., 2002, "Contact problem with partial slip for the inclined punch with rounded edges," International Journal of Fatigue, 24(11), pp. 1191-1201. [95] Mason, J. C., and Handscomb, D. C., 2003, Chebyshev Polynomials, CRC Press LLC, USA. [96] Gil, A., Segura, J., and Temme, N., 2007, Numerical Methods for Special Functions, Society for Industrial and Applied Mathematics. [97] Johnson, K. L., 1985, Contact Mechanics, Cambridge University Press, Cambridge, England. [98] Yu, W., and Blanchard, J. P., 1996, "An elastic-plastic indentation model and its solutions," Journal of Materials Research, 11(9), pp. 2358-2367. [99] Tabor, D., 2000, The Hardness of Metals, Oxford University Press, Oxford, England. [100] Johnson, K. L., 1970, "The correlation of indentation experiments," Journal of the Mechanics and Physics of Solids, 18(2), pp. 115-126. [101] Jackson, R. L., and Green, I., 2005, "A Finite Element Study of Elasto-Plastic Hemispherical Contact Against a Rigid Flat," Journal of Tribology, 127(2), pp. 343-354. [102] Eynian, M., and Altintas, Y., 2010, "Analytical Chatter Stability of Milling With Rotating Cutter Dynamics at Process Damping Speeds," Journal of Manufacturing Science and Engineering, 132(2), pp. 021012-021012-021014. [103] de Oliveira, F. B., Rodrigues, A. R., Coelho, R. T., and de Souza, A. F., 2015, "Size effect and minimum chip thickness in micromilling," International Journal of Machine Tools and Manufacture, 89, pp. 39-54. 188  [104] Grzesik, W., 1996, "A revised model for predicting surface roughness in turning," Wear, 194(1), pp. 143-148. [105] Liu, X., DeVor, R. E., and Kapoor, S. G., 2005, "An Analytical Model for the Prediction of Minimum Chip Thickness in Micromachining," Journal of Manufacturing Science and Engineering, 128(2), pp. 474-481. [106] Woon, K. S., Rahman, M., Neo, K. S., and Liu, K., 2008, "The effect of tool edge radius on the contact phenomenon of tool-based micromachining," International Journal of Machine Tools and Manufacture, 48(12), pp. 1395-1407. [107] Yuan, Z. J., Zhou, M., and Dong, S., 1996, "Effect of diamond tool sharpness on minimum cutting thickness and cutting surface integrity in ultraprecision machining," Journal of Materials Processing Technology, 62(4), pp. 327-330. [108] Islam, C., Lazoglu, I., and Altintas, Y., 2015, "A Three-Dimensional Transient Thermal Model for Machining," Journal of Manufacturing Science and Engineering, 138(2), p. 021003. [109] Rezaei, H., Sadeghi, M. H., and Budak, E., 2018, "Determination of minimum uncut chip thickness under various machining conditions during micro-milling of Ti-6Al-4V," The International Journal of Advanced Manufacturing Technology, 95(5), pp. 1617-1634. [110] Yen, Y.-C., Jain, A., and Altan, T., 2004, "A finite element analysis of orthogonal machining using different tool edge geometries," Journal of Materials Processing Technology, 146(1), pp. 72-81. [111] Engin, S., and Altintas, Y., 2001, "Mechanics and dynamics of general milling cutters.: Part I: helical end mills," International Journal of Machine Tools and Manufacture, 41(15), pp. 2195-2212. [112] Lee, P., and Altintas, Y., 1996, "Prediction of ball-end milling forces from orthogonal cutting data," International Journal of Machine Tools and Manufacture, 36(9), pp. 1059-1072. [113] Altintas, Y., 2012, Manufacturing Automation, Cambridge University Press, New York, NY, USA. [114] Budak, E., Altintas, Y., and Armarego, E. J. A., 1996, "Prediction of Milling Force Coefficients From Orthogonal Cutting Data," Journal of Manufacturing Science and Engineering, 118(2), pp. 216-224. [115] Eksioglu, C., Kilic, Z. M., and Altintas, Y., 2012, "Discrete-Time Prediction of Chatter Stability, Cutting Forces, and Surface Location Errors in Flexible Milling Systems," Journal of Manufacturing Science and Engineering, 134(6), p. 061006. [116] Merdol, S. D., 2008, "Virtual three-axis milling process simulation and optimization," PhD Thesis, The University of British Columbia, Vancouver, BC. [117] Eynian, M., 2010, "Chatter Stability of Turning and Milling with Process Damping," PhD Thesis, The University of British Columbia, Vancouver, BC. [118] Tunc, L. T., 2010, "Selection of Cutting Strategy And Parameters in Multiaxis Machining Operations for Improved Productivity," PhD Thesis, Sabanci University, Istanbul, Turkey. [119] Malekian, M., Mostofa, M. G., Park, S. S., and Jun, M. B. G., 2012, "Modeling of minimum uncut chip thickness in micro machining of aluminum," Journal of Materials Processing Technology, 212(3), pp. 553-559. [120] Ahmadi, K., 2011, "Machining chatter in flank milling and investigation of process damping in surface generation," PhD Thesis, The University of Waterloo, Waterloo, ON, Canada. 189  [121] Li, Z.-L., Tuysuz, O., Zhu, L.-M., and Altintas, Y., 2018, "Surface form error prediction in five-axis flank milling of thin-walled parts," International Journal of Machine Tools and Manufacture, 128, pp. 21-32. [122] Kilic, Z. M., 2015, "Generalized modeling of flexible machining system with arbitrary tool geometry," PhD Thesis, The University of Brisith Columbia, Vancouver, BC, Canada. [123] Brezinski, C., 2000, "Convergence acceleration during the 20th century," Journal of Computational and Applied Mathematics, 122(1–2), pp. 1-21. [124] Roberts, D. E., 1996, "On the convergence of rows of vector Padé approximants," Journal of Computational and Applied Mathematics, 70(1), pp. 95-109.   190  Appendices Appendix A  Epsilon Algorithm Perturbation of the mode shapes of the in-process workpiece provides a significant computational saving as a replacement of repetitively solved Eq.(4.9) for the varying dynamic characteristics with material removal. However, the truncation in Eq.(4.11) after the second order term might cause perturbation to become inaccurate especially for large material removals. Therefore, the Epsilon algorithm, which accelerates the convergence of linear sequences [123], is used as follows.  Consider the matrix Γ  is formed column-wise by the thmj  mode shape of the in-process workpiece 1iB   and its perturbation terms as,  11, 2, q ,, , ,...,i i i im m m p mB A A Aj j j j      Γ       (A.1) The partial summation of the first z columns of the matrix Γ  can be expressed as,    1zz ddΓ Γ   (A.2) where dΓ  is the thd  column of Γ , and zΓ  is the thz  column of the partial summation matrix Γ . The vectors ε  are recursively constructed to form a table of vectors until a single vector, which is equal to imBj  in Eq.(4.11), is obtained at the last calculation step using the columns of  Γ  as,   1 ,zOz = 1... N +2 0ε   (A.3)   2 ,z zOz = 1... N +1ε Γ   (A.4)  1 1 11 1 ( ) , ;z z z ze e e e z = 1,2,... e= 2,3,...     ε ε ε ε   (A.5) 191  where upper and lower indices of the vector ε  denote column and row numbers of the epsilon table, respectively. The vector inverse in Eq.(A.5) is expressed as [124],  11 11 1( )( )( ) ( )z zz z e ee e z z T z ze e e e    ε εε εε ε ε ε  (A.6)   192  Appendix B  Directional Coefficients for Dynamic Cutting and Process Damping Forces In Section 6.2.3, the time variable in Eq.(6.44) is replaced with the rotation angle of the jth flute using the following change of variables,  ,, ,,,,(t, z)2 2(t, z)(z) 2 (rad)(z) (rad)Then,t 0 (t, z) 0t (z) (t, z) (rad)jT j s jp j p j sT j js j p jjj j p jdt t dtTTT                 (B.1) The directional coefficients for the dynamic cutting and process damping forces are obtained using Fourier Series expansion in Eq.(6.46) and the average cutting force coefficients from Eq.(6.47) for the jth flute at axial position z. The general equations for the directional coefficient matrices are,    ,,, , ,2, , , ,,, , ,,2,,1 1(t, z)(z)1 1(t, z)(z)exp j jstexp j jstr r rxx j xy j xz jircd cd r r rr j j j yx j yy j yz jk ks j p jr r rzx j zy j zz jkrxx j xyirpd pdr j j jk ks j p ja a ae d a a aTa a ab be dT                            A AB B, ,, , ,, , ,r rj xz jr r ryx j yy j yz jr r rzx j zy j zz jkbb b bb b b        (B.2) The individual entries in Eq.(B.2) can be evaluated as follows, 193     1 , ,,21 ,21sin 2 (t, z) (t, z) (z)cos (z) (t, z) (t, z) (z)sin (z)cos 2 (t, z) (t, z)cos (z) (t, z)sin (z) (t, z) (z)(t, z)cos (z) (t, z)sin (z)4j ac p j j tc rc p j jrxx jj ac j rc j tc p jac j rc jK irK KairK irK KK K ir                   (B.3) where, ,2 (t,z)(z)1jp jire        and  2 2 22 ,4 (z)p j r    .   1 ,,21 , ,21sin 2 (t, z) (t, z)cos (z) (t, z)sin (z) (t, z) (z) )cos 2 (t, z) (t, z) (z)cos (z) (t, z) (t, z) (z)sin (z)(t, z)4j ac j rc j tc p jrxy jj ac p j j tc rc p j jtcirK irK K iaK irK KK ir                (B.4)   ,,3 3, ,2 2 2,(t, z) (z)cos (z) cos (t, z)2 (t, z)cos (t, z)(t, z) (z) sin (t, z)cos (z) (sin + cos ) + (t, z) (z) sin (z) cos (t, z)sin (z) (z) 42 (t, z)sinac p j j jtc jtc p j jjrxz j rc p j j jj p jrciKrKiKia iKrrK                (z)sin (t, z)2 (t, z)cos (z) sin (t, z)   j jac j jrK                (B.5) where, 3,2 (t, z)(z)jp jr  . 194    1 ,,21 , ,21sin 2 (t, z) (t, z)cos (z) (t, z)sin (z) (t, z) (z)cos 2 (t, z) (t, z) (z)cos (z) (t, z) (t, z) (z)sin (z)(t, z)4j ac j rc j tc p jryx jj ac p j j tc rc p j jtcirK irK KaK irK KiKr                (B.6)    1 , ,yy,21 ,21sin 2 (t, z) (t, z) (z)cos (z) (t, z) (t, z) (z)sin (z)cos 2 (t, z) (t, z)cos (z) (t, z)sin (z) (t, z) (z)(t, z)cos (z) (t, z)sin (z)4j ac p j j tc rc p j jrjj ac j rc j tc p jac j rc jK irK Kai rK irK KK K ir                   (B.7)  ,,,3 3yz, 2 2 2,(t, z) (z)cos (t, z) 2 (t, z)sin (t, z)(t, z) (z)cos (z)sin (t, z)+ (t, z) (z)sin (z)sin (t, z)cos (z) (sin + cos ) +2 (t, z)cos (zsin (z) (z) 4tc p j j tc jac p j j jrc p j j jjrjac jj p jiK rKiKiKiarKr                ) cos (t, z)+2 (t, z)sin (z) cos (t, z)   jrc j jrK              (B.8)   1 ,, 2 2 2,(z)cos (t, z) 2 sin (t, z) (t, z)cos (z) (t, z)sin (z)(z) 4p j j j rc j ac jrzx jp jir K Kar          (B.9)   1 ,, 2 2 2,(z)sin (t, z) 2 cos (t, z) (t, z)cos (z) (t, z)sin (z)(z) 4p j j j rc j ac jrzy jp jir K Kar          (B.10)   21,(t, z)cos (z) (t, z)cos (z) / sin (z)2ac j rc j jrzz ji K Kar      (B.11)   195   2 2 2 2,1, 2 2 2 2 2 2, , ,sin (z) cos2 (t, z)sin (z) (z)cos2 (t, z)4 (z) (z)sin2 (t, z)sin (z) sin2 (t, z) (z)sin (z)j j j p j jrxx jp j p j j j j p j ji r i r rbr r r i r i                            (B.12)  2 2 2 2,1, 2 2 2 2 2 2, , ,sin2 (t, z)sin (z) (z)sin2 (t, z)4 (z) (z)cos 2 (t, z)sin (z) cos2 (t, z) (z)j j p j jrxy jp j p j j j j p ji r i r rbr r r i r i                           (B.13)   , ,3 3, 2 2 2,(z)cos (t, z)sin (z) (z)sin (t, z)cos (z) cos sin2 cos (t, z) 2 sin (t, z)sin (z)sin (z) (z) 4p j j j p j jjrxz jj j jj p jibi r i rr                      (B.14)   2 2,1, 2 2 2 2 2 2 2 2, , ,(z)cos 2 (t, z)sin (z) sin 2 (t, z)sin (z)4 (z) (z)sin 2 (t, z) cos 2 (t, z) (z)p j j j j jryx jp j p j j j p jr i rbr r r i r i r i                          (B.15)  2 2,2 21yy, ,2 2 2, 2 2 2,(z)cos 2 (t, z) cos 2 (t, z)sin (z)sin (z) (z)sin 2 (t, z)sin (z)4 (z)sin 2 (t, z) (z)sin (z)p j j j jrj j p j j jp jj p j jr i rb i r rr ri r i                            (B.16)   , ,3 3yz, 2 2 2,(z) cos (t, z) (z)sin (t, z)sin (z)cos (z) cos sin2 sin (t, z) 2 cos (t, z)sin (z)sin (z) (z) 4p j j p j j jjrjj j jj p jibi r i rr                     (B.17)   1 ,, 2 2 2,cos (z) (z)cos (t, z) 2 sin (t, z)(z) 4j p j j jrzx jp ji rbr         (B.18)   1 ,, 2 2 2,cos (z) (z)sin (t, z) 2 cos (t, z)(z) 4j p j j jrzy jp ji rbr        (B.19)  21,cos (z)2 sin (z)jrzz jjibr     (B.20) 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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-0373471/manifest

Comment

Related Items