ARBITRARY LAGRANGIAN-EULERIAN FORMULATION FOR QUASI-STATIC AND DYNAMIC M E T A L FORMING SIMULATION By Hassan N . Bayoumi B.Sc, M.Sc. Cairo University, Egypt, ,1992, 1995 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF DOCTOR OF PHILOSOPHY in T H E F A C U L T Y OF G R A D U A T E STUDIES D E P A R T M E N T OF M E C H A N I C A L ENGINEERING We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA May 2001 © Hassan N . Bayoumi, 2001 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: ^ L X V * . ^ , ZJ2>0\ Abstract Many engineering problems involve large material deformation, large boundary motion and continuous changes in boundary conditions. The Arbitrary Lagrangian Eulerian (ALE) formulation has emerged in recent years as a technique that can alleviate many of the shortcomings of the traditional Lagrangian and Eulerian formulations in handling these types of problems. Using the A L E formulation the computational grid need not adhere to the material (Lagrangian) nor be fixed in space (Eulerian) but can be moved arbitrarily. Two distinct techniques are being used to implement the A L E formulation, namely the operator split approach and the fully coupled approach. A survey of the A L E literature shows that the majority of A L E implementations for quasi-static and dynamic analyses are based on the computationally convenient operator split technique. In addition, all previous dynamic A L E formulations are based on explicit time integration where no linearization is needed. This thesis presents a fully coupled implicit A L E formulation for the simulation of quasi-static and dynamic large deformation and metal forming problems. A L E virtual work equations are derived from the basic principles of continuum mechanics. A new method for the treatment of convective terms that sidesteps the computation of the spatial gradients of stresses is used in the derivation. The A L E virtual work equations are discretized using isoparametric finite elements. Full expression for the resulting A L E finite element matrices and vectors are given. A new relation that relates grid displacements with material displacements is introduced. ii The A L E finite element equations are implemented into a 2-D computer code for plane stress, plane strain and axisymmetric problems. The transfinite mapping method is used as the mesh motion scheme for internal nodes. A new treatment for mesh motion on material boundaries is introduced and implemented. Implicit, explicit and mixed implicitexplicit time integration schemes are implemented in the code. A line search technique is employed to accelerate the convergence of implicit calculations. Several quasi-static and dynamic large deformation applications are solved using the developed code. Experimental analysis of a simple V-bending process is conducted for comparison. Comparison of A L E predictions for deformed shapes, equivalent stress and plastic strain distributions and loading curves with analytical, numerical and experimental results are presented. A L E results are in good agreement with other methods of analysis. A L E is shown to prevent mesh distortion and eliminate the need for special contact treatments for problems with corner contact. 111 Table of Contents Abstract ii List of Tables ix List of Figures x Nomenclature xiii Acknowledgment xviii 1 Introduction 1 1.1 Motivation 1 1.2 Problem Description 2 1.2.1 Traditional Large Strain Formulations 2 1.2.2 The Arbitrary Lagrangian Eulerian (ALE) Formulation 6 1.3 Scope of Work 10 2 Background 13 2.1 A L E Literature Review 13 2.2 Dynamic Effects 19 2.2.1 19 Fast Transient Metal Forming iv 2.2.2 2.3 Solution of Nonlinear Equilibrium Equations in Finite Elements Summary 23 3 Derivation of A L E Governing Equations 3.1 3.2 3.3 20 25 Preliminaries 25 3.1.1 Notations 25 3.1.2 Kinematics 26 3.1.3 Continuity 28 Quasi-Static Analysis 28 3.2.1 Principle of Virtual Displacements 29 3.2.2 Incremental Decompositions 30 3.2.3 Linearization 32 3.2.4 Treatment of Convective Terms 33 3.2.5 Fully Coupled A L E Equilibrium Equation 35 Dynamic Analysis 37 3.3.1 Virtual Work Done by Inertia Forces 37 3.3.2 Decomposition of Velocities and Accelerations 38 3.3.3 Linearization of the Referential Inertia Term 38 3.3.4 Linearization of the Convective Inertia Term 39 3.3.5 Fully Coupled A L E Equation of Motion 41 v 4 Finite Element Discretization 43 4.1 Isoparametric Finite Elements 43 4.2 Discretiztion of the Quasi-Static A L E Equation 45 4.2.1 Discretization of the Lagrangian Internal Force Term 45 4.2.2 Discretization of the Lagrangian Material Stiffness Term 47 4.2.3 Discretization of the Lagrangian Geometric Stiffness Term 48 4.2.4 Discretization of the First Convective Stiffness Term 49 4.2.5 Discretization of the Second Convective Stiffness Term 50 4.2.6 Discretized Quasi-Static A L E Equation 52 4.3 Discretiztion of the Dynamic A L E Equation 52 4.3.1 Discretization of the Incremental Mass Term 54 4.3.2 Discretization of the First Velocity-Stiffness Term 54 4.3.3 Discretization of the Second Velocity-Stiffness Term 55 4.3.4 Discretization of the Inertia Force Term 56 4.3.5 Discretization of the First Convective Inertia Force Term 56 4.3.6 Discretization of the Second Convective Inertia Force Term 57 4.3.7 Discretization of the Third Convective Inertia Force Term 57 4.3.8 Discretization of the Fourth Convective Inertia Force Term 58 4.3.9 Discretized Dynamic A L E Equation 59 vi 5 Implementation 5.1 5.2 5.3 5.4 60 Mesh Motion 60 5.1.1 Grid Displacement 60 5.1.2 Mesh Motion for Interior Nodes 61 5.1.3 Mesh Motion on Free Material Boundaries 62 Solution of Nonlinear Equilibrium Equations 66 5.2.1 Quasi-static Analysis 66 5.2.2 Dynamic Analysis 67 5.2.3 Elimination of Grid Displacements on the Element Level 70 5.2.4 Frontal Solver 71 5.2.5 Line Search 72 Contact Analysis 73 5.3.1 Lagrange Multipliers 73 5.3.2 Penalty Function 74 5.3.3 Direct Constraint Method 75 Program Structure 76 5.4.1 Updated Lagrangian Routines 77 5.4.2 A L E Routines 81 5.4.3 Contact Routines 84 vii •V * 87 6 Applications 6.1 One-Dimensional Stress Wave Propagation 87 6.2 Bar Impact 90 6.3 Sheet Metal Extrusion 94 6.4 Quasi-static and Dynamic Coining 97 6.5 V-Bending 105 6.6 Plate Indentation 121 6.7 Quasi-static and Dynamic Brake Bending 128 7 Conclusions 133 7.1 Summary 133 7.2 Future Work 135 Bibliography 1 3 9 Appendix A 1 4 6 vm List o f Tables Table 6.1 Comparison of results for bar impact 91 Table 6.2 Comparison of stresses at punch corner after the first load step 99 ix List of Figures Figure 1.1 Simulation of punch indentation using the Lagrangian approach Figure 5.1 Transfinite mapping of a distorted mesh region 63 Figure 5.2 Mesh motion on free material boundaries 63 Figure 5.3 Flow chart of the Updated Lagrangian calculations 80 Figure 5.4 Flow chart of the A L E calculations 83 Figure 5.5 Flow chart with contact routines Figure 6.1 Mesh for the one-dimensional wave propagation problem 88 Figure 6.2 Stress wave amplitude and duration 88 Figure 6.3 Longitudinal stress distribution comparison, elastic case 89 Figure 6.4 Longitudinal stress distribution comparison, elastic-plastic case 89 Figure 6.5 Comparison of the Lagrangian and A L E solutions for bar impact 92 Figure 6.6 Equivalent stress and plastic strain contours for bar impact 93 Figure 6.7 Geometry and mesh for extrusion process 95 Figure 6.8 Plastic strain contours and deformed shape for extrusion problem 96 Figure 6.9 Longitudinal stress at different lateral positions 96 Figure 6.10 Geometry and initial mesh for coining process 100 Figure 6.11 ANSYS solution with no contact elements 101 Figure 6.12 A N S Y S solution with default contact parameters 101 x 4 85-86 Figure 6.13 ANSYS solution with extreme contact parameters 101 Figure 6.14 Evolution of the deformed shape during the coining process 102 Figure 6.15 Punch load-displacement curve comparison 103 Figure 6.16 Final plastic strain distribution at different punch velocities 104 Figure 6.17 Isometric sketch for the setting of the V-bending experiment , 107 Figure 6.18 Tinus Olsen U T M used in the V-bending experiment 108 Figure 6.19 Finite element mesh used in the simulation of the V-bending process.. 109 Figure 6.20 Development of the deformed shape during the experiment 110-112 Figure 6.21 Development of the deformed shape using M A R C 112(-115 Figure 6.22 Development of the deformed shape using A L E 116-118 Figure 6.23 Deformed shape using M A R C showing corner penetration 119 Figure 6.24 Comparison of load displacement curves for V-bending 120 Figure 6.25 Comparison of loading curve for the flat-faced indenter 123 Figure 6.26 Comparison of loading curve for the hemispherically-tipped indenter...124 Figure 6.27 Plate deformed shape under flat-faced indenter 125 Figure 6.28 Plate deformed shape under hemispherically-tipped indenter 125 Figure 6.29 Plastic strain distribution for flat-faced indenter 126 Figure 6.30 Plastic strain distribution for hemispherically-tipped indenter 126 Figure 6.31 Comparison of the loading curve for flat and hemispherical indenters.. .127 Figure 6.32 Tooling geometry for the brake bending process xi 129 Figure 6.33 Quasi-static load-displacement curve for the brake bending process ... 130 Figure 6.34 Dynamic load-displacement curve for the brake bending process 131 Figure 6.35 Plastic strain distribution for the brake bending process 132 xii Nomenclature quantity occurs at time t '() quantity referred to time t 'a element nodal referential material acceleration vector a element nodal incremental referential material acceleration vector 'a components of referential material acceleration vector a components of incremental referential material acceleration vector j i B element shape function derivative matrix related to first A L E stiffness B element shape function derivative matrix related to second A L E stiffness Ai t t A2 ,B i l element shape function derivative matrix related to Lagrangian material stiffness ,B i 2 element shape function derivative matrix related to Lagrangian geometric stiffness 'C element elastic-plastic material constitutive matrix 'C element first convective velocity-stiffness matrix due to A L E 'C element second convective velocity-stiffness matrix due to A L E M A2 'C element third convective velocity-stiffness matrix due to A L E 'C element fourth convective velocity-stiffness matrix due to A L E A3 A4 'C ijkl components of fourth order material constitutive tensor 'D kl components of rate of deformation tensor 'dV elemental volume xm 'dS elemental surface area E Young's modulus EP plastic modulus ,e element strain vector ,e components of strain tensor f* element equivalent nodal force vector 'f element internal force vector '/ arbitrary function '/ material derivative of arbitrary function '/' grid derivative of arbitrary function 'f components of body force vector per unit mass ' ff components of surface traction vector H element interpolation matrix I identity matrix/tensor tj 8 h element shape function at node a ' K* element equivalent stiffness matrix 'K element convective stiffness matrix due to A L E a A 'K^ 'K A2 'K L 1 element first convective stiffness matrix due to A L E element second convective stiffness matrix due to A L E element Lagrangian stiffness matrix xiv 'K element Lagrangian material stiffness matrix 'K element Lagrangian geometric stiffness matrix u 12 'M element convective mass matrix due to A L E A 'M element Lagrangian mass matrix N number of element nodal points L 'n components of unit normal to boundary surface 'S element stress matrix related to first A L E stiffness j A> 'S i 2 element stress matrix related to Lagrangian stiffness s line search parameter t time u element incremental material displacement vector u element incremental grid displacement vector u element nodal incremental material displacement vector 8 u element nodal incremental grid displacement vector w, components of incremental material displacement vector uf components of incremental grid displacement vector u components of nodal incremental material displacement vector s jk uf components of nodal incremental grid displacement vector v element nodal incremental material velocity vector k \ 8 'v element nodal incremental grid velocity vector element nodal material velocity vector xv 'v element nodal grid velocity vector g v. components of incremental material velocity vector vf components of incremental grid velocity vector 'v,. components of material point velocity vector 'vf components of grid point velocity vector 'v,. components of material acceleration vector 'v'j components of material referential acceleration vector ( 'W external work 'x element coordinate vector 'x element nodal coordinate vector X"' components of material point initial coordinates Xf components of grid point initial coordinates 'Xj components of material and grid points position vector 'x'" components of material point mapping function 'xf components of grid point mapping function 'x jk components of nodal coordinate vector 'y ijk components of stress-velocity product tensor a mesh motion parameter vector a. components of mesh motion parameter vector exl xvi B mesh motion parameter matrix P Newmark integration parameter R components of mesh motion parameter vector Ar time increment 5 virtual or variation (j) s p eq equivalent strain fa boundary curve y Newmark integration parameter v Poisson's ratio 6 angle of inclination of x' and y' axes to the global x and y axes 'p material density 'o element stress vector a equivalent stress a initial yield stress 'ay components of Cauchy stress tensor 'cjy components of Truesdell stress rate tensor Lr nodal stress components at node a y ija xvn Acknowledgement Sincere thanks and appreciation are expressed to my supervisor, Professor Mohamed Gadala, for suggesting the area of this research and for his invaluable advice and continuous guidance. His knowledge and experience greatly influenced this thesis. I am very glad for having the opportunity to study at the University of British Columbia. I would like to acknowledge the funding of my supervisor as well as of the Faculty of Graduate Studies. I am honored to be a recipient of the University Graduate Fellowship for four years. I would also like to thank my professors and previous supervisors. Special thanks are due to Professors Salah Bayoumi, Mohsen El-Mahdi, Ahmed Mansour and Aly E l Shafei. I am greatly indebted to them all. I am grateful to my parents, my brother and my friends for their steadfast support. Thanks are due to my friend Ahmed El-Bindari for being a constant source of help and encouragement. Finally, a very special thank you goes to my wife, Ingy, for her patience, love and dedication. xviii Chapter 1 INTRODUCTION 1.1 MOTIVATION The history of metal forming is based on a highly empirical, nearly artisan form of technology. Perhaps one of the most expensive and time-consuming stages of the metal forming design and manufacturing cycle is prototype development. Prior to 1980's, prototype development was mainly based on trial and error. This is in a large part due to the fact that metal forming problems were too complex to model using the available tools. The decade of the 1980's resulted in revolutionary technological advancements for metal forming industries. The vast increase in computational capabilities together with the recent advances in numerical techniques made it possible to simulate metal forming processes in the design stage. The use of the finite element method in the metal forming design process is practiced by a wide range of industries including aerospace, automotive, medical, defense and other diverse production companies. Metal forming simulation can predict possible forming defects, potential operating problems and provide optimum tooling design and running conditions leading to the reduction or elimination of the costly tool prototyping. As a by-product of the simulation, a better understanding of the mechanical properties of forming processes can be acquired. Prediction of forming loads, final workpiece geometry, potential workpiece cracking or damage, incomplete die 1 Chapter 1. 2 INTRODUCTION filling, residual stresses, excessive work hardening, sheet metal wrinkle, tearing and spring-back, optimum tooling curvatures, tool wear, and blank dimensions are examples of the areas that can be addressed by process simulation. Though metal forming simulation has significantly improved in recent years, it is noted, however, that the accuracy and utility of metal forming simulation are still questionable. 1.2 PROBLEM DESCRIPTION Metal forming processes are complex problems that generally involve large deformation, nonlinear material behavior, large boundary motion and interaction of the workpiece with the forming tools and dies. Finite element simulation of metal forming problems necessitates the use of a large strain formulation that can handle material and geometric nonlinearities. 1.2.1 Traditional large strain formulations Foundations of large strain analysis of elastic-plastic solids can be traced back to the early work of Hill [1]. The first finite element formulation for large strain problems, known as the Total Lagrangian formulation, was introduced by Hibbitt et al. [2]. Later on, McMeeking and Rice [3] pioneered the use of the Updated Lagrangian formulation. The conceptual difference between the two formulations is the reference configuration that is used for the linearization of the incremental equations of motion. In the Total Lagrangian formulation the initial configuration is used as a reference, whereas in the Updated Lagrangian formulation, the reference configuration corresponds to the last Chapter 1. INTRODUCTION 3 calculated configuration. Being referential in nature, both formulations should give the same results [4], provided that the two formulations are consistently driven. Finite element computer codes that are currently being used in metal forming simulation can be classified into two categories: general-purpose finite element codes [510] and codes that are developed specifically for metal forming simulation [11-16]. Both types of codes mainly depend on the Lagrangian (referential) formulation to simulate the large deformation behavior that is inherent in metal forming problems. In the Lagrangian formulation, the finite element mesh, or reference configuration, is fixed to the material points of the deforming body and the grid moves with the same velocity as the material. In the case of large strain problems, this can lead to excessive distortion of the finite element mesh and presents a major drawback of the formulation. Distorted meshes reduce the accuracy of numerical integrations and may ultimately result in singular matrices and computation termination. Another drawback of the Lagrangian approach is the difficulty to model non-material-associated boundary conditions. Figure 1.1 shows a punch indentation problem that has been analyzed using a general-purpose finite element package, NISA [7], which is based on the Lagrangian formulation. Because of symmetry, only half of the geometry is being modeled. The punch movement has been described as a downward constant velocity applied at the nodes right below the punch. It is required to estimate the forming loads for a 60 % reduction in the height of the workpiece. The figure shows that as the material moves under the punch load, the finite element mesh, being attached to the material, also moves with the effect of increasing the punch diameter. It is also clear from the figure that starting at approximately 50 % reduction in Chapter 1. INTRODUCTION 4 Punch Initial i itton Actual punch size Enlarged punch 1 Excessive ' mesh distortion r 40 %reductionin height 55 % reduction in height Figure 1.1. Simulation of punch indentation using the Lagrangian approach (NIS A). Chapter 1. INTRODUCTION 5 height, the finite element mesh was excessively distorted and the run was terminated at a 55 % height reduction. A more detailed description of the problem is given in Chapter 6. Using general-purpose codes, the only remedy for the mesh distortion problem is to remesh the material domain with a new finite element mesh whenever needed [17]. The material properties for the new mesh are found from the corresponding old material points by interpolation between the old and new Gauss points, where certain approximations have to be introduced. Computations are then restarted from the last converged step. Remeshing, besides being very time consuming, sometimes requires user intervention. The simulation process may need frequent remeshing to reach the required deformation level. The analyst must either run the code interactively or have enough knowledge of the solution to write a detailed set of remeshing instructions ahead of time. Only a few general-purpose codes have automatic remeshing capabilities. On the other hand, metal forming codes are usually equipped by the auto-remesh capability [11]. Automatic remeshing, though still very time consuming, can be adjusted to take place either at a specified number of incremental load steps or when excessive mesh distortion occurs. However, the problem of incorrect modeling of non-material-associated boundary conditions will persist between remeshes in addition to the approximations of the assumed interpolation. Meanwhile, several attempts have been aiming at the adaptation of Eulerian (spatial) formulation, widely used in fluid flow simulations, to large strain and metal forming problems [19]. Some special metal forming codes use the Eulerian formulation for the simulation of steady state forming processes such as extrusion and drawing [20]. In the Chapter 1. INTRODUCTION 6 Eulerian description, the finite element mesh is fixed in space. The obvious drawbacks of the Eulerian description are the difficulty to model changing material boundaries and the difficulty to model material-associated boundary conditions. This limits the applicability of the Eulerian description to the analysis of a small class of metal forming problems of the steady-state type. A n extensive bibliography of the application of finite element techniques in metal forming simulation is given in [21]. 1.2.2 The Arbitrary Lagrangian Eulerian (ALE) formulation The above discussion indicates the need for a formulation that is more suited for the simulation of large deformation problems with large boundary motions. The Arbitrary Lagrangian Eulerian (ALE) formulation has emerged in recent years as a technique that can alleviate many of the drawbacks of the traditional Lagrangian and Eulerian formulations [22, 23]. In the A L E formulation, the finite element mesh, or reference configuration, need not adhere to the material nor be fixed in space but can move arbitrarily. As the material deforms, the finite element mesh is continuously moved to meet any preset criterion (e.g. optimize elements shape) and the simulation should be completed without user intervention. Combining the merits of both the Lagrangian and Eulerian formulations, A L E can easily describe all types of boundary conditions and prevent mesh distortion. A proper A L E formulation should reduce to both the Lagrangian and Eulerian formulations as necessary. Different approaches have been adopted to develop the A L E formulation. Since the pure Lagrangian and the pure Eulerian viewpoints have complementary virtues, Chapter 1. INTRODUCTION 7 researchers have tried to produce a formulation that can combine the best of both. The differences among the available A L E formulations depend on the intended application, assumptions made in deriving the A L E equations and details of implementation. These factors determine the limitations and accuracy of the resulting formulation. A survey of the main features of the various A L E formulations available in the literature is given in Chapter 2. A L E is usually termed a coupled formulation since material deformation and convective effects are coupled in the same equations. However, two distinct techniques are being used to implement the A L E equations. The first technique is referred to as an 'operator split' or a 'fractional step' approach. Virtually all A L E analyses, with very few exceptions, are based on this strategy. In this approach, material deformation and convective effects are treated separately. Thus each time step may be divided into two steps: a regular Lagrangian step followed by an Eulerian step. In the Lagrangian step, the grid moves with the material, whereas in the Eulerian step, the Lagrangian solution is mapped to the reference grid and stresses are updated using convective effects. Time advances only during the Lagrangian step and there is no time associated with the Eulerian step. Thus it is not necessary to perform an Eulerian step in every time step. The alternative A L E approach, which has been used by fewer researchers, is known as the 'fully coupled' or the 'unsplit' approach. In this approach, the governing A L E equations are implemented and solved without disruption. A solution algorithm that can handle convective effects simultaneously with material deformation must be used. The main advantage of the operator split over the fully coupled approach is the Chapter 1. INTRODUCTION 8 reduction in the cost of implementation of A L E into current Lagrangian codes as the Lagrangian step is unchanged and only the Eulerian step algorithm needs to be added. Moreover, the decoupling of the Lagrangian and Eulerian steps results in simpler equations to be solved and simpler algorithms to be used. However, from the theoretical point of view, the fully coupled A L E approach represents a true kinematic description that employs a more rigorous scheme in considering equilibrium in each step relative to a moving reference configuration. Therefore, the operator split solution, which stems from computational convenience, is expected to be less accurate than the fully coupled solution. Although most of the A L E literature is based on the operator split methods, the author believes that as computers become faster, unsplit methods will probably dominate because of their theoretical accuracy advantage. As shall be inferred from the A L E literature review section, no attempt has been made to develop and implement a true fully coupled A L E formulation for general quasistatic and dynamic solid mechanics analyses. In addition, most of the available A L E formulations for solid mechanics applications are based on a generalization of the Eulerian formulation in which velocities are used as the independent variables. This treatment is not consistent with the traditional Lagrangian formulation that is more suited for solid mechanics applications and in which displacements are more natural to use as independent variables. The author believes that this treatment is one of the reasons for the misconception that A L E is more suited for fluid mechanics applications and for the delay in the wide use of A L E in solid mechanics applications. In deriving the A L E equations, the relationship between the material time derivative, Chapter 1. INTRODUCTION 9 grid time derivative and spatial derivatives of arbitrary quantities is substituted into the governing equations. This substitution gives rise to convective terms in the ALE equations which account for the transport of material through the grid. Convective terms are the terms that involve the spatial gradients of quantities such as stresses. Since the values of these quantities are available at the integration points, and not at the nodal points connecting the elements, these quantities are generally discontinuous across element edges and their gradients cannot be reliably computed on the element level. The problem of evaluating convective terms presents a difficult task in the ALE formulation and different numerical assumptions and treatments are being used. As a result, the treatment of convective terms is a key factor that distinguishes the different ALE formulations. Details of the derivation of the equations governing the ALE formulation are given in Chapter 3. Using the ALE formulation, the finite element mesh can be moved arbitrarily to maintain a homogeneous mesh and properly represent boundary conditions throughout the deformation process. A mesh motion scheme is necessary to continuously adapt the positions of internal nodes. On the boundaries, however, mesh motion must satisfy the boundary constraint that prevents the relative motion between the material points and mesh points in the direction normal to the boundary. This constraint ensures that the material and mesh configurations have the same boundary at all times. Implementation of the boundary constraint is not a straightforward task and some ALE developers have sidestepped its implementation completely by simply assuming Lagrangian boundaries and remeshing between load increments. Chapter 1. INTRODUCTION 10 The main drawback of the ALE technique lies in the fact that the number of unknowns is increased, being both the material and mesh displacements. In addition, ALE matrices are generally unsymmetric and require an unsymmetric equation solver. Fully coupled ALE computations are expected to be more time consuming than the corresponding Lagrangian ones especially for practical metal forming problems with large number of elements and complex contact conditions. This indicates the necessity of incorporating an efficient solution algorithm when using ALE. As indicated in Chapter 2, in a quasi-static ALE formulation, solution of equations is limited to implicit techniques in which iterations and convergence checks must be employed. However, using a dynamic ALE formulation, computations may be performed using implicit, explicit or mixed implicit-explicit calculations and solution is always much faster and easier to achieve. 1.3 SCOPE OF WORK The objective of this research is to develop and employ a true fully coupled ALE formulation in the simulation of large deformation metal forming problems. The formulation is aimed to be applicable to both quasi-static and high-speed metal forming processes. The scope of work may be summarized in the following: - Derivation of fully coupled ALE virtual work equations from the basic principles of continuum mechanics for both quasi-static and dynamic analyses. Using displacements as the independent variables, a consistent ALE formulation is sought Chapter 1. INTRODUCTION 11 that can be easily related to standard Lagrangian formulations for solid mechanics applications. Introduction of a new method for the treatment of convective terms in the equilibrium equations that avoids the resort to unjustified assumptions in calculating spatial gradients of stresses. This treatment aims to maintain the consistency and generality of the developed formulation. Discretization of the A L E virtual work equations using isoparametric finite elements. A n effort is required to put the different virtual work terms in an easy-to-compute matrix form consistent with standard nonlinear finite element calculations. Enhancement of the mesh motion scheme especially for boundary nodes. Boundary nodes should be allowed to move in the tangential direction to material boundaries to maintain a uniform distribution while satisfying the A L E boundary constraint. Implementation of the A L E formulation into a 2 - D computer code for plane stress, plane strain and axisymmetric problems. The implementation is intended to be modular and convenient for future developments. Implementation of efficient time integration schemes that allow for implicit, explicit and mixed implicit-explicit calculations. Implementation of a line search technique to accelerate the convergence of implicit calculations. Implementation of a simple contact algorithm to allow the simulation of metal forming problems involving workpiece-tool interactions. Chapter 1. INTRODUCTION 12 Application of the developed code in the simulation of several large deformation solid mechanics and metal forming problems. Comparison of results with experimental measurements or with other well established numerical techniques. Chapter 2 BACKGROUND 2.1 A L E LITERATURE REVIEW The concept of A L E was first proposed for fluid mechanics applications using the finite difference analysis [22, 23]. A L E was later introduced into the finite element analysis of fluid flow problems [24]. In this work, Hughes et al. elaborated on one of the basic concepts in A L E , which is the relationship between the material time derivative, the grid time derivative and the spatial gradients of any arbitrary quantity. This relationship is very important in deriving the governing A L E equations as indicated in Chapter 3. The A L E technique was also applied influid-structureinteraction problems to model the fluid domain while the structure domain was handled using the usual Lagrangian description [25]. The reason for the delay in the application of A L E to solid mechanics problems has been primarily due to the complexities of updating history dependent properties at integration points for arbitrary moving meshes. Huetink [26] introduced the first A L E formulation for the analysis of solid mechanics applications. Huetink used the convenience of the operator split technique in his finite element simulations of quasi-static metal forming problems. He proposed a method for updating variables at integration points. He first calculated the variables at each nodal point by averaging the integration point values from all elements that share this point. A continuous field for each variable is obtained by interpolating nodal point values using 13 Chapter 2. BACKGROUND 14 element shape functions. The gradients may then be computed using the derivatives of shape functions. This method was later refined by applying local and global smoothing of these variables to avoid numerical instabilities [27]. This technique is highlighted in Chapter 3. Schreurs and co-workers [28] discussed a similar A L E procedure. In this work, only fundamental ideas of A L E were discussed and formulation details were not given. The same authors [29] presented algorithms for the control of mesh quality, an important aspect of any A L E formulation. In controlling mesh quality, attention is mainly focussed on optimizing the element shape, while the topology of the element mesh is not changed. This is the reason that this process is often called 'mesh motion' as opposed to 'remeshing' which is commonly used in Lagrangian analyses and which usually involves optimization of element shape as well as element size. Haber [30] presented an uncommon form of the A L E formulation termed Eulerian Lagrangian Description (ELD). In this description, the total deformation in each increment is divided into separate Eulerian incremental displacements and Lagrangian incremental displacements. A n Eulerian deformation gradient defines the mapping from the initial configuration to the reference configuration while a Lagrangian deformation gradient describes the mapping from the reference configuration to the current configuration. The total deformation gradient from the initial to the current configurations is then the product of the Eulerian and Lagrangian deformation gradients. The E L D was successfully extended to include dynamic effects [31] and applied to analyze dynamic crack propagation problems [32]. However, the introduction of two sets of displacements Chapter 2. BACKGROUND 15 in this A L E formulation makes it difficult to relate to other formulations and difficult to implement into existing Lagrangian codes. Liu, Belytschko and Chang [33] presented an A L E formulation for path-dependent materials using explicit time integration. A new method for the treatment of the convective term was proposed. Although they discussed A L E equations for implicit calculations, their implementation was explicit. Explicit calculations do not require the decompositions and linearization involved in implicit calculations. The method of finding a continuous stress field by interpolation developed by Huetink [26] results in implicit calculations that are not compatible with the explicit time integration used. By defining a stress-velocity product, the computation of stress gradients was circumvented. However, this method requires interpolation of the stress-velocity product using shape functions. Thus, mixed finite elements, as opposed to the more common displacement based finite elements, must be used. This method is briefly described in Chapter 3. Liu et al. [34, 35] derived a fully coupled implicit quasi-static A L E formulation in rate form. The independent variables in this work are velocities as opposed to displacements. The stress-velocity product technique, previously developed by the same research group [33], was used to handle the convective effects. The resulting A L E equations are difficult to relate to the incremental form of the displacement based Lagrangian formulation commonly used in solid mechanics simulations. Benson [36] proposed a Simple A L E (SALE) formulation with the aim of reducing the cost of analysis and implementation into current Lagrangian codes. The S A L E formulation is an A L E formulation that is limited to a single material in each element as Chapter 2. BACKGROUND 16 well as using Lagrangian boundaries. He used the operator split technique with the central difference explicit time integration scheme. Huerta & Casadei [37] presented several A L E dynamic applications using explicit time integration. Since their calculations are explicit, they employed the same stressvelocity product technique developed by Liu et al. [33] to update stresses. They showed that A L E is a strong competitor to the classical Lagrangian formulation for fast-transient dynamic applications such as impact and coining. Extended from fluid mechanics, an implicit A L E formulation for solid mechanics problems was given by Ghosh and Kikuchi [38]. They indicated that the majority of A L E formulations available in the literature were based on explicit time integration. Explicit methods suffer from the lack of generality of application due to the stringent stability conditions, which usually necessitate the use of very small time steps. Consequently, an implicit A L E formulation was sought. In this work, however, it was assumed that the motion of the material points is quasi-static with respect to grid points, i.e. the time rate of change of the material velocity for a fixed grid point may be neglected. This simplified the derivation of the A L E equilibrium equations by neglecting terms that partially contribute to the material point acceleration. This assumption is quite common in steadystate analyses based on the Eulerian formulation in which all partial derivatives with respect to time are neglected. However, in A L E , the grid is not fixed and its motion is arbitrary. A solution for the A L E equations based on steady-state assumptions is physically meaningless if the grid velocity is to be arbitrarily modified. For quasi-static problems, it is the total material point acceleration that may be neglected. This Chapter 2. BACKGROUND 17 formulation is, therefore, incorrect and may not be warranted for quasi-static nor transient applications. An A L E formulation for quasi-static applications was developed by Wang [39]. Fully coupled A L E equations were derived in a fashion similar to that of Liu et al. [35]. Starting with the usual virtual work, integrations over spatial coordinates were first transformed into integrations over referential coordinates. Then by taking the rate of change of the virtual work expression with the reference coordinates held constant and transforming it back to the spatial coordinates, the rate form of the A L E expression of virtual work, i.e. virtual power, was obtained. Many of the basic procedures related to A L E , such as the mesh motion scheme and the need for unsymmetric solvers, were addressed in this work. A n in-house finite element code was developed based on this formulation. Several quasi-static metal forming problems were successfully simulated showing the effectiveness of A L E . However, the following comments could be made on this work: Although fully coupled A L E equations were derived, the implementation was not a strictly fully coupled one. In this work, stress calculation followed a Lagrangian procedure and the extra convective terms related to A L E were not handled within the iterations of each load increment. A separate stress updating routine was used to update material associated properties, such as stresses, strains, and current yield surface after iterations converge and before the start of a new load increment. The reason for this undertaking was mainly to avoid convergence problems. Since fully Chapter 2. BACKGROUND 18 coupled stiffness equations were solved, this formulation is a mix between the fully coupled and the operator split techniques. The treatment of convective terms in the fully coupled equations was based on the assumption of a continuous stress field by shape function interpolation developed by Huetink [26]. The same assumption was used to update stresses in the separate stress updating routine. Similar to the A L E formulations developed for fluid mechanics applications, the A L E equilibrium equations used in this work were derived from the principle of virtual power in which velocities are used as independent variables. This fact, together with assumptions used in the treatment of convective terms, resulted in A L E equations that can not be easily correlated with the standard Lagrangian formulation. The equations developed in this work were limited to quasi-static applications. The absence of inertia effects precludes the simulation of fast transient metal forming processes. In addition, quasi-static formulations necessitate the use of an implicit solver in which convergence problems are always expected. - The A L E boundary constraint was not implemented in this work. Within each load increment, boundary nodes were assumed to be pure Lagrangian. Although this treatment is theoretically correct, it does not ensure that boundary nodes are uniformly spaced on the boundaries at the end of each load increment. To overcome this problem, a remeshing routine was used to remesh material boundaries with new nodes between increments. Material properties for the new nodes were obtained from the old ones by interpolation. The remeshing of material boundaries between Chapter 2. BACKGROUND 19 increments is another reason for not considering this work a true fully coupled A L E implementation. The finite element code developed in this work had limited capabilities. Only twodimensional 4-node plane strain isoparametric elements were incorporated. No contact analysis capabilities were implemented. 2.2 DYNAMIC EFFECTS Dynamic effects play an important role in metal forming simulation. Including dynamic effects extends simulation capabilities to include fast transient metal forming processes. In addition, dynamic effects enhance the characteristics of the nonlinear finite element solution. 2.2.1 Fast transient metal forming The need for a higher production rate is a worldwide trend that gave rise to increased metal forming speeds [40-43]. Steinmann et al. [44] indicated that in sheet metal forming processes, high-speed deformations are very likely to occur. The interest is therefore to simulate the short-term transient response to these loading conditions and thus a solution strategy that can efficiently trace dynamic deformation behavior is required. Baillet et al. [15] developed a finite element program, based on the dynamic explicit method, for the simulation of the ironing process. It is shown that the higher the punch speed, the more the sheet lifts away from the die, thus changing the contact conditions and strain distribution. Similarly, in a study of the static and dynamic forming of circular Chapter 2. BACKGROUND 20 plates [45], it is shown that, while in static deformation the contact between the ram and the plate is maintained throughout the deformation history, in the dynamic case the contact between the plate and the ram changes. Fontane and Gelin [46] simulated several high speed metal forming processes, namely rolling, upsetting and wire drawing. In the rolling process simulation, the roll pressure was found to increase with the increase in rolling speed. In the wire drawing process, the optimal die angle and drawing stress increased with the increase in drawing speed. In the upsetting test, it is seen that the shape of the deformed mesh strongly depends on the upsetting speed. They concluded that the inclusion of dynamic effects in such calculations leads to modifications in the plastic flow and load distribution during forming processes and thus it is of the utmost importance. Kapinsiki studied the influence of punch velocity on the material deformation during deep drawing [47]. He noticed that punch velocity significantly influences the stress and strain distributions and that the quasi-static assumption causes large discrepancies between theoretical and experimental results. He concluded that it is necessary to take inertial effects into account to be in agreement with experimental observations. The above examples indicate that it is sometimes necessary to retain inertial terms in developing an A L E formulation for general metal forming simulation. 2.2.2 Solution of nonlinear equilibrium equations in finite elements A common solution approach for the nonlinear finite element equilibrium equations is the step-by-step incremental procedure in which it is assumed that solution for all Chapter 2. BACKGROUND 21 equilibrium positions for all time steps from time 0 to time t have been obtained, and the solution for the next equilibrium position at time t + At is required. This process is applied repetitively until the complete solution path has been solved for. A n implicit solution scheme is adopted if the solution for equilibrium at time t + At is based on using equilibrium conditions at time t + At, whereas an explicit scheme is used if equilibrium at time t + At is based on using equilibrium conditions at time t. In quasi-static analyses, solution of equilibrium equations is always based on the implicit scheme. On the other hand, it is possible to employ an implicit, an explicit or a mixed implicit-explicit solution scheme for dynamic problems. Using an implicit solution scheme, it is necessary to employ iterations and check convergence for every time step. Bathe [48] indicated that in nonlinear static analysis, iteration might not converge under certain loading conditions, whereas convergence is always achieved in dynamic analysis provided that the time step is sufficiently small. In addition, for a dynamic analysis the inertia of the system renders its dynamic response more smooth than its static one, and convergence is in general more rapid. The same observation was confirmed by Choudhry and Lee [41]. They indicated that in metal forming processes, like sheet-metal bending and frictionless deep-drawing, large rigidbody motions might be involved, leading to convergence problems. Contributions from the dynamic terms help to stabilize the computations. As an example for a sheet-metal bending process with large rotations and relatively small strains, a flat sheet of 1 mm thickness is deformed with a cylindrical punch into a channel. The sheet deformation is primarily due to bending in the regions of contact with punch and die. The maximum Chapter 2. 22 BACKGROUND punch displacement attained in modelling the bending process quasi-statically was 35 mm. However, a displacement of 53 mm was attainable when the dynamic effects were considered. They concluded that it was important to retain dynamic effects in the simulation of this problem. This problem is discussed in more details in Chapter 6. The solution of dynamic equations of motion using explicit time integration simply corresponds to a forward marching in time without iterations. Explicit integration schemes do not require a factorization of the global stiffness matrix. Moreover, if the global mass matrix is diagonal, no matrix factorization is required at all. In this case it is not necessary to assemble the global mass and stiffness matrices and the solution may be carried out on the element level and a little high-speed storage is achieved. The main shortcoming of this scheme is that it is conditionally stable. The time step At has to be smaller than a critical value At . For typical sheet metal forming operations, a time step cr in the order of one millisecond is normally sufficient [10, 49]. This may impose a stringent limitation on the efficiency of this method for the simulation of slow metal forming processes which occur in several seconds. Comparative investigations into implicit, explicit, and iterative implicit-explicit schemes for the simulation of sheet metal forming processes have been the topic for extensive recent research [49-55]. It is shown that the difficulties associated with the convergence of the implicit methods have led to a renewed interest in the use of dynamic explicit methods, even when the forming problem is essentially quasi-static [54]. These difficulties are more pronounced in the case of complex problems with large number of elements, large rigid body motions and variable contact conditions. Chapter 2. BACKGROUND 23 Although the implicit approach employs a more reliable and rigorous scheme in considering equilibrium at each step, convergence problems may be expected. If the explicit approach is used for all time steps, computation time for slow metal forming processes with fine meshes will become excessive. Thus, it might be more efficient to shift between the two techniques during the deformation process. Therefore, implicit, explicit and mixed implicit-explicit solution techniques will be implemented in the current work. 2.3 SUMMARY Survey of the literature on A L E leads to the following conclusions: Almost all of the previous A L E analyses are based on the easy-to-implement operator split technique as opposed to the theoretically more accurate fully coupled approach. Most of the previous A L E formulations where developed in a rate form similar to the Eulerian formulation, with velocities as independent variables, and can not be easily related to the Lagrangian formulation normally used for solid mechanics applications. - A l l dynamic A L E formulations available in the literature are based on the explicit solution scheme, which suffer from the lack of generality of application due to the stringent stability conditions and which usually necessitate the use of very small time steps. - No attempt has been made to develop a fully coupled A L E formulation for both quasi-static and dynamic solid mechanics applications based on the more reliable implicit solution scheme and consistent with the standard Lagrangian formulation. Chapter 2. BACKGROUND 24 The two common treatments of convective terms, namely the assumption of a continuous stress field and the use of the stress-velocity product, introduce some sort of approximations and limitations. Chapter 3 DERIVATION OF ALE GOVERNING EQUATIONS 3.1 PRELIMINARIES An implicit time-stepping approach will be used in deriving the governing A L E virtual work equilibrium equations from the basic principles of continuum mechanics [56, 57]. In this approach, we assume that the solutions for all time steps from time 0 to time t have been solved, and that the solution for time t + At is required next. This process is applied repetitively until the complete solution path has been solved for. Linearized virtual work equations for quasi-static and dynamic applications are derived. A new method for handling convective terms in the equilibrium equations is presented [58]. 3.1.1 Notations Throughout the derivation, standard indicial notations are adopted; right subscripts denote the components of a tensor and repeated subscripts imply summation. In addition, time and configuration notations similar to those used by Bathe [4] are adopted. Left superscripts indicate the configuration in which the quantity occurs whereas left subscripts indicate the configuration to which the quantity is referred. Left subscripts may be omitted be if the quantity occurs in the same configuration in which it is measured. A quantity with no left superscripts or subscripts indicates an incremental quantity from time t to t + At. 25 Chapter 3. DERIVATION OF ALE GOVERNING EQUATIONS 3.1.2 Kinematics In the A L E description, the material configuration at any time t refers to the set of material particles, whereas the reference configuration consists of a set of arbitrarily moving grid points sharing a common boundary with the set of material particles. The material configuration is identified by a set of material point coordinates Xf while the reference, or grid, configuration is identified by an independent set of grid point coordinates Xf. Let xf{XJ,t) ! and xf(X ,t) l 8 be the vector functions or the mappings that characterize the motion of the material point XJ and the grid point X 8 in space, respectively. The position of XJ at time t is given by < ='xJ(XJ,t) (3.1) Xi The set of material particles is related to the set of grid points by requiring that the two configurations share the same space at all times. Any point within the common boundary is occupied by elements of the two sets. Thus, the position of the grid point X 8 that occupies the same point in space at time t as XJ is also given by 'x as t VWJ.O (3-2) The A L E formulation requires that the inverse of (3.1) and (3.2) exist to ensure a oneto-one mapping between the two configurations. The material velocity 'v . and the grid ( point velocity 'vf at time t are given by 'v.. dx" a (3.3) Chapter 3. DERIVATION OF ALE GOVERNING EQUATIONS v° (3.4) a The boundary constraint, which ensures that the material and grid configurations have the same boundary at all times, can be expressed as ('v.-'vf)'n,. on the boundary =0 (3.5) where n is the unit normal to the boundary surface. 1 i The governing A L E equations involve the material time derivative of several quantities. The material derivative of an arbitrary function ' / is denoted by a superposed dot and is defined to be the rate of change of the function holding the material particle XT fixed df (3.6) a However, the grid configuration is the computational configuration that tracks the history of all quantities. Thus, it is convenient to define a grid time derivative, which is the time derivative of the function ' / holding the grid point Xf fixed, and denote it by a superposed prime (3.7) a The relation between the two time derivatives is given by [24] 7='/'-KV'v,')^f ox, (3.8) Chapter 3. DERIVATION OF ALE GOVERNING EQUATIONS The A L E formulation will be discretized using the isoparametric displacement based finite element method. We denote the incremental material displacements from time t to time t + At by w. and the corresponding incremental grid displacements by uf. We have ( the following relation ' + A ' + u f (3.9) where ' 'x,. is the position of the grid point in the configuration at time t + At. +A 3.1.3 Continuity The local form of conservation of mass, continuity, at time t is given by d'v 'P=-'P^T- (- ) 3 10 ox t where 'p is the material density. Using (3.8), the continuity equation with respect to an arbitrary moving grid point can be expressed as y=-v>^-(V'vf)|^ OX: 3.2 c- ) 3 11 OX; QUASI-STATIC ANALYSIS For quasi-static analysis, such as in the case of low-speed metal forming processes, inertia effects may be neglected. Employing an implicit incremental approach, the governing equilibrium equations for A L E are established for the configuration at time t + At. Since the configuration of the body at time t + At is yet unknown, an approximate Chapter 3. DERIVATION OF ALE GOVERNING EQUATIONS solution can be obtained by referring all variables to the grid configuration at time t and linearizing the resulting equations. The solution is then refined by iterations. 3.2.1 Principle of virtual displacements Since the earliest A L E formulations were developed for fluid mechanics applications, the principle of virtual power, with velocities as independent variables, found wide appeal. A L E researchers [39, 59-63] continued to use the principle of virtual power for solid mechanics applications in which displacements are more natural to use. In this work, a procedure analogous to that used by Bathe [4] in obtaining the virtual work equations for the Updated Lagrangian formulation is chosen. Thus the principle of virtual displacements is employed to express the equilibrium of the body at time t + At. It can be written in the form J* ' ^ijS js,j ' 'dV =S 'W +A, +A l+A (3.12) exl l+/ ,+A, v where ' 'o" - are the components of the Cauchy stress tensor at time t + At and e +A ;> t+Al y is the conjugate strain tensor defined by , i iL A.) = ( + (3 13) The external virtual work, S' 'W , is given by +A ext j'^' ^Su dV+ l+At S ' + * f V e c t = p t+Al y in which t+Al ff and '^'ff i j'^ffSu^dS (3.14) t+Al g are the components of the body force per unit mass and the Chapter 3. DERIVATION OF ALE GOVERNING EQUATIONS surface traction at time t + At, respectively. 3.2.2 Incremental decompositions In referring variables to the grid configuration, variables at time t + At are assumed to be composed of their respective values at time t plus an increment given by the grid time derivative of the variable multiplied by the time increment Ar. This is in contrast to the Updated Lagrangian formulation where the incremental quantities are given by the material time derivative of the variable multiplied by the time increment At. Material density at time t + At can be decomposed into 'p='p+'p'At (3.15) ,+A which, upon substituting with (3.11), gives t+M t t„f\_ /-„ ,.g\ d'p P=p-pj^-(u -ui)^n k dx dx k (3.16) k Stress components at time t + At can be expressed in terms of the stresses at time t for the same grid point plus a stress increment o^.Ar ^a^^+'^At (3.17) and using (3.8), we get ' + A V,='a,+'cT,A -K-uft-£f (3.18) The material rate of Cauchy stresses '& is calculated from the material constitutive y relation which is usually given in terms of an objective stress rate tensor such as the Chapter 3. DERIVATION OF ALE GOVERNING EQUATIONS Truesdell stress rate tensor defined by < a 9 r = ,J J ^ , d'x d\ l J ± , dx & ' ,J k k ( ] k 3 y A 9 ) } k The material constitutive relation in terms of the Truesdell stress rate is given by '<='<VA, (3-20) where 'Dy is the rate of deformation tensor given by 1 d'v. d'v, D»=-( r-+—r ) 2 d'xj d% £ (3-21) L v and 'C ijkl is the fourth order material constitutive tensor. It should be noted that the above stress and strain measures were chosen to allow the A L E formulation to be effective in large strain situations. The variation in the strain components at time t + At can be decomposed as =Sfi„ + £ , 4 ' (3-22) A in which 8 p'y is the grid time derivative of 8 p given by [Appendix A] i} *4=-I(ftfi f ! ^ ) (3,3) + 2 dx k aXj dx k dx { Substitution in (3.22) gives 8 e =8e < W „ *Pv d'x { 2 k d' Xj d'x d'x/ K (3 24) } k Incremental decomposition of elemental volume at time t + At in terms of the elemental volume at time t is given by [64] Chapter 3. DERIVATION OF ALE GOVERNING EQUATIONS 'dV='dV+'dV'At = (\ + ^-)'dV t+A (3.25) Similarly, incremental decomposition of elemental surface area is given by [64] 'dS='dS+'dS'At = [l + ^--(^f+ -^-)'n JnJdS d'x 2 c7'x„ &xj y (3.26) k where '«„, is the unit outward normal to the surface at time t. 3.2.3 Linearization Linearization is accomplished by expanding (3.12) using the previous incremental decompositions and neglecting higher orders in all incremental quantities. Substituting by (3.18), (3.24) and (3.25), the internal virtual work can be linearized as l + My ly ly 'y ' \ ^ ' - / p d V •v i - \ k dx u { k X - u t ^ ^ d V -v dx ^k d x (3.27) k Considering the external virtual work on the RHS of (3.12), the body force term can be referred to the grid configuration by using (3.16) and (3.25) to get ,+toy ly ly ^Xj. 3 -l^f "(u -u' )^-Su dV (3.28) , l l l l <V k X Similarly, using (3.26), the traction force term of the external virtual work is expressed as {'^ffSur'dS = \' 'f?[l +& • +^ - U ^ d S X k 2 S x „ + ^) nJnJSu<dS ! d X m (3.29) Chapter 3. DERIVATION OF ALE GOVERNING EQUATIONS 3.2A Treatment of convective terms Convective terms, such as the last integral on the RHS of (3.27), involve the computation of the spatial derivatives of stresses. Since the stress values are computed at the integration points, and not at the nodal points, the stress field is generally discontinuous across element edges. Thus, stress gradients may not be reliably computed on the element level when evaluating element matrices. A method of finding a continuous stress field by interpolation was developed by Huetink [26] and used in the Eulerian step of the operator split technique. Some researchers used the same method to handle convective terms in the coupled equilibrium equations [39]. In this method, integration point stresses are first extrapolated by a least square approximation to get the nodal stresses. Nodal stresses computed from each element are then averaged. A continuous stress field is then assumed in the form '^=ZV^« (3-30) a=\ where h a is the element shape function evaluated at node a, a ija are the nodal stress components at node a and N is the number of element nodal points. Finally, the spatial derivatives of stresses can be computed on the element level as ax k a = i ax k This method is very popular in the A L E literature despite the least square approximation and assumption of a continuous stress field. Another method for treating convective terms was proposed by Liu et al. [33]. A stress-velocity product is defined in the form Chapter 3. DERIVATION OF ALE GOVERNING 34 EQUATIONS KV'v/)'^ % (3.32) Differentiating (3.32) with respect to space gives the convective term as (VX)^ =^ - ( ^ - ^ ) V „ ax dx k dx k (3.33) dx k k Equation (3.33) circumvents the computation of the stress gradients by computing the gradients of the stress-velocity product instead. However, this method also requires interpolation of the stress-velocity product using shape functions in a fashion similar to (3.30) and thus mixed finite elements, rather than the more common displacement based finite elements, must be used. In this work, a new method for the treatment of the convective term that sidesteps the computation of the spatial gradients of stresses is used. Based on fundamental A L E relations, the new method offers an accurate treatment of convective terms while maintaining the convenience of using displacement based finite elements. This method involves a transformation from volume integrals to surface integrals as offered by the divergence theorem. Use is also made of the boundary constraint in (3.5). The last integral on the RHS of (3.27) can be rewritten as •v vx k , ox v , k v - k j t - l t ' ' < , vx v x y v k s » ' i v ox k ( 3 3 4 ) k No new assumptions are introduced in the calculation of the second and third integrals on the RHS of (3.34) since the spatial derivatives of displacements and strains can be computed on the element level from shape functions derivatives. Applying the divergence Chapter 3. DERIVATION OF ALE GOVERNING EQUATIONS theorem to the first integral on the RHS of (3.34) and using (3.5), we get f ^ - f ' ^ ^ ' d V = J(« -uD'a^WdS °k 'V =0 4 (3.35) 's x Substituting in (3.27), the internal virtual work becomes i'^S^ dV=j'a S^dV+ [<J,At5^dV M v iy + 'y dV j C X 'y -\%'^jt' 'V fas^p-'dV 'k + ^ - < ' ^ ' °k ) 'v k M < 3 3 6 ) X The same method can be applied to the convective body force term to avoid the computation of the spatial derivatives of density. The last integral on the RHS of (3.28) can be treated in the same manner as in (3.34) and (3.35), to give ~\t+AtrB \' p 'f Su ' dV +A, ,+A B i j' ?-Jj-( - *)Su.'dV = \ p % 8u 'dV+ +At t i ,+ B p i l+My ly 'y " Uk u k X + j'A'(« -«;)j ^ i i <-) 33 7 3.2.5 Fully coupled A L E equilibrium equation Substituting (3.36), (3.37) and (3.29) into (3.12), the principle of virtual displacements at time t + At referred to time t can be written as [cr At8^dV+ v •y 'y =S W 'l+M where ex [a.Se'dV 'V faS^'dV k OX (3-38) Chapter 3. DERIVATION OF ALE GOVERNING EQUATIONS p<t+&trB +^^(u -u?)]8u 'dV+ k 8"*W* = J >pr% >v B + ['^f^ , J S k refill f K i dx +^ - h ^ d'x 2 dx k n >p' % (u -<y^dV k + B k dx + ^ynJnJSu/dS d'x "' (3.39) m The constitutive relations in (3.19) to (3.21) can now be introduced into the first integral in (3.38) to give i'C e 8^dV+\'a 8 ri 'dV 'V 'V uklt I kl u " k t 0 " dx k );d' d'x/ Xj d'x ,J k =8 W -j'a 8f dV 'V ,+At ext (3.40) t iJ iJ where l l , J 2d%d'x y j ' Equation (3.40) is linear in the incremental displacements u and uf. Thefirsttwo l integrals on the LHS of (3.40) correspond to the Lagrangian material and geometric stiffness matrices and are exactly the same as those obtained using an Updated Lagrangian formulation. The last two integrals on the LHS establish the convective stiffness matrices due to A L E . The last term on the RHS of (3.40) corresponds to the Lagrangian internal force vector. Equation (3.40) represents the fully coupled A L E equilibrium equation. This equation can reduce to the Updated Lagrangian formulation if we choose to attach the grid to the material, i.e. uf =u n and to the Eulerian formulation if we choose to fix the grid in space, i.e. uf = 0, as limiting cases. The A L E equilibrium equation derived in this work Chapter 3. DERIVATION OF ALE GOVERNING EQUATIONS shows that A L E can be considered as a logical extension to the Lagrangian formulation and the modifications to the equilibrium equation of current Updated Lagrangian codes are clearly identified. 3.3 DYNAMIC ANALYSIS In dynamic analyses, inertia effects are included in the balance of momentum at time t + At. Inertia forces involve the material time derivative of material velocities, i.e. material accelerations 'v . In the A L E formulation, we follow the grid point in its ,+A j motion as our reference configuration. Therefore, the referential material acceleration, which is the grid time derivative of the material velocity clarity, we will denote the referential material acceleration 'v., should be used. For l+A v' by ' 'a,.. l+At +A 3.3.1 Virtual work done by inertia forces Using the relation between the two time derivatives in (3.8), the virtual work done by inertia forces can be expanded as J" / + A y + A ( v>, dv ,+A, = j p \su, t+At ,+ dv t+A, + j ^' C ' -^)^^Su A P r+My v O t + i A l dV (3.42) Xj Equation (3.42) is considered as an extra virtual work term due to inertia effects to be added to the LHS, or subtracted from the RHS, of the A L E virtual work equation, (3.40). Chapter 3. DERIVATION OF ALE GOVERNING EQUATIONS The first term on the RHS of (3.42) can be referred to as the referential inertia term whereas the second term is referred to as the convective inertia term. 3.3.2 Decomposition of velocities and accelerations The velocities and accelerations at time t + At can be related to their respective values at time t using the relations ' 'a,='a.+a,. (3.43) ' 'v,.='v,.+v,. (3.44) vf='vf+vf (3.45) +A +A ,+t, The incremental quantities a , v. and v f depend on the implicit time integration scheme i ( to be employed. To maintain a linear increment, higher orders in a , v ., v? and A? will i ( be neglected. 3.3.3 Linearization of the referential inertia term Incremental decomposition of the variables in the referential inertia term, in a manner similar to quasi-static analysis, and linearization of the result gives j '^'p'^'a.Su/^dV = j" 'p \ou/dV l+ - J - { ^ ( ^ - ' v ^ ' a ^ A t ' d V •r k dx The last integral on the RHS of (3.46) can be rewritten as 'pi^-^-Y^'a^At'dV (3-46) Chapter 3. DERIVATION OF ALE GOVERNING EQUATIONS • SX ,| V d'x k •v d x d* k k -\'p(v -'vl)^^At'dV (3.47) k •v d x k Applying the divergence theorem to the first integral above and using the boundary constraint in (3.5), we get jS['p(\ <V 'vp'^'^te'dv x k = j"'p(' -'v*y "a Su n At'dV , + Vk l l k =0 (3.48) s Substituting in (3.46), we get J '^'p'^Su'^'dV = J 'p'^Su/dV i+toy + | 'p(\-'vf) ' '? di +A >y 'y At'dV idUi) " (3.49) k X Using (3.43), the referential inertia term can be written as | ' V 'a,<V '^ = J 'p'aM'dV + J + l + Cly +A +A ly 'pa^'dV ly + \'p(\- v!)^^-At'dV ,y d'X (3.50) t V 3.3.4 Linearization of the convective inertia term The convective inertia term can be expanded into Chapter 3. DERIVATION OF ALE GOVERNING f , + J A p C l ^ A t v v J v EQUATIONS - L S u . , + A , d v r)' Y J +Al y V = l'pC 'v-'+ >v )^Su<dV A A j •v d { ^d'x x J &x k } ) k ,yOX ffxj 1 OXj k - j'p^ "vj- +"v>)^Su£t'dr •y vj v + (3.51) , x k X As before the third integral above, which involves the spatial gradients of density, can be treated using the divergence theorem and the boundary constraint, to get j l + Af ' pc\- v*)^dur'dv +A, ,+A, y O = Xj 'pC^-^D—^ou/dv d Xj dX, fy d\ ly d'Xj O X^ Cst u Xj g - j ' p ^ r ' v - ' ^ ^ S u A t ' d V iy & j X 0 k X Using (3.44) and (3.45), the convective inertia term can be written as (3.52) Chapter 3. DERIVATION OF ALE GOVERNING EQUATIONS 3.3.5 Fully coupled A L E equation of motion Combining (3.40), (3.50) and (3.53), the A L E equation of motion can be written as + J 'C e S^dV 'V + J 'pafr'dV 'V + J pi'v -'^)p-&i, dV+ , , j'pivj J 'V 'a.5^dV +J 'V kl mt -vDp-Su/dV , OXj y + hu -ut)>o- ^pL<dV k IJ U^-^yof-^'dV + =S' W - J 'cjySfr/dV - j" 'p'ofa'dV +Al exl •v 'V -JV(\-^)^^A^F-JV(V,-VJ)|^^,.VF <y ~ Ok <y X Vj X r, , , d'v, d'v" d'v. f'p^-'vt^-l-^^du^t'dV K dx dx dxj k k - [ pCv.-'vftv-'v^^ip-du^t'dV •y OX k O Xj (3.54) Chapter 3. DERIVATION OF ALE GOVERNING EQUATIONS Equation (3.54) is a linear equation in the incremental displacements «. and uf, incremental velocities v,. and vf, and incremental accelerations a and it represents the t fully coupled A L E equation of motion. The first three terms on both the LHS and RHS of (3.54) are exactly the same as in the Updated Lagrangian formulation. The first term on the LHS corresponds to the Lagrangian mass matrix. The second, third, sixth and seventh terms on the LHS and the first and second on the RHS were defined for quasi-static analysis. The fourth and fifth terms give rise to convective velocity stiffness matrices due to A L E . The third term on the RHS of (3.54) corresponds to the Lagrangian inertia force vector whereas the last four terms on the RHS are convective inertia force vectors due to ALE. The fully coupled A L E equation of motion derived in this work is consistent with the standard Lagrangian formulation. No assumptions or approximations, other than those for the purposes of linearization, have been made in the derivation. Incremental velocities and accelerations are retained in the equations for later time integration. Chapter 4 FINITE E L E M E N T DISCRETIZATION 4.1 I S O P A R A M E T R I C FINITE E L E M E N T S In this chapter, the virtual work equations for both quasi-static analysis and dynamic analysis will be discretized using two-dimensional isoparametric finite elements. Details of computer implementation of the finite element matrix equations are also discussed [65]. In the isoparametric finite element discretization, element coordinates 'x and j incremental displacements w,. and uf are interpolated using '*,=2>*% (4-1) k=i N ",=Z^A (4-2) k=\ uf=f h uf j k (4.3) k k=\ where 'x , ik u jk and uf are the nodal coordinates and incremental material and grid k displacements of degree of freedom / of nodal point k at time t, h is the element shape k function corresponding to nodal point k, and N is the number of element nodal points. In two-dimensions, equation (4.1) can be expanded as 43 Chapter 4. FINITE ELEMENT 44 DISCRETIZATION (4.4) In matrix form o K 0 y (4.5) h,. or (4.6) x=H'x where ' x is the element coordinate vector given by (4.7) x= y H is the element interpolation matrix in the form H = 0 (4.8) K 2x2 N and ' x is the nodal coordinate vector at time t given by •Y-{.. y k •••};1x2 N (4.9) Similarly, equations (4.2) and (4.3) can be written in the form (4.10) u =Hu u g where u , u , u and u g g = HiT ? (4.11) are the incremental element and nodal material and grid displacement vectors, respectively. We also have Su = RSu (4.12) Chapter 4. FINITE ELEMENT DISCRETIZATION 45 where c5u is the virtual nodal displacement vector. 4.2 DISCRETIZATION OF THE QUASI-STATIC A L E EQUATION The A L E virtual work equilibrium equation for quasi-static analysis was derived in Chapter 3 in the form 'V I 'V k k ' S'x J k l^'xj d'x/ 3'x ,J k =8 *W -\<cj 8 e ;dV t+ (4.13) ex, ij t i 'V As indicated before, the first and second terms on the LHS of (4.13) correspond to the Lagrangian material and geometric stiffness matrices, respectively, whereas the third and fourth terms establish the convective stiffness matrices due to A L E . The last term on the RHS of (4.13) corresponds to the Lagrangian internal force vector. 4.2.1 Discretization of the Lagrangian internal force term The internal force vector obtained in A L E is the same as in a Lagrangian analysis. Thus it will be discretized exactly in the same manner [48]. We start by defining a stress vector in the form cr yy We also define the incremental strain vector in the form (4.14) Chapter 4. FINITE ELEMENT XX y e t < d'x du e t (4.15) yy 2 ,e t ^ t where CT^ and e zz 46 DISCRETIZATION •+- xy d'y e d'x are the hoop stress and incremental strain in the case of axisymmetric problems. Using (4.2) fdh _ t xk v fdh _ t d'x 0 •Uyk _ dh _ , ,e = 4^M > (—— u , H k N I 4=1 k K N — — uA d'y K _ uxk v 0 dK d'y dr\ d'x *xk (4.16) 'yk 0 7=1 or (4.17) ,e=B-« where d\ d'x 0 B" d\ d'y h, 0 d\ d'y d\ d'x (4.18) Chapter 4. FINITE ELEMENT DISCRETIZATION Al We also have SF=,B 8a (4.19) Ll Using (4.14) and (4.19), the internal force term can thus be written as l'<T 5 i dV=l(Sfi) '* dV , 9 T , f 0 •v •v = j( B Su) 'a'dV 'V Ll T t = (du) ^(^f'o'dV T 'v = (<5u) 'f (4.20) r where ' f is the internal force vector given by 'f =\( B ) 'a'dV u (4.21) T l 4.2.2 Discretization of the Lagrangian material stiffness term The Lagrangian material stiffness virtual work term can be rewritten in matrix form as \'C e S dV 'V = l(S,e) 'C e<dV 'V ! ijkll kl T a l = JLB'^S&y'CCB u)'tiF 'V L 1 = (c5u) J(,B ) 'C,B "(/Fti 'V r il r i = (<5u) 'K u r (4.22) £, where ' C is the elastic-plastic material constitutive matrix and 'K Ll is the incremental material stiffness matrix given by [48] 'K 1 1 = \{,B ) C B dV LX Tt LU t (4.23) Chapter 4. FINITE ELEMENT DISCRETIZATION 48 4.2.3 Discretization of the Lagrangian geometric stiffness term The Lagrangian geometric stiffness virtual work term can be rewritten as • f W ^ - f ' a ,„ | ^ r f (4.24) d'xj d'x, By expanding the summations over i,j and k, we get r, _ ,, r„ du d5u . du dSu . du dSu . du ddu a 5,n. dV - ( cr — -+a — -+cr — -+ a — 'V 'V " d'x d'x J " 8'y 8'x d'x d'y d'y d'y du., ddu. du„y ddu, du„ ddu. du„ ddu, y .i_ y y .(_ y y +'CT -+ cT„. —: : 1- cT —• •y .h/ _<7„ y d'x e'x d'y d'x d'x d'y d'y d'y Tr 11 r r r 1 r r xy r yx r r yy (4.25) x x Using (4.10) and (4.12), we get \'a d r 'dV iJ t Jo (dn) $( B ) 'S B 'dVu 'v = L2 T t 'v T L2 L2 t (4.26) = (<5uy'K u i2 where d_\ d'x d_\ 0 0 d'y B L2 = 0 0 h, d'x d'y (4.27) Chapter 4. FINITE ELEMENT DISCRETIZATION 'a XV 'S L2 = 49 0 0 0 0 0 0 0 \y 0 0 0 '** 0 0 0 (4.28) 0 xx ><T 0 0 5x5 and the incremental nonlinear geometric stiffness matrix 'K 12 is given by [48] (4.29) ' K " = J(,B") 'S",B"'</K r 4.2.4 Discretization of the first convective stiffness term The first convective stiffness virtual work term due to A L E can be expanded as , •y ^k X d Su 3L.2 yyd'xd'y st.. at. Xit d 8u 2 t . „, , r ,d 5u 2 t ox d 5u x ox 2 x d Su 2 x y l ddu 2 x dxdy , ddu 2 x x 'V , d Su _ , , 8u . xX x ddu l ddu ,, T (4.30) It should be noted that this term involves second derivatives of displacements and consequently second derivatives of shape functions. Substituting by (4.10) and (4.12) and defining the matrices Chapter 4. FINITE ELEMENT 50 DISCRETIZATION d'x d% 2 d'y 2 _d\_ d'xd'y 0 d'x 8% 2 0 (4.31) d'y 2 d%_ d'xd'y 1 ( Sh h k k xd'x 'x K, 2) \_dr\ 'x d'y J8x2W and °x: 0 xy xy cr yy 0 *y cr (4.32) 2x8 we get J ( « * - u ! ) ' a ^ ' d V = (<5u) •y VX , •y l(,B f'S n'dV(u-u») r Al y M y k = (bv) 'K \vi-vi ) T A (4.33) g where the first convective stiffness matrix due to A L E is given by >K = A1 ji^Y's^u'dv (4.34) 4.2.5 Discretization of the second convective stiffness term The second convective stiffness virtual work term due to A L E can be expanded as Chapter 4. FINITE ELEMENT DISCRETIZATION ,{X i^'Xj ¥ d' dx ) d'x j Xj J J d'x VK k - d'x K v ¥ « 51 .du + (— dy yx s x V t x A xy yy x dx + y yx dSu^ v d'y st. . st.. > st.. y d'y d'y 1 > du ., dSu , ddu -)( a -+a -) d'y d'x d'x ' dSu , dSu ) d'x dy " dy x v d'x dSu^ c yy st.. ' ~d'y yy ^(u -u' )'a &i ydV 7 x x B x (4.35) Substituting by (4.10) and (4.12) and defining the matrix ! dh,, d'x o : 0 d'x d\ o ! d'y B = A1 0 d'y I h„ I (4.36) ! _ 0 | N 5x2 N we get j ( | - ^ ) ' ^ ^ ' ^ = W fCB^ ) 'S"B"'rfK(ll-u') i• dx, dx, dx f " "j " "j t v r k 2 r v = (6u) 'K (u-u ) T A2 g (4.37) where the second convective stiffness matrix due to A L E is given by >K A2 = \( B ) 'S B 'dV A2 t T L2 L2 (4.38) Chapter 4. FINITE ELEMENT 52 DISCRETIZATION 4.2.6 Discretized Quasi-Static A L E Equation Substituting (4.20), (4.22), (4.26), (4.33) and (4.37) into (4.13), we get, for a single element or a group of elements ('K -i-'K")u-i-('K +'K ' )(u-li*)=' 'f '-'f (4.39) 'K^u+'K^u-uO^'f^-'f (4.40) 'K ='K +'K (4.41) il >,1 / 2 +A er or where L L 1 i 2 'K ='K +'K A 'K M A2 (4.42) and 'IC* are the Lagrangian and A L E stiffness matrices, respectively. Equation L (4.40) represents the finite element matrix equation for quasi-static A L E analysis. 4.3 DISCRETIZATION OF THE DYNAMIC A L E EQUATION The A L E virtual work equation of motion for dynamic analysis was derived in Chapter 3 in the form Chapter 4. FINITE ELEMENT + j 'pafa'dV 'V DISCRETIZATION kl +\ P( v 'v^au/dV , + =S' W ' - J 'o-ySft/dV - I ex 'V <y ~ 'cj.S^dV lp(v -v^Su/dV t r +M +J 'V + \ 'C e S^'dV 'V ukh 53 J 'p'afa'dV 'V V k <v X V j X 'p(v -'vl)(-^-2-^)^bu£t'dV dx dx 8xj k f r k k -J 'pCv.-'vOCv-'v^J-ipLSu^At'dV • o x O Xj v (4.43) k As indicated in Chapter 3, the first term on the LHS of (4.43) correspond to the Lagrangian mass matrix, whereas the fourth and fifth terms give rise to convective velocity stiffness matrices due to A L E . The third term on the RHS of (4.43) corresponds to the Lagrangian inertia force vector whereas the last four terms on the RHS are convective inertia force vectors due to A L E . The discretization of velocities and accelerations is similar to that of displacements and coordinates given by (4.6), (4.10) and (4.11). Thus we can write, v=H'v (4.44) 'v =H'v ? \ g (4.45) a=H'a (4.46) v = Hv (4.47) 8 = Hv g (4.48) Chapter 4. FINITE ELEMENT DISCRETIZATION 54 a = Ha where ' \ , ' \ 8 (4.49) and 'a are the nodal material velocity, grid velocity and referential material acceleration at time t whereas v, v and a are the nodal incremental quantities g from time t to t + At. 4.3.1 Discretization of the incremental mass term The incremental Lagrangian mass term can be discretized using (4.49) to give J 'pciidu/dV = J 'V 'V 'p(HSu) (Ha) dV T = (Su) l J'pH H'JFa 'V = (Suy IVPa T r (4.50) where ' M ' is the Lagrangian mass matrix given by 'M L (4.51) = j'pB. B.'dV •v T 4.3.2 Discretization of the first convective velocity-stiffness term The first convective velocity-stiffness virtual work term can be discretized using (4.44), (4.45) and (4.47) and expressed in the form [pCv-'v^^bu'dV 8'xj •v where 'C Al = (Su) 'C ^ T A (4.52) is given by , A\ C j I^A\ tpA\ \ j 2/-l,2y-l ^2i-\,2j i i ts~<A\ *-'2/,2y-l -'2i,2 j | C = (4.53) ( 2Nx2N Chapter 4. FINITE ELEMENT DISCRETIZATION 55 in which / and j indicate node numbers from 1 to N, and ? dh dh N K •y O s XK ~t XK' ° k=\ X I s-tA\ _t(~iA\ N y k=\ _ r\ (4.55) 4.3.3 Discretization of the second convective velocity-stiffness term The second convective velocity-stiffness virtual work term is discretized using (4.44), (4.47) and (4.48) and expressed in the form \'p(Vj-v])p-to/dV <v j = m T t C \v-v ) A g (4.56) dx where C is given by 'C A2 *-'2i-l,2 v-l = C L 2 ; , 2 7-1 2H,2 j (4.57) *^2(,2 j 2Nx2N in which i and j indicate node numbers from 1 to N, and (4.58) 'V iy k=\ (4.59) ffy (4.60) •y CZ t J = k=l U X \'ph h±^v k dV k=\ °y t i i v y (4.61) Chapter 4. FINITE ELEMENT DISCRETIZATION 56 4.3.4 Discretization of the inertia force term The inertia force virtual work term is discretized using (4.46) to give l p'a du/dV l'p(HSu) (a'a)'dV = l i •v T 'V = (Su) J'pH H'JF~'a T r 'v (4.62) = (<5u) 'M "a r where ' M L is the mass matrix given by (4.51). 4.3.5 Discretization of the first convective inertia force term The first inertia force virtual work term can be discretized using (4.44), (4.45) and (4.46) and expressed in the form ''p(\-'v!)^^At'dv ox. 'V = (mY'M^ (4.63) where M is given by 'M 'M i 2i-l,2j-\ 2i-\,2j | 'M 'M i A A 1Y1 (4.64) 1V1 A A _ 2'.2y-l _ _2i,2y { 2Nx2N in which i and j indicate node numbers from 1 to N, and Oh; 'M . _^'M . A A lt2j . dh. N 2j • u x u k=i x v dh flh N < dy TT J #>£ dyt! * + (/? +H H "£ V ' )]AT DV (4.65) Chapter 4. FINITE ELEMENT DISCRETIZATION 'M ='M A 1 V 1 =0 A 2,-1,2 j 1 V 1 57 2i,2 j-l (4.66) U 4.3.6 Discretization of the second convective inertia force term The second convective inertia force virtual work term can be discretized using (4.44) and (4.45) and expressed in the form (4.67) •v where 'C j tix is given by (4.53) to (4.55). Al 4.3.7 Discretization of the third convective inertia force term The third convective inertia force virtual work term can be discretized using (4.44) and (4.45) and expressed in the form ax ox k •v k (4.68) d Xj where ' C ^ is given by 3 < A1 C j = C 2,-l,2y-l 'C 2*,2 y - l ts~iAl ^2M,2 j (4.69) ts-iAl ^21,27 2Nx2N in which i and j indicate node numbers from 1 to N, and Chapter 4. FINITE ELEMENT DISCRETIZATION T 58 'dxf^dx v dh N r)h N dh N Fih N fftT^<%-2%)]^(% ~%)WdV + (4-70) k '^='0^=0 (4.71) 4.3.8 Discretization of the fourth convective inertia force term The fourth convective inertia force virtual work term can be discretized using (4.44) and (4.45) and expressed in the form j 'pCv.-'vfX'v-'v^^i^Su^At'dV dx dxj 'v = (Su) 'C 'v T (4.72) A4 J k where C is given by ( -"2M,2 j-\ t "2i-l,2 j (4.73) ts~iA4 'C j _JhlJzL 2i,2j 2Nx2N in which i and j indicate node numbers from 1 to N, and J + + dx ox W/a<j,a'* + dx ^ffxff/h & ir < ** * * - +h dydy v dy ^ —^ k )]2 v xk v * ] A t W k ) h k k y k y k ) (4.74) Chapter 4. FINITE ELEMENT DISCRETIZATION 59 It should be noted that this term involves second order derivatives of shape functions. 4.3.9 Discretized Dynamic A L E Equation Substituting into (4.43), we get, for a single element or a group of elements ' M a + ' C \ + ' C " (v - x ) +(' K 2 t At _ = + fext 8 L 1 +' K " )u +(' K ^ +'K )(u - u ) 1 _('M+'M^)'a-('C^Vc" +'C^)'v 3 lf A2 8 (4.76) Equation (4.76) represents the finite element matrix equation for dynamic A L E analysis. Chapter 5 IMPLEMENTATION 5.1 MESH MOTION Using the A L E formulation, the finite element grid points can be moved arbitrarily to maintain a homogeneous mesh and to properly represent boundary conditions throughout the deformation process. In this research, grid displacements are first related to material displacements through a set of arbitrary mesh motion parameters. This method is very efficient in controlling the motion of the grid in different parts of the domain. Specification of pure Lagrangian and pure Eulerian degrees of freedom using this method is a simple task. It is also very efficient in reducing the number of solution variables by condensing out grid displacements prior to solution. The choice of the arbitrary mesh motion parameters for interior degrees of freedom is handled by a special mesh motion scheme. Special treatment for mesh motion on free material boundaries is necessary to satisfy the A L E boundary constraint. 5.1.1 Grid displacement In [39], grid displacements are related to material displacements using a relation of the form U* =a fJ u j+ U) 60 U) I no summation on j (5.1) Chapter 5. where a } IMPLEMENTATION and 61 are two vectors of mesh motion parameters, and the brackets in the subscript (J) indicate no summation on j. Although this relation is quite efficient in many cases, it couples grid and material displacements at the same degree of freedom only, thus restricting the implementation of the A L E boundary constraint on material boundaries. In this work, a better control over the mesh motion, especially on free material boundaries, is achieved by associating grid displacements with material displacements by the more general relation u*=a + Bu (5.2) where a and B are a vector and a matrix of mesh motion parameters, respectively. Vector a consists of appropriate grid displacements given by the mesh motion scheme while matrix B consists of factors that allow the coupling of grid and material displacements. B is usually chosen to be a diagonal matrix, i.e. grid and material displacements are coupled only at the same degree of freedom. On free material boundaries, however, it is sometimes necessary that all degrees of freedom of grid and material displacements be coupled at the same node. For two-dimensional problems, this would result in B being a tridiagonal matrix. As special cases to the general A L E motion, pure Lagrangian degrees of freedom can be obtained by setting a = 0 and B = I, whereas pure Eulerian degrees of freedom are obtained by using a = 0 and B = 0 . 5.1.2 Mesh motion for interior nodes In this work, the transfinite mapping method [66, 67] is used as the mesh motion scheme for the degrees of freedom interior to any mesh region with four known boundary Chapter 5. IMPLEMENTATION 62 curves. This method provides a homogeneous mesh and matches the boundary of a given region at an infinite number of points. Another distinct advantage of the' transfinite mapping method is that it allows the discrete representation of boundary curves, i.e. the coordinates and displacements of boundary nodes can be used to find the optimum position of the nodes internal to the region. It also allows for discontinuities in slope of boundary curves. The transfinite mapping algorithm starts by partitioning the initial mesh into a number of regions of simpler forms. Considering the mapping of a typical distorted mesh region bounded by four curves $(r,0), $(r,l), ^,(0,s) and ^(l,s) as shown in Figure 5.1, the new mesh coordinates can be obtained by mapping the region onto a unit square to give ' ' x,. (r,s) = (l- s)^ (r,0) + +A (r,l) + (1 - r)#. (0, s) + rfa (1, s) - (1 - r)(l - JM- (0,0) - (1 - r)s^ (0,1) - rsfr (1,1) - r(l - s)4> (1,0) t (5.3) where 0 < r < 1 and 0 < s < 1 are the normalized coordinates over the region. Thus, the mesh motion parameter a for degree of freedom / internal to a region can be given by the t transfinite mapping method as ct.=' 'x,.-'x,. +A (5.4) 5.1.3 Mesh motion on free material boundaries For the transfinite mapping method to give a good quality mesh within a certain region, grid points on the boundaries of this region must be uniformly distributed. Consider point k located on a free material boundary as shown in Figure 5.2. Because of its location on a free material boundary, the A L E motion of point k is given by (5.2) as Chapter 5. IMPLEMENTATION Figure 5.2. Mesh motion onfreematerial boundaries. 63 Chapter 5. IMPLEMENTATION 64 <=cc +B u +B u x xx x xy (5.5) y u =a +B u +B u (5.6) g y yx x yy y where coupling between the x and y degrees of freedom is introduced in order that the grid motion of point k satisfies the boundary constraint given by equation (3.5). The question in this section is to find the mesh motion parameters a , B , B , a , B , and B x xx xy y yx yy for node k which satisfy the boundary constraint [68]. We define a set of local axes x' and y' at grid point k such that x' is tangent to the boundary and y' is its normal. Assume that the local x' and y' axes at point k are inclined to the global x and y axes by an angle 6. The components of the incremental material and grid displacement vectors in the global and local coordinate systems are related by u' — u cos0 + u sin0 x x u' - -u sint? + u cosO (5.8) =u cose + u sine (5.9) v u' g u' x y s s x v = -u g sin6+ u cosO g v (5.7) v (5.10) g Meanwhile, the mesh motion equations referred to the local coordinate system can be written as u' =a' B' u' B' u' (5.11) g x u' g where a , B'^ and B' x xy x+ xx x+ xy y =a' +B' u' +B' u' y yx x yy (5.12) y are the A L E mesh motion parameters in the direction tangent to the boundary and which may be arbitrarily chosen, while a' , B' y yx and B' yy are in the Chapter 5. IMPLEMENTATION 65 direction normal to the boundary controlled by the boundary constraint. In this work, cubic spline interpolation is used to find new coordinates for point k such that all grid points on this free material boundary are uniformly distributed. The difference between the old and new coordinates for point k establishes the mesh motion parameter a' . In this x case B' and B' xx can be set to zero. The boundary constraint dictates that u' =u' thus g xy y giving a' = 0, B' y letting B' xx yx = B' xy = 0 and B' = 1. Substituting (5.7), (5.8) and (5.9) into (5.11) and yy = 0, we get u cos6 + u sin0 = a' g (5.13) g x Substituting (5.7), (5.8) and (5.10) into (5.12), and applying the boundary constraint a' = y Q,B' = 0 and B' = 1, we get yx yy — u sinf5 + u cos6 = —u sinO + u cos6 g g x x Solving (5.13) and (5.14) for u and u , we get =cc' cos9 + u sin 0-u 2 x (5.14) g x u y x x u =oc' sind-u sintfcostf + z^ cos # g x sin#cosf? 2 x (5.15) (5.16) Comparing equations (5.15) and (5.16) with equations (5.5) and (5.6), we get ct =ct' cos6 (5.17) 73 =sin 6> (5.18) x x 2 xv B xy =-sin<9cos<9 a =a y B vx x sin 0 =-sin<9cos<9 (5.19) (5.20) (5.21) Chapter 5. IMPLEMENTATION 66 B =cos 6 (5.22) 2 yy The above A L E mesh motion parameters ensure that mesh motion on free material boundaries is consistent with the boundary constraint given by equation (3.5). 5.2 SOLUTION OF NONLINEAR EQUILIBRIUM EQUATIONS The A L E formulation derived in the previous chapters has been implemented into a nonlinear finite element code. Details of the computer solution of the nonlinear equilibrium equations for quasi-static and dynamic analysis are outlined in this section. 5.2.1 Quasi-static analysis Due to the approximations involved in the linearization of equation (4.40), an iterative procedure must be used to ensure equilibrium. Employing the Newton-Raphson iterative scheme, equation (4.40) can be rewritten as t+Mj^w-\) (i) u where ' and u + A g ( , ) 'K I ( M ) and ' + A 'K I+MK '- (u /f( 1) _ «(')) ' 'f '_' 'fO'-o (/) +A + U CT +A = are the tangent stiffness matrices in iteration / ( < M ) (5 23) u ( 0 are the corrections to the incremental material and grid displacement vectors in iteration / and which are used to update displacements according to ' 'u ' =' V''- +u +A ' + A 'u ( g ( 0 ) +A =' + A 1) 'u g ( M ) (i) +u g ( 0 (5.24) (5.25) Equation (5.23) corresponds to the full Newton iteration scheme. The calculation and factorization of the tangent stiffness matrices represents a major computational cost per Chapter 5. IMPLEMENTATION 67 iteration. The modified Newton iteration involves fewer stiffness reformations than the full Newton scheme. Using the modified Newton iteration, equation (5.23) is written as ' K V ' V K ^ U * 0 -ii* ' )=' ( ) +A, f '-' 0[ +A/ f '(, (5.26) ,) However, convergence in the modified Newton iteration is in general slower. In this work, both methods have been implemented and the modified Newton iteration is found to be computationally more efficient. Rearranging equation (5.26) we get ('K +'K )u L A -'K u 0) A =' s(i) +Al f ext _ ' * f <'-'> (5.27) + 5.2.2 Dynamic analysis Using the modified Newton iteration, equation (4.76) for dynamic A L E analysis can be written in the form 'M a i ( 0 +'C\ ( i ) +'C' (v -v 2 in which a , v , v (,) ( , ) s ( , ) , 0 j W )+('K +'K )u i l i 2 ( i ) +('K^+'K^ )(u 2 (0 -u g < 0 ) are the corrections to the incremental material acceleration, material velocity and grid velocity vectors whereas ' ' a +A (,_1) and ' + A / ( , _ 1 ) V are the material acceleration and material velocity vectors at time t + At, iteration / - 1 . Next material acceleration, material velocity and grid velocity approximations are obtained using ' 'a =' 'a - +a +A (0 <+A ; V '+A' g(o y ) = +A 1) , A, (,--I + t+& = (/ v v-i) yg Equation (5.28) can be rewritten in the form + ) + v (0 (O gd) y (5.29) ( 5 _ 3 (5.31) 0 ) Chapter 5. IMPLEMENTATION ' M V° 68 +('C^+'C^)v '' -'C^V < ) ( / ) +('KVK V ° - ' K V ( 0 _(']yi '4.']y[ y ' ('- ) _'tQA\ tQAI tQA4y+AtyU-V z _t+M^ext_t+&tfd-\) /, +A 1 a + + (5 32) In this work, two methods are used to integrate the dynamic A L E equilibrium equations: implicit time integration using the Newmark scheme and mixed implicitexplicit time integration using the predictor-corrector algorithm. It is noted that the latter scheme might be advantageous in some cases to avoid convergence problems. (i) Implicit time integration using the Newmark scheme Using Newmark implicit time integration, the following assumptions are used [48] P)' a + At p 'a u ='u + At'v + At (~ ( 0 2 2 <' +A ( / ) v =' + At(l - y)' a + Aty ' a ,+A ,+A (5.33) (i) (5.34) (0 v where /? and y are parameters that control the accuracy and stability of integration. Rearranging (5.33) we get, At p 2 2 H Using (5.24) <' +A ( 0 a =_^[' + A At p 2 _ ' _J_ t A, ( M ) U u [ + (,) (,-i) _ , _At,v u At p 2 = ^ a L ( , - i At ( - - p)' a] + u -vi-At'y- 2 2 _ 2 }__ pya ]+ At + 1 At p 2 Comparing (5.36) with (5.29) we get u ( 0 l_ (0 U { x ) H 2 r ' J At p 2 ( 5 . 3 6 ) Chapter 5. IMPLEMENTATION 69 a ( 0 =—^u (5.37) < 0 At p 2 Substituting (5.36) into (5.34) ' <v +A ( 0 =' v + At(\ - y)' a + Aty ' a !+A = '+A/ (/-i v ) + _7_ + ( M ) Aty—^u (0 </) ( 5 3 8 ) Comparing (5.38) with (5.30) v =-^-u AtB (,) (5.39) (/) Similarly (5.40) *('-> _ L _ « ( ' ) v = U AtB Substituting (5.37), (5.39) and (5.40) into (5.32) we get 1 _ / i _L_Cc "+'C )+('K -i-'K' )]u ' y t Ar/? M + yl2 t < (, A?/? _t+Atfext_t+&tf(i-\) _/t ) -[^'C^ +'K^]u* 2 Atp j y j i _|_/ jy|/( y+A/ (/-l) _ ( / ^ ^ l / ^ . / ( 3 _ _ / ^ ^ 4 y + A / ( ; - l ) a + (/) ) v (5 41) (ii) I m p l i c i t - e x p l i c i t p r e d i c t o r - c o r r e c t o r i n t e g r a t i o n s c h e m e The implicit-explicit predictor-corrector integration scheme associated with the Newmark algorithm can allow the finite element mesh to contain two groups of elements: an implicit group and an explicit group [69, 70]. In addition, elements' integrations could be shifted between the two schemes as necessary. In the predictor phase, we set i = 1 and use the predictor values: Chapter 5. IMPLEMENTATION t + A t 70 u " =' ' u =' u + Ar'v + Ar ( 1 - /?)' a ( +A ' 'v +A p (5.42) 2 =' V =' v + Af (1 - y)' a (1) (5.43) +A ' 'a =0 +A (5.44) (1) Iterations are then performed to satisfy the equilibrium equation in (5.41). In the corrector phase of each iteration, displacements, accelerations and velocities are updated according to: ' 'u =' 'u ''- +u +A ' 'a +A (/) (0 +A ( 1) =C u - 'u )/(At B) A, (i) ,+A p (0 (5.46) 2 ' 'v =' V+A/y 'a +A (5.45) (/) +A +A (0 (5-47) 5.2.3 Elimination of grid displacements on the element level Finite element equilibrium equations, equation (5.27) for quasi-static analysis and equation (5.41) for dynamic analysis, can both be written in the general form 'Ku where ' K and 'K g respectively, while f ( 0 -'K u g g ( 0 =f are equivalent stiffness matrices corresponding to u ( , ) (5.48) (i) ( , ) and u g ( , ) , is the incremental load vector for iteration /. The relation between the material and grid displacements in (5.2) is considered as a supplementary constraint equation to the finite element equilibrium equation. By introducing this constraint on the element level, grid displacements can be condensed out of element equilibrium equations prior to solution. Substituting (5.2) into (5.48), we get Chapter 5. IMPLEMENTATION 71 ('K-'K B)u g ( / ) =f 'VK a ( g (5.49) The only limitation to this procedure is that grid and material displacements may only be coupled at degrees of freedom within one element. Equation (5.49) can be rewritten in the form 'K*u =f* (/) (5.50) where ' K * and f * are the equivalent stiffness matrix and nodal force vector respectively. Conventional finite element assembly and elimination techniques may now be applied directly to solve for the material displacements. 5.2.4 Frontal solver A frontal or a wavefront solver [71, 72] that can handle asymmetric matrices is implemented to solve the simultaneous linear equations in (5.50). The number of equations that are active after any element has been processed during solution is called the wavefront. The wavefront is determined by the sequence in which the elements are arranged or ordered. The computer time required for solution is proportional to the square of the mean wavefront. Therefore, it is crucial to be able to minimize the wavefront. A separate module was developed to reorder elements and minimize the wavefront. Elements are reordered such that the element for which each degree of freedom is first mentioned is as close as possible in sequence to the element for which it is mentioned last. Once the element sequence has been optimized, wavefront solution starts by scanning all elements to determine which element is the last to use each degree of freedom. Chapter 5. IMPLEMENTATION 72 Solution proceeds by adding the equations for degrees offreedomrelated to the current element and which occurs for the first time. The equation for a degree offreedomthat occurs for the last time is algebraically solved in terms of the remaining unknowns and eliminated from the assembled matrix by Gauss elimination. The equation is then written to a file for later back-substitution and the remaining equations are modified. The assembled matrix expands and contracts as degrees offreedommake their first and last appearance in the element definitions. Another feature that has been developed for the solver is the ability to handle multipoint constraint equations. These equations specify relations between two or more degrees offreedomduring solution. In addition, a direct matrix transformation method is used to impose coupling of degrees offreedomat the same node. 5.2.5 Line search A line search algorithm [8] has also been implemented to speed up the convergence. In some situations, the use of the full u (,) as obtained by the solver leads to solution instabilities. The line search algorithm attempts to improve the Newton-Raphson solution u (,) by scaling the solution vector by a scalar value s termed the line search parameter. Thus, equation (5.24) can be rewritten in the form + SU (0 (5.51) where 0.05 <s< 1.0 (5.52) Chapter 5. IMPLEMENTATION 73 The line search parameter s is determined by minimizing the norm of the residual force vector through iterations. 5.3 CONTACT ANALYSIS Simulation of many metal forming problems necessitates the ability to model the contact phenomena that occurs between the forming tools and the workpiece. Contact analysis is a complex problem because of the requirement to track the motion of the bodies involved and to accurately represent the friction between their surfaces. The numerical objectives are to detect the contact of the bodies and apply enough constraints and boundary conditions to simulate the frictional behavior, avoid penetration and allow for separation when necessary. Several numerical techniques have been developed to perform these objectives. These methods include the Lagrange multipliers method, the penalty function method and the direct constraint method. 5.3.1 Lagrange multipliers The Lagrange multipliers method is one of the most common techniques to treat the contact problem in the literature [18]. In performing contact analysis, we are basically solving a constrained minimization problem where the constraint is the no penetration constraint. In this method, the finite element equilibrium equations are augmented using the constraint equations by introducing an array of Lagrange multipliers as additional degrees of freedom. In contact analysis, Lagrange multipliers signify the contact pressure. Chapter 5. IMPLEMENTATION 74 Lagrange multipliers technique has often been implemented in contact procedures using special interface elements known as gap elements. If the constraints are properly written, this method ensures that penetration does not occur. Unfortunately, the introduction of Lagrange multipliers leads to numerical difficulties as their inclusion results in a nonpositive definite mathematical system. This requires additional operations with high computational costs. Another problem with this method is that there is no mass matrix associated with the Lagrange multipliers degrees of freedom. This results in a global mass matrix that cannot be inverted. This precludes the use of the Lagrange multipliers technique in explicit dynamic simulations. In addition, the use of interface elements requires a prior knowledge of where contact occurs and puts a restriction on the amount of relative motion that can occur. This may not be feasible in the simulation of many manufacturing problems. 5.3.2 Penalty function The penalty function method is an alternative procedure to numerically implement the contact constraints in an approximate way. Using this approach, some penetration occurs with the amount being determined by the penalty constant or function. The penalty approach can be considered as analogous to a nonlinear spring between the two bodies. The penalty method is relatively easy to implement and has been extensively used in explicit dynamic analysis. However, the choice of the penalty constant has a detrimental effect on the numerical stability and accuracy of the method. It is also possible to combine the penalty method with the Lagrange multipliers method in a hybrid fashion. In Chapter 5. IMPLEMENTATION 75 the hybrid methods, the contact element is derived from a complimentary energy principle by introducing the continuity on the contact surface as a constraint and treating the contact forces as additional elements [8]. 5.3.3 Direct constraint method In this method, the motion of the contact bodies is tracked, and when contact occurs, direct constraints are applied on the bodies as boundary conditions. Both kinematic constraints on transformed degrees of freedom and nodal forces may be applied. This procedure can be very accurate if the program can predict when contact occurs. The direct constraint method is the method that has been implemented in the developed finite element program. No special interface elements are required and complex changing contact conditions can be simulated since no knowledge of where contact occurs is required prior to the analysis. The procedure has been implemented as a direct node-tosurface contact only for the case of contact between the workpiece and one or more rigid bodies [62]. The procedure starts by reading the input file that defines the perimeters of the tools as sets of geometrical entities connected together, namely straight line and arc segments. In each incremental load step, the contact algorithm tracks the motion of boundary nodes of the workpiece to detect any penetration into any of the tools. Upon detection of such penetration, the algorithm determines the magnitude and direction of the penetration of the nodes and applies them as corrective prescribed displacements to reposition the nodes on the tool surface. Once the nodes are on the tool surface, appropriate contact boundary Chapter 5. IMPLEMENTATION 76 conditions (e.g. sticking, slipping with friction, etc.) are applied on the contact nodes. Both Coulomb and shear friction laws can be applied. A node is considered in contact as long as the nodal force normal to the contact surface is compressive. If the normal force becomes tensile, the node is separated from the surface and its nodal forces are set to zero. In each step, iterations are performed until no change of contact conditions is detected. The contact algorithm implementation as it stands now suffers from major drawbacks and causes many numerical problems. Among these drawbacks is that the corrective prescribed displacements are applied on the nodes that have penetrated the contact surface. This intermediate position does not satisfy the contact conditions and any corrective action based on this position will be approximate. Another approach to overcome this problem is to define a small tolerance region on both sides of the contact surface in which the nodes are assumed to lie on the surface when penetration is detected. If the penetration of the node is within the tolerance region, contact boundary conditions are applied to the node in the next load step without any corrective action. If a node penetrates the entire tolerance region in one load step, then a smaller load step size should be used and the computation is resumed based on the equilibrium position prior to penetration. 5.4 PROGRAM STRUCTURE The developed A L E formulation has been implemented into a finite element computer program. The program has been designed with emphasis on modularity and use Chapter 5. IMPLEMENTATION 11 of standard Updated Lagrangian calculations whenever possible. The program is divided into subroutines each with a specific function while the main program consists mainly of subroutine call statements. As discussed in Chapters 2 and 3, the version of the A L E formulation developed in this work can be considered an extension to the Updated Lagrangian formulation. Thus, the first step in the computer development in this work was to implement an Updated Lagrangian finite element formulation and test it. A l l the necessary routines related to the nonlinear large deformation Lagrangian mass, stiffness, and stress calculations were developed. In addition, the necessary solution schemes for the static and dynamic matrix equations were implemented and tested. The next step was to implement the necessary routines for the new A L E formulation. These routines were comprised mainly of calculations for the additional mass and stiffness matrices as well as mesh motion routines. The last step of development was to implement a simple contact algorithm to enable the program to simulate a wider class of metal forming and large deformation applications. In the following, the main computer routines called in the main program and the main flow charts developed in each of the program building stages are highlighted. 5.4.1 Updated Lagrangian routines Figure 5.3 shows a flow chart of the main program for the Lagrangian calculations. The flow chart only shows the main routines needed for the Updated Lagrangian code and the function of each routine is briefly stated. The additional subroutines that are called from the main routines are not listed here for brevity. The program starts by Chapter 5. IMPLEMENTATION 78 reading the control variables that control the size of the program arrays. These variables are needed beforehand for dynamic memory allocation purposes. Next, all variables and arrays used in the program are initialized. Problem input data, such as nodal coordinates, element connectivity, applied loads and material properties is then read and stored in the appropriate arrays. Distributed loads are then converted into their equivalent nodal values and the total applied loads are divided into incremental load steps. Newton equilibrium iteration starts by setting the boundary conditions for the current load step and current iteration. The predictor phase of the time integration algorithm is first preformed and elements lumped mass matrices are calculated (for dynamic analyses only). Elements stiffness matrices are then calculated and stored on the hard disc for later assembly. The stiffness matrices in the Lagrangian formulation include both nonlinear material and large deformation effects. The contribution of all elements to the global effective load vector is assembled prior to solution. The global effective load vector also includes out of balance internal forces from previous iterations. The frontal solution scheme is then used to solve for the displacements. Assembly of elements equivalent stiffness occurs during the solution process. The line search algorithm is then called to accelerate convergence. The corrector phase of the time integration scheme is then performed in which velocities and accelerations are computed (for dynamic analyses only). Stress integration is then performed to calculate the stresses at the Gauss (integration) points. Convergence of iterations is then checked. Three measures are used to check the convergence. These are: iterative displacements, residual (out of balance) forces and residual energy. Iterations are performed until convergence is achieved. Chapter 5. IMPLEMENTATION 79 At the end of each load increment, extrapolation of Gauss point stresses to obtain the nodal stresses takes place. Finally the program outputs the requested results for the solved load increment and moves on to the next step. The modular structure of the program is apparent from Figure 5.3. Chapter 5. 80 IMPLEMENTATION C START I CONTRO I INTTIA INPUTD I LOADPL PREFIX PREDIC ~T~ ELMASS ESTIFF I ESTIKJ I EFLOAD ~~rFRONTS ~~T~ Next Next Step Iteration LNSRCH I ~r~ CORREC STRESS } Read Control Variables Initialize Variables & Arrays Read Input Data Calculate Consistent Nodal Loads Set Boundary Conditions Integration Predictor Phase Calculate Elements Lumped Mass Matrices Calculate Elements Stiffness Matrices Calculate Extra Jaumann Stiffness Matrices Assemble Elements Effective Load Vector Frontal Solver Line Search Algorithm Integration Corrector Phase Calculate Gauss Point Stresses I CONVER I ^ ' a c u ' a t e Convergence Parameters Calculate Nodal Stresses OUTPUT Write Output for Current Step Figure 5.3. Flow chart of the Updated Lagrangian calculations. Chapter 5. IMPLEMENTATION 81 5.4.2 A L E routines The flow chart of the main program for the ALE calculations is shown in Figure 5.4. Comparing Figures 5 . 3 and 5.4, one can easily see the necessary ALE additions and changes to a Lagrangian code. One of the major changes to the Lagrangian code is the introduction of the mesh motion routines. The mesh motion routines decide how the mesh should be moved in the current iteration and choose the mesh motion parameters a and B accordingly. The user may specify how the mesh should behave on some parts of the boundaries. Boundaries can be specified to be either pure Lagrangian or Lagrangian in the normal direction and ALE (or Eulerian) in the tangential direction. After the motion of the boundaries has been set for the current iteration, the transfinite mapping method determines how the interior regions of the mesh should be moved to preserve its uniformity. Among the important additional ALE routines are those that calculate the extra stiffness and mass matrices and load vectors. Assembly of those extra terms takes place in the solver routines according to the ALE equilibrium equations in which the mesh motion parameters are introduced. A final ALE change is required to the stress calculations routines to compute the additional stress terms due to convective effects. The consistency and simplicity of the approach used in deriving and implementing this version of the ALE formulation are highly appreciated by highlighting its implementation advantages in comparison with previous implementations [ 3 9 ] : - No updating of material properties is required between load steps since the stresses calculated within the iterations already include convective effects. Chapter 5. IMPLEMENTATION - 82 No remeshing or boundary nodes adjustment takes place between increments. This is due to new mesh motion scheme on free material boundaries that allow nodes to be Lagrangian in the normal direction and to move freely in the tangential direction. - Nodal stresses are typically calculated for postprocessing purposes only and are not involved in any equilibrium calculations. This is due to the new treatment of convective effects in the equilibrium equations. - The modularity of the implementation also helped in retaining the simplicity of the program structure and allowed for clear identification of the necessary A L E changes to Lagrangian codes. Chapter 5 . 83 IMPLEMENTATION Set Mesh Motion Parameters Calculate Extra ALE Mass Matrices Calculate Extra ALE Stiffness Matrices Add Extra ALE Load Terms Solve ALE Equilibrium Equations Add Extra ALE Stress Terms Figure 5 . 4 . Flow chart of the A L E calculations. Chapter 5. IMPLEMENTATION 84 5.4.3 Contact routines As stated earlier, the last stage in the code development was to include contact analysis. Although A L E can alleviate the need for contact capabilities in some cases, however, contact analysis is essential in many metal forming and large deformation applications. In this work, contact analysis is based on the direct constraint method [62]. Figure 5.5 shows a flow chart of the A L E code with the contact routines included. The first group of contact routines read the contact related data consisting of tool geometries, list of potential contact nodes on workpiece surface and list of nodes that are initially in contact. Another routine is called in the beginning of each load step to apply the necessary displacement boundary conditions on contact tools by updating the coordinates of the tools geometry. Within the contact and equilibrium iterations, a new routine is required to correctly set the contact-related and non-contact-related boundary conditions. After equilibrium iteration converges, the contact algorithm is called to check for any possible changes (penetrations or separations) in the contact conditions. If any changes are detected, corrective action is taken and another contact iteration is performed. This process is repeated until no change in contact condition is detected for the current load step after which a new load step is applied next. It is worth noting that the implemented contact algorithm suffers from convergence problems and future modifications are required. Chapter 5. IMPLEMENTATION 85 C START I CONTRO — I — INITIA H I INPUTD I LOADPL ~~r~ CONDAT CONDSP Read Contact Data A p p l y Displacements to Contact Surfaces MOTION I CONFIX PREDIC ~ T ~ ELMASS ~~T~ EMCALE ~ T ~ ESTIFF zr: ESTIKJ ESTALE OQ© 15 Set Boundary and Contact Conditions Chapter 5. IMPLEMENTATION 86 Contact Algorithm N 'Contact* ^Convergence^ .Check. NODSTR OUTPUT Q END ) Figure 5.5. Flow chart with contact routines. Chapter 6 APPLICATIONS 6.1 ONE-DIMENSIONAL STRESS WAVE PROPAGATION The first application simulated using the developed finite element program is a wave propagation problem in a one-dimensional infinitely long elastic-plastic rod and is used to test the code with transient effects included. The same problem was reported in [33] and [37]. It should be noted that this problem doesn't require an A L E analysis and was selected because of the availability of analytical and numerical solutions. The infinitely long rod is discretized using 400 elements with a mesh size of 0.1 units as shown in Figure 6.1. The material properties of the rod are assumed to be: density p = 10000.0, Young's modulus E = 10000.0, plastic modulus E = 3333.33, yield stress o = 75.0 and p y Poisson's ratio v = 0.0. The rod is subjected to a compressive stress wave 100.0 in amplitude and 4.5 in width. The stress wave and the time interval under consideration are depicted in Figure 6.2. Any consistent system of units can be used to define the data for this problem. Figures 6.3 and 6.4 compare the longitudinal stress distribution obtained using the presented A L E formulation with the analytical solution for both the elastic and elasticplastic cases. A L E results are in good agreement with the analytical solution and with those of [33] and [37] (not shown here). Deviations from analytical solution are attributed to the difficulties associated with the numerical representation of step functions. 87 Chapter 6. APPLICATIONS 88 A o(0 1 2 400 3 0.1 V 40.0 <; Figure 6.1. Mesh for the one-dimensional wave propagation problem. 100 - 0 4.5 16.0 Figure 6.2. Stress wave amplitude and duration. t Chapter 6. 50 APPLICATIONS 89 Elastic-Plastic 0 </> Ui ALE Analytical 2> -50 4-1 C/> -100 -150 0 10 20 30 40 Spatial Coordinate Figure 6.4. Longitudinal stress distribution comparison, elastic-plastic case. Chapter 6. 6.2 90 APPLICATIONS BAR IMPACT The bar impact problem has been investigated by many researchers [35, 73-75] and is considered a standard benchmark test for transient dynamic computer codes. In this problem, a cylindrical copper bar of initial radius 3.2 mm and length 32.4 mm strikes a rigid frictionless surface. The impact velocity is 227 m/s. The material is assumed to be elasto-plastic with Young's modulus E = 117 GPa, plastic modulus EP = 100 MPa, Poisson's ratio V- 0.35, initial yield stress a = 400 MPa, and density p = 8930 kg/m . A 3 y von Mises yield surface with linear isotropic hardening is assumed. Contact conditions are imposed by simply constraining the nodes in contact with the rigid surface. An axisymmetric mesh of 250 second order 8-noded elements is used. Computations are carried out up to a time of t = 80 ps. Figure 6.5 compares the deformed shape and finite element mesh obtained using both the Lagrangian and the A L E formulations at different stages during the deformation process. In the A L E solution, boundary nodes are allowed to move in the tangential directions to the boundaries to maintain their uniform distribution while satisfying the A L E boundary constraint. The Lagrangian solution is obtained as a special case from the developed A L E formulation. It is clear that while the Lagrangian solution suffered severe mesh distortion, A L E was able to maintain a uniform mesh. It is worth mentioning, however, that the computation time for this problem using A L E is approximately twice that of the Lagrangian solution. This is mainly due to the fact that the A L E matrices are unsymmetric. Table 6.1 compares the final bar length and mushroom radius obtained using the developed code with numerical results obtained by other researchers. It is clear Chapter 6. APPLICATIONS 91 that the results obtained using the fully coupled implicit dynamic A L E formulation are in agreement with other well established numerical techniques and codes. Figure 6.6 shows the distribution of equivalent stress and equivalent plastic strain contours. Table 6.1. Comparison of results for bar impact. Reference (Code/Method) Final length (mm) Final mushroom radius (mm) 21.66 21.47 21.47 21.47 7.02 7.13 7.03 7.07 21.26-21.49 6.97-7.18 21.42-21.44 7.21-7.24 21.42 21.53 7.15 6.87 21.48 21.69 7.22 6.90 Kamoulakos [73] MARC DYNA2D DYNA3D NIKE2D Zhu & Cescotto [74] Lagrangian (different element types) Camacho & Ortiz [75] Lagrangian (different remeshing schemes) Liu et al. [35] Lagrangian A L E (explicit) Current Work Lagrangian A L E (implicit) Chapter 6. APPLICATIONS 92 t=0 t=20 us / = 40 us t= Original mesh Lagrangian ALE Lagrangian ALE 80 us Lagrangian A L E Figure 6.5. Comparison of the Lagrangian and A L E solutions for bar impact. Chapter 6. APPLICATIONS 5.13644E+08 4.80256E+08 4.46868E+08 4.1348E+08 3.80092E+08 3.46703E+08 3.13315E+08 2.79927E+08 2.46539E+08 2.13151E+08 1.79763E+08 1.46375E+08 1.12987E+08 7.95988E+07 4.62107E+07 Figure 6.6. Equivalent stress and plastic strain contours for bar impact. 93 eq 3.00 2.79 2.57 2.36 2.14 1.93 1.71 1.50 1.29 1.07 0.86 0.64 0.43 0.21 0.00 Chapter 6. 6.3 94 APPLICATIONS SHEET M E T A L EXTRUSION Extrusion is a typical metal forming problem in which large strains are expected and in which remeshing and boundary condition updating are required if the traditional Lagrangian formulation is employed. In this problem, the extrusion die produces a 25% reduction in sheet thickness and is shaped in the form of a 5 order polynomial with zero th curvature and slope at both ends. An aluminum billet of thickness 2a is forced through the die by a rigid piston pressing against the rear face of the billet and moving with a prescribed displacement. The length of the die is taken as 1.2a. The same problem was previously solved using the Updated Lagrangian approach [76]. The material properties for the aluminum billet are taken as: Young's modulus = 10 ksi (68.95 GPa), Poisson's 4 ratio = 0.3, initial yield stress = 57 ksi (393 MPa) and hardening parameter = 165 ksi (1.14 GPa). Figure 6.7 shows the initial geometry and mesh used in the simulation. Because of symmetry, only the upper half of the billet was analyzed. Using the A L E formulation, all the nodes confined to the die area are set to be Eulerian during the course of deformation. The nodes on the boundaries are set as Lagrangian. The motion of all the other nodes is controlled by the A L E mesh motion scheme. Figure 6.8 shows the plastic strain distribution and the deformed shape after a piston displacement of 2a units. Contact between the billet and the die was set as boundary constraint equations. The A L E approach was able to eliminate the need for remeshing or boundary condition updating and the desired deformation level was reached without any user intervention or special contact treatments. Variations in the longitudinal stress component in the die region at different lateral positions from the mid-plane of the Chapter 6. APPLICATIONS 95 initial configuration are shown in Figure 6.9. The obtained stresses are in good agreement with those previously published [76]. Figure 6.7. Geometry and mesh for extrusion process. Chapter 6. APPLICATIONS 96 i 502172 468694 435216 401738 368259 334781 301303 267825 234347 200869 167391 133913 100434 066956 033478 Figure 6.8. Plastic strain contours and deformed shape for extrusion problem. Gxx (ksi) Figure 6.9. Longitudinal stress at different lateral positions. Chapter 6. 6.4 APPLICATIONS 97 QUASI-STATIC AND DYNAMIC COINING In this example, a coining process, also known as punch indentation, is simulated to show the effectiveness of the A L E formulation in handling contact boundary conditions and in preventing mesh distortion. This process is simulated using both quasi-static and dynamic analyses to investigate the significance of dynamic effects. The problem is also solved using A N S Y S [8] for comparison. Figure 6.10 shows the geometry and initial mesh of the body. The workpiece is placed between two rigid frictionless tools moving with constant velocity under plane strain conditions. Because of symmetry, only one quarter of the domain is modeled. The deformation process is continued up to a 60% reduction of the original workpiece height. The material is assumed to be elastic-plastic with a Young's modulus of 200 GPa, a Poisson's ratio of 0.3, an initial yield stress of 250 MPa and a hardening parameter of 1 GPa. In the Lagrangian ANSYS simulation, it is noted that applying the downward motion of the punch as prescribed boundary conditions on the workpiece nodes initially in contact with punch would effectively cause the punch size to increase during the course of deformation as shown in Figure 6.11. When the Lagrangian formulation is employed, non-material-associated boundary conditions, such as the contact between the punch and the workpiece, is best handled using contact elements. However, contact elements are difficult to handle and may cause difficulties in achieving convergence. ANSYS, by default, uses a mixed penalty + Lagrange multipliers method in the formulation of its surface to surface contact elements. Two main parameters are used to define the characteristics of these elements. The first is the normal contact stiffness factor in the Chapter 6. APPLICATIONS 98 range of 0.001 to 100.0, with a default of 1.0. A smaller value of the normal contact stiffness factor provides for easier convergence but more penetration. The second parameter is a tolerance factor that determines if penetration compatibility is satisfied. Penetration compatibility is satisfied if penetration is within a clearance of the value of tolerance factor times the depth of the underlying element. The range of this factor is less than 1.0, with a default of 0.1. Figure 6.12 shows the deformed shape of workpiece at the final stage of the process using the default settings of the contact elements in ANSYS. It is noticed that large penetration of the punch inside the workpiece was allowed when using the default values of the normal contact stiffness and tolerance parameters. Figure 6.13 shows the deformed shape when a contact stiffness of 100.0 and a tolerance of 0.00001 are used. It is noticed that the penetration is quite smaller in this case, though not completely eliminated. Reducing the tolerance parameter beyond this value did not affect results significantly. Convergence was more difficult to achieve in the case of these extreme contact parameters. Figure 6.14 shows the evolution of the deformed shape obtained using the developed A L E formulation. The figure shows that the contact condition between the punch and the workpiece is accurately satisfied. This was easily achieved by allowing the degrees of freedom of the nodes directly under the punch to be Lagrangian in the vertical direction (to satisfy the boundary constraint) and Eulerian in the horizontal direction (to maintain the same punch size under deformation). No special contact algorithm was needed to handle the contact conditions. In addition, the transfinite mapping method was able to Chapter 6. 99 APPLICATIONS maintain a homogeneous mesh throughout the deformation history. In this analysis, the 60% reduction in height is solved using 750 incremental steps. Table 6.2 compares the stresses at the punch corner after the first load step of the total 750 steps of the analysis in the three cases: (a) ANSYS with default contact parameters; (b) ANSYS with extreme contact parameters; (c) A L E . It is clear that the stresses obtained by defining extreme contact parameters in ANSYS match the A L E analysis. However, stresses obtained by using the default values for the contact parameters, which provide for easier convergence, are inaccurate even at such early stages of the analysis. Table 6.2 Comparison of stresses at punch corner after the first load step. (a) A N S Y S with default contact parameters. (b) A N S Y S with extreme contact parameters. (c) Present A L E formulation. Sxv(MPa) Sw (MPa) Case S x (MPa) -39.50 9.31 -3.05 (a) -70.12 -234.18 23.35 (b) 23.65 -70.43 -234.76 (c) X Szz (MPa) -12.76 -91.30 -91.56 Figure 6.15 compares the load-displacement curves obtained by ANSYS and the developed A L E formulation. It is clear that the erroneous load fluctuations associated with contact analysis in the Lagrangian solution can be avoided when using A L E . The significance of dynamic effects is investigated by examining the indentation problem at different punch velocities. Figure 6.16 compares the final deformed shape and equivalent plastic strain distributions for the quasi-static case and for four different punch velocities. Dynamic effects are not significant for punch velocities less than 1 m/s since the deformed shape and the plastic strain distribution are quite similar to those of the Chapter 6. 100 APPLICATIONS quasi-static simulation. One can also notice that low punch velocities allow the workpiece to flow horizontally away from the punch, while at high velocities the workpiece tends to back extrude. 24 mm Punch 20 mm i T 60 mm Figure 6.10. Geometry and initial mesh for coining process. Chapter 6. APPLICATIONS Punch Figure 6.13. ANSYS solution with extreme contact parameters. 101 102 Chapter 6. APPLICATIONS Punch Original Shape 20 % Reduction 40 % Reduction 60 % Reduction Figure 6.14. Evolution of the deformed shape during the coining process. Chapter 6. APPLICATIONS Figure 6.15. Punch load-displacement curve comparison. 103 Chapter 6. APPLICATIONS Chapter 6. 6.5 105 APPLICATIONS V-BENDING Although sheet metal bending processes, such as the V-bending process, might seem to be simple processes, precise numerical simulation using elastic-plastic finite elements with contact analysis is not an easy task. This example serves two purposes. The first objective is to elaborate on the power of A L E in overcoming the numerical difficulties associated with the simulation of metal forming processes with corner contact. The second objective is the ability to perform simple experiments for this type of process and thus verify the A L E predictions against experimental results [77]. Figure 6.17 shows an isometric sketch for the setting of the V-bending experiment. The punch and die are made of mild steel 1018 while the sheet material is Aluminum 5052-H32 with a Young's modulus of 70 GPa, Poisson's ratio of 0.33, yield strength of 195 MPa, tensile strength of 230 MPa and break elongation of 12 %. The sheet metal is 114.3 mm (4Vi in.) long, 63.5 mm (2 A in.) wide and 6.35 mm (14 in.) thick. The punch l and die both have 90° V angles and the punch has a nose radius of 9.525 mm (0.375 in.). The bending experiment was conducted using a Tinius Olsen U T M at a constant downward speed of 0.2 in/min. Figure 6.18 shows the actual experimental setup. The V-bending process is simulated using the developed A L E code as well as using the Lagrangian F E A program M A R C [18]. Accurate treatment of the tool-workpiece contact constitutes a major problem for simulations based on the Lagrangian formulation. Huang and his coworkers [78, 79] showed that it is very difficult to numerically satisfy the sharp corner contact condition that occurs between the die edge and the workpiece without penetration. The material behavior used in the A L E simulation is assumed to be Chapter 6. APPLICATIONS 106 bilinear elastic-plastic with kinematic work hardening. Figure 6.19 shows the finite element mesh used in the simulation. Figures 6.20, 6.21 and 6.22 show the development of the deformed shape for the V bending process at different stages obtained using the experiment, M A R C and A L E simulations, respectively. Figure 6.23 shows the large corner penetration that occurs at an intermediate stage as obtained by M A R C even though a relatively fine mesh was used. Unlike A N S Y S , in which surface to surface contact can be defined, the contact algorithm in M A R C is based on node to surface type of contact. In the A L E simulation, contact analysis is only introduced between the punch and the sheet. The sharp die corner contact with the sheet was handled in A L E by using Eulerian degrees of freedom for the node at the die corner [77]. Figure 6.24 compares the load-deflection curves obtained from the experiment with the simulations. It is clear from the figure that the sharp corner contact in M A R C produces large fluctuations in the load prediction. It is also clear that the A L E predictions are smoother. The fluctuations in the A L E predictions are attributed to the contact analysis between the nose of the punch and the sheet that A L E could not eliminate. Chapter 6. APPLICATIONS Figure 6.17. Isometric sketch for the setting of the V-bending experiment. 107 Chapter 6. APPLICATIONS Figure 6.18. Tinus Olsen U T M used in the V-bending experiment. 108 Figure 6.19. Finite element mesh used in the simulation of the V-bending process. Chapter 6. APPLICATIONS (a) Figure 6.20 (a-b). Development of the deformed shape during the experiment. 110 Figure 6.20 (c-d). Development of the deformed shape during the experiment. Figure 6.20 (e-f). Development of the deformed shape during the experiment. Chapter 6. APPLICATIONS Figure 6.21 (a-b). Development of the deformed shape using M A R C . 113 Figure 6.21 (c-d). Development of the deformed shape using M A R C Chapter 6. APPLICATIONS Figure 6.21 (e-f). Development of the deformed shape using M A R C 115 Figure 6.22 (a-b). Development of the deformed shape using A L E . Chapter 6. APPLICATIONS Figure 6.22 (c-d). Development of the deformed shape using A L E . 117 Figure 6.22 (e-f). Development of the deformed shape using A L E . Chapter 6. APPLICATIONS Figure 6.23. Deformed shape using M A R C showing corner penetration. 119 Chapter 6. APPLICATIONS ' 0 1 1 1 1 120 I i 5 i i i I i 10 i i i I i 15 i 1 1 I 20 Punch Displacement (mm) Figure 6.24. Comparison of load displacement curves for V-bending. Chapter 6. 6.6 APPLICATIONS 121 PLATE INDENTATION Local loading of circular steel plates resting on a ring support is simulated using A L E and compared to experimental results [80]. The objective of this example is to compare the plate response when loaded using flat-faced and hemispherically-tipped indenters. In the case of the flat-faced indenter, A L E can easily handle the sharp corner contact condition whereas the contact algorithm will handle the case of the hemisphericallytipped indenter. The plate is 6 mm thick, 300 mm in diameter and is simply supported at a diameter of 248 mm. It is transversely loaded at its center at a constant rate of 3 mm/min using flat-faced as well as hemispherically-tipped cylindrical indenters of 12.7 mm diameter. The plate is made of mild steel with a Young's modulus of 235 GPa, yield strength of 305 MPa and ultimate tensile strength of 446 MPa. The material behavior is assumed to be bilinear elastic-plastic with a plastic modulus of 1 GPa. An axisymmetric finite element model of 300 elements was used in the simulation. Figures 6.25 and 6.26 compare the load displacement curves obtained from the experiment and the A L E simulation for the two cases of the flat-faced and the hemispherically-tipped indenters, respectively. Experiments show that the plate fails by perforation at approximately 20 mm under the flat-faced indenter and at approximately 30 mm under the hemispherically-tipped indenter. Accordingly, utmost displacements of 20 mm and 30 mm were sought in the simulations, respectively. The figures show a very good agreement between the A L E predictions and experimental results. Deviations from the experimental results are more pronounced at higher deformation levels and may be Chapter 6. 122 APPLICATIONS attributed to discrepancies between the material model and actual material behavior at higher loads. Figures 6.27 and 6.28 show the predicted deformed shape for the plate when loaded using the flat-faced and the hemispherically-tipped indenters, respectively. Figures 6.29 and 6.30 compare the deformed shape and equivalent plastic strain distribution in the indenter region. The figures show that the hemispherically-tipped indenter caused a greater degree of local indentation and bending around the indenter than the flat-faced indenter. For the flat-faced indenter, A L E eliminated the need for contact analysis between the sheet and the indenter. For the hemispherically-tipped indenter, A L E prevented mesh distortion at high levels of local deformation in the indenter region. The deformed shapes are in good agreement with those of the experiments [80]. Figure 6.31 shows that the flat-faced indenter produced a steeper load-displacement curve but with a smaller load and displacement to failure than the hemispherically-tipped indenter. The indenter profile plays an important role in the plate response. The degree of local indentation and the onset of perforation are highly dependent on the indenter profile. It is seen that the hemispherically-tipped indenter caused greater degree of local indentation and bending around the indenter than the flat-faced indenter. Although the local indentation caused plate thinning, the accompanied stretching and bending served to decrease shearing effects, thus delaying the onset of perforation. Thus, plates loaded by hemispherically-tipped indenters are able to withstand greater loads at higher deformations before perforation occurs than those loaded by flat-faced indenters. Chapter 6. APPLICATIONS Figure 6.25. Comparison of loading curve for the flat-faced indenter. 123 Chapter 6. APPLICATIONS 124 Figure 6.26. Comparison of loading curve for the hemispherically-tipped indenter. Chapter 6. APPLICATIONS Figure 6.28. Plate deformed shape under hemispherically-tipped indenter. 125 Chapter 6. APPLICATIONS 126 0.362625 0.338446 0.314268 0.290089 0.26591 0.241731 0.217552 0.193374 0.169195 0.145016 0.120837 0.096658 0.062479 0.048301 0.024122 Figure 6.29. Plastic strain distribution for flat-faced indenter. 0.598763 0.558291 0.517818 0.477346 0.436874 0.396402 0.355929 0.315457 0.274985 0.234512 0.19404 0.153568 0.113096 0.072623 0.032151 Figure 6.30. Plastic strain distribution for hemispherically-tipped indenter. Chapter 6. APPLICATIONS 127 100 Flat Indenter Hemispherical Indenter 90 / 80 70 z 6 0 ^ 50 o 40 30 h 20 / F 10h// 0 t J L J 10 1 I 20 Displacement (mm) i ' ' ' I 30 Figure 6.31. Comparison of the loading curve for flat and hemispherical indenters. Chapter 6. 6.7 APPLICATIONS 128 QUASI-STATIC AND DYNAMIC BRAKE BENDING The aim of this example is to demonstrate the significance of the dynamic formulation by comparing the results of simulation of a plane-strain brake bending process (also known as air bending) using quasi-static and dynamic analyses. Choudhry and Lee [41] pointed out that although the tooling and the process for this problem are rather simple, they simulated this problem using a dynamic Lagrangian formulation to prove the effectiveness of dynamic analysis even for sheet metal forming processes which are essentially quasi-static. It also provides a challenging test for the numerical stability of the solution algorithm due to the small strains and large rotations involved. The tooling geometry for the process is shown in Figure 6.32. The punch and the die are both 10 mm in radius. The sheet is 1 mm thick. The material behavior of the sheet is assumed to be bilinear elastic-plastic with a Young's modulus of 210 GPa, Poisson's ratio of 0.3, yield strength of 300 MPa, and a plastic modulus of 1 GPa. Free boundary conditions are assumed at the sheet edge. Figures 6.33 and 6.34 show the punch load against the punch displacement for both the quasi-static and dynamic simulations, respectively. While the quasi-static formulation gives a rather smooth curve for the load-displacement curve, the dynamic formulation exhibits oscillations that agree with the static solution only in the average sense. The maximum punch displacement attained when modeling the process as a quasi-static process was 37 mm at which point the tangent stiffness matrix became extremely illconditioned and equilibrium iterations ceased to converge. However, a displacement of Chapter 6. APPLICATIONS 129 52.5 mm was obtained when dynamic effects were considered. The results are in agreement with those previously reported [41]. Figure 6.35 shows the distribution of the equivalent plastic strain at different punch displacements. It is shown that the sheet deformation is primarily due to bending in the region in contact with the punch and virtually no strain occurs away from the punch. One is able to conclude that although inertia effects are negligible in this process, contributions from the dynamic terms seem to stabilize the convergence problems resulting from the large rigid body motions with little straining. Virtually identical strain distributions are obtained for both quasi-static and dynamic formulations when a convergent solution is obtained. Figure 6.32. Tooling geometry for the brake bending process. Chapter 6. APPLICATIONS Figure 6.33. Quasi-static load-displacement curve for the brake bending process. 130 Figure 6.34. Dynamic load-displacement curve for the brake bending process. Chapter 6. APPLICATIONS Chapter 7 CONCLUSIONS 7.1 SUMMARY The accomplishments and conclusions of the present work may be summarized in the following: Derivation of fully coupled A L E virtual work equations from the basic principles of continuum mechanics for quasi-static and dynamic analyses. The fully coupled A L E virtual work equations derived in this work are unique as they are based on a new treatment for convective terms and because displacements are used as the independent variables as opposed to velocities which is a common practice in the A L E literature. The derivation introduces implicit fully coupled dynamic A L E calculations for the first time. Introduction of a new method for the treatment of convective terms in the A L E equilibrium equations that avoids the use of unjustified assumptions in calculating spatial gradients of stresses. This new method results in a consistent formulation and maintains its generality and accuracy. Discretization of the A L E virtual work equations using isoparametric finite elements. A n effort is made to present the different A L E virtual work terms in matrix forms that are easy to compute while maintaining the resemblance of the calculations with the 133 Chapter 7. CONCLUSIONS 134 Lagrangian formulation. Implementation of the A L E finite element equations into a 2-D computer code for plane stress, plane strain and axisymmetric problems with 4, 8 and 9 node elements. Special attention is made to the modularity of the code structure. Introduction of a new general relation that relates grid displacements with material displacements. Implementation of this relation on the element level during the assembly of element matrices eliminates grid displacements from the equations and offers a substantial saving in computation time. New treatment for mesh motion on material boundaries that allows boundary nodes to move in the tangential direction to the material boundaries in order to maintain a uniform distribution while satisfying the boundary constraint. - Implementation of time integration schemes that allow for implicit, explicit and mixed implicit-explicit calculations. Implementation of a line search technique to accelerate the convergence of implicit calculations. Employment of a simple contact algorithm based on the direct constraint method. Including contact capabilities into the code allowed the simulation of more applications. - Establishing a clear connection between the developed A L E formulation and the Lagrangian formulation. A L E is shown to be a logical extension to the Updated Lagrangian formulation and the necessary modifications to current Lagrangian codes are clearly identified. The author believes that a Lagrangian-consistent A L E Chapter 7. CONCLUSIONS 135 formulation will help spread the use of A L E in general-purpose commercial finite element codes in large deformation analyses. Simulation of several large deformation quasi-static and dynamic metal forming applications using the developed A L E code and other commercial Lagrangian-based codes. Experimental analysis was carried out for one of the simulated applications. Comparison of A L E predictions with analytical, numerical and experimental results is presented. A L E results are in good agreement with the other methods of analysis. The various simulations clearly show the power of the A L E formulation in preventing mesh distortion under high degrees of deformation. A L E was also employed to eliminate the need for contact analysis and to avoid corner penetration in some applications. 7.2 FUTURE WORK The focus of the work in this thesis has been on the derivation of a consistent A L E formulation for quasi-static and dynamic solid mechanics applications, implementation into a 2-D finite element code and application in the simulation of large deformation and metal forming problems. Future developments may be directed towards the enhancement and extension of the formulation capabilities, improvement of the current code implementation and the simulation of more applications. From the theoretical point of view, the A L E formulation presented in this work may be complemented in several areas: Chapter 7. CONCLUSIONS 136 Investigation of more powerful mesh motion schemes. The transfinite mapping method, employed in this work as the mesh motion scheme, optimizes the mesh based solely on the geometry of the boundaries. Although this scheme is shown to be adequate for many applications, a better mesh can be generally produced if stress localization and strain gradient effects are taken into account when moving the mesh. This is essential if severe stress gradients are present. Isoparametric mapping can also be employed as an efficient mesh motion scheme especially for 3-D analyses. Extension of the formulation into the fluid mechanics area by deriving the virtual work equations for fluid elements. Material models for fluid elements can be introduced into the A L E governing equations following a similar procedure as that used for the constitutive equations of elastic-plastic materials. On the implementation level, desirable improvements to the current code may include: Extension of the analysis dimensions to 3-D. The developed A L E equations are general and are directly applicable in 3-D. However, this extension will require the addition of 3-D brick and shell elements. It will also necessitate the implementation of a 3-D mesh motion scheme and a surface remeshing procedure, which are not simple tasks. Modification of the current contact algorithm or implementation of a more reliable one, possibly based on Lagrange multipliers and/or penalty function. The current contact algorithm causes lots of numerical problems in achieving convergence. It is Chapter 7. CONCLUSIONS 137 essential to use a superior and more stable contact algorithm especially if 3-D calculations are introduced. On the applications side, the use of the developed A L E formulation and code can produce an impact in the simulation of several applications: Fracture mechanics is an application that is well suited for the A L E calculations. A L E offers an effective method for the simulation of crack tip advancement and control of mesh quality near the crack tip compared to the more approximate nodal separation and remeshing techniques used in the Lagrangian simulations. Using the dynamic A L E formulation developed in this thesis, dynamic fracture mechanics, in which high speed loading and high speed crack propagation take place, can be easily simulated. However, additional features to the current code implementation, such as a method of handling of newly created material surfaces and calculation of the J integral and strain energy density, are necessary. Fluid-structure interaction is another class of applications that can be efficiently simulated using A L E and the implemented solution schemes. Fluid mechanics problems are characterized with fine meshes and large material motions while material motion in structural problems is generally more restricted and require a relatively lower mesh density. A L E , with its flexible mesh motion scheme, can efficiently handle the distinct mesh motion requirements that occur in this kind of interaction problems. In addition, since mesh resolutions can be quite different at different regions in the problem, explicit integration for the relatively less dense mesh region (with larger critical time step size) and implicit integration for the more fine Chapter 7. CONCLUSIONS 138 mesh region (with restrictively small critical time step size) can be combined in the analysis leading to savings in computational time. Bibliography [I] R. Hill, "Some Basic Principles in the Mechanics of Solids Without a Natural Time", Journal of the Mechanics and Physics of Solids, Vol. 7, pp. 209-225, 1959. [2] H.D. Hibbitt, P.V. Marcal and J.R. Rice, "A Finite Element Formulation for Problems of Large Strain and Large Displacement", International Journal of Solids and Structures, Vol. 6, pp. 1069-1086, 1970. [3] R . M . McMeeking and J.R. Rice, "Finite Element Formulations for Problems of Large Elastic-Plastic Deformation", International Journal of Solids and Structures, Vol. 11, pp. 601-616, 1975. [4] K.J. Bathe, E. Ramm and E.L. Wilson, "Finite Element Formulations for Large Deformation Dynamic Analysis", International Journal for Numerical Methods in Engineering, Vol. 9, pp. 353-386, 1975. [5] J.C. Nagtegaal and N . Rebelo, "On the Development of a General Purpose Finite Element Program for Analysis of Forming Processes", International Journal for Numerical Methods in Engineering, Vol. 25, pp. 113-131, 1988. [6] ABAQUS User's Manual, Hibbitt, Karlsson, & Sorensen Inc., 1992. [7] NISA II User's Manual, Engineering Mechanics Research Corporation, 1992. [8] ANSYS User's Manual, Swanson Analysis Systems Inc., 1989. [9] P.C. Galbraith and J.O. Hallquist, "Shell Element Formulations in LS-DYNA3D: Their Use in the Modelling of Sheet Metal Forming", Journal of Materials Processing Technology, Vol. 50, pp. 158-167, 1995. [10] M.J. Finn, P.C. Galbraith, L. Wu, J.O. Hallquist, L . Lum and T.L. Lin, "Use of a Coupled Explicit-Implicit Solver for Calculating Springback in Automotive Body Panels", Journal of Materials Processing Technology, Vol. 50, pp. 395-409, 1995. [II] J. Tang, W.T. Wu and J. Walters, "Recent Development and Applications of Finite Element Method in Metal Forming", Journal of Materials Processing Technology, Vol.46, pp. 117-126, 1994. [12] A . Makinouchi and M . Kawka, "Process Simulation in Sheet Metal Forming", Journal of Materials Processing Technology, Vol. 46, pp. 291-307, 1994. [13] M . Karima, "Practical Application of Process Simulation in Stamping", Journal of 139 Bibliography 140 Materials Processing Technology, Vol. 46, pp. 309-320, 1994. [14] W. Kubli and J. Reissner, "Optimization of Sheet Metal Forming Processes Using the Special Purpose Program A U T O F O R M " , Journal of Materials Processing Technology, Vol. 50, pp. 292-305, 1995. [15] L . Baillet, M . Brunet and Y . Berthier, "Experimental and Numerical Dynamic Modelling of Ironing Process", Journal of Materials Processing Technology, Vol. 60, pp. 677-684, 1996. [16] M . Anon, "Numerical Simulation of Sheet Metal Forming Processes", Sheet Metal Industries, Vol. 32, No. 11, pp. 24-30, 1995. [17] J. Cheng and N . Kikuchi, "A Mesh Rezoning Technique for Finite Element Simulations of Metal Forming Processes", International Journal for Numerical Methods in Engineering, Vol. 23, pp. 219-228, 1986. [18] MSC.MARC User's Guide, M S C Software Corporation, 2000. [19] M.S. Gadala, G.A. Oravas and M . A . Dokainish, "A Consistent Eulerian Formulation of Large Deformation Problems in Statics and Dynamics", International Journal of Non-linear Mechanics,Vol. 18, No. 1, pp. 21-35, 1983. [20] K . A . Derbalian, E . H . Lee, R.L. Mallet and R . M . McMeeking, "Finite Element Metal Forming Analysis with Spatially Fixed Mesh", Applications of Numerical Methods to Forming Processes, A S M E , A M D - V o l . 28, pp. 39-47, 1978. [21] N . Brannberg and J. Mackerle, "Finite Element Methods and Material Processing Technology", Engineering Computations, Vol. 11, pp. 413-455, 1994. [22] W.F. Noh, "A Time Dependent Two-Space-Dimensional Coupled EulerianLagrangian Code", Methods in Computational Physics, Vol. 3, pp. 117-179, 1964. [23] C W . Hirt, A . A . Amsden and J.L. Cook, "An Arbitrary Lagrangian-Eulerian Computing Method for All Flow Speeds", Journal of Computational Physics, Vol. 14, pp. 227-253, 1974. [24] T.J.R. Hughes, W.K. Liu and T.K. Zimmermann, "Lagrangian-Eulerian Finite Element Formulation for Incompressible Viscous Flows", Computer Methods in Applied Mechanics and Engineering, Vol. 29, pp. 329-349, 1981. [25] J.M. Kennedy and T.B. Belytschko, "Theory and Application of a Finite Element Method for Arbitrary Lagrangian-Eulerian Fluids and Structures", Nuclear Engineering and Design, Vol. 68, pp. 129-146, 1981. Bibliography 141 [26] J. Huetink, "Analysis of Metal Forming Processes Based on a Combined EulerianLagrangian Finite Element Formulation", International Conference on Numerical Methods in Industrial Forming Processes, pp. 501-509, 1982. [27] J. Huetink, P.T. Vreede and J. van der Lugt, "Progess in Mixed EulerianLagrangian Finite Element Simulation of Forming Processes", International Journal for Numerical Methods in Engineering, Vol. 30, pp. 1441-1457, 1990. [28] P.J.G. Schreurs, F.E. Veldpaus and W . A . M . Brekelmans, "An Arbitrary-EulerianLagrangian Finite Element Model for the Simulation of Geometrical Non-Linear Hyper-Elastic and Elasto-Plastic Deformation Processes", International Conference on Numerical Methods in Industrial Forming Processes, pp. 491-500, 1982. [29] P.J.G. Schreurs, F.E. Veldpaus and W . A . M . Brekelmans, "Simulation of Forming Processes Using the Arbitrary Eulerian-Lagrangian Formulation", Computer Methods in Applied Mechanics and Engineering, Vol. 58, pp. 19-36, 1986. [30] R.B. Haber, "A Mixed Eulerian-Lagrangian Displacement Model for LargeDeformation Analysis in Solid Mechanics", Computer Methods in Applied Mechanics and Engineering, Vol. 43, pp. 277-292, 1984. [31] H . M . Koh and R.B. Haber, "Elastodynamic Formulation of the Eulerian-Lagrangian Kinematic Description", Journal of Applied Mechanics, Vol. 53, pp. 839-845, 1986. [32] H . M . Koh, H.S. Lee and R.B. Haber, "Dynamic Crack Propagation Using EulerianLagrangian Kinematic Descriptions", Computational Mechanics, Vol. 3, pp. 141155, 1988. [33] W.K. Liu, T. Belytschko and H . Chang, "An Arbitrary Lagrangian-Eulerian Finite Element Method for Path-Dependent Materials", Computer Methods in Applied Mechanics and Engineering, Vol. 58, pp. 227-245, 1986. [34] W.K. Liu, H . Chang, T. Belytschko, and J.S. Chen, "Arbitrary Lagrangian-Eulerian Stress Update Procedures for Forming Simulations", Advances in Inelastic Analysis, A S M E - A M D , Vol. 88, pp. 153-175, 1987. [35] W.K. Liu, H . Chang, J.S. Chen and T. Belytschko, "Arbitrary Lagrangian-Eulerian Petrov-Galerkin Finite Elements for Nonlinear Continua", Computer Methods in Applied Mechanics and Engineering, Vol. 68, pp. 259-310, 1988. [36] D.J. Benson, "An Efficient, Accurate, Simple A L E Method for Nonlinear Finite Element Programs", Computer Methods in Applied Mechanics and Engineering, Vol. 72, pp. 305-350, 1989. [37] A . Huerta and F. Casadei, "New A L E Applications in Nonlinear Fast Transient Solid Dynamics", Engineering Computations, Vol. 11, pp. 317-345, 1994. Bibliography 142 [38] S. Ghosh and N . Kikuchi, "An Arbitrary Lagrangian-Eulerian Finite Element Method for Large Deformation Analysis of Elastic-Viscoplastic Solids", Computer Methods in Applied Mechanics and Engineering, Vol. 86, pp. 127-188, 1991. [39] J. Wang, Arbitrary Lagrangian-Eulerian Method and Its Applications in Solid Mechanics, PhD Dissertation, The University of British Columbia, January 1998. [40] J. Tirosh and D. Iddan, "The Dynamics of Fast Metal Forming Processes", Journal of Mechanics and Physics of Solids, Vol. 42, pp. 611-628, 1994. [41] S. Choudhry and J. K. Lee, "Dynamic Plane-Strain Finite Element Simulation of Industrial Sheet-Metal Forming Processes", International Journal of Mechanical Sciences, Vol. 36, No. 3, pp. 189-207, 1994. [42] D.A. Schoch, "Critical Factors for Achievement of Successful High Speed Metal Forming Production", Journal of Materials Processing Technology, Vol. 46, pp. 409-414, 1994. [43] J.C. Gelin, "Dynamic Loading, Viscoplasticity and Temperature Effects on the Evolution of Damage in Metal Forming Processes", Journal of Materials Processing Technology, Vol. 32, pp. 169-178, 1992. [44] P. Steinmann, C. Miehe and E. Stein, "Fast Transient Dynamic Plane Stress Analysis of Orthotropic Hill-Type Solids at Finite Elastoplastic Strains", International Journal of Solids and Structures, Vol. 33, No. 11, pp. 1543-1562, 1996. [45] K. Kormi, R.A. Efheridge and D.C. Webb, " F E M Simulation of the Static and Dynamic Forming of Circular Plates", Journal of Materials Processing Technology, Vol. 42, pp. 451-462, 1994. [46] J.F. Fontane and J.C. Gelin, "A Finite Element Analysis of High-Speed MetalForming Processes", Annals of the CIRP, Vol. 40, pp. 277-280, 1991. [47] S. Kapinski, "Influence of the Punch Velocity on Deformation of the Material in Deep-drawn Flange", Journal of Materials Processing Technology, Vol. 34, pp. 419-424, 1992. [48] K.J. Bathe, Finite Element Procedures in Engineering Analysis, Prentice-Hall, 1982. [49] A . Makinouchi, "Sheet Metal Forming Simulation in Industry", Journal Materials Processing Technology, Vol. 60, pp. 19-26, 1996. of [50] T. Huo and E. Nakamachi, "Evaluation of the Dynamic Explicit/Elastoviscoplastic 143 Bibliography Finite Element Method in Sheet-Forming Simulation", Journal Processing Technology, Vol. 50, pp. 180-196, 1995. of Materials [51] D. Zhou and R.H. Wagoner, "Development and Application of Sheet-Forming Simulation", Journal of Materials Processing Technology, Vol. 50, pp. 1-16, 1995. [52] N . Rebelo, J.C. Nagtegaal, L . M . Taylor and R. Passmann, "Industrial Application of Implicit and Explicit Finite Element Methods to Forming Processes", Numerical Methods for Simulation of Industrial Metal Forming Processes, A S M E , CED-Vol. 5, pp. 67-76, 1992. [53] R.G. Whirley, B.E. Engelmann and R.W. Logan, "Some Aspects of Sheet Forming Simulation Using Explicit Finite Element Techniques", Numerical Methods for Simulation of Industrial Metal Forming Processes, A S M E , CED-Vol. 5, pp. 77-83, 1992. [54] D.Y. Yang, D.W. Jung, LS. Song, D.J. Yoo and J.H. Lee, "Comparative Investigation into Implicit, Explicit and Iterative Implicit/Explicit Schemes for the Simulation of Sheet-Metal Forming Processes", Journal of Materials Processing Technology, Vol. 50, pp. 39-53, 1995. [55] L. Taylor, J. Cao, A.P. Karafillis and M . C . Boyce, "Numerical Simulations of Sheet Metal Forming", Journal of Materials Processing Technology, Vol. 50, pp. 168179, 1995. [56] H.N. Bayoumi, M.S. Gadala and J. Wang, "Application of the A L E Formulation to Metal Forming Problems", 3 EUROMECH Solid Mechanics Conference, Stockholm, Sweden, August 1997. rd [57] H.N. Bayoumi, M.S. Gadala and J. Wang, "Numerical Simulation of Metal Forming Processes", Proceedings of the 6 International Conference on Numerical Methods in Industrial Forming Processes - NUMIFORM'98, Enschede, The Netherlands, June 1998. th [58] H.N. Bayoumi and M.S. Gadala, "An Implicit Arbitrary Lagrangian Eulerian Formulation for Quasi-static and Dynamic Solid Mechanics Applications", Computer Methods in Applied Mechanics and Engineering, submitted for publication, June 2000. [59] J. Wang and M.S. Gadala, "Formulation and Survey of A L E Method in Nonlinear Solid Mechanics", Finite Elements in Analysis and Design, Vol. 24, pp. 253-269, 1997. [60] M.S. Gadala and J. Wang, "Simulation of Metal Forming Processes with Finite Element Methods", International Journal for Numerical Methods in Engineering, Vol. 44, pp. 1397-1428, 1999. 144 Bibliography [61] M.S. Gadala and J. Wang, "Elasto-Plastic Finite Element Simulation of Rolling and Compression between Wedge-Shaped Dies", Journal of Materials Processing Technology, Vol. 97, pp. 132-147, 2000. [62] M.R. Movahhedy, ALE Simulation of Chip Formation in Orthogonal Metal Cutting Process, PhD Dissertation, The University of British Columbia, January 2000. [63] M.R. Movahhedy, M.S. Gadala and Y. Altintas, "Simulation of Orthogonal Metal Cutting Process Using an Arbitrary Lagrangian-Eulerian Finite Element Method", Journal of Materials Processing Technology, Vol. 103, pp. 267-275, 2000. [64] L . E . Malvern, Introduction to the Mechanics of a Continuous Medium, Prentice Hall, 1969. [65] H.N. Bayoumi and M.S. Gadala, "Finite Element Implementation of the Fully Coupled A L E Formulation", International Journal for Numerical Methods in Engineering, submitted for publication, December 2000. [66] R. Haber, M.S. Shepard, J.F. Abel, R.H. Gallagher and D.P. Greenberg, "A General Two-Dimensional, Graphical Finite Element Preprocessor Utilizing Discrete Transfinite Mappings", International Journal for Numerical Methods in Engineering, Vol. 17, pp. 1015-1044, 1981. [67] H.N. Bayoumi and M.S. Gadala, "Finite Element Analysis of Large Strain Solid Mechanics Problems", Proceedings of the 1 Canadian Conference on Nonlinear Solid Mechanics - CanCNSM'99, Victoria, BC, Canada, pp. 375-384, June 1999. st [68] H.N. Bayoumi and M.S. Gadala, "Simulation of Large Deformation Problems Using the Arbitrary Lagrangian Eulerian Formulation", Proceedings of the European Conference on Computation Mechanics - ECCM'99, Munich, Germany, pp. 1-14, September 1999. [69] T.J.R. Hughes and W.K. Liu, "Implicit-Explicit Finite Elements in Transient Analysis: Stability Theory", Journal of Applied Mechanics, Vol. 45, pp. 371-374, 1978. [70] T.J.R. Hughes and W.K. Liu, "Implicit-Explicit Finite Elements in Transient Analysis: Implementation and Numerical Examples", Journal of Applied Mechanics, Vol. 45, pp. 375-378, 1978. [71] B. Irons and S. Ahmad, Techniques of Finite Elements, J. Wiley & Sons, 1980. [72] H.N. Bayoumi, M.S. Gadala and J. Wang, "Finite Element Frontal Solver", Technical Report, Submitted to Forming Technologies Inc., Oakville, Ontario, Canada, July 1996. Bibliography 145 [73] A. Kamoulakos, "A Simple Benchmark for Impact", Bench Mark, pp. 31-35, 1990. [74] Y . Y . Zhu and S. Cescotto, "Unified and Mixed Formulation of the 4-Node Quadrilateral Elements by Assumed Strain Method: Application to Thermomechanical Problems", International Journal for Numerical Methods in Engineering, Vol. 38, pp. 685-716, 1995. [75] G.T. Camacho and M . Ortiz, "Adaptive Lagrangian Modelling of Ballistic Penetration of Metallic Targets", Computer Methods in Applied Mechanics and Engineering, Vol. 142, pp. 269-301, 1997. [76] E.H. Lee, R.L. Mallet and W.H. Yang, "Stress and Deformation Analysis of the Metal Extrusion Process", Computer Methods in Applied Mechanics and Engineering, Vol. 10, pp. 339-353, 1977. [77] S. Ko, Numerical Comparison of Various Finite Element Formulations in Nonlinear and Plasticity Problems, Master of Engineering Thesis, The University of British Columbia, in preparation. [78] Y . Huang and Y . Lu, "Elasto-Plastic Finite Element Analysis of V-Shape Sheet Bending", Journal of Materials Processing Technology, Vol. 35, pp. 129-150, 1992. [79] Y . Huang and D. Leu, "Finite Element Analysis of Contact Problems for a Sheet Metal Bending Process", Computers and Structures, Vol. 57, pp. 15-27, 1995. [80] G . G . Corbett and S.R. Reid, "Quasi-static and Dynamic Local Loading of Monolithic Simply-Supported Steel Plate", International Journal of Impact Engineering, Vol. 13, pp. 423-441, 1993. Appendix A GRID TIME DERIVATIVE OF VIRTUAL STRAIN The grid time derivative of dp.. is given by d 1 .ddu, ddu, —( H dt 2 d'x, -) L d'x. s d dSuy d d5u j dt d'xj X X +•a d% For clarity, we will use vector notation to derive (A.1) X s s ddu, We can write ——- as d'x. d ddu i a d'xj ddu. ddu d'x, d'x X* (A.2) Defining I to be the identity tensor, we have ddu f ddu' (A.3) d'x Kd'xj Taking the grid time derivative of both sides of (A.3) and rearranging gives d ddu a d'x (ddu\ X s _1 _ ddu d (ddu\ (AA) d'x a {d'x) [d'x) x e Multiplying both sides of (A.4) by —— from the right gives dx d ddu a d'x f in which — ddu' a Kd'xj ddu d (ddu\ X s d'x a {d'xj d'x X can be calculated as X ddu s 146 s (A.5) Appendix A . GRID TIME DERIVATIVE OF VIRTUAL STRAIN d (d8*X a _d U'xJ 147 d'x a ddu X s d d'x i a x g dSu d'y g dSu d^ d'x d'x dSu d'y dSu d'x d'x (A.6) Substituting (A.6) into (A.5) gives d dSu dSu d'y' dSu dSu a d'x d'x d'x d'x d'x dSu d'y g (A.7) d'x d'x Rewriting (A.7) in indicial notation gives d dSiij ddu d'v a d' d'x d'xj g (A.8) j k Xj X s Therefore = F'j 1 dSu^'vf o ^JC 2 ^ d'x d'xj at.. k d5u.d'v g d'x a L . k d'x. a L , ' (A.9)
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Arbitrary Lagrangian-Eulerian formulation for quasi-static...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Arbitrary Lagrangian-Eulerian formulation for quasi-static and dynamic metal forming simulation Bayoumi, Hassan N. 2001
pdf
Page Metadata
Item Metadata
Title | Arbitrary Lagrangian-Eulerian formulation for quasi-static and dynamic metal forming simulation |
Creator |
Bayoumi, Hassan N. |
Date Issued | 2001 |
Description | Many engineering problems involve large material deformation, large boundary motion and continuous changes in boundary conditions. The Arbitrary Lagrangian Eulerian (ALE) formulation has emerged in recent years as a technique that can alleviate many of the shortcomings of the traditional Lagrangian and Eulerian formulations in handling these types of problems. Using the ALE formulation the computational grid need not adhere to the material (Lagrangian) nor be fixed in space (Eulerian) but can be moved arbitrarily. Two distinct techniques are being used to implement the ALE formulation, namely the operator split approach and the fully coupled approach. A survey of the ALE literature shows that the majority of ALE implementations for quasi-static and dynamic analyses are based on the computationally convenient operator split technique. In addition, all previous dynamic ALE formulations are based on explicit time integration where no linearization is needed. This thesis presents a fully coupled implicit ALE formulation for the simulation of quasi-static and dynamic large deformation and metal forming problems. ALE virtual work equations are derived from the basic principles of continuum mechanics. A new method for the treatment of convective terms that sidesteps the computation of the spatial gradients of stresses is used in the derivation. The ALE virtual work equations are discretized using isoparametric finite elements. Full expression for the resulting ALE finite element matrices and vectors are given. A new relation that relates grid displacements with material displacements is introduced. The ALE finite element equations are implemented into a 2-D computer code for plane stress, plane strain and axisymmetric problems. The transfinite mapping method is used as the mesh motion scheme for internal nodes. A new treatment for mesh motion on material boundaries is introduced and implemented. Implicit, explicit and mixed implicitexplicit time integration schemes are implemented in the code. A line search technique is employed to accelerate the convergence of implicit calculations. Several quasi-static and dynamic large deformation applications are solved using the developed code. Experimental analysis of a simple V-bending process is conducted for comparison. Comparison of ALE predictions for deformed shapes, equivalent stress and plastic strain distributions and loading curves with analytical, numerical and experimental results are presented. ALE results are in good agreement with other methods of analysis. ALE is shown to prevent mesh distortion and eliminate the need for special contact treatments for problems with corner contact. |
Extent | 13432969 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-10-05 |
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. |
IsShownAt | 10.14288/1.0080952 |
URI | http://hdl.handle.net/2429/13598 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2001-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2001-714349.pdf [ 12.81MB ]
- Metadata
- JSON: 831-1.0080952.json
- JSON-LD: 831-1.0080952-ld.json
- RDF/XML (Pretty): 831-1.0080952-rdf.xml
- RDF/JSON: 831-1.0080952-rdf.json
- Turtle: 831-1.0080952-turtle.txt
- N-Triples: 831-1.0080952-rdf-ntriples.txt
- Original Record: 831-1.0080952-source.json
- Full Text
- 831-1.0080952-fulltext.txt
- Citation
- 831-1.0080952.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-0080952/manifest