A L E SIMULATION OF CHIP FORMATION IN ORTHOGONAL M E T A L C U T T I N G PROCESS By Mohammad R. Movahhedy B. Sc. University of Tehran, Iran M. Sc. University of Waterloo, Canada A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCOTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF MECHANICAL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA April 2000 Â© Mohammad R. Movahhedy, 2000 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mechanical Engineering The University of British Columbia 2324 Main Mall Vancouver, Canada V6T 1Z4 Date: y4prtl 2g , ZOOo Abstract This thesis deals with the application of the Arbitrary Lagrangian-Eulerian (ALE) finite element analysis in simulation of chip formation in orthogonal metal cutting process. A critical review of the literature in thisfieldshows that due to the very complex set of conditions present in a cutting process, the application of conventional Lagrangian and Eulerian methods for this problem is inefficient and entails numerical difficulties. In particular, the pertinent problems in the node separation technique or the remeshing approach are discussed. In contrast, the adaptivity of the mesh in an ALE analysis provides the possibility of combining the strengths of both Lagrangian and Eulerian methods in a single analysis. In this approach, the chip formation occurs as a result of plasticflowof the work material around the tool edge on the one hand, and unconstrainedflowof the material on free surfaces of the chip, on the other. Due to high deformation speed, strain rate and temperature play a significant role in the chip formation process. In this work, the ALE formulation originally proposed by Wang and Gadala [72] is extended to include rate and thermal effects. A heat transfer module is included that updates the temperaturefieldin the cutting zone at each step of the analysis. Contact algorithms are developed which are able to detect emerging contact conditions and apply contact constraints at interfaces betweenflexible-flexibleor rigid-flexible pairs. An efficient ALE mesh motion is designed that prevents element distortion in the deformation zone and at the same time facilitates the evolution of the chip size at free boundaries. Furthermore, General guidelines for designing a mesh motion strategy are presented, an algorithm for mesh sliding on free boundaries is introduced, and transfinite and isoparametric mapping techniques are adopted for moving the mesh in the interior of the body, so that the mesh remains optimal throughout the analysis. n The large deformation, rate-dependent, thermo-mechanical ALE finite element code is used to simulate chip formation in orthogonal metal cutting processes. The results of simulation of cutting low carbon steel with a carbide tool are presented as a benchmark problem, and a parametric study is conducted to investigate the influence of cutting conditions on the chip formation process. The results of these studies are verified through comparison with available experimental data. The fairly good qualitative and quantitative agreement between predicted and experimental results confirms that ALE simulation provides a realistic representation of the actual process. Finally, the influence of cutting edge geometry on the chip formation process is investigated through simulation of high speed cutting of hardened steel alloys with chamfered or worn tools of carbide or CBN type. This study shows that changing the chamfer angle does not affect the chip significantly, because the dead metal zone that is formed under the chamfer acts as the main cutting edge of the tool. However, cutting and thrust forces increase with increase in the chamfer angle. The predictions on the effects of chamfer angle are in agreement with experimental observations. A study of the influence of cutting speed on the deformation process for chamfered CBN tools shows that at higher cutting speeds, the maximum temperature on the rake face increases substantially, signifying the importance of diffusion wear at such speeds. in Table of Contents Abstract ii List of Tables vii List of Figures viii Nomenclature xii Acknowledgement xvi Dedication xvii 1 INTRODUCTION 1 1.1 MOTIVATION 1 1.2 DESCRIPTION OF ORTHOGONAL METAL CUTTING PROCESS 2 1.3 SCOPE OF THIS RESEARCH 4 2 MODELING OF C U T T I N G PROCESS IN T H E L I T E R A T U R E 2.1 ANALYTICAL AND EMPIRICAL MODELS OF CUTTING PROCESS 2.1.1 2.2 2.3 6 ... Friction and Temperature on the Rake Face REVIEW OF NUMERICAL SIMULATIONS OF CUTTING PROCESS 6 8 ... 11 2.2.1 Review of Lagrangian Cutting Analyses 14 2.2.2 Sample Lagrangian Simulations 20 2.2.3 Eulerian Simulation of Metal Cutting Process 25 ARBITRARY LAGRANGIAN-EULERIAN FORMULATION 27 2.3.1 27 History of ALE Approach iv 2.3.2 ALE Applications in Metal Cutting Simulation 30 2.3.3 The ALE Work in This Research 32 3 T H E O R Y OF A L E FINITE E L E M E N T M E T H O D 3.1 3.2 4 33 ALE FORMULATION 33 3.1.1 Governing Equations in ALE Analysis 36 3.1.2 Virtual Work Equation and Its Linearization 38 STRAIN RATE AND TEMPERATURE EFFECTS 43 3.2.1 Treatment of Rate Effects in the Literature 43 3.2.2 Formulation of the Constitutive Equation 44 3.3 HEAT TRANSFER ANALYSIS 47 3.4 CONTACT AND FRICTION 52 3.4.1 Treatment of Contact Problems in the Literature 53 3.4.2 Contact Between Deformable Bodies 55 3.4.3 Contact Between Rigid-Deformable Bodies 63 3.4.4 Contact Problems in ALE Analysis 64 N U M E R I C A L PROCEDURES A N D IMPLEMENTATION 4.1 69 MESH MOTION IN ALE ANALYSIS 69 4.1.1 Solution Procedure for Mesh Motion 70 4.1.2 Computation of Mesh Velocities 72 4.1.3 Numerical Examples of ALE Mesh Motion 84 4.2 UPDATING SOLUTION VARIABLES 90 4.3 PROGRAM STRUCTURE 94 4.3.1 Main Program 95 4.3.2 Deformation Module 95 4.3.3 Contact Module 102 v 4.3.4 Heat Transfer Module 106 5 N U M E R I C A L RESULTS A N D DISCUSSION 5.1 108 ON THE DETERMINATION OF MATERIAL PROPERTIES FOR CUTTING SIMULATION 108 5.2 CUTTING MODEL AND MESH MOTION SCHEME Ill 5.3 THE BENCHMARK CUTTING PROBLEM 118 5.3.1 Input Data 118 5.3.2 Simulation Results 119 5.4 THE EFFECTS OF CUTTING CONDITIONS ON CHIP FORMATION PROCESS 5.5 5.6 126 5.4.1 Experimental Data 126 5.4.2 Numerical Model 127 5.4.3 Numerical Results and Discussion 128 5.4.4 Influence of Cutting Conditions on Process Variables 134 ANALYSIS OF CUTTING WITH CHAMFERED AND WORN EDGE TOOLS 139 5.5.1 Finite Element Model 143 5.5.2 Experimental Work and Input Data 143 5.5.3 Influence of Chamfer Angle on the Cutting Process 5.5.4 Influence of Cutting Speed on Cutting with Chamfered Tools ORTHOGONAL CUTTING OF A TITANIUM ALLOY 6 CONCLUSIONS A N D F U T U R E WORK 145 156 159 162 6.1 DISCUSSION AND CONCLUSIONS 162 6.2 FUTURE WORK 167 Bibliography 169 vi List of Tables 5.1 material properties for the low carbon free cutting steel .119 5.2 Experimental and simulation results. The columns represent rake angle, cutting speed, cutting force, thrust force, shear angle, contact length, and maximum temperature, respectively. Note that cutting conditions in simulations (2) and (3) are different from those in experiment and simulation (1). The uncut chip thickness in all cases is 0.1 mm 123 5.3 Cutting conditions used in the parametric study 128 5.4 Tool edge geometry in cutting tests with carbide tools 5.5 Material properties for P20 and carbide tool [128] 5.6 Material properties for titanium alloy [31] 5.7 Forces in the cutting and thrust directions and chip thicknesses from simulation and experiment [145] for cutting titanium alloy TiSAUV vii 144 145 159 161 List of Figures 1.1 Schematic view of orthogonal metal cutting process 2 2.1 Force and velocity components in Merchant's model 7 2.2 Distribution of normal and frictional stresses on the rake face, left) Zorev's Model [17], right) results of Barrow, et al. [19] 9 2.3 Node separation in Lagrangian analysis 16 2.4 Lagrangian analysis with node separation; cutting with a 30Â° rake angle tool . . 22 2.5 Lagrangian analysis with node separation; cutting with a 10Â° rake angle tool . . 22 2.6 Forces in Lagrangian analysis with node separation 2.7 Lagrangian analyses with indentation using DEFORM-2D code 24 3.1 Material and mesh motion in ALE analysis 34 3.2 Bodies in contact 3.3 Corresponding nodes of tool and workpiece in contact 58 3.4 Pseudo-elements used for detecting penetration 62 3.5 Lagrangian and ALE analysis of punch indentation problem, a) initial config- 23 55 uration b) punch force in Lagrangian and ALE analyses, c) deformed shape in Lagrangian analysis, d) deformed shape in ALE analysis 3.6 4.1 65 Example of a) Lagrangian, and b) ALE analyses of cutting with restricted contact length tool 67 Transfinite mapping 79 viii 4.2 Transfinite mesh motion a) distorted mesh, b) mesh after applying mesh motion scheme, c) mesh motion with the effect of boundary motion 81 4.3 re-mapping the mesh using isoparametric method 83 4.4 Mode-I crack extension in center-cracked plate, a) initial configuration and boundary conditions, b,c) configuration and stress field after 50% and 100% crack extension, respectively, d) results from NISA 87 4.5 Forging example a) mesh in Lagrangian analysis, b) mesh in ALE analysis . . . 4.6 Flowchart of the program structure 96 4.7 Flowchart of the deformation module 97 5.1 Schematic view of the ALE model illustrating the mesh motion scheme . . . . 113 5.2 Initial model and various stages of chip formation in cutting simulations (continued onfigure5.3) 5.3 115 Various stages of chip formation in cutting simulations (continued from figure 5.2) 5.4 116 Overlaying the initial and the steady state meshes shows that the chip thickness is independent of the initial assumption of its shape 5.5 89 117 Contours of process variables in simulation (1) [a = 0Â°, V â€” 150 m/min]; 0 w (a) effective stress, (b) effective plastic strain, (c) effective plastic strain rate, and (d) temperature distribution 5.6 121 Cutting and thrust forces predicted versus cutting length in simulation (1) in comparison with experimental values [130] 5.7 distribution of predicted normal and frictional stresses and temperature on the chip-tool interface in comparison with experimental results [130] 5.8 122 124 Contours of strain rate and temperature for high speed cutting simulation (2) [a = 0Â°, V = 450 m/min] 0 125 w ix 5.9 Chip formation and temperature distribution in simulation (3) [a = â€”5Â°, 0 V = 150 m/min] 125 w 5.10 Velocityfieldin the chip and the workpiece 129 5.11 Contours of distribution of solution variables in the cutting zones 130 5.12 Distribution of (a) equivalent plastic strain rate, and (b) maximum shear stress on the shear plane. L is the length of shear plane 131 s 5.13 Temperature distribution on (a) the shear plane, (b) chip-tool interface. L and s L are length of shear plane and contact length, respectively 132 c 5.14 distribution of normal and frictional stresses and temperature on the shear plane for [a = 0Â°, V = 100 m/min, and / = 0.2mm] in comparison with 0 w experimental data [131] 135 5.15 The effect of rake angle on cutting and thrust forces, shear angle and maximum temperature on the rake face 136 5.16 Contact length and chip curling for 0Â°,6Â°, and 12Â° rake angles, from left to right, respectively. Temperature contours are also shown in the figure 136 5.17 The effect of cutting speed on forces, shear angle and temperature 137 5.18 The effect of feed rate on forces, shear angle and temperature 139 5.19 Schematic view of cutting with chamfered edge tool 141 5.20 Velocity field on the chamfered edge of tools with different chamfer angle shows the presence of a stagnated metal zone; a) â€”10Â°, b) â€”25Â°, c) â€”35Â°, and d) worn tool of radius 0.075 mm. All velocities are in (mm/s) 147 5.21 Distributions of effective stress (top) and plastic strain rate (bottom) for different chamfer angles 148 5.22 Distribution of a xx for sharp and chamfered tools of various chamfer angles . . 151 5.23 Distribution of a for sharp and chamfered tools of various chamfer angles . . 152 yy x 5.24 Cutting and thrust forces versus chamfer angle. Numerical predictions are compared with experimental and analytical results of Ren, et. al [15] 153 5.25 temperature distribution in the chip and the tool for the case of â€”35Â° chamfer angle 154 5.26 Distributions of (a) a , (b) a , (c) effective plastic strain rate, and (d) temperxx yy ature for the worn tool 155 5.27 Distributions of effective plastic strain rate (left), and temperature (right) for the CBN tools at cutting speeds of V = 240, 600, and 1000 m/min, from top w to bottom, respectively 157 5.28 effect of cutting speed on the cutting and thrust forces. Numerical predictions are compared with experimental and analytical results of Ren, et. al [15] . . . . 158 5.29 (a) Distribution of normal and frictional stresses on the rake face [a = 2Â°, V = w 60m/min, h = 0.06mm], (b), temperature on the rake face for two different cutting speeds [a = 2Â°,h â€” 0.06mm] xi 161 Nomenclature coefficients in mesh motion scheme A area of an element bi components of body force vector bf length of chamfer cutting edge C thermal capacitance matrix Gijki material constitutive tensor c specific heat of work material Ci components of contact force c D components of rate of deformation tensor B *j elastic part of rate of deformation tensor D plastic part of rate of deformation tensor i5 e v 13 thermal part of rate of deformation tensor D effective plastic strain rate E Young's modulus e specific thermal energy p.. components of deformation Gradient tensor F shear force acting along shear plane in cutting process F cutting force in cutting process F thrust force in cutting process f loading vector in FE equation fi surface tractions h uncut chip thickness p a c t xii h chip thickness h thermal convection coefficient c H ,h shape functions J Jacobian matrix K stiffness matrix in FE equations a a K ,K m c parts of stiffness matrix related to material and mesh motion, respectively k shearflowstress of material lc contact length along the tool-chip interface Is length of shear plane m strain hardening index n strain rate hardening index components of normal vector Q thermal load vector components of heatfluxvector components of rotation tensor R resultant force on the chip and the tool r, 5 normalized parametric coordinates of an isoparametric element T transformation matrix T temperature T -*â€¢ room t,At initial temperature of work material time, time increment Ui components of material displacement vector Ui components of mesh displacement vector V volume of an element V cutting speed Xlll v. v shear velocity chip velocity c components of material velocity vector Vi components of mesh velocity vector v normal component of velocity vector v tangential component of velocity vector Wi =Vi â€” Vi, convective velocity n t position vector in material coordinate system current yield stress of the material Y a, a main rake angle of cutting tool ot\ chamfer angle (second rake angle) 0 friction angle at the tool-chip interface Xi position vector in referential coordinate system $ij Kronecker's delta components of strain tensor e effective strain e effective strain rate <f> shear angle r pseudo-load vector due to rate and thermal dependence V percentage of work converted to heat K thermal conductivity coefficient A proportionality factor in flow rule coefficient of friction V Poisson's ratio e angle for coupled displacements xiv p density (Tij components of Cauchy stress tensor <x effective stress,flowstress cr normal stress n <r Jaumann rate of Cauchy stress a time rate of Cauchy stress r friction stress '{} the quantity in current time increment * *{} +A m e quantity in next time increments xv Acknowledgement IN THE N A M E OF A L L A H , THE BENEFICENT, THE MERCIFUL "He, who teaches me a word, makes me his servant" A l i Ibn Abi-Taleb (600-661 A . D . ) I wish to express my sincere gratitude to my co-supervisors, Dr. Yusuf Altintas and Dr. Mohamed S. Gadala for their patience, encouragement, invaluable advice and continuous support throughout the course of this thesis. It has been a great opportunity for me to be able to learn from their knowledge and experience. I also would like to thank the members of my Ph.D. advisory committee, Drs. D. Anderson, H. Vaughan, and I. Yellowley for their helpful discussions and suggestions at various stages of this work. I am indebted forever to my parents for their lasting love and prayerful support. I also would like to express my deepest gratitude to my wife, Fatemeh, and my children, Hossein and Kosar for their loyalty, patience, encouragements and all the help they offered me throughout the completion of this work. This thesis is dedicated to all of them. Finally, I would like to gratefully acknowledge the financial support of Ministry of Higher Education of Islamic Republic of Iran. This research is also partially supported by natural sciences and engineering research council of Canada (NSERC) under grants #11R81878 and #11R81510, as well as NSERC grant #11R80193 with General Motors of Canada and Boeing Company. These contributions are thankfully acknowledged. xvi Dedicated, As the Smallest Token of Devotion, To: T H E GREATEST OF THE MANKIND, T H E PINNACLE OF A L L THAT IS HUMANITY, T H E ULTIMATE GOAL OF THE CREATION, T H E LAST MESSENGER OF A L L A H ; MUHAMMAD, T H E PROPHET OF ISLAM (PEACE B E UPON HIM). xvii Chapter 1 INTRODUCTION 1.1 MOTIVATION Machining is one of the most common manufacturing processes. Nearly every mechanical component in use has undergone a machining operation at some stage of its manufacturing process. Therefore, the economics of metal cutting process significantly affects the overall cost offinalproducts, and there is a strong drive to reduce the time and the cost of machining operations. In the past century, a great deal of research has been devoted toward understanding the mechanics of metal cutting, with the objective of obtaining more effective cutting tools and more efficient manufacturing process plans. Traditionally, these objectives have been achieved by experimentation and prototyping. In spite of extensive research in thisfield,the basic mechanics of the cutting process and the interplay of many factors which leads to its great variety is not yet totally understood, and the search for more effective models has continued on analytical, experimental and numerical fronts. Metal cutting is a complex process in which several mechanisms are at work simultaneously, and interact with each other. This process is greatly affected by material properties, cutting conditions, tool geometry and machine-tool dynamics. In practical machining operations, such as various turning and milling operations, the geometry and kinematics of the tool is very complex. For this reason, traditionally, a much simplified model of cutting, the orthogonal cutting, is used in the fundamental study of mechanics of this process. In orthogonal cutting, the cutting velocity is perpendicular to the tool cutting edge, and if the ratio of uncut chip 1 Chapter 1. INTRODUCTION 2 width is sufficiently small, it can be modeled as a 2-D plane strain problem. A 3-D counterpart of this process is called the oblique cutting, in which the cutting direction is inclined with respect to cutting edge. However, even for this most basic process of orthogonal cutting there are still many challenging issues to deal with. 1.2 DESCRIPTION OF ORTHOGONAL M E T A L C U T T I N G PROCESS In the cutting process, shown schematically in figure 1.1, a hard tool is pressed into a softer workpiece material, producing a chip and a machined surface. The material undergoes severe Figure 1.1: Schematic view of orthogonal metal cutting process plastic deformation in passing through a highly localized shear zone extending from the tip of the tool to the free surfaces at the juncture of chip and undeformed workpiece, called primary deformation zone. Depending on the cutting speed, the deformation may occur at very high strain rates. When the cutting tool moves through the workpiece material, the chip is.separated from the workpiece and slides over the rake face of the tool. The highly pressurized friction between the chip and the rake face causes further straining of the material in the area around the Chapter 1. INTRODUCTION 3 tool surface called secondary deformation zone. In this region, the chip initially sticks to the tool, but with the reduction of pressure further up on the rake face, it slides over the tool face elastically until it separates from the tool and curls away. The plastic deformation and friction generate large amount of thermal energy in the cutting zones, which raises the temperature in the cutting zones significantly, sometimes beyond the re-crystallization temperature of the metal. Such temperature rise has the effect of softening of the material, in contrast to the hardening effects of large strains and high strain rates. The process variables are dependent on cutting conditions and tool geometry which bring about large variety to the process. A successful cutting model should be able to incorporate such effects. A reliable model of material behavior at such extreme conditions is an essential input to cutting simulation. However, the strain, strain rate and temperature in the cutting process are orders of magnitude higher than what could be measured by most of existing measurement equipment, and thus, direct data about material behavior at conditions similar to what occurs in cutting process is hardly available. Many other phenomena may occur in a cutting process such as various chip types obtained in different cutting conditions, creation of built-up edges, tool wear and dynamic effects such as chatter. However, the basic mechanism of formation of the chip in the deformation zones and the interplay of various factors in this process remain at the core of the fundamental research on mechanics of cutting. In the last century, understanding of cutting mechanics has caught the attention of many researchers, and many analytical, empirical and numerical models have been contributed to the knowledge in this field. While all these models have advanced the common knowledge about the cutting process, the ultimate usefulness of a cutting model is in its ability to predict the process variables for a given set of cutting conditions with reasonable accuracy and cost. In particular, numerical modeling of this process has attracted the attention of many researchers in recent years, because of the better insight it can provide into the mechanics of chip formation by avoiding many of the simplifications needed in other approaches. There is a prospect that one Chapter 1. INTRODUCTION 4 day, such simulations may replace costly machining tests needed for tool design and process planning. However, the reliability of a numerical cutting model is dependent on two factors; the reliability of the input data in terms of material and frictional characteristics, and the correctness and efficiency of the numerical approach, formulation and procedures used. The latter issue is the subject of this thesis. 1.3 S C O P E O F THIS R E S E A R C H Traditionally, Lagrangian and Eulerian finite element approaches have been used for simulation of metal cutting process. In this thesis, it is argued that in spite of their many advantages, each of these approaches suffer from drawbacks that make it inadequate for such purpose. Therefore, an alternative approach which can combine the strengths of both approaches and avoid their drawbacks would be desirable for simulation of metal cutting process. It is shown here that the arbitrary Lagrangian-Eulerian (ALE) formulation offers such opportunity. In ALE, the FE mesh is neither attached to the material nor fixed in space, but can have an independent and arbitrary motion assigned by the analyst. Therefore, it is a more general approach which reduces to Lagrangian and Eulerian formulations as its special cases. With a proper assignment of mesh motion, the advantages of both Lagrangian and Eulerian approaches can be combined in a single analysis. In this work, an ALE approach is used forfiniteelement simulation of orthogonal metal cutting process. Chapter 2 provides a critical review of literature on numerical simulation of the cutting process, and presents the ALE as an alternative approach. Chapter 3 describes a large deformation ALE formulation which includes the effects of strain rate and temperature rise. After linearization and discretization of the equation of virtual work, an incremental finite element is created with particular emphasis on application to metal forming and cutting processes. A contact algorithm is designed to detect contact between the tool and the workpiece Chapter 1. INTRODUCTION 5 and apply frictional conditions. The deformation process is coupled with a heat transfer module which calculates the temperaturefieldin every increment of the numerical process. Chapter 4 discusses the numerical procedures used for assigning arbitrary motion of the mesh and for mapping material associated variables to the ALE mesh. The structure of the ALEfiniteelement code is also discussed in this chapter. Numerical results of cutting simulation under various cutting conditions are presented in chapter 5. By designing a robust mesh motion scheme for the workpiece, the process of chip formation in orthogonal metal cutting process with sharp and chamfered-edge tools is simulated. A parametric study of the effects of cutting conditions on the distribution offieldvariables and on tool forces is also presented. The numerical results are verified through comparison with experimental data of cutting forces and shear angles obtained under similar cutting conditions. Finally, chapter 6 presents conclusions and outlines the future work in this field. Chapter 2 M O D E L I N G O F C U T T I N G P R O C E S S IN T H E L I T E R A T U R E 2.1 A N A L Y T I C A L A N D EMPIRICAL M O D E L S OF C U T T I N G PROCESS Analytical investigation of the cutting process made much progress in the first half of the century. In 1937, Piispanen [1] presented the card model which depicts the shearing of the material as a deck of cards inclined to the free surface. In the early 1940's Ernst and Merchant [2] developed their classic shear plane model based on free body diagram of the chip which is held in equilibrium by two equal, opposite and collinear resultant forces, one acting on shear plane, and the other on the rake face of the tool. Figure 2.1 depicts the famous Merchant's circle which shows the relationship between various force and velocity components acting on the chip. In thisfigure,V , V , and V are chip velocity, shearing velocity and workpiece c s w velocity, respectively, R is the resultant force on the tool with its projections in shear and s in frictional directions, a , (j>, and 3 are rake, shear and friction angles, respectively, L is 0 c the contact length, and h and h are uncut chip thickness and chip thickness, respectively. c The shear plane model assumes that the deformation of the material ahead of the tool occurs instantaneously on a thin plane of concentrated shear, the shear plane, and that shear and normal stresses are uniform on this plane. Across the shear plane, velocity of the workpiece material is instantaneously changed to the chip velocity. They obtained the angle between the shear plane and cutting direction, the shear angle, by making the shear plane a direction of maximum shear stress; <j>=--?- * 4 6 2 + ^ 2 V (2.1) ' Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 7 Figure 2.1: Force and velocity components in Merchant's model Later, Merchant [3] showed that the same relationship could be obtained by making the shear plane as the direction of minimum energy for a rigid-perfectly plastic material. Despite showing relatively poor agreement with experimental results in some cases, Merchant's model formed the basis for much of the research that followed, and is still used as a quick check tool by machinists. Lee and Shaffer [4] applied the theory of slip-line field of ideal plasticity to the problem of machining and estimated forces, chip thickness and shear angle in terms of tool geometry, the coefficient of friction and the yield stress of the material. L i k e Merchant's model, the deformation is assumed to occur on the shear plane, but the plastic field is extended above this plane to form a triangular region. The shear angle predicted by this theory is 4>= L-f3 + a 7 0 (2.2) Many other slip-line field and upper-bound solutions have followed the above work, in which researchers have striven to improve the shear angle prediction by including realistic factors Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 8 which affect the cutting process, such as hardening, rate and temperature effects, stickingsliding friction and more complex tool geometry. The works of Shaw, Cook and Finnie [5], Johnson and Mellor [6], Oxley and Coworkers [7], and Armarego [8] are noteworthy in this area. Palmer and Oxley [9] used experimental flowfieldsobtained by tracing the path of individual grains in low speed cutting to determine the extent of the deformation zone and the shear angle. They showed that the plastic deformation zone has substantial thickness and the streamlines form smooth curves from the workpiece into the chip. They developed slip-line field solutions from the experimental results and calculated the stress distribution in the deformation zone. Oxley and Welsh .[10] developed a parallel-sided shear zone theory which allows for strain hardening of the material. In this model, the shearflowstress changes from initial yield stress at the lower boundary of the shear zone to a higher value at its upper boundary. Later, Stevenson and Oxley [11] developed more comprehensive models by using a power law flow stress equation whose coefficients were dependent on strain rate and temperature. These effects were combined using the concept of velocity-modified temperature, and it was suggested that theflowstress for a particular material at a given strain can be assumed to be a unique function of this parameter. Oxley [7] has laid out a predictive theory of orthogonal machining based on the above chip formation model. There has been attempts tofindshear angle solutions for 3-D oblique cutting process [12, 13]. Seethaler and Yellowley [14] presented an upper bound solution for oblique cutting tools with a nose radius. Ren and Altintas [15] obtained a slip-line field solution for cutting with chamfered edge tools. Jawahir and coworkers [16] have studied chip curling and breaking using analytical and numerical methods. 2.1.1 Friction and Temperature on the Rake Face The frictional behavior between the chip and the tool is an important aspect of the mechanics of chip formation. This behavior is determined from the nature of, and relationship between, Chapter 2. MODELING OF CUTTING Contact Length PROCESS IN THE LITERATURE 9 Contact Length Figure 2.2: Distribution of normal and frictional stresses on the rake face, left) Zorev's Model [17], right) results of Barrow, et al. [19] the normal and shear stress distributions on the rake face of the tool. Considerable amount of theoretical and experimental work has been done to determine the stress distributions along the tool rake face. The stress distribution obtained from photoelastic studies by Zorev [17] in 1963 have been widely referred to by researchers. This distribution, shown in figure 2.2, suggests that the normal stress decreases exponentially from a maximum value at the tool edge to zero at the point of separation of the chip. It has been widely accepted that the contact length on the rake face is divided into two regions. In the region starting from the tool tip, the normal stress is so high that the chip material sticks to the tool. In this region, known as the sticking zone, the shear stress of the material is equal to its shear strength and has a uniform distribution. In the second region which starts from the end of the sticking zone and ends at the separation point, the normal stress is sufficiently reduced so that the proportion of the real area to the apparent area of contact is relatively small. In this region, known as the sliding region, the chip slides over the tool face and the shear stress gradually reduces to zero. In the work of Usui and Takeyama [18] the normal stress peaked near the tool tip, became almost constant in the middle part, and then gradually reduced to zero at the point of separation. Barrow, et al. [19] used split Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 10 tool dynamometry to determine rake face stress distributions in cutting of steel under a range of cutting speeds and uncut chip thicknesses. From their tests, they concluded that both normal and tangential stresses are constant in the sticking zone, after which they decrease exponentially to zero (seefigure2.2). However, their method could not determine the distribution in the first 20% of the contact length close to the tool tip. Temperature distribution on the rake face is another important aspect of the cutting process, specially since it is directly correlated with the tool wear. Various temperature models have been proposed and experimental techniques such as infrared camera and thermocouple methods have been used to estimate the temperature at the rake face as well as on the shear plane. Stephenson [20] has reviewed and assessed the predictions of these models against experimental data. It is generally accepted that the temperature has a non-uniform profile on the rake face, with its peaks at some distance from the cutting edge. Also, the temperature in the primary shear zone increases progressively from the free boundary to the rake face, and depending on the thermal properties of the tool and the workpiece materials, there could be a large gradient of temperature across the thickness of the chip. The maximum temperature generally has an increasing trend with the cutting speed, and may contribute significantly to the shear banding in the chip in some materials. The analytical and empirical models of chip formation provide valuable insight into the cutting process. Nevertheless, their use in prediction of process variables in the deformation zone is somewhat limited due to many simplification which has to be made and the need for experimental data as inputs to the model. With the surge in computational power in recent decades, numerical simulation has become an attractive alternative to such models. The advantages of numerical simulation are in that they provide detailed information about the deformation zone, they usually require no experimental cutting data as input, sophisticated material and friction models may be easily implemented, and nonlinear boundary conditions can be modeled. However, reliable material and friction models are vital to the success of a Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 11 numerical process. Many orthogonal metal cutting simulation analyses have appeared in the literature in recent years. In spite of the many promises of numerical simulation, there are aspects of cutting process which make such endeavor particularly challenging. These aspects include very large and localized deformation resulting in severe mesh distortion, large unconstrainedflowover vast free boundaries, and changing boundary conditions. But more importantly, and unlike many forming processes, there is cutting action in which the bulk of the material is branched to form the chip. This action is in a sense similar to what occurs in an elastic-plastic crack propagation problem. Simulating such cutting action is perhaps one of the most challenging aspects of this analysis. In the rest of this chapter, these aspects of numerical simulation of cutting process will be discussed through critical review of the literature in this field. 2.2 R E V I E W OF NUMERICAL SIMULATIONS OF C U T T I N G PROCESS Finite element method has been used as an alternative to the analytical approach for modeling the cutting process since early 1970's. One of the first attempts at using the FEM in metal cutting research goes back to 1973 when Klamecki [21] simulated the incipient stage of the chip formation process. He used a 3D model to show a transition from plane strain condition at the center of the chip along its width to plane stress condition at the chip surface. During the 1970's, Tay, et al. [22, 23] used the FEM to predict temperature distribution in steady state orthogonal machining. They calculated plastic and frictional work at the cutting zone from experimentally obtained cutting forces and shear angle data, and used them as source terms for the heat transfer analysis. In 1982, Usui and Shirakashi [24] made the first comprehensive simulation of the cutting process which included steady state chip formation. They used a sophisticatedflowstress model, obtained from material tests as a function of strain, strain rate and temperature. For friction on the rake face, they used an exponential function which gives Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 12 the frictional stress in terms of normal stress and shear strength of the material, and provides a smooth transition between sticking and sliding regions (see equation 4.26). In this simulation, the analysis started from assumed initial estimates of shear angle and chip shape, and a small crack at the cutting edge was inserted. The crack were extended by an element length when the tool deemed sufficiently close to the crack front node. The strain, strain rate and temperature field were calculated, the latter through afinitedifference analysis, and the chip geometry and field variables were updated iteratively until reasonable solution was obtained. Iwata, et al. [25] pursued a similar approach with a rigid-plastic material model in their simulation, but a ductile fracture criterion was used for chip formation and frictional stress was based on the Vickers hardness of the material. They compared the results of their simulation with cutting experiments carried out in scanning electron microscope (SEM) at very low speed and reported good agreement between results. In 1985, Strenskowski and Carroll [26] developed the first cutting simulation which started from incipient stage and continued until steady state was achieved, in a Lagrangian framework. They ignored the rate effects on material model and their thermal model did not include conduction in the chip. Simple Coulomb friction law was used and a crack inserted in front of the tool tip, which was extended when the effective plastic strain at the neighboring node reached a critical value. Carroll and Strenkowski [27] later developed another FE model which was based on an Eulerian frame of reference, and only included steady state cutting. A preassumed chip geometry had to be defined, based on trial and error. Later, Strenkowski and Moon [28] improved the latter model by automating the procedure for predicting and updating of chip geometry, based on the premise that the chip velocity should not have a component normal to the free boundary of the chip. Many other simulations have appeared in the literature in the last fifteen years. They vary in the numerical formulation they use, the type of chip they model, the material and frictional behavior they consider, and the features of cutting process that they include. The material Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 13 responses used in the literature range from simple rigid-plastic or elastic-plastic behavior to more complex models such as thermo-viscoplastic ones. Ueda and Manabe [29] have used the famous Gurson's yield function which includes the effect of hydrostatic pressure on the yield surface. Most of the numerical analyses have only modeled continuous chip formation. However, some attempts at modeling discontinuous chip formation have also appeared in the recent literature. Hashemi, et al. [30] modeled segmented chip formation by developing a fracture algorithm which extended the crack by splitting nodes in the direction dictated by the maximum principal stress. They used very coarse mesh and an elastic-plastic material model. Obikawa and Usui [31] modeled formation of serrated chips due to ductile cracks emanating from free surface of the chip in the machining of titanium alloys. They used a fracture criterion based on a critical value of effective plastic strain, and extended the crack in the direction of maximum shear stress by node-splitting technique. Marusich and Ortiz [32] modeled all types of chips by incorporating criteria for both brittle and ductile fracture, mesh adaptation technique, and an algorithm for propagation of the crack in an arbitrary direction as dictated by the process. Xie, et al. [33] analyzed the formation of shear localized chips in machining by developing a flow localization parameter. Furthermore, In an attempt to model the effect of built up edge on the cutting process, Komvopolous and Erpenbeck [34] add a small triangular nose to the tool. From another respect, almost all the cutting simulations have focused on the orthogonal cutting process, but recently there has been a few attempts at modeling 3-D cutting process [35]. Hushimura, et al. [36] also simulated burr formation in oblique cutting. Another aspect of the various metal cutting simulations in the literature is the type of numerical formulation used in these analyses. It is this aspect of cutting simulation which is the focus of this research, because of the great implications it has in the success of such analyses. From this point of view, two distinct formulations have been used for modeling this process; the Lagrangian approach and the Eulerian approach. In the Lagrangian formulation, the natural Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 14 approach for solid mechanics problems, thefiniteelement mesh consists of elements which are attached to the material and deform with its deformation. This approach is particularly efficient when unconstrained flow of material is involved, because the evolution of material boundaries is simulated accurately by the deformation of the mesh. Also, the formulation and implementation of this approach is easy and straight-forward. The Eulerian approach, on the other hand, is mostly used forfluidflowproblems involving a control volume. In this method, the mesh consists of elements which arefixedin space and cover the control volume, and the material properties are calculated at fixed spatial locations as materialflowsthrough the mesh. This approach has also been used to model large deformation of solids, mostly in steady state metal forming analyses. Due to the fixed nature of the mesh, it is most suitable for cases where the boundaries of the material region are known a priori, and deforming free boundaries are minimal. Examples of these cases are closed die forging and extrusion problems. Lagrangian formulation has been extensively used in simulation of cutting process, and a number of Eulerian analyses have also appeared in the literature. 2.2.1 Review of Lagrangian Cutting Analyses One of the early Lagrangian developments was introduced by Hibbit, Marcel and Rice [37] for finite strain problems. They introduced a formulation which is commonly known as Total Lagrangian Formulation (TLF) in which the variables are referred to the initial configuration of the material. Later, McMeeking and Rice [38] pioneered the use of Updated Lagrangian Formulation (ULF) which refers the variables at time t + Ai to the configuration at time t. Both methods have been widely used for plasticity analysis, and are implemented in most of the available commercial FE programs. If strains are relatively small, the Lagrangian formulation is easy to implement, efficient, and fast converging. However, difficulties arise when it is used for problems involving large deformation, nonlinear boundary conditions which change in the Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 15 course of deformation, and when mesh distortion and element entanglement are critical factors in the analysis [39]. Unfortunately, most of the metal forming processes fall in the above category, with metal cutting being one of the most extreme cases. In simulation of metal cutting process, what makes Lagrangian approach particularly attractive is the fact that it can easily model the large unconstrained flow of material and natural formation of chip into its steady state shape. In other words, no a priori assumption about the shape of the chip is required, because it naturally evolves to its proper thickness. All stages of cutting from indentation to incipient stage to steady state can be modeled in this approach. Nevertheless, there are major drawbacks in applying Lagrangian approach in a cutting simulation. One such drawback is that due to the highly localized nature of the cutting process, severe distortion of the FE mesh in front of the tool tip occurs, requiring frequent remeshing. But, perhaps the most important drawback of Lagrangian approach is in the simulation of cutting action at the tip of the tool. The formation of chip in this approach has been traditionally modeled by separation of nodes in front of the tool tip, similar to the procedure commonly used for self-similar crack propagation in fracture mechanics. In this technique, based on a separation criterion, the node in front of the tool tip is successively split into two nodes, one moving up on the rake face to form the chip, and the other moving beneath the tool to form the machined surface. A schematic view of such process is shown in figure 2.3. Despite its simplicity, the node separation method suffers from many problems which significantly affect the results of simulation. Firstly, imperative to using a node splitting method is that nodes are lined up in the path of tool tip along a pre-defined parting line that will be "unzipped" as the tool advances. Beside the fact that this requires a priori assumption, it sometimes creates problems as the parting line may be pushed out of position during the course of deformation and thus the target nodes will not be in front of the tool tip at the time of separation. In addition, the use of parting line effectively restricts the analysis only to cutting with sharp edge tools. Secondly, the successive separation of nodes often creates a crack which Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 16 7M \ V\ / ' - - Figure 2.3: Node separation in Lagrangian analysis moves ahead of the tool. Such a crack is unrealistic because it is not observed in actual cutting of ductile materials. Due tofinitelength of elements in the FE mesh, the continuous process of cutting will be modeled as a discrete process of successive jump of crack tip from one node to the next one. As a result, deformation history of the material at the vicinity of the tool tip is incorrect. Furthermore, significant oscillation in the value of cutting forces will occur, because every time a node is separated, there is a sudden reduction of the forces. The gap also shifts the points of maximum stress from the tip of the tool to a point higher on the rake face, affecting the predicted shear angle. Another drawback of node separation method is the creation of unbalanced forces when a node is released. Furthermore, the application of this technique in cutting with large negative rake angles or with chamfered or blunt tools is very doubtful. If the tool edge is not sharp, there is no definite parting line along which nodes may be separated. Another important aspect of the node separation method in a Lagrangian approach is the necessity for defining a criterion based on which a node in front of the tool tip is allowed to separate. Many criteria have been used in the literature ranging from simple geometric or strain criteria to more complex fracture mechanics criteria. Even if the type of criterion is properly chosen, there is no physical indication as to what criterion value would be appropriate. All of this makes the separation criterion more of an arbitrary nature, with great effects on the results. Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 17 A critical look at the literature shows that while most researchers have used this approach, there is no consensus on the appropriate criterion. Usui and Shirakashi [24] use a geometrical criterion which splits the node in front of tip into two nodes when it is sufficiently close to the tool tip. Komvopoulos and Erpenbeck [34] state that a plausible criterion should be based on either a distance tolerance or a critical deformation parameter indicative of the extent of plasticity. In view of uncertainties in the value of a strain criterion, they use a distance criterion. A parting line consisting of pairs of superimposed nodes is defined. Each pair is constrained to have identical displacements, but when the distance of the pair to the tool tip is smaller than a specified tolerance, the nodes are allowed to move independently. They point out that the magnitude of tolerance is crucial, because if the tolerance is too small, convergence problems will arise due to excessive distortion of elements closest to the tool tip, but on the other hand the tolerance should be small enough to minimize the gap at the tool tip and the corresponding error. The value of tolerance is in essence dependent on the element size in the tip area. They use a fixed value of 35/Lim in their simulations. Zhang and Bagchi [40] use 2-node conditional link elements along the parting line, which produce zero stiffness upon separation based on a distance criterion. They have adopted a distance tolerance of (0.1 â€” 0.3)L where L is the length of solid elements in front of the tool tip. In a later paper [41], they propose a criterion which uses the ratio of separation distance to uncut chip thickness. Shih [42] too uses a distance criterion, but with different values chosen arbitrarily. The geometric criterion is not based on any physical condition of the deformation in the cutting zone. Therefore, alternative criteria based on physical measures are proposed in the literature. Many researches have used a critical value of effective plastic strain at the node closest to the tool tip as the separation criterion. Carol and Strenkowski [27] employ a range of strain values between 0.25 and 1, and conclude that while the criterion value does not affect the chip shape and forces significantly, it has a profound effect on the residual stresses in the workpiece. In view of this, Strenkowski and Mitchum [43] attempted to obtain a physical Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 18 justification for the value of strain criterion. They simulated indentation and incipient stage of the cutting and reasoned that this value should be equal to the value of effective plastic strain at the time when the indentation transforms to incipient stage. They obtained a range of strain criterion values of [0.4 â€” 0.7] for a range of uncut chip thicknesses for 70/30 Brass. Other works which have used a strain criterion include Xie, et al. [33] who use a critical value of 0.5, and Hashemi, et al. [30] who adopt a range of 0.6 to 1.5. Another alternative for node separation criterion is the use of a measure of fracture or damage in the material. Iwata et al., [25] employ a stress based ductile fracture criterion which takes stress history into account as chip formation criterion. Lin and Lin [44] argue that the chip formation criterion should be a material constant, independent of cutting conditions. They propose a criterion based on strain energy density which is obtained from tensile test. However, as Chen and Black [45] point out, the critical strain energy density determined from uniaxial tensile test cannot be representative of conditions in a cutting operation. Huang and Black [46] conduct simulations with different types and magnitudes of chip separation criteria, including criteria based on distance and shear stress in the element ahead of the tool. They conclude that while the type of criterion is not crucial in determining the chip shape and process variables, its magnitude has a great effect on the chip separation process and distribution of stress and strain in the chip and machined surface. They also propose that a combinations of geometrical and physical measures be used for modeling both incipient and steady-state cutting stages. In view of the disagreement among the researchers about the type and magnitude of separation criterion, and in spite of attempts at justifying the choices in this regard, it is perhaps fair to say that all of these criteria are arbitrary and artificial. In the opinion of this author, this problem stems from the inherent inadequacy of the separation technique for cutting process, because it is not representative of the real cutting action which occurs in front of the tool tip. It is noted that formation of chip is essentially defined by what happens at the tip of the tool, and there seems to be a singularfieldof strain and strain rate right at that point. The gap which is Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 19 created by the node separation at the tool tip practically obscures this singularfield,and thus, is essentially unrepresentative of the cutting mechanism. The fact that the node separation technique cannot simulate cutting with large negative rake angles or with chamfered edge tools is further indication of the inadequacy of this approach. Of course, some of the aforementioned problems may be eased by increasing the density of the mesh around the tool tip, but a finer mesh gets distorted more quickly, and thus, the possibility of numerical instabilities due to mesh distortion increases. The difficulties with node separation method have prompted some researchers to look for alternative methods for simulating the cutting action within a Lagrangian framework. Ceretti, et al., [47] use a scheme that deletes the elements when they reach a critical value of accumulated damage. Such scheme entails partial loss of volume which violates the law of continuity. However, this can be minimized by using a very dense mesh around the tool tip, which in turn, increases the frequency of mesh distortion and remeshing. Besides, sometimes elements at positions other than those at the tool tip may be wrongly targeted for deletion, creating gaps in the material. More recently, the difficulties with node separation technique are overcome by simulating the cutting action as one of continuous indentation and flow of the material around the tool tip. Based on some experimental evidence, Madhavan, et al. [48] argue that machining of ductile materials is indeed a special type of indentation with an asymmetric wedge, and that the material "separation" is due to plasticflowaround the tool rather than tensile rupture implied by node separation technique. Based on the indentation analogy, cutting as a process of pure shear may be simulated by forcing the tool into the material in small increments, causing flow of the material around the tool tip. The difficulty in using this approach in a Lagrangian analysis is the large and frequent distortion of FE mesh due to very large deformation of material in front of the tool tip. In fact this seems to be the reason why researchers have resorted to node separation technique in the first place. With the surge in computational power in recent years, Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 20 it has been made possible to remedy the mesh distortion problem by frequent and automatic remeshing of the analysis domain. Madhavan et al. [48] use a scheme that triggers a remeshing when effective plastic strain exceeds a value of 0.25 in any point in the workpiece. Sekhon and Chenot [49] use similar approach in which the FE mesh is checked frequently and the mesh is automatically regenerated if a distorted element is detected. Marusich and Ortiz [32] present a thorough analysis of chip formation which also includes fracture criteria for segmented chip as well as modeling of shear localization. They use a mesh adaptation technique which refines the elements when their plastic power contents exceed a prescribed tolerance. Ng, et al. [50] define a geometric function called "local penetration" rate for triggering adaptive remeshing. This function is indicative of the amount of penetration of tool into the workpiece. Although the indentation approach overcomes many of the problems pertinent to node separation approach, it comes at a high computational cost. First, the mesh has to be very fine around the tool tip, or else the material will overlap the tool in between the nodal points, and subsequent remeshing of the workpiece will amount to loss of overlapped material. Second, due to highly localized deformation and the small size of elements in afinemesh , the mesh has to be regenerated quite frequently. It should be noted that following each remeshing, the process variables should be transferred from old mesh to new mesh, often by interpolation. Such variable transfer will inevitably introduce errors which may accumulate rapidly, and deteriorate the results if many remeshing steps are required in a single simulation run. Overall, the computational cost of this approach is sometimes prohibitive [32]. 2.2.2 Sample Lagrangian Simulations i) Lagrangian approach with node separation To demonstrate some of the problems pertinent to the Lagrangian approach with the node separation method, a series of Lagrangian simulations has been carried out. The material Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 21 is assumed to have an elastic-plastic behavior with linear isotropic strain hardening, and the effects of strain rate and temperature were neglected in this example. The Young's modulus, hardening modulus, and initial yield stress of the material are 200 GPa, 1 GPa, and 414 MP a, respectively. A shear friction law was used on the rake face which saturates the frictional stress at local shear strength of the material in the region of sticking friction. Two different tool rake angles of 10Â° and 30Â° were simulated and a small uncut chip thickness of 0.025 mm was considered. Nodes were lined up along a parting line in front of the tool tip. A critical distance criterion was used and the criterion value wasfinetuned in trial runs tofindthe best trade off between having a smaller crack ahead of the tool and avoiding the distortion of elements and excessive shifting of the parting line. The typical value of distance tolerance was less than half of the length of elements along the parting line. To compensate for the change in the equilibrium resulting from node separation, the nodal forces immediately before the separation were added to the two nodes created by splitting the released node, and equilibrium equation is recalculated. The added forces were then gradually reduced in the course of a few increments until the nodes on the chip side came into contact with the rake face. Figure 2.4 shows a simulation with a 30Â° rake angle tool. Contours of effective plastic strain are also shown in the figure. It can be seen that a crack is moving ahead of the tip, and that due to shifting of the point of maximum stress further up on the tool, the shear plane does not emanate from the tool tip. The gap created at the tip of the tool essentially obscures the singularity of strain at the tip of the tool, and the maximum value of strain calculated is less than 2. Figure 2.5 illustrates the results of a simulation with a 10Â° rake angle tool. Although the critical distance is fine-tuned to get the best result, it is seen that the machined surface is unrealistically wavy due to shifting of the parting line during the process. Choosing a larger distance tolerance will decrease this waviness at the expense of a larger gap. Figure 2.6 shows the cutting and thrust forces (i.e., forces in the direction of cutting velocity and its perpendicular direction, respectively) on the tool in a typical Lagrangian simulation. The oscillation in the forces is due to successive separation of Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 22 Figure 2.5: Lagrangian analysis with node separation; cutting with a 10Â° rake angle tool Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE Cutting O.OS 0.10 Tool 0.15 23 Force 0.20 0.25 0.30 Displacement(mm) Figure 2.6: Forces in Lagrangian analysis with node separation nodes in front of the tool tip. Every time a node is split the force equilibrium in the workpiece is changed and a sudden drop in the forces occurs. Then, the forces gradually increase until the next node is released. ii) Lagrangian simulation with indentation approach The performance of Lagrangian simulation with the indentation approach has been studied in another series of cutting simulations, carried out using the commercialfiniteelement code DEF0RM-2D [52]. This software has an automatic remeshing module which regenerates the mesh based on a given criterion. In these simulations, the tool was advanced into the workpiece at a constant speed to produce the chip. However, the entire mesh had to be regenerated frequently, because the quality of elements in front of the tool deteriorated due to excessive distortion. The fast and frequent deterioration of the mesh was mostly due to penetration of the tool tip in between the nodes of neighboring elements of the workpiece. In each remeshing step, a part of the workpiece material which had overlapped the tool due to such penetration had to be cut out. In addition, process variables such as stresses and strains had to be interpolated from the old mesh to the new mesh. The workpiece material in these examples was AISI1045 cold steel for which the flow stress is available in the software database. Thermal effects in Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE X Load vs Straka H A W f 1 J J.MD Y Load va Straka 1-000' I f too. . .l A -0.4C0XO ^-j 1-2BO L! r ii 1 LjL.l. m sal. 1 Ii i ' Figure 2.7: Lagrangian analyses with indentation using D E F O R M - 2 D code Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 25 the chip was neglected. Figure 2.7 shows two simulations, one with a rake angle of 30Â° and uncut chip thickness of 0.5 mm, and the other with a rake angle of 0Â° and uncut chip thickness of 0.2 mm, and both with 120 m/min cutting speed. The automatic remeshing module is set up in a way that the mesh around the tool tip isfinerthan other regions. It can be seen in the figure that the tool penetrates the workpiece in between its nodes, and the mesh is distorted in this region. In the course of analysis, the entire workpiece had to be remeshed in intervals of around 10-20 incremental steps, and still considerable volume of material had to be cut out due to overlapping with the tool in every remeshing. Here again, there is a trade-off between havingfinermesh to reduce penetration and volume loss and having coarser mesh to decrease the frequency of automatic remeshing steps required. The problem is more acute for a sharp tool edge as it penetrates between the nodes easier. A typical graph of forces in the cutting and thrust directions is also shown in figure 2.7. The large oscillation in the forces corresponds to occurrences of remeshing which changes the balance of forces due to material cut out and variable transfer. 2.2.3 Eulerian Simulation of Metal Cutting Process In Eulerian approach, on the other hand, none of the problems pertinent to Lagrangian formulation exist in modeling the cutting action, because this action is modeled as flow of the material around the tool. In this respect, it is similar to the indentation approach. However, since in this analysis the FE mesh isfixedspatially, no mesh distortion occurs, and no remeshing is required. In addition, because the density of the mesh in Eulerian approach is dictated only by the expected gradients of stress and strain, this approach is more efficient computationally and performs well in the region around the tool tip, at least for ductile materials. This approach, however, is not suitable for modeling the unconstrained flow of material at free boundaries. The fact that the mesh is spatiallyfixedmeans that a priori assumption should be made about Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 26 the shape of the chip and about contact conditions. In other words, because of the fixed mesh, the contact length and the free boundaries of the chip cannot develop as natural results of deformation. To overcome this difficulty, researchers have resorted to either trial and error to find the optimal assumption for chip size, or to adjustment of the boundaries outside the solution domain. In the latter approach, following Zienkiewicz et al. [53], an iterative procedure is commonly used which checks the velocities at free boundary nodes to ensure that their components normal to the boundary are zero and the boundary is a streamline. This requires starting the process with an initial assumption about the chip size, followed by progressive adjustment of the mesh, and interpolation of variables. Strenkowski and Moon [28] have automated the process of mesh adjustment. They report that the convergence of the chip geometry is fairly rapid. Wu, et al. [54] and Kim and Sin [55] have used the same technique for updating chip geometry with a thermo-viscoplastic material model. It is noted that the Eulerian approach can only model the steady state process, and modeling discontinuous chip formation in this approach may not be possible. Another issue in using an Eulerian formulation is the presence of convective terms which makes material time derivatives much more difficult to handle. Material points have to be traced back to thefixedmesh points and history dependent variables such as stresses and strains, calculated for material points, have to be interpolated in an approximate way to grid points. Many attempts have been made to achieve this. For example, Abo-Elkhier, et al. [56] use a method in which a hypothetical mesh is generated using the displaced material points and an interpolation using element shape functions is performed between the hypothetical mesh and the fixed Eulerian mesh. Such process is computationally intensive, introduces interpolation errors and creates convergence problems. In addition, most of Eulerian cutting analyses appeared in literature use theflowapproach, introduced initially by Zienkiewicz and Godbole [57], in which the material is modeled as rigid-viscoplsatic. This means that residual stresses in machined surface cannot be predicted. The above-mentioned problems make the Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 27 implementation of an Eulerian method more difficult than Lagrangian method. It may be concluded that the Eulerian approach is more suitable for processes whose external boundaries are known a priori, where the concept of a control volume is warranted. For a process such as metal cutting which involves large unconstrained boundaries, the shape and size of the chip and the contact length are important aspects of the solution, and adjustment of the boundaries outside the finite element solution seems to be an undesirable intrusion in the solution procedure. 2.3 ARBITRARY LAGRANGIAN-EULERIAN FORMULATION 2.3.1 H i s t o r y of A L E A p p r o a c h It is concluded from the above discussion on the strengths and drawbacks of the two classical numerical approaches that neither the Eulerian, nor the Lagrangian approach is well-suited for modeling processes involving large deformation and unconstrained flow of material. In fact, the strong points of each method are the weak points of the other, and the two methods complement each other in this respect. Since the deformation of material in forming processes exhibits the characteristics of both solid and fluidflows,an approach which can combine the advantages of both methods in a single analysis would be an ideal approach for modeling such processes. The arbitrary Lagrangian-Eulerian (ALE) description provides such analysis tool, and has acquired wide acceptance for being able to perform many difficult analyses. The arbitrary Lagrangian-Eulerian (ALE) analysis is a general formulation in which nodal points of the FE mesh are neither attached to the material points nor are they fixed in space. Instead, a typical mesh point can in general have an arbitrary motion as desired by the analyst. With this definition, the Lagrangian and Eulerian motions will be special cases of such arbitrary motion. Therefore, a typical ALE mesh may include degrees of freedom which are Eulerian, Lagrangian, or arbitrary (ALE) type. Theflexibilityof ALE description in adaptation of the finite element mesh provides a powerful tool for performing many difficult analyses especially Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 28 inflowproblems. The concept of ALE analysis has been born out of the need for suchflexibilityin the fluidflowproblems. This concept was first proposed by Noh [58] in the context of a 2-D finite difference method for hydrodynamics problems with movingfluidboundaries. It was later introduced into finite element method in the seventies, originally for problems involving fluid-structure interaction arising in nuclear safety analysis. Kennedy and Belytschko [59] and Donea, et al. [60] developed ALE formulations for inviscid compressiblefluidsin which the fluid domain was modeled using ALE while the solid domain was handled using the usual Lagrangian description. Hughes, et al [61] elaborated on basic concepts of the ALE analysis and presented the general kinematic theory of this description, in particular the relationship between the rate of change of physical variables in the material and referential domains. They applied their formulation to viscous incompressible and free surfacefluidflowproblems. The ALE concept was adopted for solid mechanics applications in the early 1980's. Some of the early attempts in this respect are directed toward linear-path independent material deformations. Haber and Abel [62] introduced a displacement-based "coupled Lagrangian-Eulerian" formulation for solving elastic contact problems in which the mesh adaptivity of ALE is used to accurately represent the slip compatibility conditions betweenflexiblebodies in contact. Haber and Koh [63] applied this approach to the problem of elastic crack propagation, in which the ALE mesh is conveniently moved to accommodate the moving crack tip. In the last decade, many attempts have been made to use ALE method in the analysis of history-dependent materials and in particular to metal forming processes. The suitability of ALE analysis for forming operations stems from its ability to side step problems of mesh distortion and entanglement, to handle complex boundary conditions, and to model unconstrained transient flow. However, when ALE is applied to history-dependent materials, the fact that the mesh points do not represent the same material point throughout the analysis presents a challenging issue in the implementation of this approach. In this respect, the ALE analysis is similar to Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 29 Eulerian analysis. In this context, the migration of material points with different histories through thefiniteelement mesh should be traced, and material associated variables such as stresses, strains and yield surface parameters should be updated from material points to mesh points. Liu, et al. [64] developed a general nonlinear ALE formulation which included both material and geometrical nonlinearities, and applied it to modeling of elastic-plastic wave propagation. Huetink, et al. [65] presented an early attempt to apply ALE to metal forming operations. They termed their approach "combined Eulerian-Lagrangian formulation" in which the material velocities were obtainedfirstfrom ALE equilibrium equations, followed by mesh movement and updating of variables. Liu, et al. [66] presented an ALE approach for forming processes with emphasis on the rate of external work on frictional interface, and simulated metal rolling problems. Other applications of ALE in forming processes can be found in the works of Schures, et al. [67], Benson [68], Ghosh and Kikuchi [69] and Huerta [70]. A review of ALE formulations in forming applications may be found in reference [72]. From the point of view of handling the convective terms arising from migration of material points through the mesh in the ALE analysis, ALE formulations in the literature may be cast into two different categories. First, an operator-split approach is used [60, 65, 68] in the sense that the deformation in every step of the analysis is decoupled into a Lagrangian step followed by an ALE step. In the second category, the fully coupled equations of motion involving both material and mesh velocities are solved [61,64,69,72]. The obvious advantage of the operatorsplit approach is in its simplicity and computational convenience, but this comes at the expense of accuracy. In effect, the operator-split approach is equivalent to an updated Lagrangian formulation with continuous remeshing. In this respect, this approach may be called pseudoALE, because the material motion and mesh motion do not occur simultaneously. In fact, some authors suggest that the convection part is performed only once every few steps [68, 74] which is a clear indication of the link between this approach and the Lagrangian approach with frequent remeshing. In contrast, the coupled solution is the real ALE formulation, because it Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 30 represents true kinematic description of equilibrium at each step, and thus, it is expected to be more accurate. It is appropriate at this point to outline the required characteristics of a rigorous ALE formulation; â€¢ since the ALE description of motion is general, it should be reducible to Lagrangian and Eulerian descriptions as special cases. â€¢ the ALE formulation involves two sets of nodal velocities; the material and the mesh velocities. Only one set, the material velocities, is obtained from equilibrium equations, whereas the mesh velocities are defined by the analyst through an automated scheme which decides the motion of each degree of freedom of the mesh at the start of each incremental step. Such scheme should be robust to generate the best mesh, and at the same time computationally simple and efficient to reduce the cost of an ALE step in comparison to a similar Lagrangian step, so that the use of ALE is justifiable. â€¢ a consistent, stable and efficient updating scheme for updating material associated properties must be included in the ALE procedure. 2.3.2 A L E Applications in Metal Cutting Simulation In recent years, a number of ALE models of metal cutting process have appeared in the literature. Hashemi and Stefan [73] have used an ALE hydrodynamic finite difference code, MESA, to model the chip formation in orthogonal cutting process from transient to steady state. An operator split approach divides the ALE step into Lagrangian and advection phases, with the latter phase being carried out once every four time steps. The effect of variation of rake angle on chip formation is studied. Olovsson, et al.[74] also use an operator-split ALE approach in modeling steady state orthogonal cutting problem. Similar to reference [73], the Lagrangian step is carried out Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 31 separately, and the advection phase is performed only every few steps. Furthermore, they insert a crack in front of the tool edge and simulate the chip formation process by continuously moving the node at the crack tip. In order to minimize the loss of material due to crack tip motion, they introduce a special crack tip element which conforms the element to the tool rake face. It may be noted that neither the approach used is truly ALE, nor the adaptivity of ALE is used in modeling theflowof material around the tool. It is not clear why an unrealistic crack has to be inserted in front of the tool within an ALE analysis, since it will suffer from some of the drawbacks of the node-separation technique. In their work, the difficulties in re-mapping of the solution variables from Lagrangian step to the ALE mesh is discussed, and a linear least square interpolation is used for this purpose. An arbitrary damping factor is introduced for variable re-mapping which compromises accuracy in the interest of keeping the monotonicity and improving convergence. In contrast, Rokotomalala and his coworkers [75,76,77] use a truly coupled ALE approach and make use of adaptivity of the mesh in ALE to simulate chip formation in the cutting analysis. They present the ALE form of the conservation laws and solve the mass and energy using finite volume method, while the momentum equation is solved using FEM. As their mesh motion scheme, they use an ad-hoc rezoning algorithm which uses geometric criteria to minimize the distortion of elements. The procedure for re-mapping solution variables is not discussed in these articles. The cutting simulation is started from a well-developed chip shape of arbitrary thickness, and the thickness and contact length is updated in the ALE process. However, the growth and curling of the chip is not obvious in the model, and the contact algorithm for updating contact surface is not discussed. Finally, in reference [71], the merits of using ALE for metal cutting simulation are discussed in contrast to Lagrangian and Eulerian approaches. Chapter 2. MODELING OF CUTTING PROCESS IN THE LITERATURE 2.3.3 32 The A L E Work in This Research In the light of advantages and drawbacks of both Lagrangian and Eulerian formulations, an ALE approach is used to simulate the chip formation process under various cutting conditions in orthogonal metal cutting applications. The original formulation presented by Wang and Gadala [72] is extended to include strain rate and thermal effects. Frictional load terms are added to the equation of virtual work and algorithms for implementation of contact constraints for soft-soft and rigid-soft multi-body contacts are developed. A heat transfer analysis is coupled with the stress analysis to calculate the temperature field in the tool and the workpiece, and update material properties at every step of the analysis. The transfinite and isoparametric mesh generation methods are used as the mesh motion scheme, and a procedure for applying tangential mesh motion on free boundaries is described. The mesh velocities in the finite element equations are handled on element level, using a static condensation scheme, which reduces computational cost. Algorithms for updating material-associated process variables are also presented. Chapter 3 THEORY OF A L E FINITE E L E M E N T METHOD In this chapter, the theory of arbitrary Lagrangian-Eulerian formulation is discussed and various algorithms for its implementation in the analysis of rate and temperature dependent metalworking processes are presented. 3.1 ALE FORMULATION The arbitrary Lagrangian-Eulerian description of material motion is an extension of the two commonly used analysis viewpoints, the Lagrangian and the Eulerian, which are used to describe the motion of the material from a given reference time to the time under consideration. In a Lagrangian viewpoint, the deformation of specific material particles of the continuum is traced between the configuration at the reference time and the configuration at the current time, whereas in the Eulerian viewpoint, attention isfixedon particular spatial points occupied by the continuum at any given time. When the material undergoes deformation from initial time to time t, the configuration of the continuum changes from Â°R to R. These two configurations t are usually referred to as material and spatial domains, respectively. Throughout this section, notation similar to reference [78] is used where a left superscript indicates the configuration at which the quantity occurs, whereas a left subscript indicates the configuration with respect to which the quantity is measured. If the quantity is measured with respect to same configuration at which it occurs, the left subscript is dropped. A variable with no super- or subscript indicates an increment of the quantity from time t to time t + At. Also, standard indicial notation is used in which a right subscript denotes components of a tensor and repeated indices represent 33 Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD summation over the admissible range. Furthermore, the notation 34 stands for A material particle occupying point P is identified by position vector Â°Xi = (Â°xi, Â°x , Â°x ) 2 s in the material coordinate system. In the deformed configuration at time t, the same particle is located at point Q and has the position vector *Xi in the spatial coordinate system (see figure 3.1). The motion of the material between points P and Q is defined by the velocity vector Configuration at time <+A/ Configuration at timet Figure 3.1: Material and mesh motion in ALE analysis d Xj l (3.1) Vi In a Lagrangian description, the initial configuration of a material point is mapped onto its current configuration and the motion of the particles is described by x = f {Â°x ,t) (3.2) t i i j whereas in the Eulerian description, a particle occupying position Xi in the current configuration l Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 35 is traced back to its original position in the initial configuration, i.e. Â°x = g ( x ,t) i i i j (3.3) When the continuum is discretized intofiniteelements, the Lagrangian description corresponds to a mesh which is attached to the material and deforms with it, while the Eulerian description corresponds to afixedmesh through which the material flows. In the arbitrary LagrangianEulerian description, the computational mesh is used as a reference system independent of the material deformation. In other words, in ALE description, the motion of the particles is referred to a computational reference system which is neither adhered to the material norfixedin space. Rather, the reference system may assume an arbitrarily assigned motion of its own. Let the referential coordinate Xi indicate the marker of a mesh point, e.g. afiniteelement node. Then, t the particle motion in space is described by x = x { x ,t). t t i i j (3.4) Then, the components of the mesh velocity are expressed as d xl vi = ~ | (3.5) Equation ( 3 . 4 ) maps the grid configuration onto the spatial configuration. Although the motion of the mesh in ALE description is arbitrary, there should be a one-to-one mapping between the material and referential domains, i.e. there should exist the inverse mapping Xi = Xi[%,t) t (3-6) The necessary and sufficient condition for the existence of such inverse mapping is that the determinant of the Jacobian of the transformation is nonzero, i.e., Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 36 This condition should be guaranteed by the scheme which governs the motion of the mesh. Furthermore, the boundaries of the two domains should coincide with each other at all times. This requires that no material should cross the boundaries of the referential domain. Mathematically, this condition may be expressed as (vi â€” rn = 0 (3.8) wherera;are the components of the boundary normal vector. Equation (3.8) states that no convective velocity occurs in the direction normal to the boundary, i.e., the boundary should be Lagrangian in the normal direction, while there is no restriction on the motion of the mesh tangent to the boundary. It is possible to derive the governing equations of the motion with respect to the referential system. Since in the conservation laws of the continua, the material time rate of a function f(x, t) is needed, it is necessary to express this quantity with respect to the referential system. Denoting the time rates with respect to material and referential systems with "dot" and "prime" superscripts, respectively, the relationship between the two derivatives is given by [64]; (3.9) where Wi = v; â€” v~i is the convective velocity of the material particles relative to the mesh in the referential system. Equation (3.9) gives the relationship between the material associated quantities and mesh point associated quantities and makes it possible to trace the history of material deformation in the ALE description. 3.1.1 Governing Equations in A L E Analysis The governing equations of the continuum are expressed in terms of laws of conservation and the constitutive equation. The conservation laws of mass, momentum and energy are expressed in strong form in the spatial coordinate system, respectively, as: Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 37 i> = ~ p f | (3.10) pvi = -~ + pbi (3.11) OXj dvi dqi pe = try- h pr - â€” OXj OXj (3.12) where p is the density, o-ij are the components of the Cauchy stress tensor, b{ are the components of body force vector, e = e(p, T) is the specific internal energy, r is the internal distributed heat source, qi are the components of heatfluxvector, and T is the temperature. The constitutive equation is usually expressed in terms of rate of the Cauchy stress tensor. In large deformation analyses, an objective stress rate should be used to maintain the objectivity of the constitutive relations. One such rate is the Jaumann stress rate denoted here by &ij and defined as CTij = ^<Tik{vj,k - Vk,j) - ~ <7jk(v k CTij - it - (3.13) Vk,i) where the terms in parentheses are components of the spin tensor. The constitutive relation can be written in the rate form as Vij = CijkiDki (3-14) where DM = (vk,i + vi k)/2 are the components of the rate of deformation tensor, and Cijki are t the components of the fourth order stress-strain constitutive tensor. The conservation laws may be expressed with respect to the referential coordinate system by substituting the rate terms from equation (3.9): 9 p V i + dx~- = ~ dx~- >dx- = d^ +Wi pW p e + p W j dx- = p + dx- aiJ p b + ' ( 3 i ( 3 -dVj pr ' ( 3 1 5 ) 1 6 ) - 1 7 ) Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 38 where ()' indicates the time rate of a variable with referential configuration held constant. Also, the time rate of Cauchy stress in equation (3.13) may be expressed in the referential system as crij = a' + w -^ {j 3.1.2 (3.18) k Virtual Work Equation and Its Linearization Many forming operations may be treated as quasi-static processes, in which the inertia plays no significant role. Also, since in forming processes the plastic strains are significantly larger than their elastic counterparts, the density may be assumed to remain constant due to plastic incompressibility. As a result, the continuity equation is satisfied automatically. With these assumptions, the weak form of equation of conservation of momentum may be expressed as the equation of virtual work. Linearization of this equation is carried out by adopting an incremental approach. The derivation in this section follows that of Wang [39]. In an increment from time t to time t + At, the equilibrium of the body is expressed in the form of principle of virtual work, / Jt+Aty t + A t ^ - ^ h ^ O Xj V Sup %d V+ = / t+A Jt+At t+At [ t+At Jt+At r V Su^fid^S J (3.19) ' v S where /; and hi are the components of body force and surface traction vectors, respectively, and Sui is the variation of an admissible displacement field. In order to facilitate the solution of equation (3.19), all quantities may be transformed into the known computational reference configuration at time t. This configuration is chosen to coincide with the material configuration at time t. Therefore, the transformation is only due to the change in or motion of computational reference domain in the increment At. The relationships between the volume, area, and partial derivatives in the two configurations arefirstestablished. For the volume, we have [79] d*+ V = (TV \t +{d*V)'At = cfV \t M x x +(|^)Ai (3.20) Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 39 The term involving spatial derivatives may be expressed as dSui and since t + A t x = x + u~i, where t k k dSiii d Xk t (3.21) is the displacement of a mesh point in from time t to Ui time t + At, by substitution in equation (3.21) and letting At â€”v 0, we obtain .dSui.. dSuid^k Therefore, we have dSui dSui ~ Wx~ (dSuA' V x + k dSui \Â¥x~ ) ~Â¥x~ k V x k dSuid^k . ~Â¥x~ 'Â¥x~i , , ( 3 , 2 3 ) k The change in surface area is represented by the following equation o? S +At = (TS \t +(d*S)'At (3.24) x It can be shown that [79] : <^>' = <^S-'"'S'"<) - (3 25) where rii are the components of surface unit normal at time t. The Cauchy stress in equation t (3.19) is expressed in terms of material time rate of stress; = V V +%At = ^ y x |. + i^an - At x (3.26) Similarly, the body forces and surface tractions may be transformed to configuration at time t, f = fi + 7/At = 7* + {^h - At (3.27) % = % + %At = % + K - V k ^ - ) At (3.28) t+At t+A l Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 40 After substitution of equations (3.20-3.28) into equation (3.19) and neglecting all terms involving orders of Ai higher than one, the equation of virtual work in the ALE description may be written as: Atd*V AiJ 5(3.29) t Thefirstterm on the left hand side and thefirstand second terms on the right hand side are cancelled out on the account of the equilibrium at time t. Then, thefinalALE equation may be written as (3.30) This equation is general in the sense that it may be applied to variousfinitedeformation problems with generalized types of loads and boundary conditions. The rate form of the equation makes it easy to implement rate dependent material constitutive relations. Other characteristics of this formulation are: â€¢ The equation is written explicitly in terms of both material and mesh velocities. Therefore, the mesh motion and material deformation occur simultaneously, but independently, at each increment. In this sense, this equation represents a fully coupled ALE formulation in contrast to the operator-split techniques. This makes the formulation flexible and Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 41 adaptive in handling various types of boundary conditions and in modeling evolution of free surfaces. â€¢ To incorporate the constitutive equation in the rate form, the Jaumann rate of Cauchy stress may be introduced into equation (3.30) to substitute the time rate of stress. Then the equation takes the following form; (3.31) â€¢ The formulation reduces to either updated Lagrangian or Eulerian formulations as special cases. In the Lagrangian case, the material and mesh velocity coincide, Vi = and thus the convective velocity Wi vanishes. Then, the equation takes the form (3.32) This equation is similar to the updated Lagrangian formulation given by McMeeking and Rice [38] except for the load correction terms in the last two integrals. The load correction terms are useful in handling of deformation dependent loads. On the other hand, if mesh velocity is set to zero V{ = 0, then u>i = V{. Upon substitution of these conditions into equation (3.30), it reduces to the Eulerian formulation; Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 42 d fi d*S fi ~ v dx ) l (TV+ [ Sm k l (3.33) k â€¢ Since two sets of velocities are defined for the motion, the mesh velocity and the material velocity, one set of these variables must be specified prior to solving the equation in every incremental step. The mesh velocities have to be assigned in a scheme outlining the motion of the mesh. This gives the user the flexibility to design algorithms which maintain the mesh in optimal shape throughout the analysis. â€¢ The drawback of ALE method in contrast to its Lagrangian counterpart is that, similar to Eulerian formulation, it involves convective terms which should be handled carefully. Furthermore, since variables such as strain and stress are associated with material particles, and are obtained at material points, an updating scheme is necessary to transfer these values to the mesh points at the end of each incremental step. The ALEfiniteelement equations are obtained from discretization of equation (3.31) [39]. The final form of these equations may be cast in the standard form as [K + K L where K and K L N L N L ]v+K v= f G (3.34) are the linear and nonlinear stiffness matrices related to material defor- mation, K is the stiffness term resulting from grid motion, and f is the load load vector. G This equation yields an unsymmetrical stiffness matrix oi N x N equilibrium equations for 2N variables, where N is the total number of degrees of freedom in the model. Therefore, supplementary equations are required for solving this equation. Such equations are provided by a mesh motion scheme that describes the desired motion of the computational mesh within an incremental step. Often, such scheme is designed in a way that mesh distortion is minimized and the regular shape of the elements are maintained throughout the analysis. The details on design and implementation of mesh motion scheme are described in chapter 4. Chapter 3. THEORY OF ALE FINITE ELEMENT 3.2 METHOD 43 STRAIN R A T E A N D T E M P E R A T U R E EFFECTS The inelastic deformation of most metals is known to be rate dependent even at room temperature, but these effects generally become more significant at elevated temperatures. The properties of materials also vary considerably with temperature. Since significant heat generation occurs during most metalworking operations, it is necessary to include rate and temperature effects in the numerical modeling of these processes. 3.2.1 Treatment of Rate Effects in the Literature The rate effect on the material deformation is usually considered in the form of viscoplasticity. The classical treatment of viscoplasticity is based on the 'over-stress' model of Perzyna [81], in which the rate of viscoplastic strain is determined by the instantaneous difference between the effective stress and the current yield stress of the material. In a certain class of large deformation analysis, known as the flow approach, the deformation of a solid is modeled as the flow of a non-Newtonian viscous fluid [57]. In this approach, the elastic deformation of material is neglected in comparison to its plastic deformation, and thus the metal deforms as a fluid shearing under deviatoric stress. As a result, the solution procedure becomes similar to that of nonlinear small strain elasticity problems [82]. A variation of the flow approach which includes elastic strains through an iterative procedure was presented by Dawson and Thompson [83]. In contrast, in the solid approach, the elastic strains in the material are the basis for calculating the stress, and commonly, an incremental approach is used in which the distribution of variables at time t is known, and those distributions at time t + At is sought [84]. The over-stress model has been widely used in both solid and flow approaches. In another class of viscoplastic solutions, constitutive equations are based on explicit definition of internal state variables that represent the underlying microscopic state of material. For example, Hart [85] used two state variables called the anelastic strain tensor and Hardness, Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 44 and Anand [86] introduced internal variable called deformation resistance which represented the resistance of the material to inelasticflow.Evolution equations are required for such internal variables. Finally, the rate effects on plastic deformation of material may also be considered by direct extension of rate-independent plasticity to include strain rate -and temperature- effects. In this approach, the size of yield surface is dependent on these variables [87, 88]. This approach is simple in formulation and implementation, and does not require defining any material dependent parameter other than those used inflowstress of the material. This approach is used in this work, and is described in the following. 3.2.2 F o r m u l a t i o n of the C o n s t i t u t i v e E q u a t i o n The rate form of the constitutive equation may be written by replacing the strain increment with the rate of deformation tensor Dij and the stress increment with an objective stress rate, such as Jaumann stress rate <T;J (equation 3.14). Infinitestrain plasticity, the deformation gradient tensor F may be decomposed into elastic and plastic parts F = F F . This multiplicative e p decomposition may be approximated by an additive decomposition of the rate of deformation tensor if the elastic strains are small, as in most metals, and the rate of total strain is defined as the rate of deformation [89]. Then, the stress may be obtained as a function of elastic strain. Extending this decomposition to include thermal strains, the rate of deformation tensor can be decomposed into elastic, plastic, and thermal parts; Dii = Dlj + 1% + Dl (3.35) where it is assumed that the elastic part is independent of rate and temperature effects. For an isotropic elastic material response, the rate form of generalized Hook's law is written as *X = C; DU jU = C ? . â€ž (D - Du - D ) T U kl (3.36) Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 45 where vE C? jkl = + ( 1 â€ž + ) ( 1 _ 2 l / *Â» s ) s is the fourth order elastic tensor, E is Young's modulus, v is Poisson's ratio, and G is the shear modulus of the material. The plastic part of rate of deformation tensor is governed by the normality rule of plastic law; â€¢ OF Dl = A â€” (3.37) in which A is the proportionality factor, F is a plastic potential function which is the same as the yield function in associated plasticity, and may be given by F{a Y) = f(a )-Y ih = 0 ij (3.38) and Y is the current uniaxial yield strength of the material, which may be dependent on variables such as effective plastic strain, effective plastic strain rate, and temperature; Y = R(e ,D ,T) p (3.39) p in which the effective plastic rate of deformation D = {^D -D ^j is the measure of strain p p p 2 rate. Using von-Mises yield function, /(cr ) = a, the plastic consistency condition F = 0 i;) results in df . i ai do- i:i _dRde OR dD p 11 dRdT p ~ a^^JT + am n + Â«nâ€ž deP dt 8DP dt dT dt (6AO) substituting equations (3.36) and (3.37) in equation (3.40) yields j^QW (D i - A ^ - - D^j = HD p k + MD V + NT (3.41) where the strain hardening, strain rate hardening and temperature softening parameters are denoted by H = Jp, M = and N = |f, respectively, and associated plasticity is assumed. By introducing the rate of plastic work density as aD = o - y D ? , p (3.42) Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 46 and after substituting in equation (3.41), and some algebraic manipulation, it can be shown that i ^-Ct (D -DZ)-{Mf iu u df ne. df d<r , mnrs 9<T n l U<Tm df 9 n r m NT) + " V ' / d<T â€ž m Upon substitution of equations (3.37) and (3.43) in equation (3.36), we find = Ct? (D - D ) + ^ (3.44) T kl where U / kl C 9/ df ne ijkl - \^ijkl _dÂ±_ne df TTg df \ da ^rnnrs da . ~ a damn e , mn mn r and Ch i^( MbP df da 'â€¢v - k n e m n r s mn + ) Nf df , TLX Q~rnn df â€¢ d a , ~ S da r mn Thefirstterm on the right hand side of equation (3.44) corresponds to usual rate independent plasticity, and the second term is the modification due to rate and temperature dependence of the material response. When equation (3.44) is substituted in the ALE equation (3.31), the second term will act as a pseudo-load term and may be added to the right hand side of the FE equations. Finally, thermal strains in equation (3.44) are calculated from = aSi/f (3.45) where a is the coefficient of thermal expansion. In order to implement the above constitutive equation in thefiniteelement procedure, it is â€”P -p necessary to define variables H, M, N, D and T. Theflowstress function R(e?,D ,T) is provided for the material from experimental tests. The partial derivatives of this function with respect to its variables are approximated as R (*+ *e ,*5 , T) - R A p p t (^//JP/T) - *gp H â€” 4 R /3VT) t+At Â£)p M -- R * â€” = ( TÂ§P tB , T) p t+M - R Â»*T-*T " ( 3 - 4 6 ) Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 47 where estimates of variables at time t + At are obtained from forward or central difference approximations; t+At- P t P tf)PAt = Z t+AtÂ£jp _ + 2 D â€” ~ T) t t+Atrp p 2*7"* t At p tâ€”Atrp ^2 47) -p Similarly, D and T may be estimated by using a Taylor series in which the second order term is approached by thefirstorder term, resulting in: t^j _ __v )_ _ t-Atjy At )(t'Ti t-Atrp\ t T = tr A t ~ (3.48) It is noted that at each increment of the analysis, the rate of deformation tensor is calculated from converged velocities. The estimate of variables at time t + At could be updated during iterations within each loading increment. However, our experience showed that such scheme may induce oscillations in the value of plastic strain rate. Instead, the above variables are not updated in each iteration. To compensate for possible deviations from flow stress curve due to approximations in equations (3.46), theflowstress value for each integration point is updated with converged variables at the end of each increment, and the stresses are scaled back on to the yield surface at time t -\- At. Any imbalance created by modification of stresses is then calculated and added to the residual forces in the next incremental step. This proved to be very stable, and with moderate size of incremental steps, minimal scaling of the stresses were required, and the imbalance forces turned out to be small. 3.3 H E A T T R A N S F E R ANALYSIS In a metalworking operation, substantial amount of work is dissipated due to plastic deformation and friction. Most of the dissipated work is converted into heat and raises the temperature of Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 48 the workpiece and the tool. Since the deformation of most materials and their properties are temperature dependent, the assumption of isothermality is not appropriate for a process involving significant heat generation, and it is essential to determine the extent of heat generation and heat transfer within the material and the tool. The thermal and stress analyses are in fact coupled. It is possible to obtain fully coupled thermo-mechanical equations by account of energy balance and dissipation from the first and second laws of thermodynamics [87]. Alternatively, In an incremental analysis, the coupling may be achieved in a staggered fashion. In this approach, at each increment of loading, an isothermal deformation analysis is carried out with known temperature distribution at the start of the step. Then, the work dissipated in the current step is used as the source of heat generation for an iso-stress thermal analysis which gives the temperature field at the end of current step. The thermal and mechanical variables could be updated by iterating between the two processes until convergence is achieved [90]. However, if the time step is sufficiently small and dependence of material properties on temperature is weak, such iteration may not be necessary. The time step used in thermal process may be different from that of stress analysis, especially if an explicit time integration is used in thermal process, because such approach is only conditionally stable. The general form of energy equation for a moving object with respect to a fixed set of coordinates is pcT + pcviT^ + qi,i = q (3.49) in which v is velocity, T and T are temperature and its rate of change, p is density, c is specific heat, q is heatfluxper unit area and q is the rate of heat generation per unit volume. Using Fourier's law for isotropic medium qi = â€”fcT,,where k is the conductivity coefficient, equation (3.49) can be written as pet + pcviTi - {kT ) i = q . ti t (3.50) Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 49 The boundary conditions for this equation are qini = -q, \ +h (T. - T ) \s +ea (T - T ) | , 4 Sl f 2 and 4 s Ss T \ =T Si s which represent specified flux, convection, radiation and specified temperature on respective parts of boundary, Si to 5 . Here, h is the convection coefficient, Tf is thefluidtemperature, 4 T is the temperature at the walls of enclosing space, e is the surface emissivity, a is the w Stefan-Boltzman constant, T is the surface temperature, and n is the surface normal vector. s Equation (3.50) is a combined convection-diffusion equation in which thefirstterm represents transient effects, the second term is the velocity convection or mass transport conductivity, and the third term is due to heat conduction. Applying the divergence theorem on the third term to include boundary conditions (ignoring radiation here), and discretizing the equation using a weighted residual method yields [91] CT + K T = Q (3.51) where C = J cWN dV K = K + K + K = / B k B d V + / pcWv BdV + / h W N d V T P T c v h w Jv T Jv Jv Q = Q + Q + Q = I WqdV + / q WdV + [ h T W d V g s h s f JSi Jv JS 2 in which W and N are the weight function and shape functions, and matrices B and B w involve derivatives of weight and shape functions, respectively. Utilizing a one-parameter time stepping scheme t+M ^t T T + _ y ^ a f + t+Atf^ a A ( t 3 5 2 thefinalequation for transient heat analysis is written in the following incremental form; (A7 c + a K ) A T = ( ~ )*Q + * 1 a a + A t Q - * K T (- ) 3 53 ) Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 50 where A T = *+ T â€” *T is the temperature rise in the current step. At It should be noted that the velocity convection term in equation (3.50) is an unsymmetrical term even if W = N, and thus requires an asymmetric equation solver. The presence of this term in the equation has considerable implications from the numerical point of view. It is known [53] that such terms in a conventional FE analysis may introduce large spurious oscillations in the solution which sometimes blur out the actual temperature results. Though such oscillations may be avoided by using afinemesh, often a special technique called the "upwinding scheme" is used which eliminates the oscillations in the solution. Many upwinding methods are available in the literature [92]. For example, a weighting function of the form W = N + BF may be used, where 3 is the matrix of upwind parameters and F is the vector of upwind test functions. A proper set of upwind parameters removes the oscillation and results in a stable solution. However, the drawback of the upwinding scheme is in the introduction of artificial diffusion terms which entails some inaccuracies [91]. In solid mechanics applications, the velocity convection term is a feature of an Eulerian formulation in which the mesh is spatially fixed. Since this formulation is usually used for studying steady state problems, thefirstterm of equation (3.50) involving transient effects is dropped. On the other hand, in a Lagrangian formulation the coordinate system is attached and moves with material points, so the velocity convection term vanishes. In the more general ALE formulation, both transient and convective terms could be present, but the velocity convection term should be modified to act on the difference between mesh and material velocities, i.e., pcT' + pc(vi - Vi)T - {kTi) ti ti = q (3.54) This equation is equivalent to the ALE energy equation (3.17) if the internal volume source r is neglected, e = cT, and the rate of heat generation is defined as dvi dxj (3.55) Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 51 where use has been made of the symmetry of the stress tensor. In practice, not all the dissipated work is converted to heat, and the dissipation is only due to the plastic part of the rate of deformation tensor. Therefore, q = rjaD (3.56) p where n is a factor reflecting the percentage of work converted to heat, usually assumed to be around 90%. The heat generated due to friction at the contact surface between the tool and the workpiece is incorporated by defining q = s (3.57) n\o-fVf\ where 07 and Vf stand for frictional stress and tangential sliding velocity at the contact surface. The heat flux due to friction is shared between the tool and the workpiece. Following [49], the heat flux may be divided between the two bodies based on their thermal properties I 9Â« â€” Q I t \toolâ€” a +a where a = (k p c ) t t t t 1/2 , and a = w Q â„¢ / q q$ \workpieceâ€” ?Â») t _ 9Â« E O \ ^O.OoJ a + a w t w [k p c ) l l 2 w w w The above arguments about the need for upwinding schemes also applies to the ALE heat equation due to the presence of velocity convection terms. However, an alternative approach may be adopted in which one can take advantage of the flexibility of the ALE method to solve the heat transfer problem in a Lagrangian frame of reference, and thus avoid the need for upwinding schemes and their side effects. This is carried out as follows; at the end of stress analysis in each step, the new material configuration is known. The heat transfer equation is solved for this configuration to obtain the change in temperature distribution resulting from heat generation in the mechanical part of the step. Since a material reference frame is used for the thermal part, it will be a pure Lagrangian process, and no convective term enters the heat transfer equation. Once the temperaturefieldis calculated for the material mesh, it can be transferred to the computational mesh in the same way the other material associated variables Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 52 are mapped. This will ensure the consistency of mapping scheme for all mechanical and thermal variables. It should be noted that while this scheme is in effect an operator split between the mechanical and thermal processes, it is not expected to introduce extra approximations beyond those adopted in a staggered approach, because the split between the two processes has already been assumed when this approach is used. This scheme is easy to implement, and the tests performed proved it to be stable and accurate. 3.4 C O N T A C T A N D FRICTION A common feature of most metalworking operations is the presence of highly-pressured areas of contact between a hard tool and the workpiece. Depending on the degree of lubrication, the contact may be characterized as frictional or frictionless. Quite often, the extent of the contact area is not generally known a priori, and must be calculated as part of the solution. In many practical problems, the contact area between the tool and the workpiece evolves in the course of deformation, and the frictional conditions may change from sticking friction to sliding friction. As a result, a contact algorithm is needed that is able to â€¢ detect new areas of contact between the two bodies, â€¢ apply appropriate contact constraints and frictional conditions at the interface, depending on whether sticking or sliding conditions exist at a given point, â€¢ detect any loss of contact at some interfacial areas, and ensures that the contact forces on released areas are zero. In addition, finding the direction of application of frictional forces is often a numerical challenge. If the frictional force is applied in the wrong direction, it will assist tangential motion instead of resisting it, yielding erroneous results. Other challenges of contact problems include the treatment of inter-penetration between contacting objects, and distinguishing between Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 53 sticking and sliding frictions. The Coulomb friction law is frequently used in modeling sliding friction, although more rigorous laws have also been introduced based on micro-structural analyses [93]. 3.4.1 Treatment of Contact Problems in the Literature Many numerical contact algorithms have been proposed in the literature for applying the constraints arising from contact between deformable bodies. They include Lagrange multiplier method, penalty method, gap and friction element, and direct method. The Lagrange multiplier method is one of the most common approaches and is implemented in some commercial codes such as ADINA [94]. In this method, the constraint equations in the form Gu + d = 0, where u is the vector of displacements on the interface, is enforced using an array of Lagrange multipliers A, and the stiffness and load terms in the finite element equation Ku = f are augmented in the following form K G G T 0 ' fu] \-\ i f 1 -d The physical significance of Lagrange multipliers is the contact pressure. Therefore, in this approach both displacements and contact pressures are treated as independent variables [95, 96, 97], and the constraint equations are strictly satisfied. However, in this approach, the size of the stiffness matrix varies as the number of active contact nodes changes. Furthermore, the stiffness matrix contains zero entries on its diagonal. In the penalty function method, on the other hand, the constraint equations are satisfied in an approximate way, by introducing a large penalty number a [98, 99, 100]. Only, in the limit of a â€” > oo, these equations are satisfied exactly. In this approach, unknown vector of variables contains only displacements (or velocities), and the contact pressures are obtained Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 54 after displacement solution is obtained, using A = a (Gu + d) p As a result, the penalty method is numerically convenient, and is easy to implement into existing FE codes. This approach has been implemented in many FE codes such as NIKE2D[98], and DEFORM2D [52]. However, the accuracy of the solution in this method is strongly dependent on the penalty number, while there is no clear rule to choose this parameter. Another common treatment of contact problems is through insertion of specially designed gap or friction elements on the interface between the two bodies. Typically, these elements connect the nodes of contacting bodies, and become active only when contact is established. The stiffness of the element in the normal and the tangential directions is usually defined by the user. The difficulty in the node-to-node gap elements is that such element can be greatly distorted if large relative displacement occurs between the two contacting bodies, resulting in poorly conditioned stiffness matrices. Furthermore, parameters of the element has to be defined by the user, while the convergence and accuracy of the solution is greatly dependent on these parameters. Various types of gap and friction elements are used in general purpose FE codes such as ANSYS [101]. Finally, if one of the contacting bodies is significantly harder that the other and its deformation is not of primary interest, a direct approach may be adopted in which the harder body is assumed to be rigid, and the contact constraints on the softer material are enforced through the application of force and displacement boundary conditions. This approach is used widely in analysis of metal forming processes [102,103]. In such simulations, the forming tool is commonly modeled as a rigid material. However, if elastic deformation and the stressfieldin the tool is of interest, this approach is not suitable. In what follows, two distinct contact algorithms are presented for applying contact conditions on the deforming bodies, within the framework of an ALE analysis. In thefirstone, the contact is assumed to occur between pairs of deformable bodies and a rigorous scheme for Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 55 applying constraints arising from contact is developed. In the second case, a rigid-deformable pair is considered and a direct scheme is developed to apply the constraints through boundary conditions. 3.4.2 Contact Between Deformable Bodies Consider the case of two bodies A and B that are in contact over a common surface described by S , as shown infigure3.2. In line with common notation in the literature, we may call body c Figure 3.2: Bodies in contact A as the contactor and body B as the target, or equivalently for a forming process, workpiece and tool, respectively. As the contact area between the two bodies evolves, a contact algorithm should be employed to update the contact area and to apply appropriate contact constraints and frictional conditions. Here, the bodies are assumed to be deformable, though usually the tool material is much harder than the workpiece material in a forming process and its deformation may remain in the elastic range. Here, the formulation is presented for the general ALE case, which includes the possibility of mesh motion tangent to the contact boundaries. The FE region of analysis may be viewed as a single body with an internal boundary over which slippage occurs, i.e., the tangential component of velocity is discontinuous, and frictional work is done at this boundary. The normal component of velocity at this boundary is continuous, so that no Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 56 gap will be created. Also, the force vectors on the two sides of this boundary are equal, but opposite. Alternatively, the analysis region may be divided into two continuous sub-regions that share a common interface on which the above constraints are applied to ensure the compatibility of their normal displacements and equilibrium of the contact forces on the contact boundary. In deriving the finite element equations from principle of virtual work, the work of frictional forces at the sliding surface should be added to the work of other external forces for each sub-region [104]. In the following, It is assumed that the nonlinear finite element problem is solved incrementally for the increment from time t to time t -f At. If the forces due to contact at the interface S are distinguished from all the other surface c tractions, a term is added to the right hand side of ALE equation (3.31) to account for power of frictional forces. Then, the FE equations for body A will have the following closed form K W = P4 + where and (3.59) are the stiffness and load terms arising from ALE equation. The latter term includes the rate form of body forces and non-contact external forces, in addition to pseudo-load terms arising from application of mesh motion scheme as will be described in section 4.1.1. C' is the vector of rate of contact loads on body A which arises from the contact load term A which is added to the principle of virtual power in equation (3.31) in the form U= A [ Suf (cf + cfvÂ£ tk - af vÂ£ nf k tj + (vÂ£ - vÂ£) o-f. nf) dS k c (3.60) where cf is the material rate of contact tractions on body A, and all variables are referred to time t. After descretization of the ALE equations and implementing the mesh motion scheme to express grid velocities in terms of material velocities, the terms related to contact forces (equation 3.60) may be re-arranged in the form TI = K > A A + RX + C l (3.61) Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 57 such that the appropriate terms can be added to stiffness and load terms, resulting in the final equation for body A; K v A =R +C A A (3.62) A where KA = K A + K A , R A = R A + R-AÂ» CA= C A and C is the global vector of rate of contact tractions applied on body A. Similar equations A can be written for body B. Thefinalequation for the wholefiniteelement region is obtained by combining equation (3.62) with the similar equation for body B; K 0 A R/ (3.63) < 0 K RT B The stiffness and load matrices and vectors may be partitioned to distinguish between contact and non-contact nodes in the two bodies. KAA KAW 0 0 KWA Kww 0 0 v v A w < > 0 0 KXT KTB v 0 0 KBT KBB VB T = < RA 0 Rw Cw Rx CT (3.64) 0 > where the contacting nodes in body A are denoted by subscript W, while those denoted by subscript A are the rest of nodes in that body. Similarly, the contact and non-contact nodes in body B are denoted by T and B, respectively. Note that there are no non-contact external forces on contact nodes, and the terms Rw and Rx are the pseudo-load terms arising from coefficient aj in the mesh motion scheme. These terms are referred to time t and their Values are known in the solution of current increment. Also, in writing equation (3.64), we have assumed that the numbering of contact nodes is done in a way that such partition is possible. Though not required, such a form makes the numerical calculation more efficient, because the stiffness matrix remains near-banded [96]. Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 58 The velocities and contact traction rates in equation (3.64) are all unknown. Thus, If N t is the number of total velocity degrees of freedom in the model and N is the number of c components of contact forces on the two bodies, there are (N + N ) unknowns for only N t c t equations in (3.64). The remaining equations should come from the constraints applied at the interface of the two bodies due to contact. These constraint equations ensure the compatibility of the displacements and the equilibrium of contact forces at the interface of the two bodies. The constraints can be defined depending on the status of contact at any particular contact point, i.e., sticking contact or sliding contact. In the remaining part of this section we will limit the discussion to a 2-D case, although its extension to 3-D is straightforward. a) Sticking Contact The sticking contact conditions prevail when the normal force at a point is sufficiently large to prevent any sliding at that point. Under the Coulomb friction law, this condition can be mathematically expressed as \c \ < p, \c \ , where c and c are contact tractions in the normal t n n t and tangential directions, respectively, and p, is the coefficient of static friction. We also s assume that when a workpiece nodefirstcomes into contact and penetrates the tool, the sticking contact conditions exist. At a sticking point, there is no relative motion between the workpiece and the tool. When the node K of the workpiece penetrates an element of the tool on the side bounded by nodes / and J as shown infigure3.3, a consistentfiniteelement approach is used Figure 3.3: Corresponding nodes of tool and workpiece in contact Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 59 by expressing the velocity constraint for such a node in terms of the shape functions of the tool nodes, vÂ« = HW + H v J m J m +v p m (3.65) where H and H are shape functions for the line element IJ corresponding to nodes / and J, 1 3 respectively, and v^ are the components of the penetration vector. This vector is known from previous iteration in which an overlap between the bodies is detected, prompting a re-iteration to apply contact constraint at the overlap point. Note that although a two-node line element is shown above, the constraint equation can be written for a contacting element with any number of nodes on its boundary. This constraint equation ensures that the penetration is compensated, and that the workpiece node remains on the surface of the tool. Considering a set of contact nodes as in figure 3.3, there are 6 velocity variables (v^, v^v^) and 6 force variables ( , c^, ), (m = 1,2) which are unknown, for only 6 equations. Using Equation (3.65), the two variables v^ are eliminated by expressing them in terms of other variables. Corresponding to the above three nodes there are six contact force components on the right-hand side of equation (3.64). Another four of these variables can be eliminated by using a consistent approach to express the force constraints in terms of shape functions; & = c J m = H cÂ«, m=l,2 J (3.66) Therefore, the number of unknown velocity variables eliminated is equivalent to the number of unknown force variables added. If the stiffness matrix and contact force vector in equation (3.64) are rearranged to include the unknown contact forces in the vector of unknown variables on the left-hand side, a stiffness matrix of size (N x N ) is obtained, and the equations can be t t solved for unknown velocities and contact forces. b) Sliding Contact If the tangential force at an interface point is larger than frictional capacity of the bodies, Chapter 3. THEORY OF ALE FINITE ELEMENT 60 METHOD relative tangential motion between the bodies will occur, and the contact is characterized as sliding contact. According to Coulomb friction law, such condition will occur if \c \ > p \c \. t s n Once the sliding conditions prevail at a point, the relationship between normal and tangential components of the contact force is given by \c \ = p,d \c \ in which (id is the dynamic coefficient t n of friction. For a sliding node of the workpiece, there is one velocity constraint and one force constraint. While the tangential component of velocity of the sliding node is not constrained, its normal component is constrained with respect to the corresponding nodes of the tool, v* = H'vi + H v , J v =v J n n (3.67) ini Here, n; are the Cartesian components of the unit normal to the boundary at the point considered. This equation can be transformed to the global coordinate system using an appropriate transformation matrix; Tuvf = H'Tijv]*- (3.68) H Tijvj J which represents a coupling between the Cartesian components of the velocity. The transformation matrix T may be defined as T= cosO â€”sind (3.69) sinO cosO where 6 is the angle between the local and the global coordinate systems. Therefore, one velocity degree of freedom for node K can be eliminated from the equations. Similarly, the normal components of contact forces on the tool can be expressed in terms of the normal force at node K\ cl c J n = H cl = H c* l J =â€¢ T ] = ijC H'Tttf Tirf = H T f J ijC (3.70) which allows for elimination of force variables related to nodes I and J. Then, the friction law Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 61 is used to relate the components of the force vector for node K; i- i V t \v \ t When this equation is transformed to the global coordinate system, the two components of the force are coupled, and one of them can be eliminated from the equations. Therefore, the remaining unknown variables are two velocity variables for each of nodes I and J , and one velocity variable and one force variable for node K. Again, the number of equations equals the number of unknown variables and the system of equations can be solved after rearranging the equations to put the unknown force terms in the solution vector on the left-hand side. Once the system of equations is solved, the constraint equations can be used tofindthe eliminated variables, i.e., vf = H vl : + Hv, J J mj cf = Â£ H c?, a mj cf = Â£ H cf, a (3.71) where mj and mj are number of contact nodes on the workpiece which contribute to the forces on nodes / and J of the tool. The above approach in applying the compatibility and equilibrium constraints over the contact interface is more efficient numerically than the Lagrange multiplier and penalty method commonly used in the literature to apply these constraints. Compared to the Lagrange multiplier method, the size of the stiffness matrix is not increased in this approach, and no zero diagonal term is created. Also, this approach does not suffer from the difficulty in penalty method that arises from the dependence of its solution on the choice of penalty term. This approach only requires algorithms for implementation of nodal constraints and rearrangement of matrices. Thefirstalgorithm is readily available in mostfiniteelement programs, and the second can be easily developed. A similar approach is used by Voyiadjis, et. al. [105] for Lagrangian contact problems, though they neither considered the penetration, nor treated sliding contact. It is also noted that although only Coulomb friction is considered above, it is possible to introduce other friction models such as shear friction. This will be elaborated in the next section. Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 62 Next, we discuss the detection of penetration and denning of normal vectors at a contact point. The surface nodes at candidate contact areas are specified at the start of the analysis as the potential contact nodes. At the start of each incremental step of the analysis, all the contact nodes should be mapped to a proper segment of the tool contact boundary, e.g. node K in figure 3.3 is associated with segment IJ. To ensure a unique mapping for each tool node, a pseudoelement approach [101] is used in which each segment of the workpiece contact boundary is surrounded by a pseudo element, as shown infigure3.4. The thickness of the pseudo-elements - - - 1- ^ Figure 3.4: Pseudo-elements used for detecting penetration can be selected for an analysis, depending on the size of time steps. A potential contact node on the tool that falls between the two bounding line of a segment is associated with that segment. If the potential node is inside the respective pseudo-element, it is considered close to contact. Otherwise, it will not be examined further in this step. After the solution for current step based on the contact status at the start of the step has converged, the updated positions of "inside" potential nodes will be examined against their respective segments on the tool to detect any overlap of the two bodies. If a workpiece node is determined to have penetrated the workpiece, a penetration vector from the node to the closest point on the segment is generated, and the node will be given a 'sticking' status in the following re-iteration of the current step. The status of all contact nodes is examined following each solution, and updating and re-iteration of the step continues until there is no change in contact conditions in this step. The tool segments can also be used for defining a local coordinate system at a contact point. Assuming that straight-edge elements are used, the tangential and normal directions Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 63 are denned for each segment of the boundary, using the positions of the nodal points on the segment. Finally, the release of a node from contact is determined based on the sign of the force applied on a workpiece node. The contact will be lost as soon as the normal force on such a node becomes tensile. Upon detection of a node release, the contributions of this node to contact forces are set to zero, and the step is re-calculated to achieve equilibrium. Alternative approaches that base the decision for release of a node on the accumulated effects of its neighboring segments are also used in the literature [95]. 3.4.3 Contact Between Rigid-Deformable Bodies If one of the bodies in contact can be assumed as rigid, the.contact problem can be solved using a direct approach. Such assumption is common in metal forming operations, in which the tool is usually much harder than the workpiece, and its small -and often elastic- deformation is negligible compared to the much larger plastic deformation in the workpiece. In such circumstances, the tool is not included in thefiniteelement region of analysis, and the contact area becomes an external boundary of the region. Consequently, the contact conditions at the interface can be applied directly as velocity or force boundary conditions on external boundaries of the workpiece. In this work, a direct node-to-surface contact algorithm is designed for application in processes such as metal cutting. In this algorithm, the perimeter of the tool is defined as a set of geometrical entities such as line and arc segments connected together. In each incremental step of the analysis, the algorithm inspects the boundary nodes to detect any penetration into the tool area. Upon detection of such penetration, the algorithm determines the extent of the penetration, and applies a corrective velocity boundary condition reset the node on the tool surface. Once the node is on the surface of the tool, appropriate contact boundary conditions are applied on the node which constrain it to remain on the surface of the tool. Also, frictional Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 64 force on the node is calculated based on Coulomb or shear friction laws. A node is considered in contact as long as the nodal force normal to the contact interface is compressive. As soon as the normal force on a node becomes tensile, the node is released from contact and its contact forces are set to zero. At the end of equilibrium solution in each step, the contact algorithm checks the contact conditions of all potential contact nodes to detect any change in the contact status of a node. If such a change is detected, contact conditions are updated and the step is re-iterated to uphold the equilibrium conditions. Further procedural details of this algorithm are presented in chapter 4. Finally, the drawback of the rigid tool assumption in metal cutting applications is that the effects of contact forces on the tool cannot be studied, e.g. in study of tool life and wear. Furthermore, because the tool is not explicitly modeled in this approach, the heat loss through the tool cannot be directly modeled. Therefore, either an equivalent thermal boundary condition should be applied, or a separate FE mesh should be used which only becomes active in a staggered thermo-mechanical analysis. The latter approach is taken in this work. 3.4.4 Contact Problems in A L E Analysis Contact problems have been traditionally solved in a Lagrangian frame of reference. However, when a node-to-surface algorithm such as the one described above is used in a Lagrangian analysis, difficulties may arise due to inaccurate modeling of the borders of the contact region. The detection of contact in this type of algorithms is based on penetration of tool surface by a node of the workpiece. Conversely, the tool is allowed to penetrate the workpiece in between its nodal points. Therefore, depending on the density of the FE mesh, some overlaps between the tool and the workpiece, especially at sharp corners of the tool, may occur. To minimize such penetrations, a fine mesh is usually used around the contact zone. This problem is demonstrated here using examples of punch indentation and cutting with a restricted contact length. Figure Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 65 3.5-a shows the initial configuration of the mesh in the frictionless punch problem and the figure 3.5-c shows the Lagrangian configuration after 45% reduction in the workpiece height. As the punch advances, the nodes on the contact surface move to the right and successively (a) Punch T 0 10000 â€¢ i - i i Punch Displacement (mm) (c) (d) Punch T , ,V-irrrr, i i Figure 3.5: Lagrangian and A L E analysis of punch indentation problem, a) initial configuration b) punch force in Lagrangian and A L E analyses, c) deformed shape in Lagrangian analysis, d) deformed shape in A L E analysis leave the tool surface. However, once a nodal point is released from contact, there is no material point represented at the edge of the tool until the next node arrives at this point. As a result, an overlap develops between the tool and the workpiece as shown in figure 3.5-c, and the predicted punch force fluctuates, because every contact release at the edge of the tool results in a sudden drop in the force (figure 3.5-b). Furthermore, the point of maximum pressure at Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 66 the edge will continuously shift back and forth, affecting the deformation history and resulting in underestimation of the punch force. Of course, the severity of these problems will depend on the density of the mesh. However, afinermesh is also more prone to distortion resulting in convergence problems. Wang and Gadala [72], solved a similar Lagrangian punch problem with a refined mesh of around 500 element using DEFORM-2D code. However, the predicted force was still oscillatory, and an unrealistic deformed shape was obtained. Similar results were reported by Voyidjis and Frouzesh [51]. In contrast, the adaptivity of the mesh in an ALE approach can be used to avoid such drawbacks. On a contact boundary, the mesh must follow material motion in the direction normal to the boundary, but it may be moved arbitrarily in the tangential direction. This flexibility of the mesh is instrumental in accurate representation of contact conditions in an ALE analysis; â€¢ In problems where the extent of contact zone is known a priori, the nodes can be moved in a way that there is always a node at the border of the contact zone. â€¢ In contact between deformable bodies, it may be possible to design the motion of a pair of contacting nodes in a way that they always remain in contact with each other. Such a node-to-node contact avoids the approximations involved in enforcing contact constraints (3.65-3.70) through the use of shape functions. â€¢ In problems where the area of contact zone changes in the course of deformation, it is possible to insert nodes at the transition points between sticking and sliding contact zones, and between sliding and contact release zones, such that these transitions are modeled accurately. However, the position of such transition points should be predicted beforehand. As an example, Hiriandja and Haber [106] used the equations of sensitivity of contact forces to grid motion to predict these positions. Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 67 In solving the punch indentation example as an A L E analysis, the nodes under the tool may be kept stationary in the tangential direction, while they are allowed to follow deformation in the normal direction. A s a result, there w i l l always be a node at the edge of the tool, and contact conditions are simulated accurately. Figure 3.5-d shows the configuration of the mesh i n A L E analysis after 4 5 % reduction. It can be seen that in spite of the coarse mesh used, the resulting shape is realistic, and the mesh has remained regular. Figure 3.5-b compares the punch force in both approaches, showing that the punch force predicted in A L E analysis has a uniformly increasing trend as expected in practice. The example of cutting with a restricted contact length is shown in figure 3.6. Again, (b) Figure 3.6: E x a m p l e of a) Lagrangian, and b) A L E analyses of cutting w i t h restricted contact length tool Chapter 3. THEORY OF ALE FINITE ELEMENT METHOD 68 when such a node-to-surface contact algorithm is used within a Lagrangian analysis, overlaps at the tool tip and at the end of contact length occur as shown in figure 3.6-a. The same force fluctuation due to successive release of the nodes is also observed here. The overlaps at the tool tip should be treated by node separation or frequent remeshing in a Lagrangian analysis. Similar overlap was also demonstrated in the results of DEFORM-2D cutting analysis shown in figure 2.7. In contrast, in an ALE analysis, the motion of the nodes may be designed in a way that the contact at the edges of the tool is never lost, and no overlap at the tool tip occurs. Figure 3.6-b shows a typical treatment in ALE analysis, in which the tangential motion of the nodes on the contact surface is Eulerian, while their normal motion is Lagrangian. As a result, accurate solution can be obtained even with a relatively coarse mesh. Finally, it should be noted that the overlaps are the characteristics of the contact algorithm used. Thus, they can be avoided by using a more sophisticated contact algorithm in a Lagrangian analysis, at the cost of higher computation intensity. Chapter 4 N U M E R I C A L PROCEDURES A N D IMPLEMENTATION In implementing an ALE large deformation simulation program, two major numerical procedures are needed in addition to the common Lagrangian routines; a mesh motion scheme which assigns the velocities of all degrees of freedom of the mesh in each incremental step, and a variable updating scheme that updates solution variables from material points to mesh nodal points. The latter issue is common between ALE and Eulerian analyses, but the former issue is unique to ALE. In this chapter, the numerical procedures to achieve these tasks are first presented. Then, the structure of the ALE code developed in this work is discussed. On each subject, reference is also made to some of previous works in that field. 4.1 M E S H MOTION IN A L E ANALYSIS The incremental ALEfiniteelement equation may be cast in a standard form as K v + K (v-v) = f m where K m c (4.1) is the Lagrangian tangent stiffness matrix involving both linear and nonlinear terms, K is the stiffness matrix resulted from mesh motion, and f is vector of rate of incremental load. c For a model with N degrees of freedom, the above equation yields an unsymrnetric stiffness matrix with 2N unknown variables, while only N equilibrium equations are available. Thus, supplementary equations are needed to solve this set of equations. Such equations are provided by assigning the desired mesh velocities in each incremental step at the start of that step. Any desired criterion can be used in deciding the mesh motion, reflecting the 'arbitrary' and adaptive 69 Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 70 nature of the ALE. For example, the mesh motion can be designed such that mesh distortion is minimized and element shapes are optimal. This is particularly useful in areas of the FE mesh where large and localized deformation rapidly deteriorates the mesh. In assigning mesh velocities, two issues have to be addressed. First, how the velocities for each DOF of the model are found at each increment. Usually, a mesh motion scheme is employed which automates velocity assignment. Secondly, given the mesh velocities within a step, how these variables are incorporated in equation (4.1) to facilitate the solution of this equation. We discuss the latter issue first. 4.1.1 Solution Procedure for Mesh Motion Two types of treatments of the supplementary equations are reported in the literature. Some researchers have treated these equations as constraint equations that are solved simultaneously with equilibrium equations for mesh and material velocities [67]. However, this approach is not efficient computationally because the number of equations increases considerably. Alternatively, the supplementary equations are treated separately, so that the size of stiffness matrix is not increased. In this approach, the supplementary equations are solved at the start of the step (time t), and mesh velocities are eliminated from equation (4.1) before solving it. The drawback of this approach is that the mesh velocities will have to be based on the nodal positions at time t, and may not guarantee an optimal mesh at time t + At. Nevertheless, the deviation from optimal mesh is usually insignificant, as the two meshes are only different by a small perturbation. The elimination approach is efficient computationally, and has been used by many researchers [64, 68]. On the other hand, in defining the supplementary equations, it should be noted that a general ALE model might include purely Lagrangian or Eulerian degrees of freedom. As a result, the mesh motion scheme should be designed in a way that its reduction to these types of motion Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 71 is straightforward. Gadala and Wang [107] introduced a mesh motion equation that facilitates the reduction of ALE motion into its special cases. This equation expresses the relationship between material and mesh nodal velocities as Vi = a' ) + b^Vi 1 (no sum on i) (4-2) where aW and b^ are scalar coefficients for z-th degree of freedom of the node. This equation describes the mesh velocity as a combination of an arbitrary motion, defined by coefficient a, and a motion dependent on material deformation, applied through coefficient b. As a result, all types of motions, namely, Lagrangian, Eulerian and ALE, may be specified by suitable choice of the two coefficients. However, problem arises when this scheme is used for assigning mesh sliding, i.e. tangential motion of nodes on external boundaries, because the mesh velocity is only coupled with the material velocity at the same DOE In reference [107] nodes on free boundaries are assigned pure Lagrangian motions in all directions, and are only adjusted between increments and outside solution domain. In order to deal with this problem, we modify equation (4.2) to a more general form Vi = di + BijVj (4.3) where the coefficient Bij are now scalar elements of a matrix. This modified form of mesh motion equation allows a coupling between degrees of freedom of a node, and is instrumental in specification of mesh sliding on external boundaries. Using this scheme, different types of motion are realized by the following assignments; where if di = 0 and b^ = 5ij Vi = Vi if cii = 0 and b^ = 0 Vi = 0 Eulerian motion else, a,i ^ 0 Vi arbitrary general ALE motion. is the Kronecker delta. Lagrangian motion Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 72 When the mesh motion is described by equation (4.3), it is easy to handle the resulting supplementary equations in processing of stiffness matrix. Here, the mesh motion equation is employed to eliminate mesh velocities from equilibrium equations on the element level, i.e. before element stiffness matrices are assembled. This is equivalent to static condensation of mesh velocities and reduces the computational intensity of the solution considerably. By writing equation (4.3) in an appropriate vector form and substituting it in equation (4.1), we have (K m - K ) v + K ( a + Bv) = f c c (4.4) which may be re-arranged to obtain {K m - K (I - B)} v = f - K a C c (4.5) It may be easily verified that this equation reduces to Eulerian or Lagrangian formulations with proper choices of the coefficients. 4.1.2 Computation of Mesh Velocities The next issue is how mesh velocities in each incremental step may be determined such that element distortion is avoided and a high quality mesh is maintained throughout the analysis. For this purpose, various regions of the mesh may require different types of motion. Therefore, a mesh motion scheme is needed that automatically determines the mesh velocities based on some pre-defined optimality criteria. It is useful here to draw an analogy between assigning mesh velocities and mesh generation algorithms. Generating a mesh usually involves three steps. First, a model of complex geometry is divided into smaller areas of primitive shapes, and the type of elements and mesh densities for each region of the model are specified by the analyst. Next, the distribution of nodal points on external boundaries of the body is determined, and finally the interior of the body is meshed. The last two steps are usually carried out automatically by mesh generators, but the first step, Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 73 which is often a crucial one in the success of analysis, is dependent on particulars of each problem, and only general guidelines can be defined for this step. Likewise, the process of assigning mesh velocities may be divided into three consecutive steps. First, an overall strategy for motion of the mesh in different areas of the model has to be designed. Next, the motion of the mesh on boundaries is specified, and finally, the velocities of the mesh in the interior of the body are computed. Similarly, the second and third steps may be automated, while the first step is dependent on the experience of the analyst. These steps are detailed in the following sections. Step 1: Designing overall mesh motion strategy In this step of mesh motion, an overall strategy for the motion of the mesh is outlined, which determines what type of motion is best suited for each region of the mesh. It should be noted that the mesh motion scheme is to maintain the quality of the mesh in successive steps through applying small incremental motions to nodal points, and as a result, the performance of mesh motion scheme may only be evaluated at the end of analysis. Therefore, one has to predict the deformation of various parts of the model in the course of analysis and the final shape of the body. In this respect, this step is more difficult than corresponding step in mesh generation task. On the other hand, often, the initial connectivities of the mesh need not be changed in mesh motion scheme. General types of motion that might be assigned to each degree of freedom (DOF) of the mesh include pure Lagrangian motion in which the nodal DOF is allowed to move freely under load, pure Eulerian motion in which the DOF is spatially fixed, and general ALE motion in which a specific motion is assigned to the nodal DOF. Some general guidelines for designing an overall strategy for a forming problem: â€¢ In areas where large deformation of material, -and consequently, fast deterioration of the Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 74 mesh- is expected to occur, an Eulerian type of motion is often most efficient. Suitably fine Eulerian mesh may be adopted to capture the high stress and strain gradients in these areas without the drawback of quick mesh deterioration. An example of such area is the mesh around the tip of the tool in metal cutting problems. â€¢ On free boundaries, the mesh is required to be Lagrangian in the direction normal to the boundary. This requirement ensures one-to-one mapping of material and mesh domains and allows for natural evolution of free boundaries under applied loads. On the other hand, the boundary nodes may be arbitrarily moved in the tangential direction to prevent element distortion or mesh entanglement. Otherwise, it is often convenient to assign Lagrangian motion to all the boundary degrees of freedom. The same applies for interior nodal points. â€¢ On boundary nodes where displacement boundary conditions are applied, a Lagrangian motion in the direction of the applied displacement is required. However, when the body is modeled as a control volume, Eulerian motion might be assigned to such boundaries. â€¢ On contact boundaries, a motion conforming to the motion of the tool in the direction normal to the tool surface should be assigned, but the motion in the tangential direction is arbitrary. A suitable tangential motion on the nodes at contact interface is sometimes beneficial in accurate modeling of contact conditions. For example, in a flexible-flexible contact between a tool and a workpiece, the nodes of the workpiece may be moved tangentially so that they always coincide with a node on the tool surface. This will remove the approximation involved in enforcing contact constraints through the use of element shape functions. In problems where contact area changes in the course of deformation, a procedure may be put in place that changes the motion scheme for a contact node once it is released from contact. Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 75 â€¢ For the interior nodes of the body where a general ALE motion is considered, the nodes are repositioned in a way that element shapes are optimized. However, if no significant element distortion is expected in an area surrounded by free external boundaries, a pure Lagrangian motion may be more efficient. Similarly, for an interior region whose boundaries are not moved during the analysis, a pure Eulerian motion is most convenient. Step 2: Assigning motion to the boundaries In general, the boundaries of a deformation domain may be classified as one of the two categories; those whosefinalshapes are known a priori, and those whose deformed shapes are not known in advance, but are obtained as part of the solution. Examples of thefirsttype are boundaries with prescribed displacement or velocity and contact boundaries, and examples of the second type are free boundaries and boundaries on which force or pressure is applied. For all external boundaries, equation (3.8) requires that the motion of boundary nodes in the normal direction follows the motion of the material points, while there is no restriction for tangential motion of these nodes. For boundaries with pre-known deformed shape, such requirement presents no difficulty, and it is sufficient to assign a motion in the normal direction equal to the known material motion in this direction. In the tangential direction, the node may be moved as desired to achieve computational efficiency. On the other hand, for a boundary with unknown deformed shape, normal and tangential directions are not known in advance, and usually vary during the course of deformation. Therefore, it is essential to implement a motion scheme which continuously updates these directions, and ensures that equation (3.8) is satisfied. It is sometimes necessary to move boundary nodes in the tangential direction in order to avoid deterioration of the mesh. In dealing with this problem, Gadala and Wang [107] used a scheme that assigns Lagrangian motion to all external boundaries, followed by readjustment of nodal positions on the boundary outside solution domain. Huetink, et al. [65], used an Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 76 operator-split ALE method and updated material surface using spline curve fitting after each Lagrangian step. The mesh motion equation proposed above provides the possibility of coupling various components of nodal velocities through coefficient tensor. Once the tangential coordinate system at a nodal point on the boundary is known, it may be used to move the node tangent to the external boundary. If the tangential coordinate system at such a node is identified by angle ft between the tangent direction and global ^-direction in the two dimensional space, the transformation between the global and tangential coordinate systems for mesh velocities is given by e n s ft ft â€” i n ft cos â€” ssin 6 < 1f Tu 1 (4.6) < sin ft Vy 6 cos and similar relationship can be written for material velocities. Equation (3.8) requires that mesh and material velocities in the normal direction are equal, while the mesh can move arbitrarily in the tangential direction; v = v = â€” sin ft â€¢ v + cos 6 â€¢ v (4.7) v = C (4.8) n n x v t where C is any desired velocity defined for the tangential direction. When this equation is substituted in equation (4.6), the global and tangential velocities of the mesh point can be expressed in terms of their material counterparts in an expanded form of equation (4.3) as sin 2 ft â€” sinftcos ft sinftcos ft â€” cos 6 2 v a (4.9) { y v In order to assign a desired tangential motion for a boundary node, one only needs to assign the desired velocity C and define the tangential directionftfor that node. In each incremental step, the tangential direction at a nodal point corresponds to the path along which the material point at the node moves from time t to time t + At. However, the Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 77 mesh velocities are assigned before the step is calculated and so, the direction is not known, and an approximation is necessary. A reasonable approximation of the tangential direction may be obtained considering the direction of motion from time t â€” At to time t, assuming that the change in node path between successive steps is small. As a result, the calculation of tangential directions always lags behind deformation by one step. Although it is possible to update the nodal directions iteratively within a step, such iteration may be practically unnecessary, because each new mesh is a small perturbation away from the last mesh. In thefirstincrement, a pure Lagrangian motion may be assigned to such nodes, and tangential motion is started from the second increment. An important case of use of equation (4.9) is when the boundary nodes are required to maintain their tangential position on the free boundary, while moving freely in the direction normal to the surface. By setting C = 0, the appropriate mesh motion is obtained. Finally, it is often necessary to slide the boundary nodes in the tangent direction in a way that the element sizes on the boundary remains the same or keep a certain ratio. A simple numerical curve fitting, e.g. splinefitting[108], of the nodal positions may be used tofindthe desired location of the node. The desired velocity of the node can then be implemented through the coefficient C in equation (4.9). Step 3: Moving interior nodal point Once the motion of external boundaries of the model is designed, it is possible to assign the velocities of interior nodes. For models with complex shape, it is common to divide the model into smaller regions with primitive shapes, such as triangular or quadrilateral regions. The mesh motion algorithm should be simple and efficient, because it will be used repeatedly in each step. Also, because it is often not necessary to redefine mesh layout and connectivity beyond the initial mesh generation, an algorithm that can make use of the given mesh topology and Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 78 only re-locates the nodal positions to improve mesh quality is most efficient in mesh motion scheme. Various mesh generation algorithms have been employed for assigning mesh velocities in the ALE literature. Solving partial differential equations forms one class of such algorithms. Liu, et al. [64] solve Laplace differential equation which maps the node space (I, J) into the real space (x,y) together with well-arranged boundary nodes. A similar approach is used by Benson [68] who usesfinitedifference mesh relaxation stencils which are created by minimizing a functional or inverting the solution to a partial differential equations. It may be noted that while solving partial differential equations may create high quality mesh, it often amounts to prohibitive computational cost if it has to be repeated for every step. Another approach used by some researchers [77, 109] is to apply some formulae based on geometric considerations. Giuliani [109] has formulated an equation for mesh velocity, which repositions a nodal point in relation with its neighboring nodes, such that the shape of elements sharing that node remains regular. This may not guarantee, however, the global optimality of the mesh in a general problem. An alternative approach is to utilize algebraic interpolation methods and introduce a set of algebraic equations for mesh motion. These methods are generally simple and efficient, but introduce curvefittingerrors in the description of regions whose boundaries may not be exactly described by polynomials of the same order as those appearing in the interpolation functions [80]. A more viable algorithm for mesh motion is the transfinite mesh generation scheme [80] which creates a structured mesh, but is inefficient if a part of the body has significantly higher mesh density, and transition regions are needed. In this section, a mesh motion scheme that incorporates both transfinite and isoparametric mapping techniques is described. In each region of the mesh, either of the two algorithms may be used at the discretion of the analyst. The transfinite method provides a homogeneous mesh and matches the boundary of a given region at an infinite number of points, so that no curve fitting error is introduced. The isoparametric method, on the other hand, is able to re-mesh a Chapter 4. NUMERICAL PROCEDURES Potential Region AND IMPLEMENTATION 79 Modified Region Figure 4.1: Transfinite mapping region of varying density, but some curve fitting errors may be introduced. Also, a procedure is introduced which accounts for the motion of the boundary. i) Transfinite mapping method: The transfinite mesh generation method is originally designed for creating a mesh on a geometric region with specified boundaries. If the borders of a quadrilateral region are identified as (pi(r, 0), <f>i(r, 1), <j6j(0, s), and <^(1, s) as shown schematically infigure4.1, the coordinate of an interior mesh point Pi is given as: Pi{r,s) = (1 - s)(f>i{r,0) + scj>i{r,l) + (1 - r)&(0, s) + r&(l, s) -(1 - r)(l - *)&(0,0) - (1 - r)s<f>i{0,1) - r*&(l, 1) - r(l - a)&(l, 0) (4.10) where (r, 5) are normalized coordinates, i = x,y; i.e., <f> (r, 0) and </>(r, 0), etc. are the a: x y and y coordinates of boundary curves, and P (r, s), P (r, s) are nodal x and y coordinates of x y point P, respectively. In the discrete version of this algorithm, the functions <j> are expressed as series of discrete nodal points on the borders of the region. It may be noted that the position of interior nodes is totally dependent on distribution of nodal points on the borders of the region. For example, if the nodes are uniformly distributed on all borders, a homogeneous mesh will be created. The boundary curves may be represented by as many nodes as needed, Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 80 and boundaries with discontinuous slopes can be represented without introducing any error beyond initial discretization of the boundary. If the position of a node in the distorted mesh is given by are PÂ°(r, s), then the mesh velocity for this node in the next increment is given by â€ž]->,>,â€ž = * M - * ( - â€¢ Â« ) ( 4 . N ) where P{ (r, s) is the optimal coordinate calculated by equation (4.10). In applying this velocity to the nodal point through the mesh motion equation (4.3), it is suffice to set a; = and Bij = 0. Figures 4.2(a) and 4.2(b) show an example of a distorted mesh and its modification by transfinite method. It can be seen that although the boundaries of this region have discontinuous slope, the algorithm still produces a homogeneous mesh. It should be noted, however, that the mesh at time t + At would be optimal only if the boundaries of the region do not move in this increment. However, since mesh motion and loading increment occur simultaneously, the resulted mesh is based on boundaries at time t. Here, we outline a procedure that accounts for the motion of the boundaries within the increment, and optimizes the interior mesh in anticipation of final shape of the region at the end of increment. We use the same function as in equation (4.11) to interpolate the boundary velocities and find the modification term; vl \r,s) 2 = (1 - s)Zi(r,0) + &(r,l) + (1 - r)d(0,s) + r Â£ ( M ) a -(1 - r)(l - a)Â£(0,0) - (1 - r)^(0,1) - rsfc(l, 1) - r(l - s)&(l,0) (4.12) where &(r, 0), &(r, 1), &(0, s), and &(1,5) are the velocities of the boundary nodes at time t + At. However, if a border of the region is part of an external free boundary, these velocities are unknown at the start of increment. In order to avoid solving implicit equations, we apply the nodal velocities at time t â€” At to approximate the values at time t during the mesh motion. Therefore, for such boundary points; t (r,s) = - v (r,s) i t At i (4.13) Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION (a) (b) 81 (c) Figure 4.2: Transfinite mesh motion a) distorted mesh, b) mesh after applying mesh motion scheme, c) mesh motion with the effect of boundary motion where i = x,y and the mesh velocities on external boundaries are ensured to conform to material velocities in the normal direction. Again, if higher quality mesh is required, some iterations may be introduced, but such measure is practically unnecessary because the quality of the mesh created as above is good enough for most practical applications. Finally, the mesh motion for the nodes in the interior of the domains is given by: ^ = vl + v\ = cn 1] 2) (4.14) For the example shown in figure 4.2, if the distorted mesh at time t has boundaries DE and BE moving at some inhomogeneous speed, using the above mesh motion scheme, we get the mesh at time t + At as shown infigure4.2(c). So, even when boundaries move at inhomogeneous speeds, which is the usual case infiniteelement large deformation problems, the scheme may still give a homogeneous mesh. The only limitation of this algorithms is that the number of nodal points on the opposite sides of the meshed region should be equal. This restricts the method from being applied to more complex geometries where the density of the mesh has to Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 82 be different at various points due to geometrical or analysis reasons, ii) Isoparametric mapping method: The generic isoparametric mapping as a natural extension of isoparametric elements was originally proposed by Zienkiewicz and Phillips [110] for mesh generation. This method is based on mapping between curvilinear coordinates and Cartesian coordinates using interpolation functions. Here, we use this method for remeshing by considering the region to be meshed as a super-element in which the internal nodes of the patch will be mapped using the global coordinates of the boundary nodes of the region, called master nodes, and the parametric nodal coordinates. Here, the method is presented for quadrilateral regions only, but it can be easily extended to triangular regions as well. For a given mesh with specified master nodes on its boundaries, the first step is to find the parametric coordinates (r, s) of all nodal points, for a region of quadrilateral shape with straight edges, the parametric coordinates of a node may be obtained by solving the nonlinear equation resulted from inverting the shape functions of a serendipity element; i i E y^i Ax - E ^ i U + ViSiVs) i =< E ViSi (4.15) 4y - E y ; ( l + TiSirs)- where (x, y) and (r, s) are the global and local coordinates of the node considered, and (a:;, yA and (ri,Si) i = 1 to 4 are the global and local coordinates of the four corner nodes of the patch; master nodes, respectively. This equation needs to be solved only once at the start of the analysis for all nodal points in the region and the values of the local coordinates of the nodes are stored. Then, at every step of the ALE analysis, the current global coordinates of the master nodes and the local coordinates calculated above, are used to find the optimal position of the node in the next incremental step using the direct mapping = Â£ (' Ni r ) ' s Xi y = X) (' Ni r )yi' s * = 1 t o 4 Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 83 The velocity of the nodal points in the next step is calculated from the difference between new and old global coordinates of the node as in equation (4.11). The continuous mapping of the mesh keeps the mesh regular throughout the analysis, without the need to remesh the region. In order to be able to handle curved boundaries, the above analysis may be extended to higher order super-elements. For example, regions with quadratic or cubic surfaces may be mapped using respective shape functions in a super-element with 8 or 12 master nodes, though it will require solving a higher order nonlinear equation. The choice and positions of the master nodes may be used to change the density of the nodes in the region. A standard method such as Newton-Raphson method may be used to solve the nonlinear equation at the start of the analysis [111]. Once this step is completed, the direct mapping is used tofindthe velocity of the nodal points in successive steps. The effect of boundary motion may be included in a procedure similar to that in the transfinite method. Figure 4.3 shows a region with an inhomogeneous initial mesh in which the parametric coordinates of the nodes are obtained from inversion of Figure 4.3: re-mapping the mesh using isoparametric method Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 84 shape functions. These coordinates are then used to map the mesh to new regions created by dramatically altering the boundaries of the region. It is seen that the mapping has kept the regularity of the mesh. An important limitation of the isoparametric mapping is that if the boundaries of the region cannot be represented by a polynomial of the order used in mapping, the boundaries may be altered in the process. Therefore, this scheme is not suitable for regions with complex boundaries or boundaries with discontinuous slope. When this scheme is used in combination with the transfinite method, a robust mesh motion scheme is obtained. In this case, a body with complex geometry is divided into smaller regions, and the analyst can use either of the above mappings in each region at discretion. It is preferable to use transfinite for regions with complex boundaries, but with structured mesh, while the isoparametric mapping is used in regions where change of mesh density is required, such as in mesh transition regions. A note should be made here about the use of mesh motion scheme in keeping the accuracy of analysis when higher order elements are used. One drawback of using elements with mid-side nodes in large deformation analysis is that the change in position of these nodes with respect to corner nodes may degrade element performance. In an ALE analysis, a special procedure may be employed to ensure the optimal positioning of mid-side nodes, and thus improve element conditions. On external boundaries, the mid-side nodes may be moved tangent to the boundary for this purpose, as outlined earlier. 4.1.3 Numerical Examples of A L E Mesh Motion i) Elastic-Plastic Crack Extension Problem Extension of crack in a material creates new boundaries for the material domain. In the simulation of crack propagation using FEM, the topology of the FE mesh has to be changed to accommodate creation of new boundaries. In a Lagrangian framework, a common method to Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 85 deal with crack extension is the nodal release method in which the nodes are separated along a pre-defined path of the crack based on a given fracture criterion. The disadvantage of this approach is that the advancement of the crack is modeled as a discontinuous process in which the crack tip jumps from one node to the next one. Therefore, The accuracy of the analysis becomes greatly dependent on the element size. In elastic-plastic crack propagation, inaccurate deformation history and unrealistic crack profiles may result. A more elaborate modeling of crack propagation is achieved by successive remeshing of the FE mesh to accommodate changes in the material boundaries. The remeshing may be total or partial. In contrast, the ALE method can be used to model crack extension more easily and efficiently, because in this approach the crack extension is modeled as a continuous process. Haber and his coworkers [112] used a displacement-based ALE method for elastodynamic crack propagation. While using mesh motion to advance the crack tip, they remeshed the model periodically to allow the crack to get longer. Here we present the analysis of quasi-static self-similar crack extension using ALE analysis. Using theflexibilityof the mesh in ALE, a scheme is designed that extends the crack smoothly and continuously. As a result, all points on the path of crack experience the center of the singularfield,and accurate deformation history and realistic crack profile is obtained. The analysis is also computationally efficient. Figure 4.4 shows a mode-I crack extension in a center-cracked plate under plane strain conditions. The analysis is started with a crack of initial length a and the plate is subjected to a specified displacement along its horizontal edges. Here, it is assumed that under incrementally applied displacement boundary condition, the crack is extended at a known velocity, but it can also be based on a fracture criterion. At each increment, the crack is extended by a distance Ac*. Four node quadrilateral elements are used in this simulation, but quarter-node singular elements or any other element type can also be used with no difficulty. Due to symmetry, only a quarter of the plate is modeled. Figure 4.4(a) shows the FE mesh at the start of simulation. Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 86 Boundary conditions on the bottom and left edges represent symmetric boundaries, while a total displacement of 0.1 mm is applied on the top edge. The initial crack length is 6 mm and it is extended to a final length of 12 mm. The mesh motion strategy for this problem is as follows [113]. A region around the crack tip is needed that moves with the crack tip and has sufficient resolution to exhibit crack tip singularity. To achieve this, the mesh is divided into three regions, a moving region at the middle, an expanding region to the left, and a shrinking region to the right. The middle region has afinergrid and remains unchanged in the course of analysis. This region moves as a whole with the same velocity as that of the crack tip node, so that its shape is not changed. As this region moves, the region to its left is expanded and the region to its right is contracted, while the mesh motion scheme maintains the uniformity of the mesh in these regions. The boundaries of the plate are Lagrangian in normal direction. Figure 4.4(b) shows the topology of the mesh and the contour of the stresses in Y-direction after crack extension of 3 mm. It can be seen that the singularfieldhas developed around the tip and is moving with it. Figure 4.4(c) shows similar results after the crack has extended by 6 mm. It is observed that the crack is extended smoothly and without changing the connectivity of the mesh or mesh distortion. This analysis preserves the history and profile of the elastoplastic crack. It has also been verified that forces on crack faces in the normal direction are negligible, which conforms to theoretical expectation [113]. To verify these results, a similar case is analyzed using a Lagrangian approach in the commercial FE code NISA [114]. Since it is not possible to model crack extension in NISA, snap-shot analysis corresponding to various crack tip positions are carried out, and thus, some differences will be expected due to different deformation histories. Nevertheless, the results, shown for 6 mm crack length extension infigure4.4(d), shows good agreement with ALE results. Chapter 4. i i. NUMERICAL PROCEDURES AND IMPLEMENTATION t i^ I i i crack tip t 11 / i^/ 1 1 / 14 i. iH tI ii 1â€¢ J ing r&g/oi7 (a) (b) cracfcripposition -3.00E..08 (c) -1.22E-.08 5.56E.-07 2.33E-.08 * 4.11E-.08 5.89E-.08 7.67E4-08 8.44E-.08 1.12E-fOÂ» 1.30E-K) â€¢ 1.0E7) 220.4 186.0 crack tip position 153.0 137.0 121.0 (d) 100.0 83.0O 72.00 . 50.00 . 30.00 . 10.00 . -10.00 . -22.O0 . -34.00 -47.00 Figure 4.4: Mode-I crack extension in center-cracked plate, a) initial configuration and boundary conditions, b,c) configuration and stressfieldafter 50% and 100% crack extension, respectively, d) results from NISA Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 88 ii) Flat Punch Indentation Problem In Flat punch indentation problem, a rectangular plate is coined by a flat punch. This problem was previously solved in chapter 3 with a contact algorithm. However, if the die is assumed to be rigid and frictionless, the punch motion may be simulated by applying velocity boundary conditions. Such a treatment in a Lagrangian formulation will result in increasing the punch size in the course of deformation, and yields unrealistic deformed shape [72]. Furthermore, the elements under the punch are distorted and the analysis may be aborted incomplete due to lack of convergence. The convergence problem is more acute if afinermesh is used. Thus, in a Lagrangian analysis of such problem, a contact algorithm as well as periodic remeshing are required. On the other hand, it is obvious that an Eulerian analysis cannot be used for such unsteady state problem. In contrast, in an ALE analysis, these problems may be side-stepped easily with a proper mesh motion scheme as described below. The model is divided into two rectangular regions, the region under the punch and the region to its right. In thefirstregion, all the nodal points arefixedin the x-direction. The nodes under the punch are free in the y-direction so that they follow the motion of the punch, while the nodes on the left and right borders of the region are adjusted by the boundary motion scheme so that they are equidistant. Similarly, in the second region, the nodes on the free andfixedboundaries are moved in the tangent direction so that the ratio of element sizes in the horizontal direction remains the same. Then, the transfinite method is used to reposition the interior nodes in both regions such that the mesh remains optimal throughout the analysis. The resulting deformed shape is similar to what is shown in figure 3.5. Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 89 iii) Forging Problem The forging problem is presented here to demonstrate the use of the isoparametric mesh motion scheme in problems where the mesh distribution is not uniform. A curved shape die is used to deform a material block under plain strain conditions. The contact algorithm is used to apply contact conditions at the interface of the tool and the workpiece [115]. Only a quarter of the block is modeled due to symmetry, and appropriate boundary conditions are applied on the bottom and left edges to represent symmetry conditions. Figures 4.5(a) and 4.5(b) show the deformed shapes in a Lagrangian and ALE analysis, respectively. It is noted that the (b) Figure 4.5: Forging example a) mesh in Lagrangian analysis, b) mesh in ALE analysis number of nodes on opposite edges of the model is not equal, hence, the transfinite method can not be used for this model. In the Lagrangian analysis, it can be seen that the mesh condition has deteriorated around the inflection point of the die, and some penetration has Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 90 occurred. At this stage of deformation, the Lagrangian simulation is aborted unsuccessfully due to ill-conditioning. In the ALE analysis, the nodes on the interface are Lagrangian in the normal direction, but move tangentially to remain equidistant. The bottom boundary is kept Lagrangian in both directions, because large distortion is not expected there. On the left and right boundaries, nodes are free in normal direction but are adjusted tangentially to remain equidistant. The isoparametric mapping method is used to reposition interior nodal points so that the mesh remains regular throughout the analysis. In ALE analysis, the simulation may be continued beyond the stage shown infigure4.5(b), because elements are not distorted. 4.2 U P D A T I N G SOLUTION V A R I A B L E S Upon specifying mesh velocities and condensing them out of the ALE finite element equations, the solution of equation (4.5) may be obtained by conventionalfiniteelement assembly and elimination techniques, and strains, and stresses are calculated. In general, the solution variables in FEA, such as stresses and strains, are mostly material path dependent, and are obtained at the Gauss quadrature points of the material element at the end of the step. Subsequently, these variables should be updated to the new nodal points of the mesh by tracing the deformation of the material point corresponding to the new nodal point in the configuration at the start of the step. This process is critical in ALE analysis, and may become a source of large errors and degradation of solution accuracy. A Similar step is usually required in remeshing schemes in Lagrangian analyses in which the variables are transferred to the new mesh from the old mesh. In Lagrangian analysis, the remeshing is often carried out every few steps or when the old mesh is deteriorated whereas in ALE, the mesh is optimized at every step, and thus errors due to element ill-conditioning are avoided. Nevertheless, repeated updating of variables makes the solution prone to other types of errors, e.g. interpolation errors, which may accumulate gradually and degrade the solution. Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 91 Also, in Lagrangian remeshing, no physical relationship between the old and new meshes may be established, and the variable transfer is based on geometrical considerations only, while in ALE, the rate equation (3.9) provides an explicit relationship between mesh and material associated variables. Various updating algorithms are used in the literature. Schreurs, et al. [67] use an updating scheme in which a "pseudo material mesh" is defined and the solution variables are obtained by extrapolating from integration points to the nodal points of this pseudo-mesh. Next, iterations are used tofinedthe pseudo-element in which a mesh point is located and the variables are obtained at the mesh points by interpolation. Since this approach requires one extrapolation and one interpolation for each step, the accuracy of the solution may be degraded. Ghosh and Kikuchi [69] update the variables to a mid-step configuration. Another procedure is used by Derbalian, et al. [116] for an Eulerian analysis, in which shape functions conjugate to displacement shape functions are used to approximate a stressfieldwhich is continuous across inter-element boundaries. This approach requires solving a large system of equation for each mesh configuration, and may not be suitable for an ALE analysis. An alternative approach is to base variable updating on the ALE rate equation (3.9) [65, 69]. Since such scheme is based on a continuum mechanics equation rather than pure numerical interpolation, it is deemed more consistent with the ALE analysis. At the end of an incremental step, the positions of nodal and Gauss points of the mesh are known. For any physical quantity / , the change due to material motion A f is related to the x change due to mesh motion A / through the equation x A where f t+At = f( Xi,t t+At Qt+Atf x / = A . / + ^Krr(vi ~ Vi)At (4.16) + At) is the value of / at material point. If the material and computational domains coincide at time t, the total value of variable at the new grid position is Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 92 obtained as a t+At f f= t+At t /+A./ + 5 T r Q t+At f ^f+l^L^-vAM ^(iJ -t, )A*= i i (4.17) Therefore, a physical quantity may be updated from material points to grid point directly. The only difficulty in implementation of this scheme is that the gradient of the quantity should be obtained at the material integration point. Such gradients may be obtained using element shape functions within an element. For example, The equivalent strain at an integration point n may be updated to mesh point by assuming e = J2 h e where n is the number of nodes per a a a=l element, and h are the shape functions. However, these variables are computed at integration a points, and should be extrapolated to nodal points for this purpose. Usually, a local least square smoothing [117] followed by weighted nodal averaging is used for such extrapolation. In a discrete local smoothing, the quantityfieldover the element is approximated by a polynomial whose coefficients are obtained from minimizing a squared function. Some remarks should be made here with respect to the effect of local smoothing on the accuracy of local smoothing; 1. Local remeshing and averaging has the effect of numerical diffusion on the results, and the stressfieldobtained after smoothing may not satisfy equilibrium. In addition, often the current yield stress updated at a mesh integration point by the same procedure does not satisfy plastic consistency. Huetink, et al. [65] applied a weight factor to control the degree of smoothing and showed that diffusion is minimized if this factor is chosen proportional to the ratio of displacement increment and element size. In dealing with the violation of equilibrium and plastic consistency by the updating process, a return mapping scheme may be used to update the stresses back to the yield surface. Nevertheless, this requires further iterations. To avoid such iterations, a relatively crude but often effective way is to scale the stresses back to the updated yield surface and add the residual force created by the force imbalance due to scaling to the load increment in the next step. In our Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 93 numerical experiments, this method proved to be effective and the percentage of stress scaling and amount of residual forces were negligible. It may be noted, however, that in regions where high gradients of stress and strain exist, the updating scheme may perform poorly by introducing too much diffusion so that nodal values obtained from smoothing may become unreasonable. This often occurs when there is large differences between Gauss-point values of a variable within an element, resulting in poor extrapolation to nodal points, and can be improved by refining the mesh in that locality. 2. It is often desirable that the updating algorithm is monotonic, i.e. the extreme values of updated variables do not exceed the values of variables in the element and its neighboring elements. Lack of monotonicity is also due to interpolation functions used in smoothing process. A possible treatment is suggested by Olovsson et al. [74] who have introduced a damping factor which damps the coefficients of smoothing function if the extreme values of the variable is exceeded. 3. In using higher order elements, the least square smoothing may give rise to unreasonable nodal values if the stress/strain values are extrapolated by a linear function within the element. In such circumstances, one has to either use higher order interpolation functions, or use a lower order updating scheme. For example, an 8-node quadrilateral element may be treated as a 4-node bilinear element for the purpose of updating, and the values at the middle nodes may be obtained by averaging respective corner nodes. It should be noted that this will introduce further diffusion into the solution. 4. In areas of high gradients of variables, if nodal variables are averaged before the updating equation (4.17) is applied, severe instabilities may occur after sufficient number of steps. It is therefore advisable that unaveraged variables are stored and used for updating variables within an element. Chapter 4. 4.3 NUMERICAL PROCEDURES AND IMPLEMENTATION 94 P R O G R A M STRUCTURE Based on the formulation and numerical algorithms presented in this thesis, a general 2 D finite element program is developed with emphasis on applications to metalworking processes involving large deformation, rate dependency, heat transfer and frictional contact. The program is an extension of the one developed by Wang [ 3 9 ] , and consists of a main routine and three major modules; the deformation module, the contact module, and the heat transfer module. The overall flowchart of the program which describes algorithmic relationship between the main program and modules is displayed in figure 4 . 6 , and the components are described in following sections. The main features of the program may be summarized in the following: â€¢ A general ALE formulation is implemented which may be reduced to updated Lagrangian or Eulerian formulations by proper choice of mesh motion coefficients. â€¢ The mesh motion scheme of equation ( 4 . 3 ) facilitates the task of assigning mesh motion to various parts of the body, including free and contact boundaries. The user only needs to specify the type of motion required on the boundaries and interior of the body through specifying a set of coefficients. â€¢ General elastoplastic material models including strain hardening and rate and temperature dependent plastic behavior may be used. â€¢ Velocity is the primary variable in the FE solution, which is convenient for implementation of rate dependent plasticity, and for study of velocity effects in forming processes. â€¢ The solver handles unsymmetric stiffness matrices and coupled velocity boundary conditions. The mesh velocities are condensed out of the stiffness matrix at the element level. Full or modified Newton-Raphson solution methods may be used. â€¢ A variety of linear and quadratic 2 D elements is available. Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 95 â€¢ Deformation processes may be modeled as isothermal or non-isothermal. In the latter case, the heat transfer in the model is calculated using a staggered approach. Heat transfer between tool and workpiece may also be modeled. â€¢ The contact between the tool and the workpiece is modeled using an efficient direct approach, and different types of friction models and frictional boundary conditions may be applied. â€¢ Input and output data format of the program is compatible with some commercial FE and graphic codes, so that file transfer, pre- and post processing are simple. 4.3.1 Main Program The main program is the driver for various modules and handles input/output tasks, as described in figure 4.6. As input to the model, the user specifies the type of analysis, elements, boundary conditions, loading, material properties, thermal properties, convergence criteria, etc. in addition to data related to contact and heat transfer analyses. The program can import FE model created by some commercial codes such as NISA. The incremental analysis is started by initializing the variables and proceeds by calling various modules and managing the flow of data between them, until the analysis is completed. The main program also generates output files for post-processing. Among the output files are those which can be imported into the graphics software TECPLOT [118], for generating mesh configurations, contours of variables, velocity fields, load-displacement curves and animation files. Such data may also be stored in intermediate stages of the analysis at a period specified by user. 4.3.2 Deformation Module The deformation module is at the core of the FE program. This module, which is based on the original development of Wang [39], has been extensively modified in this work by including rate Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 96 Start ~~^) Input mesh, loading contact and thermal data Initialize variables * Apply load increment Solve for velocities and stresses Update contact conditions and apply frictional load Yes Contact conditions ^changed?_ No Output results Figure 4.6: Flowchart of the program structure and thermal effects, implementing general elastoplastic material models, expanding element library, and modifying mesh motion and variable updating procedures. The main function of the deformation module is to obtain velocity and stress solutions for an applied load increment. Figure 4.7 shows the flowchart of this module and flow of variables in an increment of the analysis. Some of the procedures used in this module is briefly described in the following sections. Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION i^Step n: from main program^) Determine mesh motion coefficients a, B for this step Initialize variables u, Ao Apply load increment Form element stiffness materices Eliminate mesh velocities, assemble and solve F E equations for velocities Find t + A t a from stress integration and return mapping Calculate residual forces and check equilibrium Apply residual forces as load increment Update reference configuration r Update material associated variables to mesh points i r (^Return to the main program^) Figure 4 . 7 : Flowchart of the deformation module 97 Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 98 i) S t i f f n e s s M a t r i x P r o c e s s i n g Following elimination of the mesh velocities on the element level and assembling the element stiffness matrices, the global unsymmetric ALE (K + K L stiffness matrix is obtained as N L )v = f (4.18) In this equation, the rate effects are included by substituting equation (3.44) in the linear term of stiffness matrix, i.e. K L = J B 5-dV = J B C B d V - j B (CD* - T) dV T T (4.19) T where B is the strain-displacement matrix. The second term on the right hand side is the correction term due to rate and thermal effects, which acts as a pseudo-load effect. It is noted that this matrix includes terms involving spatial derivatives of Cauchy stress tensor Such terms arise from mesh motion in ALE. It is, therefore, necessary to obtain stress gradients to form the stiffness matrix. This is a challenging issue that may give rise to instability of the solution process. A possible treatment of such term is to apply Gauss theorem to the appropriate term [119]; / Jv w <Tij 8u j dV = / w (TijUk Suij dS k ii: Js k Jv w ,k o-ijtik Suu dV k Jv w a- Suij dV k k (4.20) However, second order spatial derivatives of displacements will appear in this equation, which in turn, require that the trial functions for displacement field be C continuous. Alternatively, 1 the stress gradients may be obtained by way of interpolation within an element; Since the stresses are generally calculated at integration points, they should be first extrapolated to nodal points before equation (4.21) may be used for obtaining stress gradients. In addition, Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 99 different stress values are obtained at a node when extrapolated from each of the elements sharing that node, and a weighted averaging has to be used to smooth the stress values, which entails further approximation. Such approximations usually yield reasonably accurate and stable results in areas of low or moderate stress gradients. But, if the stress gradient is large, this scheme may result in large errors that may diminish solution accuracy after sufficient number of time steps. Since the stress gradient is an element quantity, it is advisable that unaveraged nodal stresses are used in calculation of this quantity to reduce the amount of approximation involved. In other words, the element nodal stresses are stored in each iteration to be used for stiffness calculation in the iteration that follows. Another consideration in this respect is that, similar to mixed formulation FE problems, it is preferable to use interpolation functions for stresses which are lower order than those used for interpolation of displacement variables. ii) Solution Procedures Since the stiffness matrix in ALE analysis is unsymmetric, an unsymmetric equation solver is needed. The solver used in this work is an unsymmetric front solver. Another feature of the solver is its ability to handle coupled displacement boundary conditions. A direct method for imposing coupled DOF boundary conditions is employed which constrains a node to move along a given surface. In this approach, a transformation matrix is defined which maps between global and tangential coordinate systems at a node. Using this transformation matrix, the stiffness matrix for the pertinent node is transformed to local coordinate system on the element level. Then, by substituting this equation in the global stiffness matrix, normal velocity DOF of the constrained node is eliminated from the equations. After assembling and solving the equations, a tangential velocity is obtained for the coupled node which may, in turn, be transformed to global coordinate system by an inverse transformation. Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 100 iii) Integration of Constitutive Equations The rate constitutive equation in finite deformation analysis is often expressed in terms of an objective stress rate and the rate of deformation tensor. Upon calculation of the rate of deformation tensor from velocity solution, the constitutive equation should be integrated to obtain the stresses. In an incremental solution procedure, the constitutive equation is numerically integrated over a time step to obtain the stress increment in the step, which is then added to the stress at the end of previous step to obtain total current stress. Due to the fact that the configuration of the body changes in large deformation analysis, the incremental objectivity of the stress integration should be ensured. In addition, the integration process should uphold the plastic consistency of the stress, i.e. the stress state should remain within the elastic range or on the yield surface. The incremental objectivity of stress tensor within an increment is ensured by using an objective stress rate such as Jaumann or Truesdell rates of Cauchy stress tensor, so that; CijkiDkiAt = frijAt (4.22) In this equation, the rate of deformation tensor may be referred to an intermediate deformation state a, (4.23) which depending on the choice of parameter 0 < a < 1 will result in explicit or implicit integration schemes. Other procedures for upholding objectivity of integration scheme may also be found in the literature. Hughes and Winget [121] proposed an algorithm for integration of Jaumann rate of stress which is based on pure geometrical concepts; (4.24) where t+ t Rij is the rotation tensor for the current increment. The first term in this equation ensures the objectivity of accumulated stress at time t with respect to rigid body motion. Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 101 Similar algorithm is used by Pinsky, et al. [122] for Truesdell rate based on mathematical tensor transformation. Wang and Gadala [120] give an alternative derivation of integration algorithm based on physical definition of the Cauchy stress. In their derivation, the deformation is decomposed into a pure rigid body rotation and a pure strain, resulting in t + A t where Fij t+At = a t g^ X i "a = t + t+Atp t A t F im cr F t+ At mn t nj +C D At l ijkl kl (4.25) is the deformation gradient tensor. For a pure rigid body rotation, equations (4.24) and (4.25) lead to same results. The latter equation is shown to be equivalent to the integration of Truesdell stress rate. The incremental plastic consistency of the integration algorithm is usually achieved by a return mapping scheme which retains the plastic stress state on the yield surface. Many such algorithms are available in the literature. A commonly used procedure is the operator split algorithm in which the stress state is updated in two steps; an elastic predictor followed by a plastic corrector. A particular problem in return mapping procedures is theflowdirection along which the stress state is projected onto the yield surface. This problem is discussed by Ortiz and Popov [123] who cast the various algorithms into two generalized families; the trapezoidal family and the mid-point family. Ortiz and Simo [124] describe an iterative procedure which includes successive updating of plastic states andflowdirections until the plastic stress is appropriately mapped onto the yield surface. In this work, a trapezoidal return mapping is used in which the mapping is divided into smaller steps for increasing the accuracy of mapping. iv) Convergence Criteria In an incremental FE solution, various convergence criteria may be used. In this program, the termination criterion is based on convergence of three variables; displacement, force, and energy [39]. In the displacement criterion, the ratio of the norm of total displacements in current iteration to the same norm in thefirstiteration within the increment is compared to a given Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 102 tolerance. The force criterion is based on the ratio of norm of residual forces in the current iteration to the maximum norm of residuals in previous iterations, while in the energy criterion, the work done by iterative residual forces is compared to the work of initial residual forces. It is generally accepted that a combination of any two criteria is effective for most engineering problems [125]. 4.3.3 Contact Module The contact module is called from the main program after the solution for load increment is converged in the current step. The current implementation is based on direct contact algorithm between rigid-flexible bodies (section 3.4.3). The main function of this module is to update the status of contact at contact nodes and calculate and apply frictional forces for the next increment. When the control is returned to the main program, a control parameter determines if the solution is to be re-iterated or it can proceed to the next load step. Also, this module provides initial contact configuration and tool profile at the start of the analysis. The tasks performed in this module are described below. i) Defining Tool Profile and Attributes, and Initial Contact At the start of the analysis, the geometry and attributes of the tool have to be defined. Upon calling this subroutine from the main program, a user input file containing the geometry of the tool is read and its attributes are determined. In the input file, the tool is defined as a series of geometric entities such as lines and arcs which are connected together, by specifying their geometric properties such as start, end, or center points. Only a part of the tool that may come into contact needs to be defined. Each geometric entity defined in this way is called a segment. The program internally calculates the geometric attributions of the segments, such as line and arc equations, radius, slope, normal vector, etc., which will be used in detecting penetration, Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 103 contact, or release of nodes. The order of definition of segment properties determines which side of the segment corresponds to the physical area of the tool, and the normal vectors are defined outward with respect to the tool. Any number of objects may be defined as contacting the workpiece, e.g. the upper and lower dies in a forging process. A list of boundary nodes which may come into contact with the tool is provided by the user at the start of the analysis as potential contact nodes. Throughout the analysis, these nodes are continuously checked to detect changes in their contact conditions. Also, initial contact between workpiece and tool may be defined at the outset of the analysis. ii) Checking for Penetration and Node Release In each incremental step, the contact status of each potential contact node is checked to determine if any of its contact conditions has changed. The contact status of a node is identified as one of the following three states: not in-contact, just penetrated, and in-contact. At the end of a converged step, the algorithm checks the contact status of the nodes based on their current status. If a node is currently not in contact, the algorithm determines if the node has penetrated the physical area of the tool. Such determination is based on checking the position of the node with respect to the perimeter of the tool. If it is determined that the node is within the tool area, the algorithm associates the node to its closest tool segment, and determines a penetration displacement vector which returns the node to the surface of the tool along a path normal to the associated tool segment. The contact status of the node is then changed to 'just penetrated'. Next, aflagis set for re-iteration of the step, in which the penetration vector is super-imposed on the body as a displacement/velocity boundary condition. The implicit assumption here is that the solution of current step can be obtained by way of superposition of the two solutions; one due to incremental load in this step, and the other due to velocity boundary conditions which Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 104 is applied to offset the penetration. Although such assumption may not be strictly correct for a plastic contact analysis, it is a reasonable assumption if the steps are not very large, because the penetration vector is usually small. A more elaborate approach would be to re-calculate the step by breaking the step size into sub-steps which bring the contacting nodes to the surface of the tool one at a time [96]. The drawback of this approach is that if more than one node contacts in the same step, it is not always trivial to calculate the step size that brings into contact one node at a time. Once the penetration vector is applied, the node status is changed to 'in-contact', and the node is allowed to slide on the tool segment. The nodal forces obtained from compensation of penetration are used to determine the status of the contact, and to find the direction of frictional force in the case of sliding contact. A right-hand local coordinate system for the contact node is defined based on its associated segment. Since nodal positions change continuously, the local system and associated segment of each contact node has to be updated in each step. For a node currently in contact, the algorithm checks if the node has been released from contact as a result of deformation in the current step. A node is considered released if the normal contact force is no longer compressive. The status of such a node is set to 'not in-contact', and all related variables, e.g. contact forces, associated segments, etc. are reset. iii) Applying Frictional Forces In the direct contact algorithm, a node 'in contact' is constrained to remain on the surface of its associated tool segment. This constraint is applied by coupling the nodal degrees of freedom such that the normal component of nodal velocity is zero. The frictional force on a sliding node is applied as an external tangential force in the direction opposing the sliding nodal velocity. The magnitude and direction of frictional force has to be based on nodal normal forces and velocities at the last incremental step, but a re-iteration of the step is necessary if a change in Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 105 the direction of the velocity is detected. The sticking frictional condition may be applied as an ever-slippage case, and the piecewise continuous transition between sticking and sliding can be replaced by a smooth function to avoid instability. For metal cutting application, a nonlinear equation is frequently used in the literature [24]; (4.26) where k, a , and r are shear strength, normal stress and frictional stress, respectively, and n p is the friction coefficient. Shear strength of the material at a nodal point is obtained from the local value of flow stress. Equation (4.26) provides a smooth transition between sticking and sliding types of friction. Depending on the ratio of normal stress to the shear strength, the function in the brackets varies between 1 and 0, representing extreme cases of shear friction and threshold of release, respectively. For small normal forces, the equation approaches Coulomb's friction law . This can be verified by expansion of the exponential function in Taylor series, and neglecting higher order terms. Alternatively, Coulomb and shear friction types may be directly implemented. iv) Re-Iteration of Steps Detection of any changes in contact conditions following a solution triggers a re-iteration of the current step in a loop that continues until no further change is detected. Here, small negative and positive cut-off forces are introduced to create a buffer zone between node release and repenetration. This renders the algorithm to be robust and makes it to converge in a few steps. In our numerical experiments, this buffer range is proved to be effective in preventing oscillations due to fluctuation in the force quantities that are physically insignificant. An alternative way to deal with such a problem is to base the status of contact of a node on the distributed traction on its adjacent contactor segments [95]. Chapter 4. 4.3.4 NUMERICAL PROCEDURES AND IMPLEMENTATION 106 Heat Transfer Module The thermal analysis is carried out following the mechanical analysis in an incremental step. The deformation process within the increment is solved as an isothermal process with the temperature distribution at time t. The plastic and frictional work dissipated in this step are mostly converted to heat, and will act as volume or surface heat sources for the thermal process. In the ALE analysis, the heat transfer part is performed based on material mesh and before updating mesh configuration. Once the new temperature field is obtained, it is updated to mesh points using the variable updating scheme. The general capabilities of the heat transfer module are listed below. â€¢ Transient or steady state thermal analysis may be modeled. â€¢ A one-parameter time stepping method is used, which allows for implicit or explicit solutions. â€¢ Various types of thermal boundary conditions including boundary convection, volume heat source, initial temperaturefield,specified temperature, and boundaryfluxare available. Radiation is neglected to avoid introducing nonlinearity. â€¢ The capacitance matrix may be modeled as consistent and lumped. â€¢ Time step in the transient heat analysis may be different from that of the deformation process. In other words, the deformation time-step may be divided into smaller sub-steps for the thermal analysis. This is particularly useful if an explicit thermal solution is selected. â€¢ The same mesh used for deformation analysis is also employed in thermal analysis. â€¢ For a rigid tool, the heat transfer between tool and workpiece may also be modeled. In this case, the FE mesh covering the tool is only active in the thermal part of the process. Chapter 4. NUMERICAL PROCEDURES AND IMPLEMENTATION 107 â€¢ Thermal properties of the material may be temperature dependent. These properties are updated in between incremental steps according to the equation or data table provided by the user. The thermal properties and analysis choices is obtained from user files at the start of analysis. Chapter 5 N U M E R I C A L RESULTS A N D DISCUSSION In this chapter, results of numerical simulation of metal cutting processes are presented. First, a brief discussion on the determination of material property data as input to a cutting simulation is presented. In section 5.2, thefiniteelement model and the mesh motion scheme for ALE cutting simulations are described. Next, results of numerical simulation of cutting process under various cutting conditions are presented, and a parametric study on the effects of cutting conditions on the chip formation process is carried out. In particular, the influence of tool edge geometry is studied in section 5.5, where machining with tools of sharp, chamfered and worn edges are simulated. The numerical results are verified through comparison with experimental results in the literature for similar cutting conditions. 5.1 ON T H E DETERMINATION OF MATERIAL PROPERTIES FOR C U T T I N G SIMULATION The reliability of FE simulation results is to a large degree dependent on the validity of the assumptions made and the material properties that is input to the simulation. In particular, the flow stress model of the work material and the frictional properties of the tool-workpiece pair are the two most sensitive inputs to be provided. While in many solid mechanics applications, the data obtained from standard material testing such as tension, compression, and torsion tests provide a good representation of material response, such data are generally inadequate for extreme conditions such as those that exist in an ordinary machining process. In fact, even at moderate cutting speeds, strains, strain rates and temperatures in the cutting zone 108 Chapter 5. NUMERICAL RESULTS AND DISCUSSION 109 are orders of magnitude higher than what can be measured by many of the available testing equipment. Nevertheless, due to absence of better data,flowstress models obtained from standard material tests have been largely used in the literature. Some early works have used very simplistic plasticity models such as linear hardening models, sometimes even without rate and temperature effects. More rigorous flow stress data at high strain rate and temperature have also been obtained from material tests such as impact compression testing (Hopkinson bar) or ballistic testing. Usui and Shirakashi [24] used the Hopkinson bar test to obtain flow stress at strain rates up to 2000 s _1 and temperatures up to 800Â° C. But, the rates attainable in this test are still smaller than the strain rates in practical cutting processes, and the results have to be extrapolated in order to be used for such applications. Such extrapolation may amount to gross approximation, as the material response might be very different at higher rates and temperatures. In addition, sometimes under a combination of cutting conditions, metallurgical phase transformation may occur which may not be reproducible in a material testing method." The inadequacy of standard material testing methods in obtaining material response at conditions similar to those present in a cutting process has prompted many researches to turn to the cutting process itself as a way of determining material properties in such extreme conditions. Since 1970's there has been several attempts at obtainingflowstress data directly through performing cutting tests under carefully controlled conditions. Most notable in this respect is the work of Oxley and his coworkers [7] who cast the material constitutive equation into a power law form a = c^eâ„¢, where a-^ and n are parameters dependent on a variable called "velocity modified temperature" which combines the effects of strain rate and temperature. These parameters were obtained based on experimentally measured force and chip thickness data. The strain rate was calculated from an empirical equation derived from analyzing the deformationfieldsin explosive quick stop experiments, and temperature rise in the primary and secondary deformation zones is calculated based on a model proposed by Boothroyd [126]. Stephenson [127] applied statistical methods on the data from a set of carefully designed Chapter 5. NUMERICAL RESULTS AND DISCUSSION 110 cutting and drilling tests on Al 6061 and 1018 steel and concluded that there is a weak correlation between geometric stress and strain/strain rate measures, but a strong correlation between material properties and cutting conditions. Hefittedpolynomial of several parameters to obtain stress, friction coefficient and cutting force in terms of velocity, feed rate and rake angle. Ozel and Altan [128] combined experimental data withfiniteelement simulations to determine flow stress and frictional conditions for cutting process. In their approach, average values of strain, strain rate and temperature were obtained from FE simulations in an iterative process in which theflowand friction parameters were changed until forces obtained from simulation matched those measured in cutting test under similar conditions. A least square method was then used tofitthe results for various cutting conditions, and to obtain the constitutive equation of the material in terms of strain, strain rate and temperature. The difficulty with the use of cutting process for determination offlowstress is in the fact that variables such as stress, strain and strain rate used in aflowmodel cannot be measured directly on the shear plane. Therefore, these variables are commonly obtained from measurable quantities such as tool forces, chip thickness and contact length, often based on a geometrical deformation model. Although some of these variables are known to vary by orders of magnitude over the length and width of the shear zone, only average values may be obtained from such models. The small size of the cutting zone and the high speed of the deformation makes it difficult to obtain any accurate measurement of such variables on the shear zone. As a result, a flow stress model obtained from cutting experiments involves various levels of approximation. One variable that can be measured directly and with reasonable accuracy is the average temperature on the shear plane. However, such measurement is not attempted in many of the reported works in the literature, perhaps due to the difficulty in accessing the very small shear zone. Recently, Lei, et al [129] used an infrared camera to measure average temperature at the shear zone. Combined with an analytical method for determining the strain rate, they obtained a power law model which is explicit in strain, strain rate and temperature. Chapter 5. NUMERICAL RESULTS AND DISCUSSION 111 From this brief review, it may be concluded that while cutting process itself is the best testing method for determining the material and frictional response under such extreme conditions, such an approach entails considerable degree of approximation and should be carefully verified. A more accurate determination of material response is only possible if more direct and detailed information can be obtained from the shear zones. In the absence of such data, the finite element solutions of cutting process, especially at high speeds, is often subject to errors due to inaccurate input data. 5.2 C U T T I N G MODEL A N D MESH MOTION S C H E M E In this section, the finite element model for ALE simulation of chip formation in orthogonal metal cutting process under plane strain conditions is described. In cutting of ductile materials, the chip formation process involves continuous flow of material around the tool tip as well as unconstrained evolution of free boundaries of the chip. The advantage of using ALE approach for the cutting process is that the flexibility of mesh in this approach can be used to efficiently model both of the above features in a single analysis. The role of mesh motion strategy is, therefore, crucial to the success of such analysis. The general strategy is to use an Eulerian type of motion around the tool tip so that neither node separation nor high density mesh and frequent remeshing is required, while a Lagrangian type of motion is used in modeling of unconstrained flow on the free boundaries, which defines the shape and size of the chip. A cutting process is generally started by indentation of the tool by the workpiece followed by a transient stage known as the incipient chip formation, until the steady state is reached and process variables are stabilized. Since the steady state of cutting is of interest in this work, the incipient stage of cutting is not modeled and the analysis starts from a model which includes a pre-assumed chip of arbitrary size and shape. Although a similar pre-formed chip is commonly used in Eulerian cutting simulation, such assumption in ALE analysis is different Chapter 5. NUMERICAL RESULTS AND DISCUSSION 112 in two essential respects; 1. this assumption is not an inherent requirement of the ALE, and is made for numerical convenience only. In other words, it is conceptually possible to design a mesh motion scheme for modeling all stages of the process. Nevertheless, the assumption of initial chip geometry simplifies the mesh motion scheme considerably, and is numerically convenient if only steady state cutting is of interest. 2. in ALE analysis, the initial chip shape does not affect itsfinalshape, because the deformation process automatically corrects the assumed geometry as part of the solution, whereas in Eulerian analysis the chip size cannot be modified, unless it is done outside the deformation process on a trial and error basis. Nevertheless, there will be a transient stage in the ALE analysis, in which the initial chip will evolve into its proper shape and size, and the variables reach steady state. Figure 5.1 shows a schematic view of the ALE model used in the simulations. In this model, the tool is assumed as rigid and stationary, and the cutting action is simulated by advancing the workpiece against the tool through applying velocity boundary conditions on its left and bottom surfaces. The assumption of rigid tool is common in literature on simulation of chip formation in metal cutting process, because the elastic deformation of the tool is negligible compared to the very large plastic deformation in the chip, and is not expected to affect chip formation process significantly. The stressfieldin the tool is, however, important in the study of tool life and wear. Since the current study is focused on the chip formation process, the tool is assumed as rigid. On the other hand, there is significant heat transfer between the tool and the workpiece which affects the cutting process. Therefore, the tool is included in the thermal part of the ALE simulation. The direct node to surface contact algorithm is used for modeling the interaction between tool and workpiece. Chapter 5. NUMERICAL RESULTS AND DISCUSSION 1 B 113 Â° //. A Figure 5.1: Schematic view of the A L E model illustrating the mesh motion scheme Next, the mesh motion scheme for the ALE model is described. The workpiece is divided into a number of regions as shown infigure5.1. The central region of the mesh, marked as region A, is assumed to remainfixedthroughout the analysis. This region includes the area around the tool tip and part of the interface between the tool and the workpiece, though the nodes on the contact or free surfaces are Lagrangian in the normal direction. The inflow of material is modeled in region B where a velocity boundary condition is applied at its left surface. In the course of deformation, this region gradually and steadily shrinks as the material migrates into region A. The outflow of material through the machined part is modeled by steady expansion of regions C. Similarly, the increase in the chip length and chip curling is achieved by expanding the mesh in region E. Thefinalconfiguration of the model is also shown infigure5.1. The mesh in region E is initially very dense. In regions B, C and E, all external boundaries are Chapter 5. NUMERICAL RESULTS AND DISCUSSION 114 Lagrangian in the normal direction, but are moved tangentially in a way that certain ratio of element size is maintained and the mesh is denser closer to the deformation zones. The mesh velocities in the interior of these regions are determined by transfinite method based on position of nodes on their borders. Region D is the most important region, because it is bounded by a free boundary on one side and a fixed boundary on the opposite side. The unconstrained flow on the free boundary at the back of the chip is crucial in realistic evolution of the chip to its proper size and shape, while the large deformation occurring in this region calls for a scheme that avoids mesh distortion. To achieve this, the nodes on the free boundary are marked as Lagrangian in the normal direction and Eulerian in the tangential direction. This can be easily implemented by setting C = 0 in equation 4.9. In this way, the external free boundary is developed in a natural way, and at the same time, the mesh is continuously readjusted by transfinite method to adapt to the modified shape and maintain an optimal mesh throughout the analysis. It is interesting to note that region E is treated as a quadrilateral zone in the transfinite method, which demonstrates the capability of this method in modeling boundaries with discontinuous slopes and various curvatures. Figure 5.2 and 5.3 show typical mesh configurations at initial and various stages of chip formation. It may be seen that with this mesh motion scheme, the flow of material around the tool is modeled without resorting to node separation or frequent remeshing, while on the other hand the chip evolves naturally and automatically to itsfinalshape and size. The density of the mesh around the shear plane remains unchanged throughout the analysis irrespective of large material deformation there. Apart from the tail of the chip which corresponds to the initially assumed chip thickness, the chip is seen to acquire a fairly constant thickness after sufficient length of cut. In order to demonstrate the independence of the chip size from the initial assumption, the initial and steady state chip shapes in a typical simulation are overlaid in figure 5.4, where it is seen that the steady state chip thickness is determined solely by the deformation process. Chapter 5. NUMERICAL RESULTS AND DISCUSSION 115 Figure 5.2: Initial model and various stages of chip formation in cutting simulations (continued on figure 5.3) Chapter 5. NUMERICAL RESULTS AND DISCUSSION Figure 5.3: Various stages of chip formation in cutting simulations (continued from fig 5.2) Chapter 5. NUMERICAL RESULTS AND DISCUSSION 117 Figure 54 .: Overlaying the initial and the steady state meshes shows that the chip thickness is independent of the initial assumption of its shape. Based on the above FE model, the results of cutting simulations under various cutting conditions and with different tool edge geometries are presented in the following sections. The first set of simulation results pertains to a test problem that has been shared among some leading researchers in thisfield,and is used here as a benchmark problem for the purpose of validation of numerical procedures. In the second set of results, a parametric study of the influence of cutting conditions on the solution variables is carried out. In the third set of results, the effects of geometric features of the cutting tool is investigated by simulating tools with sharp, chamfered, or worn (blunt) edge. Finally, the results for orthogonal cutting of a Titanium alloy are presented. In each case, the material and frictional characteristics are adopted from the literature, and the solution variables are verified against the available experimental results. Chapter 5. NUMERICAL RESULTS AND DISCUSSION 5.3 5.3.1 118 T H E B E N C H M A R K CUTTING P R O B L E M Input Data Van Luttervelt, et al. [130] report on a round-robin exercise in which the results of Lagrangian and Eulerian cutting simulations by several researchers for a set of material and cutting conditions are compared with experimental data. Here, the same cutting process is simulated as a benchmark problem for validation of our numerical approach. Two other simulations with same material but different cutting conditions are also presented. The workpiece material is low carbon free cutting steel cut by P20 cemented carbide tool. The input data, cutting conditions, and experimental results, reported in tables 5.1 and 5.2, as well as the flow stress model are adopted from [130]. The constitutive equation for this material as a function of strain, strain rate and temperature is obtained from Hopkinson bar tests and has the general form of Y = AT [-^â€”\ (5.1) l \iooo; k ' where for the low carbon free cutting steel, the parameters are given by A = 910exp (-0.0011T) + 120exp {-4 x 10 (T - 280) } _5 2 + 50exp {-nr (T- 260) } 5 2 m = 0.16exp(-0.0017T)+ 0.09exp{-3 x 10- (T- 370) } 5 2 n = 0.018 + 3.8 x 10~ T 5 in which Y is in MP a and T is in Â°C. The thermal properties of material are functions of temperature as shown in table 5.1, but they are assumed to remain constant within an increment, and are only updated between increments. Friction on the rake face is modeled using the nonlinear equation (4.26). The shearflowstress of the material at a nodal point is obtained from its localflowstress and the parameter p is obtained from experimental measurement of normal and frictional stresses using a split tool dynamometer. This measurement shows that the Chapter 5. NUMERICAL RESULTS AND DISCUSSION Young's Modulus Poisson Ratio Density Thermal Conductivity Specific Heat 119 200 GPa 0.3 7850 Kg/m 62 - 0.044T W/mÂ°C 450 + 0.38T J/KgÂ°C 3 Table 5.1: material properties for the low carbon free cutting steel normal and frictional stresses are almost equal when the normal stress is below 360 MP a, but the frictional stress saturates at this value at higher normal stresses. Based on this experiment, a value of p, = 1.0 is suggested in reference [130]. Equation (4.26) models both sticking and sliding types of friction on the rake face and provides a smooth transition between these two states. 5.3.2 Simulation Results The first simulation is carried out for a rake angle of 0Â°, Cutting speed of 150 m/min and uncut chip thickness of 0.1 mm/rev. The experimental data reported in table 5.2 is used for verification of the simulation results. In this simulation, the FE mesh consists of about 1200 linear triangular elements and 600 nodes. The relatively small number of elements used reflects the efficiency of ALE analysis, because the density requirement of the mesh is only dictated by expected gradients of stresses and strains, and not by other factors such as node separation or tool penetration. A total cutting length of 1.5 mm is simulated which is sufficient for the simulation variables to reach the steady state and to observe chip curling. The simulation steady state is denned as a state after which no significant change in process variables in the cutting zone is observed. At this stage, the simulation variables may be considered as representing the actual steady state cutting process. The above cutting length is achieved in 1200 incremental steps. For convenience, the total strain and strain rate variables in equation (5.1) are replaced by plastic strain and plastic strain rate, respectively. This approximation is justified because Chapter 5. NUMERICAL RESULTS AND DISCUSSION 120 the elastic strains are much smaller than their plastic counterparts in the cutting zones. The tool area is also meshed, but this mesh is only active in the thermal part of the analysis in each incremental step. Although the tool is practically considered as sharp, a finite edge radius of 0.01 mm is added to the tip, which makes it more realistic. The added edge radius facilitates the modeling of the flow of material around the tool edge, and allows for observation of the interaction of the tool with the newly machined surface, known as the "tertiary deformation zone". The steady state configuration of the model is similar to figure 5.3, where it may be seen that the mesh chip formation remains regular throughout the analysis. The distributions of solution variables including effective stress, effective plastic strain, effective plastic strain rate, and temperature in the cutting zones are shown in figure 5.5. The effective plastic strain rate peaks on the tool tip where it reaches as high as 7 x 10 s' , but takes much smaller values in 4 1 other parts of the deformation zones. It is noted that the strain rate is practically zero outside the deformation zones, and thus, the size of the primary and secondary deformation zones may be estimated from this distribution. The effective stress takes its maximum on the shear plane as expected. However, the stress values are relatively smaller on the rake face, in spite of higher strain and strain rates on the tool surface. This may be explained by considering the high temperature field on the rake face which results in significant material softening in this area. In general, thermal softening competes with strain and strain hardening effects in the cutting zones, and, depending on material properties, it may become the dominant factor on the rake face. Such phenomenon occurs in cutting of low carbon steel considered here due to the significant dependence of its flow stress on temperature. For example, at strain of 1 and strain rate of 1000 s , the flow stress drops by 400 MP a when temperature increases from -1 room temperature to 600Â° C. It is also generally known that at higher strain values, the strain has little effect on the deformation process [7]. The effective plastic strain is very high on the rake face, around the tip and on the machined surface, the latter due to ploughing action of the Chapter 5. NUMERICAL RESULTS AND DISCUSSION 121 Figure 5.5: Contours of process variables in simulation (1) [a = 0Â°, V = 150m/min]; (a) effective stress, (b) effective plastic strain, (c) effective plastic strain rate, and (d) temperature distribution 0 Chapter 5. NUMERICAL RESULTS AND DISCUSSION 122 Tool Advance (mm) Figure 5.6: Cutting and thrust forces predicted versus cutting length in simulation (1) in comparison with experimental values [130]. tool in the tertiary deformation zone. Also, the temperaturefieldis shown in both the tool and the workpiece. Temperature reaches its peak value on the rake face of the tool and is also high on the machined surface. A large temperature gradient across the width of the chip is observed, which is due to larger plastic straining and friction on the rake face. Figure 5.6 depicts cutting and thrust forces on the tool versus cutting length and compares them with experimental results. The difference between the experimental and simulated force components can be attributed to many factors, among them, the material and friction models and their parameters, as well as the size of edge radius, which is not reported for the experimental results in the reference. Despite the differences, the.forces on the rake face are reasonably close to the experimental values. Other simulation results are also reported in table 5.2. The shear angle, contact length, and Chapter 5. NUMERICAL RESULTS AND DISCUSSION c*o Case Experiment (1) Simulation (1) Simulation (2) Simulation (3) 0Â° 0Â° 0Â° -5Â° V w m/min F N/mm F N/mm 150 150 450 150 174 207 190 225 83 96 85 107 r c t 123 <P 18.8Â° 22Â° 24Â° 25Â° mm Tmax Â°C 0.60 0.55 0.48 0.40 590 571 738 570 Table 5.2: Experimental and simulation results. The columns represent rake angle, cutting speed, cutting force, thrust force, shear angle, contact length, and maximum temperature, respectively. Note that cutting conditions in simulations (2) and (3) are different from those in experiment and simulation (1). The uncut chip thickness in all cases is 0.1 mm. maximum temperature obtained from simulation are all close to the experimentally measured values, verifying the closeness of simulation to the actual process. Figure 5.7 shows the distributions of normal and frictional stresses as well as temperature on the tool-chip interface in comparison with reported experimental results [130]. For both normal and frictional stresses, the simulated distributions show similar trend to that of the experiment, though the simulated values are in general higher, which explains the larger tool forces predicted. The predicted temperatures, on the other hand, are in general lower than the experimental values, and temperature peak occurs closer to the tool tip. In addition, the predicted temperatures drop with a faster pace toward the end of contact length than what experimental distributions shows. Overall, the comparison between simulated and experimental variables verifies that the FE model provides a fairly good representation of the actual process in this case. Similar comparisons of predictions from various FE codes with experimental results for this particular case are reported in reference [130], which show a wide scatter in simulation results for each of the above variables. To study the effects of change of cutting conditions on the process, two more simulations are presented here, one with higher cutting speed, and the other with negative rake angle. In the former case, marked as simulation (2), the same cutting conditions as in simulation (1) is used Chapter 5. NUMERICAL RESULTS AND DISCUSSION 124 but the cutting speed is increased to 450 m/min. As shown in contours of figure 5.8, there is a marked difference in the values of strain rate and temperature compared to simulation (1). The maximum temperature in this case is about 150Â°C higher than that in case (1). The forces on the tool, contact length and shear angle slightly decrease for higher cutting speed. Finally, A rake angle of â€”5Â° is used in simulation (3) to examine the capability of ALE simulation in handling negative rake angle tools. Other cutting conditions are the same as those in simulation (1). Figure 5.9 shows the chip formation and temperature distribution in this simulation and 1800 O 600 5 500 E 400 0.2 0.3 C o n t a c t length (mm) 0.4 0.2 0.3 C o n t a c t length (mm) 0.4 0.5 0.6 0.6 â€” â€¢ â€” Experiment - 0.1 0.2 - -Simulation _ . _ ,0.3 â€ž , . Contact length (mm) 0.4 0.5 Figure 5.7: distribution of predicted normal and frictional stresses and temperature on the chip-tool interface in comparison with experimental results [130]. Chapter 5. NUMERICAL RESULTS AND DISCUSSION Plastic Strain Rate 1/sec 164222 153274 142326 131378 120429 109481 98533.1 87585 76636.9 65688.8 54740.6 43792.5 32844.4 21896.2 10948.1 125 Temperature Figure 5.8: Contours of strain rate and temperature for high speed cutting simulation (2) [oo = 0Â°, V = 450 m/min] TEMP 534.15 499.71 465.26 430.82 396.37 361.93 327.48 293.04 258.59 224.15 189.70 155.26 120.81 86.374 51.929 Figure 5.9: Chip formation and temperature distribution in simulation (3) [a = -5Â°, 0 V â€” 150 m/min] Chapter 5. NUMERICAL 126 RESULTS AND DISCUSSION other results are reported in table 5.2. While cutting and thrust forces, and shear angle increase in this case, contact length is reduced due to apparent faster curling of the chip. The maximum temperature is in the same range as that of simulation (1). It should be noted that experimental data are not available for comparison with simulations (2) and (3). 5.4 T H E EFFECTS OF C U T T I N G CONDITIONS ON CHIP F O R M A TION PROCESS In this section, a parametric study on the influence of cutting conditions on the chip formation process is presented in which the results of several numerical simulations are reported and verified through comparing with corresponding experimental data. . In the following,firstthe experimental data are presented, and then, the simulation results are discussed and compared with experimental data. 5.4.1 Experimental Data The experimental data and material and friction models used in this study are adopted from the work of Childs and Maekawa [131] which studies the effects of cutting variables on the cutting process. They carried out a series of turning tests on hot forged bars of BS970 708M40 low alloy steel under a range of cutting conditions. Cutting forces on the tool and chip thickness were measured and their components on the plane normal to the cutting edge were calculated from the measured data. Theflowstress of the material is obtained from Hopkinson bar tests in the general form of equation (5.1), where its temperature dependent parameters are given as A = 1500exp(-0.0018T) + 380exp{-l x 10" (T- 445) } 5 2 + m = 0.135exp(-0.0012T) + 0.07exp {-2 x 1(T (T - 465) } 5 n = 0.017 + 6.8 x 10 T _5 2 Chapter 5. NUMERICAL RESULTS AND DISCUSSION 127 The friction model used here is also similar to that of the benchmark problem. The parameter p. in equation (4.26) has been estimated by measuring normal and friction stress distribution using a split tool dynamometer, which gives a value of 1.3 for this material. A cemented carbide tool with a clearance angle of 6Â°, an approach angle of 14Â°, and a nose radius of 0.25 mm is used in the turning operations. It is noted that although these turning tests are not strictly orthogonal, they are expected to closely approximate orthogonal cutting tests on the normal plane, because the nose radius and feed rate (0.25 mm and 0.125 mm/rev, respectively) are much smaller the width of cut (1.25 mm). Other material data is given in reference [131]. The thermal properties of the tool and workpiece are not given in the reference, and typical values for this type of material are adopted from handbooks for the simulations. 5.4.2 Numerical Model The range of cutting conditions for the numerical simulations in this study are given in table 5.3. The typicalfiniteelement model and the ALE mesh motion strategy are similar to what was described in section 5.2. A total cutting length of 1.75 mm is considered for each case, in which the simulation variables reach steady state condition and chip curling is observed. This cutting length is achieved in 1500 incremental steps. The tool is assumed to have a finite edge radius of 0.01 mm. The initial configuration and the formation of the chip after 0.5,1.0, and 1.5 mm of cutting for a typical case correspond tofigures5.2 and 5.3. It is seen that the chip is steadily formed, and the size and the curling of the chip are evolved in a natural unconstrained fashion, and the final chip has a fairly constant thickness. Chapter 5. NUMERICAL RESULTS AND DISCUSSION case Rake Angle Feed Rate (mm/rev) 1 2 3 4 0Â° 6Â° 12Â° 6Â° 0.125 0.125 0.125 0.254 128 Cutting Speed (m/min) 150 150 150 150 and and and and 300 300 300 300 Table 5.3: Cutting conditions used in the parametric study 5.4.3 Numerical Results and Discussion First, the general characteristics of the chip formation process in this set of cutting simulations are described. The results presented in this section pertain mainly to the case of cutting with 0Â° rake angle tool and feed rate of 0.125 mm/rev. The results of other cases are similar in nature, though different in magnitude. i) Material flow Figure 5.10 depicts the distribution of nodal velocity vectors in the cutting region. The streamlines clearly show the flow of material around the tool tip to form the chip and the machined surface. On the lower portion of the tool rake face, the material has relatively small velocity, because the chip sticks to the tool at this area due to existence of very high pressure between the chip and the tool. With reducing pressure further up on the tool, the chip velocity gradually increases in the sliding friction zone, until the chip is separated from the tool surface. ii) Strain and strain rate Typical distributions of equivalent plastic strain and strain rate in the cutting zone are shown in figures 5.11(a) and (b). For a sharp edge tool, the tip of the tool is a singular point at which the strain is theoretically infinite, similar to a crack tip. For a tool with a finite edge radius, the strain is expected to reach very high values at the tip, which may be captured in numerical simulation if a sufficiently fine mesh is used in this area. It can be seen in figure 5.11(a) that a Chapter 5. NUMERICAL RESULTS AND DISCUSSION I 1I M M M M fill I N I t t t t t I ! 1 1 Hi Hi- 129 Cutting Tool ' in j 11 11 / /; n Figure 5.10: Velocity field in the chip and the workpiece very large value of strain is attained in the immediate neighborhood of the tip, which rapidly reduces to much smaller values at the back of the chip. The average equivalent plastic strain on the shear plane is around 2.5, whereas the average strain predicted by shear plane theory for a sharp tool gives a value of 1.75 based on experimentally obtained shear angle. The strain is also very high on the rake face where the highly pressurized friction between the tool and the chip, assisted by thermal softening of the material, causes further straining in the secondary deformation zone. Also evident in the figure is the high magnitude of strain on the machined surface, which is the result of the interaction between the nose of the tool and the workpiece at the tertiary deformation zone. The edge radius of tool on the flank side acts as a tool with a large negative rake angle, ploughing the material before it leaves the tool surface. Similar pattern is seen in the distribution of equivalent plastic strain rate in figure 5.11(b). As expected, high values of strain rate are found in the primary and secondary deformation zones. The strain rate reaches its peak in the immediate neighborhood of the tool tip, where it has a magnitude of order 10 s 5 - 1 for the cutting speed of 150 m/min. Such high rate at Chapter 5. NUMERICAL RESULTS AND DISCUSSION (a) Equivalent plastic strain contour (b) Equivalent plastic strain rate contour SEQ 1200 1114.29 1028.57 â€” 942.857 857.143 771.429 685.714 600 514.286 428.571 342.857 257.143 171.429 85.7143 0 - Cutting Tool (c) Temperature distribution (d) Effective stress contour Figure 5.11: Contours of distribution of solution variables in the cutting zones Chapter 5. NUMERICAL RESULTS AND DISCUSSION 131 this point is due to the combination of cutting action of the tool and ploughing action of the tool flank. The strain rate is also high on the rake face suggesting substantial plastic straining of material in the area of sticking friction on the rake. The strain rate contour provides an indication of the size of the two deformation zones. On the other hand, there exists a high gradient of strain rate across the thickness of chip, as the strain rate decreases substantially toward the back of the chip. The profile of equivalent plastic strain rate on the primary shear plane is plotted in Figure 5.12(a) for two cutting speeds. As expected, the magnitude of strain rate is a strong function of the cutting speed, such that doubling the speed almost doubles the strain rate. 2.0E+05 -pa- V=150 m/min 1.5E+05 V=300 m/min <B 1.0E+05 CC c '5 V=150 m/min 55 5.0E+04 0.0E+00 V=300 m/min 0 0.25 0.5 0.75 Normalized Distance, x/Ls (a) 1 0.25 0.5 0.75 1 Normalized Distance, x/Ls (b) Figure 5.12: Distribution of (a) equivalent plastic strain rate, and (b) maximum shear stress on the shear plane. L is the length of shear plane 3 iii) Temperature The distribution of temperature within the chip, workpiece, and cutting tool is shown in figure 5.11(c). The temperature reaches its highest values on the rake face. The high temperature at the primary deformation zone is mainly attributed to severe plastic deformation. At the tip of the tool, temperature reaches a very high value, around 630 Â°C in the case shown. As the Chapter 5. NUMERICAL RESULTS AND DISCUSSION 132 chip moves higher on the rake face, the temperature is further increased due to the combined effects of friction and additional plastic straining in the secondary deformation zone, until it reaches its peak above the middle of the contact zone. It then decreases as the chip moves toward the point of separation. The location of the peak seems to be related to the point of transformation of sticking friction into sliding condition. The above trends appear to conform to the experimental observations [7]. A large temperature difference of the order of a few hundred degrees is observed across the thickness of the chip, which seems to be the direct result of high strain gradient. Figure 5.13(a) shows the profile of temperature on the shear plane for two different cutting speeds. It is seen that temperature drops very quickly on this plane. The strong dependence of temperature on the cutting speed is also evident in the figure. (a) (b) Figure 5.13: Temperature distribution on (a) the shear plane, (b) chip-tool interface. L and L are length of shear plane and contact length, respectively. s c This is because as the speed increases, adiabatic thermal conditions prevail, resulting in a larger temperature gradient across the chip. The temperature profile on the rake face is plotted in figure 5.13(b), which shows its peak on the rake face, depending on the cutting speed. In the contours offigure5.11(c), the distribution of temperature in the tool is also shown. There is a Chapter 5. NUMERICAL RESULTS AND DISCUSSION 133 large heat flow into the tool which, in combination with friction, raises the temperature along its rake face to the same level as that of the chip. iv) Effective stress Figure 5.11(d) shows the distribution of effective stress in the chip and the workpiece. Large stress values are found in the primary and secondary deformation zones. The highest stress values occur at the primary deformation zone, and the stresses on the rake face are substantially smaller. This is due to the fact that the stress is a function of combined effects of strain, strain rate and temperature, and hardening effects of strain and strain rate compete with softening effects of temperature rise in the deformation zone. Although both strain and strain rate are high on the rake face, the large temperature rise seems to dominate the other factors on the rake face, resulting in considerable softening of the material and thus, reduction of stress values. On the other hand, the temperature on the free surface of the chip is significantly lower, and the effects of relatively large strain and strain rate increase the effective stress of the material in that region. The distribution of maximum shear stress on the shear plane is shown in figure 5.12(b) for two cutting speeds. The typical value of average shear stress at the primary shear zone for the given material is about 650 MPa from simulations, compared to a value of around 700 MPa predicted in [131] based on shear plane theory. As the chip moves toward releasing from the tool rake face, it is gradually unloaded and the stress reduces until it vanishes at the point of separation. v) Contact stresses on the rake face In reference [131], the distribution of normal and frictional contact stresses on the rake face is measured by split tool tests for the case of cutting with 0Â° rake angle, feed rate of 0.2 mm and cutting speed of 100 m/min. For the purpose of comparison, similar cutting conditions Chapter 5. NUMERICAL RESULTS AND DISCUSSION 134 are simulated and the distribution of stresses on the rake face are reported in figure 5.14 along with experimental results. It is seen that the measured and predicted trends of both normal and frictional stresses on the rake face are quite similar. However, there is considerable difference between the values of stresses, especially the normal stress, on the rake face. The differences may be partially attributed to the approximations involved in material and frictional models. The predicted normal stresses reduce at a more monotonic pace than the measured values. The predicted normal and friction stress distributions are somewhat similar to the distribution suggested by Barrow, et al., [19] (see figure 2.2). Temperature distribution is also shown in figure 5.14. vi) Cutting and thrust forces The predicted cutting and thrust forces under various cutting conditions are compared with experimental results [ 131 ] in figures 5.15 to 5.18. In general, good agreement between simulated and measured forces is observed, considering the many simplifying assumptions which have to be made in numerical analysis (frictional conditions, material model) as well as in extracting orthogonal cutting forces from turning operation in the above experiments. The predicted cutting forces are in general larger than the experimental values, while the predicted thrust forces are lower. One contributing factor might be the size of edge radius assumed in the simulation, since this parameter is not specified in reference [131]. The magnitude of thrust force is also strongly affected by the coefficient of friction. 5.4.4 Influence of Cutting Conditions on Process Variables i) Effects of rake angle Three different rake angles are used to study the effect of rake angle on the chip formation process. Figure 5.15(a) depicts the trend of cutting and thrust forces with change in rake angle. Chapter 5. NUMERICAL RESULTS AND DISCUSSION 0 0.1 0.2 0.3 0.4 , . â€ž0.5 0.6 0.7 135 0.8 0.9 1 700 450 h 400 -I 1 1 0 0.2 0.4 1 0.6 ULc 1 1 0.8 1 1.2 Figure 5.14: distribution of normal and frictional stresses and temperature on the shear plane for [an = 0Â°, V = 100m/min, and / = 0.2mm] in comparison with experimental data [131] w It can be seen that the forces on the rake face are slightly decreased at higher rake angles, which is similar to the trend for experimental cutting forces. The decrease in forces at higher rake angle can be attributed to the fact that as the rake angle increases, the chip remains on the rake of the tool for a shorter length, and the curling of the chip occurs faster and with a smaller curl radius. Figure 5.16 compares the configuration and contact length of the chip after same cutting length for different rake angles. The contact lengths are 0.51, 0.40, and 0.36 mm for 0Â°, 6Â° and 12Â° rake angles, respectively. The effects of rake angle on shear angle is demonstrated in figure 5.15(b), which shows that the shear angle has a rising trend with the rising rake angle, similar Chapter 5. NUMERICAL RESULTS AND DISCUSSION 136 Figure 5.15: The effect of rake angle on cutting and thrust forces, shear angle and maximum temperature on the rake face Figure 5.16: Contact length and chip curling for 0Â°,6Â°, and 12Â° rake angles, from left to right, respectively. Temperature contours are also shown in the figure Chapter 5. NUMERICAL RESULTS AND DISCUSSION 137 to what is obtained from experiment. However, the predicted shear angles are higher than their experimental counterparts at low cutting speed, but they are closer in the case of higher cutting speed. Also shown in figure 5.15(b) is the predicted change in maximum temperature on the rake face versus the rake angle. ii) Effects of cutting speed Two different cutting speeds have been simulated for each combination of rake angle and feed rate. Figure 5.17(a) shows the effect of cutting speed on the cutting and thrust forces in each case, in comparison to the experimental forces. It is seen that the forces slightly 400 300 46 850 _ o 42 750 Measured Predicted 38 Measured Predicted 0 deg. rake angle 6 deg. rake angle 12 deg. rake angle 0 0 deg. rake angle A 6 deg. rake angle 12 deg. rake angle â€¢ 34 Q. 650 E 550 250 S. o o> c 200 450 3 0 < S 26 350 .c in 150 100 150 200 250 300 Cutting Speed (m/min) (a) 350 22 -1 100 150 200 250 300 250 350 Cutting Speed (m/min) (b) Figure 5.17: The effect of cutting speed on forces, shear angle and temperature decrease with increasing the cutting speed. Same trend is also observed in experimental results, though there are some exceptions. Figure 5.17(b) shows the predicted and experimental shear angles versus cutting speed. The experimental shear angles increase markedly when cutting Chapter 5. NUMERICAL RESULTS AND DISCUSSION 138 speed is increased from 150 m/min to 225 m/min, but there is less change for the speed of 300 m/min. In contrast, the difference between predicted shear angles at different speeds are smaller, although they have the same rising trend. The cutting speed has a profound effect on the maximum temperature in the chip, as shown in thefigure.As explained earlier, prevailing of adiabatic thermal conditions at higher speed is responsible for increase in local temperatures and larger gradient of temperature across the chip thickness. This, in turn, increases the softening of the material on the rake face, resulting in lower magnitudes of stress on the tool, which seems to be the reason for lower cutting forces at higher speed. Obviously, the most profound effect of cutting speed is seen in the magnitude of strain rate in the cutting zone, such that it is essentially proportional to the cutting speed. iii) Influence of feed rate Two different feed rates of 0.125 and 0.254 mm are simulated. As expected, increasing the feed rate results in higher forces, simply because there is more material to remove. Comparing the forces at different feed rates infigure5.18(a) shows that the forces are almost doubled with doubling the feed rate. The same trend is seen for the experimental forces. The predicted and experimental shear angles are plotted infigure5.18(b) as a function of feed rate. The shear angle clearly increases at higher feed rates in both simulation and experiment. However, the predicted shear angles are higher than the experimental ones. The predicted maximum temperature versus feed rate is also plotted in the abovefigure,which shows that the feed rate has little effect on the temperature in the chip. The shape of the chip is also affected by the feed rate, because the contact length increases with the increase in feed rate, and the curl radius becomes larger. Chapter 5. NUMERICAL RESULTS AND DISCUSSION 750 139 50 â€” Measured Predicted O â€¢ V= 150 m/min V= 300 m/min 900 a 800 â€” o 45 3 Is 600 700 Â£ -O E 3 Measured Predicted 450 O â€¢ 600 | V= 150 m/min V= 300 m/min 500 35 400 S 300 o g>30 < 150 0.1 0.15 0.2 0.25 Feed Rate (mm/rev) 0.3 25 0.1 300 200 0.15 0.2 0.25 0.3 Feed Rate (mm/rev) Figure 5.18: The effect of feed rate on forces, shear angle and temperature 5.5 ANALYSIS OF C U T T I N G WITH C H A M F E R E D AND W O R N E D G E TOOLS The productivity in high speed machining is hindered mainly by chatter vibrations and tool life. Once the chatter is eliminated, the material removal rate is constrained by the tool wear and thermo-mechanical stresses in the cutting edge, which are dependent on the tool geometry, tool and work material compositions, cutting speed and chip load [132]. Two common modes of tool failures in high speed machining are tool breakage and diffusion wear. In the former, the sharp edge of the tool is chipped due to high stresses, while in the latter, the tool edge is quickly worn out due to breaking of material bonds at high temperatures. One way of preventing tool breakage is to strengthen the tool edge by creating a chamfer or a Chapter 5. NUMERICAL RESULTS AND DISCUSSION 140 large edge radius at the cutting edge, or selecting speed and chip load that keeps the thermomechanical stresses within a safe level. Also, the diffusion wear of the tool might be slowed down if the cutting speed and chip load are selected in such a way that the tool-chip interface temperature remains under the critical diffusion limits of tool binding materials. In industry, either micro grain carbides or tool materials with higher resistance to diffusion such as Cubic Boron Nitride (CBN) are used in high speed machining. Both tools are used with chamfer edges in order to absorb high thermo-mechanical stresses in high speed machining of hardened steels. Researchers from different perspectives have investigated the effect of tool edge geometry on the chip removal process. Even the "nominally sharp" tools have afiniteedge radius that influences the cutting forces. The force contribution due tofiniteedge radius, i.e. the "ploughing force", is dependent on its size, i.e. it increases with gradual wear of the cutting tool or edge radius. It is generally believed that the edge or ploughing forces do not contribute directly to the chip removal process, but their presence affects surface integrity and residual stresses in the machined surface. Furthermore, it is known that the magnitude of such forces is dependent on the size of theflankwear of the tool. Since it is very difficult to measure ploughing forces directly, they are usually extracted from total cutting forces based on an analytical or empirical model. A common approach for identification of ploughing forces is based on extrapolation of forces measured at various uncut chip thickness h to zero thickness [133]. However, some experimental observations show that this approach may amount to considerable overestimation of ploughing forces [134]. Analytical models based on plasticity theory have also been presented for identification of these forces (see for example [135]). Elanayar and Shin [136] developed a procedure for separation of ploughing forces by accounting for the change in the state of indentation of the rigid tool into the work piece. Montgomery and Altintas [137] included the contribution of ploughing forces to cutting forces in dynamic milling. Waldorf, et al. [138] investigated the identification of ploughing forces by contrasting two models of Chapter 5. NUMERICAL RESULTS AND DISCUSSION 141 material flow when machining with a small cutting edge radius. One model assumes that there exists a separation point on the cutting edge at which the oncoming material is diverted to either form the chip or the machined surface, while the other model postulates the existence of a stable built-up edge adhering to the tool edge and cutting the workpiece at its apex. They developed force prediction models based on either assumptions and compared the predictions with experimental results for both small-edge and large-edge radius tools. The comparisons were shown to be in favor of the existence of built-up edge. When a chamfer is introduced to the tool geometry, the chamfered edge acts as the primary rake of the tool with a limited length and at a relatively large negative angle a 1 ; and the main rake of the tool becomes the secondary rake at a positive, neutral or slightly negative angle a , 0 as shown in figure 5.19. The chamfer enhances the performance of the tool by strengthening Figure 5.19: Schematic view of cutting with chamfered edge tool the edge and reducing the possibility of breakage. Experimental investigations have shown that the primary mechanism of cutting with chamfered edge tools is one of creation of a dead Chapter 5. NUMERICAL RESULTS AND DISCUSSION 142 metal zone made of hardened work material which is trapped under the chamfer and almost completely fills the chamfer. Such dead zone acts as the main cutting edge of the tool, and thus protects the tool surface from wearing under heavy cutting conditions. The presence of the dead metal zone is known to be independent of cutting conditions, which distinguishes it from the built up edge phenomenon that is observed on the edge of sharp tools at certain cutting conditions. The drawback of cutting with dead metal zone on the chamfered edge are that the forces on the tool are increased, and dimensional accuracy may be compromised due to cutting with the dead zone which its size may vary during cutting. There has been a number of experimental and analytical studies on the effects of chamfer on the performance of the tool. Hirao, et al. [139] observed from experiments that in cutting with different chamfer angles, the chip thickness is not affected by the chamfer angle because the dead metal zone effectively replaces the missing nose of the chamfered edge. They also observed that the chamfer angle strongly affects the thrust force on the tool, and to a lesser degree, the cutting force. Zhang, et al. [140] concluded from their quick stop tests and microphotography of the chip formation process that the existence of the dead metal zone is independent of the cutting speed, rake angle and chamfer angle, and that this zone acts as a rough extruding die, extruding the metal in the lower layer into the workpiece. They proposed a slip-linefieldaround the edge of the cutting tool and predicted the shear angle by minimizing the energy consumed in the primary and secondary deformation zones as well as in the dead zone. More recently, Ren and Altintas [15] proposed another slip-linefiledsolution and laid out a predictive model of cutting with chamfered tools based on minimum energy and identification of materialflowstress as a function of strain, strain rate and temperature. In spite of the increasing importance of chamfered tools, relatively little research has been done toward better understanding of the mechanism of chip formation in cutting with such tools, and on the effects of tool edge geometry on cutting variables. Numerical simulation of this process can provide a more in-depth analysis of this process. Chapter 5. NUMERICAL RESULTS AND DISCUSSION 5.5.1 143 Finite Element Model In this work, the ALEfiniteelement method is used for simulation of cutting of P20 mold steel with chamfered tools of various chamfer angles. The purpose of this study is to investigate the influence of the tool edge geometry on the chip formation process. It is noteworthy that such a cutting process in which the tool is not sharp could not be simulated by a Lagrangian analysis with node separation, because for such tools, there is no clear separation line along which nodes may be separated. In the ALE analysis, on the other hand, such geometric feature poses no difficulty, because the cutting action is modeled as continuousflowof material around the tool edge. Thefiniteelement model and mesh motion scheme for the simulations reported here are similar to those for sharp tools, except that the mesh around the tool is adjusted to include the chamfer. The FE mesh representing material under the chamfer is assigned an Eulerian type of motion. Afilletof very small radius is added to theflankof the tool to account for the finite radius of this edge, and the interaction of the material with theflankof the tool may be modeled. In all cases, 6-node triangular elements with 3-point integration were used. Around 750 elements and 1500 nodes were needed to model the workpiece and the tool, and a sufficiently fine mesh was used around the chamfered edge and on the shear zone to capture the high gradients of solution variables there. A typical simulation involved around 1000 incremental steps for a cutting length of 1 mm. 5.5.2 Experimental Work and Input Data A series of orthogonal cutting tests with chamfered carbide and CBN tools were performed in the UBC manufacturing automation lab, as reported by Ren [141]. Disks made of P20 mold steel were cut in plunge mode on a CNC turning center, and cutting forces, chip thickness and contact length were measured. In the first set of tests, blank carbide tools of ISO-S10 class Chapter 5. NUMERICAL RESULTS AND DISCUSSION case 1 2 3 a 144 b f (mm) c 0 0Â° 0Â° 0Â° -10Â° -25Â° -35Â° 0.0902 0.0841 0.0863 Table 5.4: Tool edge geometry in cutting tests with carbide tools were ground to different chamfer angles a i and lengths b f as listed in table 5.4. The disks c were turned at a chip load of 0.1 mm and two cutting speeds of 240 and 600 m/min. The main rake angle of the tool was 0Â° in all cases. In the second set of tests, Mitsubishi MB820 CBN inserts with main rake angles of â€”5Â° and chamfer angle of â€”25Â° and chamfer length of 0.1 mm were used to cut P20 disks at three different cutting speeds of 240,600, and 1000 m/min. The feed rate in the second set of tests was 0.06 mm. The materialflowstress model for P20 steel was provided by the ERC net shape manufacturing group in Ohio state university [142]. This model is obtained by fitting data obtained from cutting tests to a rate and temperature dependentflowstress model based on the procedure outlined by Oxley [7]; Y = (A + Be )(l + m (5.2) C\nk)(D-ET* ) n where A = 178.5 MPa, B = 462.4 MPa, C = 0.0438, D = 1.51, and E = 0.99 are the coefficient of this equation for this material, the strain sensitivity parameter m is equal to 0.169, the variable T* is defined as rp* where T m e i t i n g T- T - - room 1 Tmelting L T room = 1480 Â°C for P20, and temperature sensitivity parameter n = 0.666. The range of validity of this equation is given as 0.9 < e < 1.5, 2 x 10 < e < 8 x 10 4 5 and 600 < T < 1200 Â°C. Nevertheless, in the absence of other data, this equation has been used in the simulations here for strain and temperature values higher than the above ranges. The other mechanical and thermal properties of the workpiece and tool material is adopted from reference [128], as reported in table 5.5. The chamfer length in all simulations with Chapter 5. NUMERICAL Material properties Young's modulus Poisson's ratio Thermal conductivity Specific heat Thermal expansion coeff. RESULTS AND DISCUSSION Workpiece (AISI P20) 260 GPa 0.3 145 Tool (Carbide grade S10) (assumed rigid) 120 W/mÂ°C 343.3 J/Kg/Â°C 51.5W/mÂ°C 470 JyKg lÂ°C (1.3 - 1.4) x l O - ^ 6 0 1 5.2 x 1 0 - C 6o _1 Table 5.5: Material properties for P20 and carbide tool [128] carbide tool is assumed to be 0.09 mm. The shear type of friction is used in the sticking zone on the rake face, where the frictional stress is assumed to be equivalent to the local shear flow stress. In the sliding zone, the Coulomb friction law is used with an average friction coefficient obtained from the ratio of experimental thrust and cutting forces. 5.5.3 Influence of Chamfer Angle on the Cutting Process The effects of the chamfer angle on the chip formation process is analyzed by comparing the results of cutting simulations under same cutting conditions but with different chamfer angles. Three chamfer angles of 10Â°, 25Â° and 35Â° are modeled and compared with the cutting process with a "sharp" tool. It should be noted that by "sharp tool" we mean a tool that has an edge radius much smaller than the uncut chip thickness. In addition, cutting by a tool with "worn" edge is also modeled. The worn tool is considered as a tool that has an edge radius of the same size order of the uncut chip thickness, and represents the common situation when the originally sharp tool has been blunted due to wearing of the edge after some length of cutting. In the simulation, the edge radius of the worn tool is assumed to be 0.075 mm for an uncut chip thickness of 0.1 mm. i) Existence of the dead metal zone Of particular interest in the study of the chamfered tools is the existence of the dead zone under the tool chamfer. Earlier experimental studies report that the dead zone is formed under the Chapter 5. NUMERICAL RESULTS AND DISCUSSION 146 chamfer independent of the cutting conditions and chamfer angle, and that it almost completely fills the missing nose and acts as the actual cutting edge for the tool. While it is not possible to show the separation of part of the work material as the dead zone in the simulation without a fracture scheme, the existence of such zone may be confirmed by studying the velocity profile of the work material around the tool edge. Figure 5.20 shows the velocityfieldfor different chamfer angles as well as for the worn tool. It is seen that in all cases, part of the material is trapped under the chamfer, and is stagnated there. The size of the dead material zone in which the material velocity is close to zero, shown in dark blue in the contour plots, is dependent on the chamfer angle. This zone almostfillsthe missing nose, except for its lower portion on the flank radius, and its shape is consistent with quick stop tests of Hirao, et al. [139]. The dead zone is also formed on the nose of the worn tool in the same way as for the chamfered tool. It is seen in thefiguresthat the work material velocity increases form stagnation around the dead zone to the workpiece velocity under the toolflank,suggesting a large degree of straining at this region, which seems to be an extended tertiary deformation zone. It is also observed that the tip of the dead zone may not be at the same level of theflankof the tool which causes considerable friction between the toolflankand the machined surface, resulting in toolflankwear, and leads to dimensional inaccuracy. In fact, the surfacefinishseems to be greatly dependent on the shape and size of the dead zone. Jacobson, et al. [144] report that various types of dead metal zones on the tool may become occasionally unstable or change in form and size, and sometimes disintegrate fully or partially resulting in deposits on the cut surface [144]. Nevertheless, the dead metal zone on the chamfered edge seems to be much more stable than a built up edge on a sharp tool, though it is possible that its size might change during cutting. The influence of chamfer angle on cutting variables is demonstrated infigure5.21 where contours of effective stress and effective plastic strain rate contours are shown for three different chamfer angles. The shear zone is evident in thesefigureswhich, consistent with experimental observations, is of triangular shape emanating from the tool tip, and widening toward the free Chapter 5. NUMERICAL RESULTS AND DISCUSSION 147 Figure 5.20: Velocity field on the chamfered edge of tools with different chamfer angle shows the presence of a stagnated metal zone; a) -10Â°, b) -25Â°, c) -35Â°, and d) worn tool of radius 0.075 mm. All velocities are in (mm/s). Chapter 5. NUMERICAL RESULTS AND DISCUSSION 148 Figure 5.21: Distributions of effective stress (top) and plastic strain rate (bottom) for different chamfer angles Chapter 5. NUMERICAL RESULTS AND DISCUSSION 149 surface at the back of the chip. The effective stress, which exhibits the combined effects of strain, strain rate and temperature, increases as the material moves toward the center of primary shear zone due to strain and strain rate hardening until it reaches a maximum at the center of the shear zone. The stress has a fairly constant value across the center of the shear zone, but have lower values close to the tool tip where significant temperature rise softens the material considerably. The strain rate distribution has very large values at the flank of the tool where the work material is extruded under the trapped dead zone. It then drops quickly to much smaller values toward the free surface of the chip, and is almost constant for a large part of the length of the shear zone. High values of strain rate is also observed on the lower portion of the secondary rake face. It is interesting to note that for large chamfer angles, this variable has small values on the tool chamfer where the dead zone is stagnated on the tool chamfer. Instead, the straining occurs around the outer perimeter of the dead zone. The stress and strain rate distributions show little dependency on the chamfer angle. The size and shape of the shear zone and the stress and strain rate values are almost similar in all cases. Also the chip thickness, and thus the shear angle, show little change in the simulated cases. Such independence of the cutting process from the chamfer angle is supported by experimental observations. Researchers [139] have reported that the chips produced with different chamfer angles are almost identical. This is attributed to the presence of the dead metal zone which fills the chamfer and makes the cutting process almost equivalent for different chamfer angles. It is also observed experimentally that the shear angle shows little dependence on the chamfer angle. Zhang, et al. [140] have reported that the shear angle for a chamfered tool is decreased by a maximum 2Â° â€” 3Â° compared to that for a sharp tool. Chapter 5. NUMERICAL RESULTS AND DISCUSSION 150 ii) Cutting and trust forces, stresses, and temperature Although with and without chamfer, the chip is formed essentially in the same way, the total cutting forces on the tool are quite different in each case. The dependency of the tool forces on the chamfer angle is evident from the distribution and magnitude of stresses in the feed and thrust directions. Figures 5.22 and 5.23 compare the distributions of stresses a xx and a yy for sharp and chamfered tools of three different chamfer angles, respectively. It is seen that the distribution of stresses in the cutting direction <r is essentially similar for all the four xx cases shown. On the other hand, the stresses in the thrust direction a yy show marked increase with increase of the chamfer angle. It is noted in figure 5.23 that the area of the concentration of high a yy values turns downward toward the center of the workpiece as the chamfer angle increases. As noted by [140], the work material is extruded under the dead zone to form the machined surface. The stress distribution is consistent with what is observed in metal extrusion processes, where the change of die angle does not affect the stresses in the extrusion direction, but significantly affects the stresses in the perpendicular direction. The a yy stresses under the chamfer support the built up nose of the tool in this area. The larger stress values at larger chamfer angles mean that the material is strained to a larger extent, and the residual stresses on the machined surface will be larger too. Figure 5.24 compares the cutting and thrust forces obtained from experiment and FE simulation. Also shown in the figure are the analytical force predictions of Ren and Altintas [15]. It is seen that both cutting and thrust forces are increased with the increase of chamfer angle. However, whereas the force in the cutting direction is only slightly affected by the chamfer angle, the increase in the thrust forces is much more pronounced. This observation is consistent with the experimental results shown in the figure, and is also reported in other experimental studies [139,140]. The cutting force obtained from simulation is in average about 10% larger than experimental values, while the thrust forces are only slightly underestimated in Chapter 5. NUMERICAL RESULTS AND DISCUSSION 151 Figure 5.23: Distribution of a yy for sharp and chamfered tools of various chamfer angles Chapter 5. NUMERICAL RESULTS AND DISCUSSION 153 300 350 300 _-B- â€”a ..A 250 ...A'' X 200 200 O ^â€” 150 Experiment 150 â€” 0 ~ â€” Experiment - -B â€” Simulation â€” -aâ€” 100 Analytical X s~ -â€¢ *- Blunt tool X 100 Simulation - Analytical Blunt tool 50 10 20 30 Chamfer Angle (Deg.) 40 10 20 30 Chamfer Angle (Deg.) 40 Figure 5.24: Cutting and thrust forces versus chamfer angle. Numerical predictions are compared with experimental and analytical results of Ren, et. al [15] the simulations. Such differences may be partially attributed to the uncertainty in the material and frictional models used in the simulation, as well as the difference between the actual edge radius during cutting test, and the one used in the simulation. The forces predicted from slip-line solution show the same trend as that of experimental andfiniteelement results, though they are considerably overestimated. The forces predicted for the worn tool are also shown infigure5.24. It is seen that while the cutting force increases only slightly, the thrust force is further increased compared to that for the tool with â€”35Â° chamfer angle. It may be concluded that the lower portion of the edge radius of the worn tool acts as a tool with a very large chamfer angle, resulting in larger thrust forces. Simulation results show that the distribution of temperature in the shear zone and at the tool-chip interface is not affected significantly by the chamfer angle. The same conclusion is also obtained by Ren [15] from analytical solution based on the model of-Lowen and Shaw. A typical distribution of temperature is shown infigure5.25 for the case of â€”35Â° chamfer angle Chapter 5. NUMERICAL RESULTS AND DISCUSSION 154 which shows the maximum temperature of around 1000 C on the rake face of the tool. It is 0 noted that due to high conductivity of the tool material, the tool reaches a very high temperature after a relatively short length of cutting. Also, there exists a large gradient of the temperature across the width of the chip. Figure 5.25: temperature distribution in the chip and the tool for the case of â€”35Â° chamfer angle Finally,figure5.26 shows the distribution of stresses, effective plastic strain rate and temperature for the worn tool. Compared to the stress distributions on chamfered edges, the maximum stress in the cutting direction for the worn tool shows slight increase, while the tangential stress follows the same increasing trend as for chamfered tools. An interesting observation in the contour of strain rate is that its maximum is considerably smaller than those observed for chamfered tools (seefigure5.21). This may be because in this case there is no Chapter 5. NUMERICAL RESULTS AND DISCUSSION 155 Chapter 5. NUMERICAL RESULTS AND DISCUSSION 156 sharp edge and the material is strained over a larger area by the worn edge at theflankof the tool. The temperature distribution for the worn edge shows that while its peak occurs on the chip-tool interface close to the point of separation, theflankof the tool is also at a very high temperature due to excessive straining and friction of the material with toolflankat this region. 5.5.4 Influence of Cutting Speed on Cutting with Chamfered Tools The 510 carbide inserts used in the experimental work reported above are worn out quickly if used for machining at speeds above 400 m/min. In order to investigate the influence of cutting speed on the cutting process with chamfered tools, CBN inserts are used with a constant chamfer angle of â€”25Â°, a chamfer length of 0.1mm and a main rake angle of â€”5Â°. Three different cutting speeds of 240, 600 and 1000 m/min are tested. Similar cutting conditions are simulated using the ALEfiniteelement program. Figure 5.27 shows the distributions of effective plastic strain rate and temperature obtained from cutting simulations under the three different cutting speeds. From the simulations, the shear angle is predicted at around 27Â°, and shows only slight dependency on the cutting speed, whereas the experimental results show a shear angle of 28Â° for V = 240 m/min , and 32Â° w for the V = 600 and 1000 m/min. The distribution of strain rate in the cutting zones is not w influenced significantly for different cutting speed, but its magnitude is immensely affected by the cutting speed. The maximum strain rate reaches close to order of 10 s' 6 1 for the cutting speed of 1000 m/min. The magnitude and distribution of temperature is also affected by the cutting speed. The maximum temperature increases form around 1000Â° C to 1500Â° C when the cutting speed increases from 240 to 1000 m/min. A similar trend is reported by Ren [15] for the chip-tool interface temperature, where a range of temperatures of 1150Â°C to 1700Â°C is predicted for the above range of speeds. On the other hand, the distribution of temperaturefieldis also affected Chapter 5. NUMERICAL RESULTS AND DISCUSSION 157 Figure 5.27: Distributions of effective plastic strain rate (left), and temperature (right) for the CBN tools at cutting speeds of V = 240, 600, and 1000 m/min, from top to bottom, respectively w Chapter 5. NUMERICAL RESULTS AND DISCUSSION 250 158 250 â€” â€¢ â€” Experiment 225 225 & ffl : 200 '175 150 m Experiment â€” -B- Simulation ? E 2 200 o o o u. to 175 â€” -B â€” Simulation - Analytical 150 Analytical 125 200 400 600 800 Cutting Speed (m/min) 1000 125 200 400 600 800 Cutting Speed (m/min) 1000 Figure 5.28: effect of cutting speed on the cutting and thrust forces. Numerical predictions are compared with experimental and analytical results of Ren, et. al [15] by speed as the region of high temperature on the chip-tool interface becomes narrower and the temperature gradient across the chip width becomes larger at higher cutting speeds. This phenomenon may be attributed to the fact that at higher speeds, there is less time for spread of thermal energy in the chip through conduction, and the thermal conditions move closer to adiabatic conditions. It is also noted in the contours of temperature in the tool that less thermal energy is conducted inside the tool compared to the case of cutting with carbide tool, which is due to lower thermal conductivity of the CBN. The areas of high temperatures are concentrated on the chip-tool interface and the flank of the tool. Figure 5.28 shows the effect of cutting speed on the cutting and thrust forces. The cutting force shows a decreasing trend with temperature, which may be attributed to the more softening of the material at higher temperatures. The difference between the experimental and simulated cutting forces is around 15%. Shown in the figure is also the forces predicted by the analytical slip-line model of Ren [15]. The thrust forces obtained from simulation are also decreased at higher speeds, and are in good agreement with the experimental results. The analytically Chapter 5. NUMERICAL RESULTS AND DISCUSSION 159 predicted thrust forces show larger differences with the experimental values. From the discussion in this chapter on the cutting with chamfered tools, it may be concluded that introducing chamfer to the tool edge enhances the performance of the tool by strengthening it against breakage. While chip formation process is not affected significantly by the presence of the chamfer, the cutting forces are increased, and larger friction of the work material with theflankof the tool may result in rapid wear and dimensional inaccuracy. The above study supports the experimental observations on the role of cutting edge in the chip formation process. Nevertheless, only a limited range of material, cutting conditions and edge geometries has been considered here. Further study of this subject may include varying chamfer length, ratio of chamfer size to uncut chip thickness, a broader range of cutting conditions, and experimental verification of residual stress prediction. 5.6 ORTHOGONAL C U T T I N G OF A TITANIUM A L L O Y In this section, as an example of cutting non-ferrous materials, results of simulations of cutting TiGAHV titanium alloy is presented. Titanium alloys are extensively used in the aerospace industry due to their mechanical properties and light weight. Due to poor thermal conduction of titanium alloys, machining of these materials are normally done at lower range of cutting speeds. High speed cutting of these alloys usually results in producing segmented chips due to adiabatic shear banding. Here, a range of cutting speeds in which continuous chip is produced is considered. Material properties Workpiece (Ti6Al4V) Shear modulus Bulk modulus Thermal conductivity Specific heat Density 41.7 GPa 101.9 GPa 6.7W/mÂ°K 542 J/Kg/Â°K 4430 Kg/m 3 Table 5.6: Material properties for titanium alloy [31] Chapter 5. NUMERICAL RESULTS AND DISCUSSION 160 The material flow stress for this material is given in the same form as equation (5.2), where material parameters are given as [143]: A = 782.7MPa, B = 498.4MPa, C = 0.028, D = 1.0, E = 1.0, m = 0.28, n = 1.0, and T melting = 1676Â°<7. Other material properties are adopted from reference [31] and are shown in table 5.6. The tool is of Carbide type and its thermal properties is given in table 5.5. The uncut chip thickness in the simulations are 0.06mm and 0.1mm, and the tool rake angle is 2Â°. Two cutting speeds of 30 m/min and 60 m/min are simulated. The above cutting conditions correspond to cutting tests by Lee and Altintas [145] who turned Ti6AlAV tubes of 3.81mm thickness in orthogonal mode for a wide range of cutting conditions. In these tests, cutting forces were measured using a force dynamometer and the chip thickness was measured by weighing a certain length of the chip. Figure 5.29(a) presents a typical distribution of stresses on the rake face of the tool. The distribution of both normal and frictional stresses on are similar the model of Barrow, as shown infigure2.2. The distribution of temperature on the rake face for two different cutting speeds are shown infigure5.29(b) where it seen that the maximum temperature has significantly risen when the cutting speed is increased from 30 m/min to 60 m/min. This large increase is attributed to the poor thermal conduction in titanium. Finally, cutting forces and chip thicknesses obtained from simulation are compared with their experimental counterparts in table 5.7. It is observed that the predicted cutting forces are reasonably close to the experimental forces. The predicted forces in cutting direction are generally higher than the experimental values, whereas the predicted forces in thrust direction are lower than their measured ones. On the other hand, the chip thicknesses predicted by the ALE simulation are in general larger than the measured values, by up to 20%. Chapter 5. NUMERICAL RESULTS AND DISCUSSION 161 -â€¢â€”Normal Stress - a â€” Frictional Stress 0.04 0.06 0.08 Distance from tip (mm) 0.12 (a) 1000 -â€¢â€” Vc= 60 m/min -aâ€”Vc=30 m/min _ 900 n 800 700 - â€ž . ^ 0.02 R . â€ž B ~ - a = f c f c 7 i f c ^ ^ 0.04 ... 0.06 â€ž 0.08 distance from tip (mm) (b) Figure 5.29: (a) Distribution of normal and frictional stresses on the rake face [a = 2Â°, V = 60m/min, h = 0.06mm], (b), temperature on the rake face for two different cutting speeds [a = 2Â°,h = 0.06mm] w Cutting Conditions h [mm) V (m/min) w 0.06 0.06 0.1 30 60 60 F (N/mm) Sim. Exp. c 140 137 220 119 120 188 F {N/mm) Sim. Exp. t 60 65 93 85 88 102 h c (mm) Sim. Exp. 0.91 0.89 0.13 0.76 0.73 0.11 Table 5.7: Forces in the cutting and thrust directions and chip thicknesses from simulation and experiment [145] for cutting titanium alloy Ti6AlAV. Chapter 6 CONCLUSIONS A N D F U T U R E WORK 6.1 DISCUSSION A N D CONCLUSIONS This thesis deals with the application of ALEfiniteelement approach to simulation of chip formation in orthogonal metal cutting process. Modeling this process is often a numerical challenge due to the very complex set of conditions present in a tiny shear zone and along the chip-tool interface, where material is severely deformed at high strain rate and temperature. A critical review of the literature on numerical simulation of this process shows that the conventional Lagrangian and Eulerian formulations used are in general inefficient for this purpose. Due to very large deformation of material in the cutting zone, problems such as mesh distortion and loadfluctuationoften arise in applying Lagrangian formulation to such analysis. Furthermore, in a Lagrangian simulation, the cutting action of the tool is often simulated by way of separating the nodes in front of the tool along a pre-defined parting line. It is argued here that such approach is inherently inadequate for cutting process, because it essentially treats the plasticflowof the material around the tool tip as a crack propagation problem, since a crack will always run ahead of the tool. Other problems pertaining to this approach include the need for a separation criterion that is often chosen arbitrarily, creation of unbalanced forces, and inability to model cutting with negative rake angles. Furthermore, the application of this approach for the cutting process is limited to tools with sharp edge, as there is no clear parting line in cutting with chamfered or blunt edge tools. Consequently, a more realistic approach that models the 162 Chapter 6. CONCLUSIONS AND FUTURE WORK 163 cutting action as one of continuous materialflowaround the tool tip has gained popularity in recent years. Within a Lagrangian framework, this approach would require frequent remeshing and transfer of variables between consecutive meshes, which is computationally expensive, and prone to convergence problems due to accumulation of errors. The Eulerian approach, on the other hand, can efficiently model theflowof material around the tool tip, because the Eulerian mesh is spatiallyfixedand no distortion occurs. However, in this approach, the boundaries of material domain should be known a priori. Therefore, the natural formation of chip as a result of unconstrainedflowon free boundaries may not be modeled directly, and iterative updating of free boundaries outside solution domain is required. In this work, it is shown that the ALE method has the potential of combining the strengths of both Eulerian and Lagrangian approaches in a single analysis. In an ALE analysis, the user has control over the motion of each degree of freedom of the mesh. Thus, if, for example, the mesh around the tool tip isfixedspatially, while the mesh on the free boundaries of the chip is allowed to follow material deformation, both features of the cutting process can be modeled efficiently without resorting to node separation, frequent remeshing, or iterative updating of free surfaces. Due to the very large deformation and high-pressure friction in a cutting operation, strain rate and temperature play a significant role in the deformation process. An ALE large deformation FE formulation originally developed by Wang and Gadala [72] has been extensively modified to include rate and thermal effects. A heat transfer module is designed which updates the temperaturefieldin the tool and the workpiece, following the deformation in each incremental step of the analysis. Thermal energy in the cutting zone is released by plastic deformation and friction at the chip-tool interface. To model the interaction between the tool and the chip, two algorithms are designed which can handle contact betweenflexiblebodies and rigid-flexible bodies, respectively. In thefirstalgorithm, the contact conditions between the bodies are established through enforcing compatibility of the displacements and equilibrium of forces on Chapter 6. CONCLUSIONS AND FUTURE WORK 164 the contact interface. In the rigid-flexible algorithm, the contact conditions are directly applied as boundary conditions on theflexiblebody. These algorithms are capable of detecting and compensating for penetration, and may be used to properly apply frictional forces, and to detect contact release. Theflexibilityof the ALE mesh can be used on the contact surface to model contact conditions more accurately. The success of an ALE simulation is dependent on the mesh motion scheme used. In this work, an efficient mesh motion algorithm is introduced which facilitates the assignment of mesh velocities and their reduction to Lagrangian or Eulerian degrees of freedom, if needed. An overview of the logical steps in mesh motion assignment is provided and general guidelines for designing a mesh motion strategy are presented. A scheme for mesh sliding on free boundaries is introduced which uses a single parameter to move a node tangent to the boundary, in accordance to ALE mapping requirement. Transfinite and Isoparametric mesh generation techniques are adopted for moving the mesh in the interior of the body, so that the mesh remains optimal throughout the analysis. Following the deformation in each incremental step, the mesh is updated to its targeted configuration, and material associated variables such as stresses and strains are transferred from material mesh to computational mesh, using the ALE rate equation. Combined with local least square smoothing, the updating scheme provides an efficient alternative to commonly used interpolation and extrapolation algorithms in Lagrangian and Eulerian formulations. It is essential here to point out some of the challenging issues in using this approach that should be carefully addressed. The presence of convective terms in the ALE equations is one of such issues that require careful treatment. Another challenging issue is the procedure for mapping material associated variables to grid points. The mapping techniques used for this purpose usually yield reasonably accurate results in areas of low stress gradients, but their behavior may deteriorate in the areas where the gradients of stress are high, such as in the primary deformation zone in orthogonal cutting process. In such cases, because of Chapter 6. CONCLUSIONS AND FUTURE WORK 165 large differences between stresses at neighboring Gauss points, interpolation may result in introduction of errors that accumulate rapidly. This difficulty may not be easily by refining the mesh in these areas, because more interpolations are required, and the need for more intensive computation partially offsets the purpose of using ALE approach. Also, The response of the mapping technique seems to be dependent on the type of element used. In our experience, for some element types, the element response became gradually stiff, resulting in erroneous chip shapes. The application of ALE finite element method to simulate the orthogonal metal cutting process is made possible by developing a large deformation rate dependent thermo-mechanical code along with designing an effective mesh motion strategy. The finite element simulation provides a thorough insight to this process, calculates detailed distributions of process variables in the cutting zones, and makes it possible to study the effects of cutting conditions and geometrical features of the tool on the chip formation process, without having to resort to costly experimentation. The results of finite element simulation of cutting low carbon steel with a carbide tool are presented as a benchmark problem shared among the leading researchers in this area, and are compared with experimental data. The comparison shows a good agreement between predicted and experimental results, and verifies the potential of ALE approach in efficient and accurate modeling of the cutting process. Also, a parametric study is conducted to investigate the influence of cutting conditions such as rake angle, cutting speed and feed rate on the chip formation process. The results of this study are compared with experimental results under similar conditions, and show fairly good agreement in most cases, confirming that the simulations provide a realistic representation of the actual process. Finally, the influence of cutting edge geometry on the chip formation process is investigated through simulation of cutting with chamfered or worn tools of carbide or CBN types. Chamfered CBN tools are frequently used in high speed machining of hard materials, because the chamfer Chapter 6. CONCLUSIONS AND FUTURE WORK 166 strengthens the edge against chipping, and the CBN material has high resistance to diffusion wear at very high temperatures present in high speed machining. The results of simulation of cutting process with different chamfer angles show that chamfer angle does not affect the chip significantly, because in all cases, a dead metal zone is formed under the chamfer and acts as the main cutting edge of the tool. A zone of stagnant material under the chamfer is evident in simulation results, with a size proportional to the size and angle of the chamfer. The study of the effects of chamfer angle on the stresses in the deformation zones and forces on the tool shows that both cutting and thrust forces increase with the increase in the chamfer angle, but the increase in the thrust force is more significant. It is observed that the region of high stresses in the thrust direction turns inward as the chamfer angle increases. Similar trend is also seen for a worn tool with a nose radius of the same order as the uncut chip thickness. The predictions on the effects of chamfer angle are in agreement with experimental observations. These results are also compared with those predicted by the slip line filed model of Ren and Altitnas [15] which show good agreement. The influence of cutting speed on the deformation process for chamfered tools is studied using CBN tools, because such tools are used in high speed machining of hardened steel alloys due to its higher resistance againstflankwear caused by diffusion. The most significant effect of the cutting speed is in the large increase in the temperature on the tool-chip interface. The simulation results show that the maximum temperature on the rake face increases from around 1050Â°C to 1500Â°C when speed is increased from 240 to 1000m/min. Considering the fact that the diffusion limit for binding of tool materials is 1300Â° C for carbide tools and 1600Â°C for CBN tools used in the experiments, the importance of using chamfered-edge CBN tool in high speed machining becomes evident. Chapter 6. CONCLUSIONS AND FUTURE WORK 6.2 167 F U T U R E WORK The application of ALE method for simulation of manufacturing processes, and in particular, metal cutting process shows the promising modeling potentials of this solution approach which are otherwise difficult to achieve using conventional methods. From theoretical point of view, the ALE formulation presented in this work may be complemented in several areas; â€¢ The algorithm for mapping of material associated variables is a particular area that needs further investigation, especially for regions where high gradients of solution variables exist. â€¢ The mesh motion scheme is another area which may be further explored. Designing an efficient mesh motion strategy is a critical and challenging aspect of an ALE analysis, and sometimes trial runs for models of complicated geometry are required. Further automation of this task is essential for making the ALE approach available to inexperienced users. â€¢ The extension of the present formulation to include dynamic effects is another avenue for research, which will make it possible to model manufacturing processes for which the dynamic effects are of primary importance. From the viewpoint of application in cutting simulations, the ALE approach may be extended in several ways, among them; â€¢ Simulation of cracked and segmented chips. While it is not difficult to implement a module for creation of cracks in the material, the challenge is in designing a mesh motion scheme that can adapt to creation of new boundary surfaces and keep the mesh regular at the same time. Chapter 6. CONCLUSIONS AND FUTURE WORK 168 â€¢ Study of tool life, tool wear, chip curling and chip breakage. â€¢ Extension of the application to 3-D cases. Such extension is straightforward from the point of view of the deformation model. However, contact algorithm and mesh motion schemes should be generalized to 3-D geometry. With this capability, oblique cutting process as well as more practical processes such as turning can be modeled. â€¢ Extension to intermittent cutting in processes such as milling. Again, the biggest challenge in this regard seems to be in designing the mesh motion scheme for such simulations. Combined with dynamic effects of the machine and the tool, the cutting forces predicted fromfiniteelement simulation may become part of a fully integrated computer aided manufacturing process. Bibliography [1] Piispanen, V., 1948, "Theory of Formation of Metal Chips", Journal of Applied Physics, 27, pp.877-881. [2] Ernst, H. and Merchant, M.E., 1941, "Chip Formation, Friction, and High Quality Machined Surfaces", Trans. ASME, 29 , pp.299- 378. [3] Merchant, M . E . , 1945, "Mechanics of the Metal Cutting Process, II. Plasticity Conditions in Orthogonal Cutting ", Journal of Applied Physics, 16, pp318-324. [4] Lee, E. H., and Shaffer, B. W., 1951, "The Theory of Plasticity Applied to Problem of Machining" Journal of Applied Mechanics, 13, pp. 405-413. [5] Shaw, M . C , Cook, N. H., and Finnie, I., 1953, " Shear Angle Relationship in Metal Cutting ", Trans. ASME, 75, pp. 273-288. [6] Johnson, W., and Mellor, P. B.,1973, Engineering Plasticity, Van Norstad Reinhold Company Ltd., London, UK. [7] Oxley, P. L. B., 1989, Mechanics of Machining, An Analytical Approach to Assessing Machinability, Ellis Horwood Ltd. [8] Armarego, E . J. A., and Brown, R.H., 1969, The Machining of Metal, PrenticeHall Inc.. [9] Palmer, W. B. and Oxley, P. L. B., 1959, "Mechanics of Metal Cutting", Proc. Institution of Mechanical Engineers, 173, pp. 623-654. [10] Oxley, P. L. B., and Welsh, M . J. M . , 1963, "Calculating the Shear Angle in Orthogonal Metal Cutting From Fundamental Stress, Strain Rate Properties of the Work Material", Proc. 4th Int. Machine Tool Research Conf., Pergamon Press, Oxford, pp. 73-86. [11] Stevenson, M. G., and Oxley, P. L. B., 1970-1971, " An Experimental Investigation of the Influence of Strain-Rate and Temperature on the Flow Stress Properties of a Low Carbon Steel Using a Machining test ", Proc. Institution of Mechanical Engineers, 185, pp. 741-754. [12] Usui, E . , Hirota, A., and Masuko, M . , 1978, "Analytical Prediction of Three Dimensional Cutting Process, Part 1: Basic Cutting Model and Energy Approach", Journal of Engineering for Industry, 100, pp. 222-228. 169 Bibliography 170 [13] Shamoto, E . , and Altintas, Y . , 1997, "Prediction of Shear Angle in Oblique Cutting with Maximum Shear Stress and Minimum Energy Principle", Journal of Manufacturing Science and Technology, 121, (1999) 399-407. [14] Seethaler, R. J . and Yellowley, I., 1997, "An Upper Bound Cutting Model for Oblique Cutting Tools with a Nose Radius", Int. Journal of Machine Tools and Manufacture, 37, pp. 119. [15] Ren, H . , and Altintas, Y . , 1999, "Mechanics of Machining with Chamfered Tools", Journal of Manufacturing Science and Engineering, to appear. [16] Jawahir, I. S. and van Luttervelt, 1993, "Recent Developments in Chip Control Research and Applications", Annals of the GIRP, 42, pp. 659-693. [17] Zorev, N . N . , 1963, " Interrelationship Between Shear Processes Occurring Along Tool Face and on Shear Plane in Metal Cutting ", Proc. International Research in Production Engineering Conf, ASME, 85, pp. 42-49. [18] Usui, E . and Takeyama, H . , 1960, "A Photoelastic Analysis of Machining Stresses", Journal of Engineering for Industry, 82, pp. 303-308. [19] Barrow, G . , Graham, W . , Kurimoto, T . , and Leong, Y . F . , 1982, "Determination of Rake Face Stress Distribution in Orthogonal Cutting ", Int. Journal of Machine Tool Design and Research, 22, pp. 75-85. [20] Stephenson, D . A . , 1991, "Assessment of Steady-State Metal Cutting Temperature Models Based on Simultaneous Infraxed and Thermocouple Data ", Journal of Engineering for Industry, 113, pp. 121-128. [21] Klamecki, B. A . , 1973, "Incipient Chip Formation in Metal Cutting-A three Dimension Finite Element Analysis", Ph.D dissertation, Univ. of Illinois. [22] Tay, A . 0., Stevenson, M . G . , and de Vahl Davis, G . , 1974, "Using the Finite Element Method to Determine Temperature Distributions in Orthogonal Machining" Proc. Institution of Mechanical Engineers, 188, pp. 627-638. [23] Tay, A . 0., Stevenson, M . G . , de Vahl Davis, G . , and oxley, P. L . B . , 1976, "A Numerical Method for Calculating Temperature Distributions in Machining From Force and Shear Angle Measurements", Int. Journal of Machine Tool Design and Research, 16, pp. 335-349. [24] Usui, E . , and Shirakashi, T . , 1982, "Mechanics of Machining- From Descriptive to Predictive Theory", On the art of Cutting Metals, 75 Years Later, A S M E - P E D 7, pp. 13-35. Bibliography 171 [25] Iwata, K . , Osakada, A . , and Terasaka, Y . , 1984, "Process Modeling of Orthogonal Cutting by Rigid-Plastic Finite Element Method", Journal of Engineering Materials and Technology, 106, pp. 132-138. [26] Strenkowski, J . S., and Carroll III, J . T . , 1985, "A finite Element Model of Orthogonal Metal Cutting", Journal of Engineering for Industry, 107, pp. 349-354. [27] Carroll III, J . T . , and Strenkowski, J.S., 1988 " Finite Element Models of Orthogonal Cutting with Application to Single Point Diamond Turning" Int. Journal of Mechanical Sciences, 30, pp. 899-920. [28] Strenkowski, J.S., and Moon, K . J . , 1990, "Finite Element Prediction of Chip Geometry and Tool/Workpiece Temperature Distribution in Orthogonal Metal Cutting", Journal of Engineering for Industry, 112 , pp. 313-318. [29] Ueda, K . , and Manabe, K . , 1993, "Elastic-Plastic F E M Simulation of Chip Formation Mechanism Based on Gurson's Yield Function", Trans. Japan Society of Mechanical Engineers, 59, pp. 1274-1279. [30] Hashemi, J . , Tseng, A . A . , and Chou, P. C , 1994, "Finite Element Modeling of Segmented Chip Formation in High-Speed Orthogonal Machining", Journal of Materials Engineering and Performance, 3, pp. 712-721. [31] Obikawa, T . and Usui, E . , 1996, "Computational Machining of Titanium Alloy, Finite Element Modeling and a Few Results", Journal of Manufacturing Science and Engineering, 118, pp. 208-215. [32] Marusich, T . D. and Ortiz, M . , 1995, "Modeling and Simulation of High-Speed Machining", Int. Journal of Numerical Methods in Engineering, 38, pp. 3675-3694. [33] Xie, J . Q. , Bayoumi, A . E . and Zbib, H . M . , 1996, "A Study on Shear Banding in Chip Formation of Orthogonal Machining", Int. Journal of Machine Tools and manufacture, 36, pp. 835-847. [34] Komvopoulos, K . and Erpenbeck, S. A . , 1991, "Finite Element Modeling of Orthogonal Cutting", Journal of Engineering for Industry, 113, pp. 253-267. [35] Ueda, K . , and Manabe, K . , 1993, "Rigid-Plastic F E M Analysis of Three Dimensional Deformation Field in Chip Formation Process", Annals of CIRP, 42, pp. 35-38. [36] Hashimura, M . , Ueda, K . , Dornfeld, D., and Manabe, K . , 1995, "Analysis of Three Dimensional Burr Formation in Oblique Cutting", Annals of CIRP, 44, pp. 27-30. Bibliography 172 [37] Hibbitt, H.D., Marcal, P. V . , and Rice, J . R., 1970, "A Finite Element Formulation for Problems of Large Strain and Large Displacement", Int. Journal of Solids and Structures, 6, pp. 1069-1086. [38] McMeeking, R. M . , and Rice, J . R. ,1975, "Finite Element Formulation for Problems of Large Elastic-Plastic Deformation", Int. Journal of Solids and Structures, 11, pp. 601-616. [39] Wang, J . , 1998, Arbitrary Lagrangian Eulerian Method and Its Applications in Solid Mechanics, PhD Dissertation, U B C . [40] Zhang, B . and Bagchi, A . , 1994, "Finite Element Simulation of Chip Formation and Comparison with Machining Experiment", Journal of Engineering for Industry, 116, pp. 289- 297. [41] Zhang, B . , and Bagchi, A . , 1994, "A Study of Chip Formation and Its Approximation in Finite Element Simulation of Continuous Chip Formation", ASME Materials Issues in Machining-II and The Physics of Machining Processes-II Stephenson, D. A . , and Stevenson, R. (Eds.), pp. 157-174. [42] Shih, A . J . , "Finite Element Simulation of Orthogonal Metal Cutting", 1995, Journal of Engineering for Industry, 117, pp. 84-93. [43] Strenkowski, J . S., and Mitchum, G . L., 19, "An Improved Finite Element Model of Orthogonal Metal Cutting", Proc. ASME Winter Annual Meeting, , pp. 506-509. [44] Lin, Z. C . and Lin, S. Y . , 1992, "A Coupled Finite Element Model of ThermoElastic-Plastic Large Deformation for Orthogonal Cutting", Journal of Engineering Materials and Technology, 114, pp. 218-226. [45] Chen, A . G . and Black, J . T . , 1994, " F E M Modeling in Metal Cutting", Manufacturing Review, 7, pp. 120-133. [46] Huang, J . M . and Black, J . T . , 1996, "An evaluation of Chip Separation Criteria for the F E M Simulation of Machining", Journal of Manufacturing Science and Engineering, 118, pp. 545-554. [47] Ceretti, J . E . , Fallbehmer, P., Wu, W. T . and Altan, T . J . , 1996, "Application of 2D F E M on Chip Formation in Orthogonal Cutting", Journal of Material Processing Technology, 59, pp. 169-181. [48] Madhavan, V . , Chandrasekar, S. and Farris, T . N . , 1993, "Mechanistic Model of Machining as an Indentation Process", Materials Issues in Machining III and The Physics of Machining Process, D. A . Stephenson and R. Stevenson (Eds.), The Mineral, Metals and Materials Society, pp. 187-209. Bibliography 173 [49] Sekhon, G . S. and Chenot, J . L . , 1993, "Numerical Simulation of Continuous Chip Formation During Non-Steady Orthogonal Cutting", Engineering Computations, 10, pp. 31-48. [50] Ng, E . G . , Aspinwall, D. K . , Brazil, D., and Monaghan, J . , 1999, "Modelling of Temperature and Forces When Orthogonally Machining Hardened Steel", Int. Journal of Machine Tools and Manufacture, 39, pp. 885-903. [51] Voyiadjis, C . Z., and Foroozesh, M . , 1991, "A Finite Strain Total Lagrangian Finite Element Solution for Metal Extrusion Problems", Computer Methods in Applied Mechanics and Engineering, 86, pp. 337-370. [52] DEFORM-2D user guide, 1994, Columbus, Ohio. Scientific Forming Technologies Corporation, [53] Zienkiewicz, O. C , 1977, The Finite Element Method, Third edition, McGraw-Hill. [54] Wu, J.-S., Dillon, J.R., and Lu, W . - Y . , 1996, "Thermo-Viscoplastic Modeling of Machining Process Using a Mixed Finite Element Method", Journal of Manufacturing Science and Engineering., 118, pp. 470-482. [55] K i m , K . W . and Sin, H - C , 1996, "Development of a Thermo-Viscoplastic Cutting Model Using Finite Element Method", Int. Journal of Machine Tools and Manufacture, 36, pp. 379- 397. [56] Abo-Elkhier, M . , Oravas, G . E . And Dokainish, M . A . , 1988, "A Consistent Eulerian Formulation For Large Deformation Analysis with Reference to MetalExtrusion Process", Int. Journal of Non-Linear Mechanics, 52, pp. 37-52. [57] Zienkiewicz, 0. C. and Godbole, P. N . , 1974, "Flow of Plastic and Visco-Plastic Solids with Special Reference to Extrusion and Forming Processes", Int. Journal of Numerical Methods in Engineering, 8, pp. 3-16. [58] Noh, W . F . , 1964, "A Time-Dependent Two-Space-Dimensional Coupled EulerianLagrangian Code", Methods in Computational Physics, Vol. 3, Alderb, B . et al. (Eds), Academic Press, New York, U S A . [59] Belytschko, T . , and Kennedy, J . M . , 1978, "Computer Methods for Subassembly Simulation", Nuclear Engineering and Design, 47, pp. 17-38. [60] Donea, J . , Fasoli-stella, D. P., Giuliani, S., 1977, " Lagrangian and Eulerian Finite Element Techniques for Transient Fluid-Structure Interaction Problems ", Transactions of SMiRT-4, Paper B l / 2 , San Francisco, U S A , pp. 1-12. Bibliography 174 [61] Hughes, T. J. R., Liu, W. K., and Zimmerman, T. K., 1981, "Lagrangian-Eulerian Finite Element Formulation for Incompressible Viscous Flows", Computer Methods in Applied Mechanics and Engineering, 29, pp. 329-349. [62] Haber, R. B., 1984, "A Mixed Eulerian-Lagrangian Displacement Model for Large Deformation Analysis in Solid Mechanics", Computer Methods in Applied Mechanics and Engineering, 43, pp. 277-292. [63] Haber, R. B., and Koh, H. M . , 1985, "Explicit Expressions for Energy Release Rates Using Virtual Crack Extensions,", Int. Journal of Numerical Methods in Engineering, 21, pp. 301-315. [64] Liu, W. K., Chang, H., Chen, J. S., and Belystchko, T., 1988, "Arbitrary Lagrangian-Eulerian Petrov-Galerkin Finite Elements for Nonlinear Continua", Computer Methods in Applied Mechanics and Engineering , 68, pp. 259-310. [65] Huetink, J., Vreede, P. T., Van der Lugt, J., 1990, "Progress in Mixed EulerianLagrangian Finite Element Simulation of Forming Processes", Int. Journal of Numerical Methods in Engineering, 30, pp. 1441-1457. [66] Liu, W. K., Chen, J. S., Belystchko, T., and Zhang, Y. F., 1991, "Adaptive ALE Finite Elements with Particular Reference to External Work Rate on Frictional Interface", Computer Methods in Applied Mechanics and Engineering, 93, pp. 189216. [67] Schures, P. J. G., Veldpaus, F. E . , and Brekelmans, W. A. M . , 1986, "Simulation of Forming Processes Using the Arbitrary Eulerian-Lagrangian Formulation", Computer Methods in Applied Mechanics and Engineering, 58, pp. 19-36. [68] Benson, D. J., 1989, "An Efficient, Accurate, Simple ALE Method for nonlinear Finite Element Programs", Computer Methods in Applied Mechanics and Engineering, 72, pp. 305- 350. [69] Ghosh, S. and Kikuchi, N., 1988, "Finite Element Formulation for the Simulation of Hot Sheet Metal Forming Process", Int. Journal of Engineering Sciences, 26, pp. 143-161. [70] Huerta, A., and Casadei, F., 1994, "New ALE Applications in Nonlinear FastTransient Solid Dynamic", Engineering Computations, 11, pp. 317-345. [71] Movahhedy, M . R., Gadala M . S., and Altintas, Y., 1999, "Simulation of Orthogonal Metal Cutting Process Using Arbitrary Lagrangian-Eulerian Finite Element Method", Materials Processing Technology, to appear. Bibliography 175 [72] Wang, J . and Gadala M . S.,1996, "Formulation and survey of A L E method in nonlinear solid mechanics,Finite Elements in Analysis and Design,24 , pp. 253269. [73] Hashemi, J . and Stefan, R., 1994, "Mesa Model of Orthogonal Cutting Process", Recent Advances in Structural Mechanics, ASME-PVP, 295, pp. 113-125. [74] Olovsson, L . , Nilsson, L . , and Simonsson, K . , 1999, "An A L E Formulation for the Solution of Two-Dimensional Metal Cutting Problems", Computers and Structures, 72, pp. 497-507. [75] Rakotomalala, R. and Joyot, P., 1993, "Arbitrary Lagrangian-Eulerian Thermomechanical Finite Element Model of Material Cutting", Communications in Nu- merical methods in Engineering, 9, pp. 975-987. [76] Pantale, 0., Rakotomalala, R., Touratier, M . , Hakem, N . , 1996, "A Three Dimensional Numerical Model of Orthogonal and Oblique Metal Cutting Processes", Engineering Systems Design and Analysis, ASME-PD, 75, pp. 199-205. [77] Joyot, P., Rakotomalala, R., Pantale, 0., Touratier, M . , and Hakem, N . , 1998, "A Numerical Simulation of Steady State Metal Cutting", Proc. Institution of Mechanical Engineers, Part C, 212, pp. 331-341. [78] Bathe, K . ,1996, Finite Element Procedures, Prentice-Hall Inc., New Jersey, U S A . [79] Truesdell, C. and Toupin, R., 1960, "The Classical Field Theories", Encyclopedia of Physics Vol. III-l, pp. 552-610, Springer, Berlin. [80] Haber, R., Shephard, M.S., Abel, J . E . , Gallagher, R . H . , and Greenberg, D.P., 1981, "A General Two Dimensional Graphic Finite Element Preprocessor Utilizing Discrete Transfinite Mapping", Int. Journal of Numerical Methods in Engineering, 17, pp. 1015-1044. [81] Perzyna, P., 1966, "Fundamental Problems in Viscoplasticity", Advances in Ap- plied Mechanics, 9, pp. 243-377. [82] Zienkiewicz, 0. C , Jain, P. C. and Onate, E . , 1978, "Flow of Solids During Forming and Extrusion: Some Aspects of Numerical Analysis", Int. Journal of Solids and Structures, 14 , pp. 15-38. [83] Dawson, P. R. and Thompson, E . G.,1978, "Finite Element Analysis of Steady State Elasto-viscoplastic Flow by the Initial Stress Rate method", Int. Journal of Numerical Methods in Engineering, 12, pp. 47-57. Bibliography 176 Hughes, T. J . R., and Taylor, R. L., 1978, "Unconditionally Stable Algorithms for Quasi-Static Elasto/Visco-Plastic Finite Element Analysis", Computers and Structures, 8, pp. 169-173. Hart, E.W., 1976, "Constitutive Relation for the Nonelastic Deformation of Metals", Journal of Engineering Materials and Technology, , pp. 193-201. Anand, L. and Zavaliangos, A., 1990, "Hot Working - Constitutive Equations and Computational Procedures", Annals of the CIRP, 39/1, pp. 235-238. Hsu, T. R., 1986, The Finite Element Method in Thermomechanics, Prentice-Hall Inc., New Jersey, USA. Huang, Y. M . , and Lu, Y. H., 1991, "An Elasto-Plastic Rate Dependent Finite Element Analysis of the Metal Forming Process", Computers and Structures, 6, pp. 615-622. Rice, J.R., 1975, "Continuum Mechanics and Thermodynamics of Plasticity in Relation to Microscale Deformation Mechanisms", Constitutive equations in plasticity, A.S. Argon (ed.), MIT press. Kobayashi, S., 1984, "Thermo-viscoplastic analysis of metal forming problems by the finite element method", Numerical Analysis of Forming Processes, Pittman, et al. [Eds.], John Wiley and Sons, pp. 45-70. Huebner, K. H., Thornton, E. A., and Byrom, T. G., 1995, The Finite Element Method for Engineers, John Wiley and Sons, Inc., New York, USA. Kardestuncer, H., 1987, Finite Element Handbook, McGraw-Hill Book Co., NY, USA. Oden, J . T., and Pires, E. B., 1983, "Nonlocal and Nonlinear Friction law and Variational Principles for Contact Problems in Elasticity", Journal of Applied Mechanics, 50, pp. 67-76. "ADINA-A Finite Element Program for Automatic Dynamic Incremental Nonlinear Analysis", Report AE81-1, ADYNA Engineering, Watertown, MA 02172. Bathe, k. J., and Chaudhary, A., 1985, "A Solution Method for Planar and Axisymmetric Contact Problems", Int. Journal of Numerical Methods in Engineering, 21, pp. 65-88. Pascoe S. K., and Motershead, J. E., 1988, "Linear Elastic Contact problems Using Curved Elements and Including Dynamic Friction", Int. Journal of Numerical Methods in Engineering, 26, pp. 1631-1643. Bibliography 177 [97] Saran, M . J. and Wagoner, R. H. , 1991, "A Consistent Implicit Formulation for Nonlinear Finite Element Modeling with Contact and Friction", Journal of Applied Mechanics, 58, pp. 499-512. [98] HaUquist, J. 0., 1979, "NIKE2D: An Implicit Finite Deformation, Finite Element Code for Analyzing the Static and Dynamic Response of Two Dimensional Solids", Report UCRL-25678, Lawrence Livermore Laboratory, University of California. [99] Oden, J. T., 1981, "Exterior Penalty Methods for Contact Problems in Elasticity", Nonlinear Finite Element Analysis in Structural Mechanics, W. Wunderlich, et al. (Eds.), Springer Verlag, New York, USA. [100] Flippa, C. A., 1978, "Iterative Procedures for Improving Penalty Function Solutions of Algebraic Systems", Int. Journal of Numerical Methods in Engineering, 12, pp. 821-836. [101] Engineering Analysis System User Manual, Swanson Analysis Systems Inc. Houston, Texas, 1992. [102] Negtegaal, J. C , and Rebelo, N., 1988, "On the Development of A General Purpose Finite Element Program for Analysis of Forming Processes", Int. Journal of Numerical Methods in Engineering, 26, pp. 113-131. [103] Germain, Y., Chung, K., and Wagoner, R. H., 1989, "A Rigid Viscoplastic Finite Element Program for Sheet Metal Forming Analysis", Int. Journal of Mechanical Sciences, 31, pp. 1-24. [104] Malvern, L. E . , 1969, Introduction to the Mechanics of a Continuous Medium, Prentice HaU Inc., NJ, USA. [105] Voyiadjis, G. Z., Poe, A. A., and Kiousis, P. D., 1986, "Finite Strain Elasto-Plastic Solution for Contact Problems", Journal of Engineering Mechanics, 112, pp. 273292. [106] Hariandja, B. H. and Haber, R.B., 1986, Adaptive Finite Element Analysis of Nonlinear Frictional Contact with Mixed Eulerian-Lagrangian Coordinates, Technical Report, University of Illinois, USA. [107] Gadala, M . S. and Wang, J., 1998, "A Practical Procedure for Mesh Motion in Arbitrary Lagrangian-Eulerian Method", Engineering with Computers, 14, pp. 223234. [108] Bayoumi, H. N. and Gadala, M . S., 1999, "Finite Element Analysis of Large Strain Solid Mechanics Problems", Proceedings of the first Canadian conference on nonlinear problems in solid mechanics Victoria, Canada, pp. 375-384. Bibliography 178 in Finite Elements in Analysis and Design. [109] Giuliani, S., 1982, "An Algorithm for Continuous Rezoning of the Hydrodynamic Grid in Arbitrary Lagrangian-Eulerian Computer Codes", Nuclear Engineering and Design, pp. 205-212. [110] Zienkiewicz, 0. C. and Philips, D. V., 1971,"An Automatic Mesh Generation Scheme for Plane and Curved Surfaces by Isoparametric Coordinates", Int. J. Numerical Methods in Engineering , 3, pp. 519-528. [Ill] Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P, 1992, Numerical Recipes in Fortran, Cambridge University Press. [112] Koh, H. M . , Lee, H. S. and Haber, R. B., 1988, "Dynamic Crack Propagation Analysis Using Eulerian-Lagrangian Kinematic Description", Computational Mechanics, 3, pp. 141-155. [113] Movahhedy, M . R., Gadala, M . S., and Altintas, Y., 1998, "Simulation of ElasticPlastic Crack Propagation Using Arbitrary Lagrangian-Eulerian Finite Element Method", Proc. 12th ASCE Engineering Mechanics Conference, Murakami, H., et al. (eds.), La Jolla, California, pp. 110-113. [114] Engineering Mechanics Research Corporation, 1997, Users Manual for NISA Systems Analysis, Michigan, USA. [115] Movahhedy, M . R., Gadala, M . S. and Altintas, Y., 2000, "FE Modeling of Chip Formation in Orthogonal Metal Cutting Process: An ALE Approach", Machining Science and Technology, 4, pp. 15-42. [116] Derbalian, K. A., Lee, E . H., MaUet, R. L., and McMeeking, R. M . , 1978, "Finite Element Metal Forming Analysis with Spatially Fixed Mesh", Application of Numerical Methods in Forming Processes, ASME-AMD, 28, pp. 39-47. [117] Hinton, E . and Campbell, J . S., 1974, "Local and Global Smoothing of Discontinuous Finite Element Functions Using a Least Square Method", Int. Journal for Numerical Methods in Engineering, 8, pp. 461-480. [118] A M T E C Engineering, Inc., 1996, TECPLOT User's Manual, Version 7, Washing- ton, USA [119] Wang J. and Gadala M . S., 1998, "ALE Formulation and Its Application in Solid Mechanics", Computer Methods in Applied Mechanics and Engineering, 167, pp. 33-55. Bibliography 179 [120] Gadala M.S. and Wang, J., 1999, "Simulation of Metal Forming Processes with Finite Element Method", Int. Journal for Numerical Methods in Engineering, 44, pp. 1397-1428. [121] Hughes, T. J. R. and Winget, J . , 1980, "Finite Rotation Effects in numerical Integration of Rate Constitutive Equations Arising in Large Deformation Analysis", Int. Journal for Numerical Methods in Engineering, 15, pp. 1862-1867. [122] Pinsky, P. M., Ortiz, M., and Pister, K. S., 1983, "Numerical Integration of Rate Constitutive Equations in finite Deformation Analysis", Computer Methods in Applied Mechanics and Engineering, 40, pp. 137-158. [123] Ortiz, M . , and Popov, E. P., 1985, "Accuracy and Stability of Integration Algorithms for Elastoplastic Constitutive Relations", Int. J. Numerical Methods in Engineering, 21, pp. 1561-1576. [124] Ortiz, M . and Simo, J. C , 1985, "An Analysis of a New Class of Integration Algorithms for Elastoplastic Constitutive Relations", Int. Journal for Numerical Methods in Engineering, 23, pp. 353-366. [125] Owen, D. J. R., and Hinton, E., 1980, Finite Elements in Plasticity, Theory and Practice, Pineridge press Ltd., Swansea, UK. [126] Boothroyd, G., 1963, "Temperatures in Orthogonal Metal Cutting", Proceedings of Institution of Mechanical Engineers, 177, pp. 789-802. [127] Stephenson, D. A., 1989, "Material Characterization for Metal Cutting Force Modeling", Journal of Engineering Materials and Technology, 111, pp.210-219. [128] Ozel, T., and Altan, T., 1998, "Determination of Workpiece Flow Stress and Friction at the Chip-Tool Contact for High Speed Machining", Int. Journal of Machine Tools and Manufacture, to appear. [129] Lei, S., Shin, Y. C , and Incropera, F. P., 1999, "Material Constitutive Modeling Under High Strain Rates and Temperatures Through Orthogonal Machining Tests", Journal Manufacturing Science and Engineering, 121, pp. 577-585. [130] Van Luttervelt, C.A., Childs, T.H.C., Jawahir, I.S., Klocke, F., and Venuvinod, P.K., 1998, "Present Situation and Future Trends in Modeling of Machining Operations", Annals of the CIRP, 47(2) , pp. 587-628. [131] Childs, T.H.C., and Maekawa, K., 1990, "Computer Aided Simulation and Experimental Studies of Chip Flow and Tool Wear in the Turning of Low Alloy Steels by Cemented Carbide Tools", Wear, 139, pp. 235-250. Bibliography 180 Altan, T., FaUbohmer, P., Rodrigues, C. A., and Ozel, T., 1998, "High Speed Cutting of Cast Iron and Alloy Steels - State of Research",Proceedings of CIRPVDI Conf. on High Performance Tools, Dusseldorf, Germany, pp. 309-330. Armarego, E.J.A., 1993, Material Removal Processes- An Intermediate Course, University of Melbourne, Australia. Arsecularatne, J.A., 1997, "On Tool-Chip Inerface Stress Distributions, Ploughing Force and Size Effect in Machining", int. Journal of Machine Tools and Manufacture, 37, pp. 885-899. Waldorf, D. J., DeVor, R. E . , and Kapoor, S. G., 1997, "A Slip-line Field for Ploughing During Orthogonal Cutting", Proceedings of 1997 ASME International Mechanical Engineering Congress and Exposition, MED 6-1, pp. 121-128. Elanayar, S. and Shin, Y.C., 1996, "Modeling of Tool Forces for Worn Tools: Flank Wear Effects", Journal of Manufacturing Science and Technology, 118, pp. 359366. Montgomery, D., Altintas, Y., 1991, "Mechanism of Cutting Force and Surface Generation in Dynamic Milling", Transactions of ASME, Journal of Engineering for Industry, 113, pp. 160-168. Waldorf, D. J., DeVor, R. E., and Kapoor, S. G., 1999, "An Evaluation of ploughing Models for Orthogonal Machining", Journal of Manufacturing Science and Engineering, 121, pp. 550-558. Hirao, M., Tlusty, R., Sowerby, R., and Chandra, G., 1982, "Chip Formation with Chamfered Tools", Journal of Engineering for Industry, 104, pp. 339-342. Zhang, H. T., Liu, P. D. and Hu, R. S., 1991, "A Three Zone Model and Solution of Shear Angle in Orthogonal Machining", Wear, 143, pp. 29-43. Ren, H., 1998, Mechanics of Machining with Chamfered Tools, M.Sc. Thesis, UBC. Engineering Research Center, Net Shape Manufacturing, Ohio State University, Private communication. Shatla, M . and Altan, T., 1998, "Use of Analytical Based Computer Modeling to Analyze 2-D and 3-D Turning and Face Milling Operations", ERC/NSM Report No. HPM/ERC/NSM-98-R-028. Jacobson, S. and Wallen, P., 1988, "A new Classification System for Dead Zone in Metal Cutting", Int. Journal of Machine Tools and Manufacture, 28, pp. 529-538. Bibliography 181 [145] Lee, P. and Altintas, Y . , 1996, "Prediction of Ball-End Milling Forces from Orthogonal Cutting Data", Int. Journal of Machine Tools and Manufacture, 36, pp. 1059-1072.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Ale simulation of chip formation in orthogonal metal...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Ale simulation of chip formation in orthogonal metal cutting process Movahhedy, Mohammad R. 2000
pdf
Page Metadata
Item Metadata
Title | Ale simulation of chip formation in orthogonal metal cutting process |
Creator |
Movahhedy, Mohammad R. |
Date Issued | 2000 |
Description | This thesis deals with the application of the Arbitrary Lagrangian-Eulerian (ALE) finite element analysis in simulation of chip formation in orthogonal metal cutting process. A critical review of the literature in this field shows that due to the very complex set of conditions present in a cutting process, the application of conventional Lagrangian and Eulerian methods for this problem is inefficient and entails numerical difficulties. In particular, the pertinent problems in the node separation technique or the remeshing approach are discussed. In contrast, the adaptivity of the mesh in an ALE analysis provides the possibility of combining the strengths of both Lagrangian and Eulerian methods in a single analysis. In this approach, the chip formation occurs as a result of plastic flow of the work material around the tool edge on the one hand, and unconstrained flow of the material on free surfaces of the chip, on the other. Due to high deformation speed, strain rate and temperature play a significant role in the chip formation process. In this work, the ALE formulation originally proposed by Wang and Gadala [72] is extended to include rate and thermal effects. A heat transfer module is included that updates the temperature field in the cutting zone at each step of the analysis. Contact algorithms are developed which are able to detect emerging contact conditions and apply contact constraints at interfaces between flexible-flexible or rigid-flexible pairs. An efficient ALE mesh motion is designed that prevents element distortion in the deformation zone and at the same time facilitates the evolution of the chip size at free boundaries. Furthermore, General guidelines for designing a mesh motion strategy are presented, an algorithm for mesh sliding on free boundaries is introduced, and transfinite and isoparametric mapping techniques are adopted for moving the mesh in the interior of the body, so that the mesh remains optimal throughout the analysis The large deformation, rate-dependent, thermo-mechanical ALE finite element code is used to simulate chip formation in orthogonal metal cutting processes. The results of simulation of cutting low carbon steel with a carbide tool are presented as a benchmark problem, and a parametric study is conducted to investigate the influence of cutting conditions on the chip formation process. The results of these studies are verified through comparison with available experimental data. The fairly good qualitative and quantitative agreement between predicted and experimental results confirms that ALE simulation provides a realistic representation of the actual process. Finally, the influence of cutting edge geometry on the chip formation process is investigated through simulation of high speed cutting of hardened steel alloys with chamfered or worn tools of carbide or CBN type. This study shows that changing the chamfer angle does not affect the chip significantly, because the dead metal zone that is formed under the chamfer acts as the main cutting edge of the tool. However, cutting and thrust forces increase with increase in the chamfer angle. The predictions on the effects of chamfer angle are in agreement with experimental observations. A study of the influence of cutting speed on the deformation process for chamfered CBN tools shows that at higher cutting speeds, the maximum temperature on the rake face increases substantially, signifying the importance of diffusion wear at such speeds. |
Extent | 13901477 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-07-27 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0089713 |
URI | http://hdl.handle.net/2429/11341 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2000-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2000-48680X.pdf [ 13.26MB ]
- Metadata
- JSON: 831-1.0089713.json
- JSON-LD: 831-1.0089713-ld.json
- RDF/XML (Pretty): 831-1.0089713-rdf.xml
- RDF/JSON: 831-1.0089713-rdf.json
- Turtle: 831-1.0089713-turtle.txt
- N-Triples: 831-1.0089713-rdf-ntriples.txt
- Original Record: 831-1.0089713-source.json
- Full Text
- 831-1.0089713-fulltext.txt
- Citation
- 831-1.0089713.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0089713/manifest