"Applied Science, Faculty of"@en . "Mechanical Engineering, Department of"@en . "DSpace"@en . "UBCV"@en . "Wang, Jin"@en . "2009-06-02T19:42:10Z"@en . "1998"@en . "Doctor of Philosophy - PhD"@en . "University of British Columbia"@en . "The applicability and accuracy of existing finite element formulation methods for finite\r\nstrain deformation and metal forming problems are investigated. It is shown that the existing\r\nformulation methods, both Lagrangian and Eulerian type, are not suitable nor efficient for\r\nlarge deformation problems especially when boundary conditions change during deformation,\r\nas is the case in most metal forming problems. This creates a need for a more general and\r\nefficient type of formulation. An Arbitrary Lagrangian-Eulerian (ALE) method is presented\r\nfor the general application in solid mechanics and large deformation problems. A consistent\r\nALE formulation is developed from the virtual work equation transformed to an arbitrary\r\ncomputational reference configuration. The formulation presents a general approach to ALE\r\nmethod in solid mechanics applications. It includes load correction terms and it is suitable\r\nfor both rate-dependent as well as rate-independent material constitutive laws. The proposed\r\nformulation reduces to both updated Lagrangian and Eulerian formulations as special cases.\r\nThe formulation is presented in a form that makes the programming an extension to existing\r\nLagrangian and Eulerian type programs.\r\nAn efficient mesh motion scheme for ALE formulation is developed with a procedure for\r\nhandling boundary motion within the scheme, which can ensure homogeneous mesh results. A\r\npractical and more efficient numerical method is presented to handle supplementary constraint\r\nequations on element level rather than on the global level. Different numerical algorithms for\r\nthe integration of the rate type constitutive equation are investigated and coupled with the return\r\nmapping algorithm to provide plastic incremental consistency. A numerical procedure for stress\r\nintegration is developed based on the physical meaning of stress. Jaumann and Truesdell rates\r\nare taken as the objective stress rates in the constitutive equation. An alternative numerical treatment for rate of deformation tensor [sup t]D[sub ij] is presented to maintain incremental objectivity\r\nof the tensor. It is shown by numerical examples that the use of Truesdell stress rate with a\r\ndeveloped numerical integration procedure gives consistently more accurate results than other\r\nprocedures presented. An algorithm for updating material associated properties is presented\r\nand applied in simulation of various metal forming problems.\r\nA 2D finite element program, ALEFE, based on the presented formulation is developed and\r\ntested. The program may reduce to an updated Lagrangian or Eulerian methods as special cases.\r\nThe mesh motion for the whole domain is controlled by the motion of the boundary nodes.\r\nThe program can handle unsymmetric stiffness matrices and coupled displacements/velocity\r\nboundary conditions. The input data is designed to be similar to available commercial finite\r\nelement codes, so that the data generation phase may be directly imported from these programs.\r\nThe output data format is designed to be compatible with general graphic simulation and data\r\nprocessing commercial softwares, so that contour, x-y and deformed mesh plots may be easily\r\ncreated from the output data of ALEFE.\r\nVarious benchmark and practical problems are simulated by the developed program. Practical\r\nsimulation cases include flat punch forging process, sheet metal extrusion process, necking\r\nbifurcation of a bar in tension, a steady strip rolling and compression between wedge-shaped\r\ndies. Numerical results are compared with analytical solutions or experimental results available\r\nin literature."@en . "https://circle.library.ubc.ca/rest/handle/2429/8616?expand=metadata"@en . "8790981 bytes"@en . "application/pdf"@en . "A R B I T R A R Y L A G R A N G I A N - E U L E R I A N M E T H O D A N D ITS A P P L I C A T I O N I N S O L I D M E C H A N I C S By Jin Wang B. A. Sc., M . A. Sc. University of Science and Technology Beijing A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF D O C T O R OF PHILOSOPHY in T H E FACULTY OF GRADUATE STUDIES DEPARTMENT OF MECHANICAL ENGINEERING We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA April 1998 \u00C2\u00A9 Jin Wang, 1998 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 The University of British Columbia Vancouver, Canada Date Ajrfd l6> DE-6 (2/88) A b s t r a c t The applicability and accuracy of existing finite element formulation methods for finite strain deformation and metal forming problems are investigated. It is shown that the existing formulation methods, both Lagrangian and Eulerian type, are not suitable nor efficient for large deformation problems especially when boundary conditions change during deformation, as is the case in most metal forming problems. This creates a need for a more general and efficient type of formulation. An Arbitrary Lagrangian-Eulerian (ALE) method is presented for the general application in solid mechanics and large deformation problems. A consistent A L E formulation is developed from the virtual work equation transformed to an arbitrary computational reference configuration. The formulation presents a general approach to A L E method in solid mechanics applications. It includes load correction terms and it is suitable for both rate-dependent as well as rate-independent material constitutive laws. The proposed formulation reduces to both updated Lagrangian and Eulerian formulations as special cases. The formulation is presented in a form that makes the programming an extension to existing Lagrangian and Eulerian type programs. An efficient mesh motion scheme for A L E formulation is developed with a procedure for handling boundary motion within the scheme, which can ensure homogeneous mesh results. A practical and more efficient numerical method is presented to handle supplementary constraint equations on element level rather than on the global level. Different numerical algorithms for the integration of the rate type constitutive equation are investigated and coupled with the return mapping algorithm to provide plastic incremental consistency. A numerical procedure for stress integration is developed based on the physical meaning of stress. Jaumann and Truesdell rates are taken as the objective stress rates in the constitutive equation. An alternative numerical ii treatment for rate of deformation tensor tDij is presented to maintain incremental objectivity of the tensor. It is shown by numerical examples that the use of Truesdell stress rate with a developed numerical integration procedure gives consistently more accurate results than other procedures presented. An algorithm for updating material associated properties is presented and applied in simulation of various metal forming problems. A 2D finite element program, ALEFE, based on the presented formulation is developed and tested. The program may reduce to an updated Lagrangian or Eulerian methods as special cases. The mesh motion for the whole domain is controlled by the motion of the boundary nodes. The program can handle unsymmetric stiffness matrices and coupled displacements/velocity boundary conditions. The input data is designed to be similar to available commercial finite element codes, so that the data generation phase may be directly imported from these programs. The output data format is designed to be compatible with general graphic simulation and data processing commercial softwares, so that contour, x-y and deformed mesh plots may be easily created from the output data of ALEFE. Various benchmark and practical problems are simulated by the developed program. Practi-cal simulation cases include flat punch forging process, sheet metal extrusion process, necking bifurcation of a bar in tension, a steady strip rolling and compression between wedge-shaped dies. Numerical results are compared with analytical solutions or experimental results available in literature. m Table of Contents Abstract ii List of Tables ix List of Figures x Nomenclature xiii Acknowledgement xix 1 INTRODUCTION 1 1.1 NONLINEAR FINITE ELEMENT METHODS 1 1.2 LAGRANGIAN FORMULATIONS 2 1.2.1 Sample Examples 2 1.2.2 Limitation of Lagrangian Type Formulation 12 1.3 EULERIAN FORMULATION 14 1.4 ARBITRARY LAGRANGIAN-EULERIAN METHOD 15 1.5 OBJECTIVES AND SCOPE 16 2 ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 18 2.1 SURVEY OF A L E METHODS 18 2.2 CONSISTENT A L E FORMULATION 24 2.2.1 Geometric and Kinematic Description 24 2.2.2 A L E Formulation 26 iv 2.3 DISCUSSION AND CHARACTERISTICS OF PROPOSED FORMULATION 31 3 N U M E R I C A L I M P L E M E N T A T I O N O F A L E 35 3.1 INTRODUCTION AND BRIEF SURVEY 35 3.1.1 Mesh Motion and Stiffness Matrix Processing 35 3.1.2 Stress Integration 37 3.1.3 Updating Material Associated Properties 38 3.2 STIFFNESS MATRIX PROCESSING 40 3.3 MESH MOTION SCHEME IN A L E 42 3.3.1 Transfinite Mapping Method 42 3.3.2 Consideration of Boundary Motion 45 3.3.3 The Motion of Nodes on Boundaries 48 3.3.4 Implementation and Discussion 49 3.4 INTEGRATION OF STRESS RATE 50 3.4.1 The Integration Algorithm 50 3.4.2 Consistency of the Integration Algorithm 53 3.4.3 Alternative Bases for Integration Schemes 54 3.5 INCREMENTAL OBJECTIVITY AND NUMERICAL TREATMENT OF *D U 56 3.5.1 General Considerations 56 3.5.2 Forward Difference Algorithm to Calculate t Dy 59 3.5.3 Forward Difference with a Modification Term 59 3.6 STRESS TRANSFORMATION AND RETURN MAPPING 61 3.6.1 General Considerations 61 3.6.2 Scheme A: Transformation after Return Mapping 62 3.6.3 Scheme B: Transformation Before Return Mapping 63 3.7 UPDATING OF MATERIAL ASSOCIATED PROPERTIES 64 v 3.7.1 Basic Updating Scheme 64 3.7.2 Material Associated Properties at Nodal Points 65 4 T W O - D I M E N S I O N A L A L E P R O G R A M ( A L E F E ) 69 4.1 PROGRAM A L E F E 69 4.1.1 Features of A L E F E 69 4.1.2 Flow Chart of the Main Program 71 4.2 MAIN SUBROUTINES 74 4.2.1 Functions of Subroutines 74 4.2.2 Flow Chart of Mesh Motion 76 4.3 CONVERGENCE CRITERIA 76 4.3.1 General Discussion 76 4.3.2 Displacement Criterion 78 4.3.3 Force Criterion V 78 4.3.4 Energy Criterion 79 4.4 COUPLED DISPLACEMENT CONDITION 80 5 N U M E R I C A L E X A M P L E S 83 5.1 OVERVIEW 83 5.2 ROTATION OF A N ELASTICALLY DEFORMED ELEMENT 83 5.3 ELASTIC-PLASTIC S M A L L DEFLECTION ANALYSIS OF A SQUARE E L E M E N T 85 5.4 PURE BENDING OF A RECTANGULAR PLATE 89 5.4.1 Problem Description 89 5.4.2 Finite Element Models 90 5.4.3 Results of Lagrangian Formulation 90 5.4.4 A L E Simulation and Result 93 vi 5.5 PUNCH FORGING PROCESS 93 5.5.1 Model and Updated Lagrangian Results 93 5.5.2 A L E Simulation 96 5.5.3 Simulations from an Inhomogeneous Mesh 101 5.6 TENSILE NECKING ANALYSIS OF A N ELASTIC-PLASTIC BAR 105 5.7 SHEET M E T A L EXTRUSION I l l 5.7.1 Extrusion through a Curved Die I l l 5.7.2 Extrusion through a Straight Die 114 5.7.3 Discussion 121 5.8 STRIP ROLLING 122 5.8.1 Particular Difficulties in Rolling Process 122 5.8.2 Deformation Conditions in Rolling 124 5.8.3 Numerical Treatments in the Simulation 125 5.8.4 Simulation of Neutral Point Location 127 5.8.5 Simulation of Rolling Pressures 131 5.8.6 Simulation of Rolling Force and Torque 133 5.8.7 Concluding Remarks 136 5.9 COMPRESSION BETWEEN WEDGE-SHAPED DIES 138 5.9.1 Deformation Conditions 138 5.9.2 Finite Element Model 140 5.9.3 Slip Line Field Solution 141 5.9.4 Comparison of Results 142 5.9.5 Development of Deformed Shapes and Plastic Zones 142 6 C O N C L U S I O N S A N D F U T U R E W O R K 148 6.1 CONCLUSIONS 148 vii 6.2 FUTURE WORK 151 Bibliography 153 Appendices 160 A Material and Referential Time Rate 160 B Updated Lagrangian Formulation 163 C Mapping Between Material and Mesh Domain 165 D Stress Increment Procedure 166 v i i i List o f Tables 5.1 Displacements of top right corner from all simulations 102 5.2 Geometric dimensions of rolling process, (Length unit: mm) 124 5.3 Specimen dimensions in compression between wedge-shaped dies 139 ix List of Figures 1.1 The meshes before and after deformation 3 1.2 Not updated boundary condition 4 1.3 Distribution of stress ax on outer surface 5 1.4 The finite element model 7 1.5 The deformed shapes by updated Lagrangian method 8 1.6 Finite element model used in DEFORM2 simulation 9 1.7 Deformed mesh from DEFORM2 10 1.8 Load-height reduction result from DEFORM2 11 2.1 Schematic diagram for domains and mapping in A L E description . 24 3.1 Transfinite mapping 44 3.2 Mesh modification 45 3.3 The mesh motion in A L E 47 3.4 The material configurations 51 3.5 Rigid body rotation of point Pl 57 4.1 Flow chart of the main program 72 4.2 The flow chart for mesh control 77 4.3 The relation between global and local velocities 80 4.4 The flow chart for frontal solution 82 5.1 Element with combined extension and rigid body rotation 84 5.2 The stress results of single element under rotation 86 x 5.3 Elastic-plastic small deflection model 87 5.4 Load versus displacement in x direction 88 5.5 Pure bending of a rectangular plate 89 5.6 Finite element models for the pure bending of a plate 91 5.7 Pure bending of a plate: Lagrangian Formulation Results 92 5.8 Pure bending of a plate: A L E Results 94 5.9 The current yield stress from A L E 97 5.10 Load-height reduction result from A L E 98 5.11 Contours of ay from A L E 99 5.12 The current yield stress without boundary node adjustment 100 5.13 Inhomogeneous mesh for punch forging simulation 101 5.14 Comparison of current yield stresses from NISA 103 5.15 The A L E result of current yield stress from inhomogeneous mesh 104 5.16 The finite element model for tensile necking 106 5.17 The necking 8 and engineering straining 7 107 5.18 Comparison of two algorithms for Truesdell rate '. 109 5.19 Contours of equivalent stress (Truesdell rate, transformation A) 110 5.20 Original geometry and mesh in extrusion through a curved die 112 5.21 Contours of plastic strain after piston moved 2a 113 5.22 Comparison of ax from A L E and reference [12] 115 5.23 Comparison of rxy from A L E and reference [12] 116 5.24 The finite element model for straight die extrusion 117 5.25 The current yield stress (plastic zone) from A L E 118 5.26 The initial and deformed mesh from NISA 119 5.27 Comparison of extrusion force 120 5.28 Pressure distribution on piston face 121 XI 5.29 The definition of geometry in rolling 123 5.30 Finite element model (Case 1: R/h0 = 12.5, Reduction 14.17%) 125 5.31 Locations of neutral points when velocity criterion is applied 128 5.32 Rolling pressure and friction stress when velocity criterion is applied 129 5.33 Locations of neutral points when force condition is applied 131 5.34 Rolling pressure and friction stress when force criterion is applied 132 5.35 Comparison of rolling pressure, small R/h0 133 5.36 Comparison of rolling pressure, large R/h0 134 5.37 Comparison of rolling force 135 5.38 Comparison of rolling torque 135 5.39 Equivalent plastic strain (Case 1: R/h0 = 12.5, Reduction = 14.17%) 136 5.40 Current yield stress at different steps (Case 5: R/hQ = 39.0, Reduction = 22.83%) 137 5.41 General shape of the specimen and die 139 5.42 The finite element model for Case 1 141 5.43 The variation of vertical load on the die with compressive strain 143 5.44 The variation of mean vertical die pressure with current strip thickness . . . . 144 5.45 Distribution of end velocities 145 5.46 Equivalent plastic strain at different steps 147 A . l The description of relation between mesh motion and material motion 161 xii Nomenclature General Variables: coefficients in relation of material motion and mesh motion area of an element b dimension of a square element Ci coordinates of a point P old coordinates of a point P Cijki material constitutive tensor dp punch displacement Da components of rate of deformation tensor D?- plastic part of Dij t+At^ij components of deformation tensor E Young's modulus f a general function or a physical quantity fc the value of / at computational referential point f? components of body forces ff components of surface tractions /- unsmoothed quantity at node a smoothed property at node a F external force deformation gradient tensor G, Gij the gradient of displacement from time t to t + At h height of the specimen in punch indentation X l l l hu values of shape function for displacements in i direction at D.O.F. I ha smoothing shape function at node a H hardening modulus i index for directions / coordinate axes (i = l,2,3ora:,y,z) / index for D.O.F. I unit matrix J ratio of material volume to the volume in referential state Jc Jacobian matrix k,k0 current and initial shear yield stresses, respectively K , ku stiffness matrix related to 'v K c , kjj stiffness matrix related to 4 v c K equivalent stiffness matrix L length of the elastic-plastic bar m modification term for Dij n number of D.O.F in an element nn number nodes per element rii components of the unit vector normal to boundary surface N total number of D.O.F in the model Ns number of sampling points in an element p normal pressure pi external nodal force at D.O.F. / p^ external nodal force in kth iteration at D.O.F. / p equivalent load vector q distributed tensile force qi equivalent nodal force at D.O.F. / xiv qj^ equivalent nodal force in kth iteration at D.O.F. / v normalised parametric coordinate 0 < r < 1 TI residual nodal force at D.O.F. / residual nodal force in kth at D.O.F. / R , Rij rotation tensor s normalised parametric coordinate 0 < s < 1 S surface of deformed body on which tractions are prescribed Su 2nd Piola Kirchhoff stress Sip 2nd Piola Kirchhoff stress at an intermediate state t time td,tf,te displacement, force and energy tolerances, respectively T transformation matrix Ur total radial displacement Ux total displacement in x direction Ui displacement in i direction ui material displacement at D.O.F. J from time t to t + At u\u00C2\u00B0j mesh point displacement at D.O.F. / from time t to t + At v velocity of punch or die v material particle velocity vector Vi material velocity in direction i vj material velocity at D.O.F. / v \u00C2\u00B0 computational reference velocity vector v 1 mesh velocity in direction i V i ( a ) material velocity of node a at direction i in global coordinates vi(a) material velocity of node a at direction i in local coordinates xv v\ mesh velocity at D.O.F. J v ' material particle velocity vector in local coordinate system lV volume of deformed body at time t w diameter of the circular punch in punch indentation W diameter of the circular specimen in punch indentation W, Wij spin tensor x position vector of material coordinate x * position vector of imaginary material coordinate xi material particle coordinate in axis i a coefficient in implicit integration algorithm f3 contact angle in rolling Sij Kronecker delta At time increment Awj 5' displacement increment in kth iteration at D.O.F. / A x / change of physical quantity due to material motion A objective function in the least square method 4>i{r,s) boundary coordinate of a geometric domain X position vector of referential coordinate Xi mesh point coordinate in axis i u> angular velocity V / gradient of / Left Superscripts: t the quantity is the value at time t t + At the quantity is the value at time t + At Left Subscripts: t the quantity is with respect to the configuration at time t Superimposed Symbols: xvii material time rate of a quantity A computational reference time rate of a quantity T transformation of a matrix xvi i i Acknowledgement I would like to express my deep appreciation to my supervisor, Dr. Mohamed S. Gadala, for his invaluable help, judicious advice and understanding throughout my graduate studies. His patient guidance and continuous encouragement effectively helped me to proceed with this work. I express my gratitude and appreciation to Dr. Henry Vaughan for his invaluable suggestions and advice. I would also like to thank Dr. Stanley G. Hutton and Dr. Reza Vaziri for their help and valuable discussions and suggestions at various stages of this work. I acknowledge the financial support of Natural Science and Engineering Research Council of Canada, of UBC-UGF award and, finally, the financial support from my supervisor. Special thanks are also due to my wife Chunmei and my parents who offered incessant moral support and encouragement. To them and to my new born son, Roy, I dedicate this thesis. x i x Chapter 1 I N T R O D U C T I O N 1.1 N O N L I N E A R F I N I T E E L E M E N T M E T H O D S Although the application of finite element methods to many nonlinear problems has been successfully carried out, there are numerous areas of concern and investigation in that regard. One such area, which we consider here, deals with the solution of large strain plasticity and metal forming problems with pseudo type boundary nonlinearities. A number of critical difficulties arise in the finite element analysis of such problems. Among these difficulties are; proper formulation, mesh distortion, proper modelling of the contact boundary conditions, incorporation of the plastic incompressibility condition and accounting for plastic anisotropy. Foundations of large strain analysis of elastic-plastic solids may be traced back to the early work of Hill [1]. It took some time, however, until Hibbitt, Marcal and Rice [2] introduced the first finite element formulation for large strain problems. In their approach they used a total Lagrangian formulation (TLF). Later on, McMeeking and Rice [3] pioneered the use of updated Lagrangian formulation (ULF) in the same area of applications. The two formulation methods have been widely used for both steady and non steady static large plastic strain problems. On the other hand, Eulerian formulation (EF) has been initially introduced for finite element applications in the fluid mechanics area. Several trials aiming at adapting this formulation to large strain and metal forming problems were attempted [4, 5, 6, 7]. Owing to the difficulties in obtaining material time derivatives in spatial reference frame, no generally accepted Eulerian formulation is available for such problems. Also, since the mesh is spatially fixed in EF, it is 1 Chapter 1. INTRODUCTION 2 not easy to simulate non-steady static or dynamic behavior. Trials have been made to relate the moving material points (in terms of F E Gauss points) to the fixed spatial mesh [6, 7]. Much work is still required, however, to refine and establish this approach. This chapter investigates the applicability and accuracy of existing formulation methods to the finite strain deformation problems. The objective is to examine some of the above aspects through the use of available FE commercial programs. The programs used in the analysis are DEFORM2 [8], ANSYS [9] and NISA [10]. It is required to provide an assessment for the existing formulation of such problems in order to propose required modifications or new formulation. The applications considered in this regard are flat punch indentation, plane strain extrusion and forging problems. 1.2 L A G R A N G I A N F O R M U L A T I O N S 1.2.1 Sample Examples (i) Flat Punch Indentation Flat punch indentation into a slab is mathematically the same as the compression of a specimen between two fiat parallel punches[ll]. In the axisymmetric compression model, w and W are the diameters of a circular punch and solid circular specimen, respectively and h is the height of the specimen. The ratios W/w = 2.7 and h/w = 1.7 are used in the simulation. The material properties are Young's modulus E - 6.9 x 104 MPa, Poisson's ratio v = 0.33, hardening modulus H = 138 MPa with initial yield stress o-yieido = 90 MPa. Axi-symmetric 8-node isoparametric elements together with 2-node (node-to-node) gap elements are used in the model. Assuming a smooth punch, an indentation is made up to a punch displacement dp = 0.15w. The simulation was carried out using both NISA and ANSYS programs. The mesh topology before and after deformation is presented in Figure 1.1. The radial displacement UT, the second principle stress cr2, and the equivalent plastic strain epeq are used for comparison of the Chapter 1. INTRODUCTION 3 results. The results obtained from the two programs showed close agreement in displacement and strains but revealed noticeable discrepancies when comparing stress values. Also, there is noticeable disagreement between the results from both programs and similar ones reported in the literature [11]. w/2 i \u00E2\u0080\u00A24-Punch Size Increasement ~~1 i / \. W/2 Figure 1.1: The meshes before and after deformation (ii) Metal Extrusion Another problem analyzed with NIS A and ANS YS is the steady state deformation process of plane-strain metal extrusion. Metal is forced into a symmetric die after sliding between smooth rigid plates. The die is approximated with a straight line and produces a 25% thickness reduction over a distance 1.2a, where a is half the thickness of the original sheet, as shown by grey lines in Figure 1.2. The material properties are Young's modulus E = 6.9 x 104 MP a, Poisson's ratio v = 0.3, strain hardening i f = 1.1 x 103 MPa and initial yield stress o-yieid0 = 4.0 x 102 MP a. The metal extrusion problem is modelled by plane strain 8-node isoparametric elements and Chapter 1. INTRODUCTION 4 Figure 1.2: Not updated boundary condition 2-node gap elements. The nodal and deformed shape after the piston moves 0.3a are shown in Figure 1.2 by grey and black lines, respectively. The distributions of displacement Ux and stress ax, both in the extrusion direction, are examined. The distributions of Ux and ax from NISA and ANSYS are very similar and the values of Ux are close to each other with a difference of less than 0.4%. The values of ax are, however, quite different; and show a difference of up to 25%. Once again, the results obtained here are quite different from those reported in the literature [12]. For frictionless die-workpiece interface, the results with and without gap elements, in both NISA and ANSYS, are significantly different. As an example, distribution of ax on the outer surface of the work-piece is shown in Figure 1.3. It is not expected that the introduction of node-to-node gap elements should influence the frictionless die case. The obtained results show, however, a noticeable effect of adding gap elements to the frictionless die case. Chapter 1. INTRODUCTION 5 Ox(MPa) 3.5E2 0.0E0 -3.5E2 -7.0E2 --10.5E2 - \u00E2\u0080\u0094 \u00E2\u0080\u00A2 \u00E2\u0080\u0094 Ox(MPa) No G a p Element \u00E2\u0080\u0094 A \u00E2\u0080\u0094 Ox(MPa) G a p Element -14.0E2 , i , . , , , ^ . . i . . , . i . , . . i . . . . i . . . . i -2.0 -1.5 -1.0 -0.5 0 0.5 1.0 \u00E2\u0084\u00A2a. Figure 1.3: Distribution of stress Figure 1.6: Fini te element model used in D E F O R M 2 simulation Chapter 1. INTRODUCTION 10 the predicted deformed boundary of the free surface near the punch corner is unusual and not rational. Another shortcoming is the fluctuations in loads predicted by Lagrangian method. Figure 1.8 shows a plot for the reduction of the specimen height versus the applied load, from DEFORM2. The result shows noticeable load fluctuation that is pertinent to the updated Lagrangian formulation. Similar fluctuation is observed in the plane strain metal extrusion example. The extrusion pressure increases steadily up to the state when the billet fills the die. As the billet exits the die, the interface nodes that are constrained to move on the die are released. This causes a drop in the strain energy of the billet which in turn causes a drop in the extrusion pressure. The studies performed by Lee [12] and Voyiadjis [13] confirm the existence of these fluctuations. (mm) 20.7 H 10.9 ~\ 1.0 -i -8.8 -i \u00E2\u0080\u0094 j . . r . r . | . . T . r . r T . . r r . y . T . . r r . \u00C2\u00A5 . T . T . r ^ -8.7 2.9 14.5 26.0 37.6 (mm) Figure 1.7: Deformed mesh from DEFORM2 Chapter 1. INTRODUCTION 11 2.200 1.760 1.320 o T3 0.440 0.00 j A T I i \ j ; f A / ,. . j . . / 1 j / r ! / ' -\u00E2\u0080\u00A2 \ 4,-1-' : 0.0 0.12 0.24 0.36 R e d u c t i o n 0.48 0.60 Figure 1.8: Load-height reduction result from DEFORM2 Chapter 1. INTRODUCTION 12 1.2.2 Limitation of Lagrangian Type Formulation Most of the available commercial FE programs are based on Lagrangian type of formulation. Such formulations are efficient and quite suitable for handling nonlinear problems in which small strains prevail, boundary condition nonlinearities do not change with the course of deformation and where mesh distortion is not a critical factor in the analysis. Unfortunately, all of the above restrictions are violated in any metal forming process. The above examples reveal some of the main shortcomings in the application of Lagrangian type of formulations to metal forming problems. It is shown that large plastic strains with contact boundary conditions are not suitable to be analyzed by such formulations. Contact boundary conditions with sharp edges or corners can not be applied precisely and load fluctuations are evident with large strains. The problem of mesh distortions and element entanglement pose a serious drawback on the use of such formulation. It is not always feasible to update the mesh manually, and even if it is possible to do so, it will involve major interaction and time involvement by the user of such schemes. On the other hand, some automatic mesh rezoning methods are developed for Lagrangian finite element analyses of metal forming problems [14,15,16]. However, these methods are not so robust or efficient to remedy the mesh distortions. Furthermore, they do not eliminate load fluctuations. These remeshing methods are based on rediscretizing the deformed configuration after certain special deformation. Obviously, the determination of the special deformation is very important and can only be based on user's experience. The remesh has to be applied on the whole domain, which is actually not necessary. In general, these methods can not change the mesh topologies, so they may not improve the overall accuracy of the finite element calculations. The locations of nodes in the new mesh are calculated by algebraic interpolation method or by solving differential equations. The algebraic interpolation will introduce curve-fitting error, especially when the boundary shape is complicated. Solving differential equation can generate Chapter 1. INTRODUCTION 13 higher quality finite element meshes, but it is time consuming and, more significantly, may still produce poorly shaped meshes near boundaries in some cases [17]. The problems of contact and friction boundary conditions pose another serious concern in the solution of metal forming applications. Contact condition, which is defined from the geometrical compatibility on the contact surface or the impenetrability condition, has been introduced by the direct method, elements with special mechanical properties, the Lagrangian multiplier method and the penalty function method [18]. Since the contact area is a prior unknown, the boundary conditions of the contact problems are determined as part of the solution. Although the contact surfaces obtained by the Lagrangian multiplier method satisfy the contact condition in the integral sense, additional unknowns are required and the total number of unknowns in the system equations increases. In addition, the associated tangent matrix may be indefinite and has zero diagonal entries that pose some difficulties in the solution steps [19]. For the penalty function method, the solution results satisfy the contact conditions only approximately. The accuracy of the approximate solution depends strongly on the penalty parameter, so the correct choice for these parameters is the essence of the algorithm. When the penalty parameter is chosen to be too large, it leads to numerical problems in the form of loss of accuracy in the solution. On the other hand, a too small choice, results in unacceptable penetration of one body into the other. The node to node or node to surface gap or contact elements have been recently introduced in commercial FE packages [9, 10]. Such an approach is quite natural to the formulation of the F E method because the overall stiffness assembly process is unchanged; the contact elements can be treated as a separate material property group such as plasticity, creep, etc. In spite of their various advantages, these contact or gap elements exhibit several important shortcomings [20]. They undergo anomalous response behavior when employed in situations where large deformation kinematics are needed to generate closure, that is reflected in improper Chapter 1. INTRODUCTION 14 stiffness characterizations, i.e., poorly conditioned Jacobians. They are also awkward to apply in situations requiring friction effects and they require significant amount of equilibrium iterations. Furthermore, parameters of the contact element have to be chosen by the user, and the convergence speed and the accuracy of the solution are quite dependent on the choice of these parameters. Sometimes, especially when friction effect is incorporated, convergence may be difficult to achieve, as indicated by the experience from solving the plane strain forging problem. 1.3 E U L E R I A N F O R M U L A T I O N The main difference between the Eulerian formulation and Lagrangian ones is that the deformation of the material moving through a fixed region in space is determined as a function of the current position and the time t instead of determining the deformation of the material element by following its motion in space [6]. The independent variables in the Eulerian description are current position x of the body-point X and time t. For the material motion, x itself becomes dependent on the time t, which complicates material time derivatives and other relations when approached strictly by Eulerian, or spatial formulation. Although many authors have treated the Lagrangian and the updated Lagrangian formulation, there has been little effort concerning the development of a consistent Eulerian formulation in solid mechanics applications. Some developments, e.g. Gadala et al. [6], give the proper material time derivative in spatial description but original geometric parameters are included in the final equilibrium equation. Therefore, such approach may not be strictly considered as purely Eulerian formulation. In the Eulerian method, it is necessary to properly determine the state of the material-associated properties (e.g., strain, stress) of the material points momentarily occupying the Chapter 1. INTRODUCTION 15 integration points at the beginning of each step. The method of updating such material-associated properties is not well developed and further investigations need to be considered. Abo-Elkhier [21] utilized an imaginary finite element mesh for updating, but no discussion nor proper assessment of the accuracy of the method is given. Derbalian [22] discussed a conjugate scheme for interpolation of the stress and applied the procedure to simulate steady state static metal extrusion. By constructing a set of conjugate shape functions (bi-orthogonal), a consistent approximation for the stress field may be obtained which is continuous across inter-element boundaries and involves less mean error. However, this method introduces a large system of equations which may be solved iteratively during each incremental step. Also, the conjugate method is originally proposed for improving the accuracy of incremental stress fields, rather than the total stress [23]. Eulerian formulation is generally more suitable for the study of flow problems in a fixed region of known shape. Therefore, the analysis of steady state metal forming processes may be achieved by Eulerian formulation. Because it introduces other difficulties like appropriate representation of free body, it is less suited for domains whose boundaries or interfaces move substantially and it is not easy to simulate non-steady static or dynamic behavior within the frame work of this formulation. 1.4 A R B I T R A R Y L A G R A N G I A N - E U L E R I A N M E T H O D From the above discussion, it is believed that a new formulation type which combines the advantages of the Lagrangian and Eulerian methods is essential for the accurate simulation of metal forming processes [24, 25, 26]. This new type of formulation is called Arbitrary Lagrangian Eulerian (ALE) formulation. The key point in differentiating the A L E formulation from Lagrangian or Eulerian type formulation is that in A L E a reference computation domain that can move arbitrarily and independently from the material is introduced. The movement Chapter 1. INTRODUCTION 16 of the reference domain is represented by a set of grid points, which may be interpreted as the movement of a finite element mesh. Therefore, in an A L E formulation, the finite element mesh need not adhere to the material or be fixed in space but may be moved arbitrarily relative to the material. A proper A L E formulation should reduce to Lagrangian formulation if we choose to use the same motion for the computational and material meshes. On the other hand, if we choose to fix the computational mesh, an A L E formulation should reduce to Eulerian formulation. Combining the merits of both Lagrangian and Eulerian formulation, A L E is more suitable to handle mesh distortion and entanglement and special boundary condition changes in metal forming problems. If the nodes on tool-workpiece interface are specified as Eulerian points in tangential direction, it may eliminate load fluctuations, may describe precisely any contact boundary conditions and make boundary condition updating no longer necessary after each incremental step. 1.5 O B J E C T I V E S A N D S C O P E A L E method has the potential to eliminate problems caused by Lagrangian or Eulerian methods in simulation of general finite deformation problems. Although the concept of the A L E method was first proposed in mid-seventies [27], its use in solid mechanics problems has been restricted mainly because of the additional effort required in satisfying the deformation history-dependent stress-strain relationship and updating complexities. Relatively little attempts have been made to generalize and present A L E formulation in a general way and to provide a clear connection between A L E and Lagrangian as well as Eulerian approach. New aspects pertinent to the arbitrary mesh motion in A L E have not been investigated intensively. The objective and steps of this research may be summarized in the following: - Survey of existing A L E schemes with emphasis on their pertinent characteristics, reason er 1. INTRODUCTION 17 for differences and their applicability to large deformation metal forming problems. Develop a complete A L E formulation that may be applied to general finite deformation cases and degenerate to either updated Lagrangian or Eulerian formulation as extreme cases [28]. Emphasise on the generality of the developed formulation to use various material laws, stress rates and on its special application to metal forming problems. Investigate an efficient motion scheme for the A L E method and a numerical algorithm to process the supplementary equations arising from mesh motion. Existing schemes tend to introduce higher degree of complexity in the implementation and require more CPU processing time. Study the effect of various stress integration methods on objectivity and plastic consis-tency of results and improve the numerical algorithms. Investigate the combination of stress updating algorithms with the usual return mapping schemes for small deformation. Also investigate approximate numerical integration procedures to ensure objectivity of stress rates and special tensors used in the formulation. Investigate a practical and more efficient method to update material associated properties for A L E approach. This is a key point in the formulation because of the arbitrary mesh motion in A L E . The scheme should be based on proper continuum mechanics equations rather than on numerical interpolation and approximation. Implement the developed formulation and methods to develop a general 2D finite element program and to simulate bench mark problems and metal forming processes to show the features of A L E method. Chapter 2 A R B I T R A R Y L A G R A N G I A N E U L E R I A N ( A L E ) F O R M U L A T I O N 2.1 S U R V E Y O F A L E M E T H O D S The concept of A L E was first proposed in the mid-seventies [27] under the name of \"coupled Eulerian-Lagrangian\", and used in conjunction with a finite difference scheme to solve two-dimensional hydrodynamics problems with moving fluid boundaries. Later, the A L E method was introduced into the finite element method [29] from fluid mechanics for modelling the solid-fluid interaction and free-surface problems, and introduced in [30] to solve nonlinear analysis of nuclear safety. Since then, A L E method has mainly been used in fluid and linear-path independent solids, where stress states are solely determined by the instantaneous displacement or velocity fields [31, 32]. Only recently the A L E method has been applied to finite strain deformation problems in solid mechanics [24, 25, 28]. Huetink [33,34] introduced a finite element method to simulate metal forming process under the name of \"combined Eulerian-Lagrangian formulation\". In his formulation, the material rate of change of the equation of virtual power, with changing integration area, was used to derive the final equation. However, in the final discretized equilibrium equations, only material velocities were included. As in the updated Lagrangian method, the material velocities have to be calculated from the final equations, then the new finite element mesh can be updated by arbitrarily moving the old finite element mesh in the material domain. On the boundaries, however, the mesh can only move in the tangential direction of the material domain. Following this, is the updating of the material-associated properties. In this process, material motion and 18 Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 19 mesh motion do not occur at the same time and the discretized equilibrium equations do not include simultaneously mesh velocities and material velocities as unknowns. In this regard, Huetink formulation may be considered as an updated Lagrangian formulation (ULF) coupled with remeshing after each incremental step. Benson [35] described a method called simple A L E . It seems that the method was developed mainly for fluid mechanics. The method was introduced to reduce the computer time of large scale analyses involving finite deformation through the use of explicit finite element schemes. The time step size of the explicit program is, therefore, limited by element mesh dimensions. In order to modify the mesh, the A L E method was implemented and the equations were derived by substituting the relationship between the material time derivative and reference configuration time derivative into the governing equations for a continuum in a Lagrangian coordinate system. An operator split was advocated to decouple the equations so that for every increment, two steps were needed. First, a Lagrangian step was performed, in which the mesh moved with the material during the step. The solution was then mapped from the Lagrangian mesh to the reference one in order to complete the A L E step. This method may be also considered as an updated Lagrangian formulation coupled with remeshing after increment. An improved version of A L E formulation was presented by Haber [36]. In Haber's work two displacement variables were considered to be primary unknowns, i.e., the Lagrangian and Eulerian type of displacements. This allowed some distinction between the material and the computational motion. The new technique was applied to large-deformation frictional contact and fracture mechanics. It was demonstrated that the mixed Eulerian-Lagrangian description (i.e., ALE) may be used to vary the crack length in a continuous way and may easily handle contact boundary condition. The paper does not, however, clearly indicate which displacement stands for mesh motion and there is no discussion about how to move the computational finite element mesh. Furthermore, one of the most demanding aspects in A L E , the updating of Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 20 material-associated properties was not addressed. It is also assumed that the formulation is limited to linear elastic material and the load is independent of displacement, therefore, the final formulation may not be directly applicable to non-linear material and deformation-dependent loads. Similar to the above work, Yamada introduced two types of displacements called the material and spatial incremental displacements and derived an A L E formulation for plane deformation of hyperelastic material [37]. For such kind of materials, the stress can be uniquely determined from the strain energy density and it is independent of strain histories. Thus, no state quantity of the particles needs to be introduced except the deformation gradient tensor. The application of formulation to only path-independent hyperelastic material models eliminated the need for rigorously addressing the problem of calculating and updating material associated properties, and makes the formulation not applicable to the general solid mechanics deformation problems. An A L E formulation specifically derived for solid mechanics was described by Schreurs [38, 39]. In [39], the difference between the CRS (Computational Reference System) derivative and the MRS (Material Reference System) derivative of a physical quantity was clearly stated. Starting from the equilibrium equation in material domain, the weak form was set up by the principle of weighted residuals. The final discretized equilibrium equations were obtained by transforming integration over material domain to a reference domain. However, the physical meaning of the reference domain was not clearly addressed. A computational mesh moving method was created by assuming the grid point as a material point of a linear isotropic fictitious body. Assuming one shape of a stress-free element being optimal, linear set of algebraic equations were set up to have the deformed mesh elements recover to the optimal shape. By solving the equations, the grid nodal displacements were determined. This method is time consuming and does not guarantee an optimal mesh because the real material is generally elastic-plastic. Follower points or cells are used to represent the material associated properties and deformation history. These properties at new integration points decided by CRS points Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 21 were updated by extrapolation and interpolations from the follower points. Hughes [40] elaborated on some basic concepts related to A L E . An important one is the relationship between the material time derivative and referential time derivative * / A of a physical quantity */. These should satisfy the following relation: '/\u00E2\u0080\u00A2 = 7A+ (2.1) (J X{ where txi is material coordinate, tvi and tv^ are material and mesh velocities individually. Liu [41] applied Equation (2.1) to the momentum equation and obtained an equation with respect to arbitrary reference volume. Petrov-Galerkin formulation was then utilized to set up the final discretized equilibrium equation. A stress updating procedure was developed to calculate the stress values at the quadrature and nodal points. The constitutive law was refor-mulated by introducing a stress-velocity product; the Petrov-Galerkin finite element method was used to set up the equivalent weak form equations of the constitutive law. These equations were solved simultaneously with the final equilibrium equations to get the solution. With this scheme, some wave propagation problems and elastic-plastic dynamic deformation processes have been simulated. In the same paper, Liu also transformed the integration extending over the material domain to referential domain and derived the referential time derivative of internal virtual work. In a later development [42], Liu considered the frictional interface of plane de-formation problems and derived equilibrium equations for such case. The Laplace differential equation and fourth order differential equation were proposed as mesh generator to manage mesh movement. The application of proposed A L E algorithm to metal rolling simulation was presented. It was shown that for small rollers, both Lagrangian and A L E meshes are feasible and the agreement between the two is quite good. For realistic roller size, a Lagrangian mesh failed to complete the simulation whereas the A L E mesh performed quite well. Extended from fluid mechanics, an A L E formulation for solid mechanics was reported by Ghosh et al. [24, 43, 44]. Different from the methods discussed previously, the Reynolds Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 22 transport theorem widely used in fluid mechanics was applied to an arbitrary moving control volume to obtain the field equations of mechanics with respect to an arbitrarily moving grid point. The weak form of the differential equations were obtained by using appropriate weighting functions and integrating over the current grid volume. In Ghosh's work it was assumed that the motion of the material with respect to the grid is quasi-static, i.e., | = 0. In the opinion of the author, this assumption is not warranted. In quasi-static problems, it is the material point acceleration, as expressed by Equation (2.2), that may be neglected instead of the motion of the material with respect to the grid points, because the term contributes only partially to the material point acceleration. dt d*vi dt \u00E2\u0080\u00A2x 3 On the other hand, the exclusive characteristics of the A L E , that the grid points can move arbitrarily, is restricted so much by this assumption. Using implicit time integration scheme, Ghosh integrated the weak form equation over the current control volume. This makes the calculation laborious because, unlike transient field problems, the integration volume here is changing with time. In [43] Ghosh described a way to update variables to nodal points of arbitrary motion. Pseudo-material elements were constructed first and material-associated variables were evaluated by interpolation. Local algebraic and elliptic mesh generator was implemented to perform mesh management in highly local deformation areas [44]. Although Ghosh indicated that the specific cases of A L E formulation should be either updated Lagrangian or Eulerian formulation, no verification is presented in his work showing that the formulation may reduce to these special cases. A particular A L E implementation; ALE3D, in metal forming simulations is recently in-troduced by Cough et al. [25]. The basic computational cycle consists of a Lagrangian step followed by an advection step. At the end of Lagrangian phase of the cycle the velocities and nodal positions are updated. At this point, the user has several options. If the user opts Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 23 to run a pure Lagrangian mode, no further action is taken and the code proceeds to the next time step. If a pure Eulerian calculation is desired, the nodes are placed back in their original positions. The user has options available to tailor the evolution of the mesh in order to maximize either efficiency or accuracy. The mesh updating scheme implemented in ALE3D is a finite element-based equipotential method. For the constitutive law, the Jaumann rate is used for the stress tensor and von Mises yield condition is applied. In the paper [25], only program features are described and no formulation or rigorous algorithm is given. An overall description of A L E is presented in [45] by Huerta and Casadei. It is clearly indicated that the most important challenge for the A L E technique lies in its extension to solid mechanics problems in general, and, in particular, to non-linear solid mechanics where path dependent material behavior is fairly common. It is pointed out that the best choice for the mesh motion or velocities and a low cost algorithm for updating the material-related properties constitute the major problems. However, no particular schemes are given or presented. Some primary governing issues, e.g., the conservation laws, constitutive equations and boundary conditions, are presented, but no complete formulation is given. The applications considered in the paper are some elasto-plastic problems, with concentration on fast-transient solid dynamics showing the effectiveness of A L E to impact problems when the explicit integration scheme is used. Although several forms of A L E schemes have been discussed in the literature, it is noted, however, that these formulations always concentrate on certain aspects of the problem which makes the outcome formulation incomplete and suitable only for specific applications. Crucial points in the formulation, e.g., the final expression of equilibrium equations in A L E frame, the relation between computational and material points, evaluation of material related data, are not clearly addressed in the literature and need further investigation. Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 24 2 . 2 CONSISTENT ALE FORMULATION 2 . 2 . 1 Geometric and Kinematic Description t Boundary of Computat ional Domain F X 3 , lX3 Figure 2.1: Schematic diagram for domains and mapping in A L E description Assuming a material particle P^xif x2* x3) in a material reference system (MRS) at time t has velocity '^(t = 1,2,3). The particle moves to Q(t+Atx1,t+At x2,t+At x3) at time t + At. Similarly, a grid point P\u00C2\u00B0(txi,t X2? X3) m the corresponding computational reference system (CRS) at time t has velocity = 1,2,3) and moves to Qc{t+AtXi,t+At X2,t+At Xs) at time t + At, as shown in Figure 2.1. In A L E description, the motion of the finite element grid need not adhere to the material particle and may be controlled by the finite element user in accordance with users judgement. However, we should ensure a one-to-one mapping between Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 2 5 the material domain and computational domain at any time t, i.e., for each: 'xi = xi(txutX2,tX3) (i = 1 , 2 , 3 ) ( 2 . 3 ) we have only one: *Xi = Xi{txi,tx2,tx3) (i = 1 , 2 , 3 ) The necessary and sufficient condition for the inverse relation to exist, is that the determinant of the Jacobian transformation, ( 2 . 4 ) | d*Xi is non-vanishing. A chosen mesh motion scheme should satisfy the above requirement. In addition, at every stage of motion, the boundary of the deformed configuration of the material body (solid line in Figure 2.1) must coincide with the boundary of the configuration of the computational reference domain (dashed line in Figure 2.1), i.e., on boundary: ('vi-'vtfm = 0 ( 2 . 5 ) where trii is the component of the unit vector normal to the boundary. The physical interpretation of Equation (2.5) is that no normal convective velocity occurs across the boundary if the surface particles remain on the surface. It is then possible to derive the governing equations with respect to the referential domain. Since txi is our computational reference coordinate system, it is necessary to express the material time derivative of a function / in terms of the time derivative with respect to tXi-Assume that the material time derivative and the computational reference time derivative are denoted by a superscripted dot \"\u00E2\u0080\u00A2\" and a cap \" A \" , respectively. The relationship between the two derivatives is given by Equation (2.1) (See Appendix A for details). Equation (2.1) is important since it makes it possible to track the material deformation history when A L E is used. Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 26 2.2.2 A L E Formulation Many problems in solid mechanics may be treated as quasi-static, so that the continuity equation can be satisfied automatically. Starting from the principle of virtual work at time t + At, the equilibrium equation may be expressed as: I \"\"vijSt+tteijf+^V = [ t+Atf?8uidt+AtV+ [ \"\"ffSuJ+^S (2.6) Jt+Aty J Jt+Aty J 1 Jt+AtS J i v ' Notation similar to those used by Bathe [46,47] are adopted. Referring to a general motion of a body in a fixed Cartesian coordinate system, as shown in Figure 2.1, the left superscripts indicate the configuration at which the quantity occurs, and the left subscripts indicate the configuration with respect to which the quantity is measured. The left subscript may not be used if the quantity under consideration occurs in the same configuration in which it is measured. A quantity with no left superscript or subscript indicates an increment from time t to t + At. Therefore, t+Ataij in Equation 2.6 is the Cartesian component of Cauchy stress referred to the configuration at time t + At. The deformation tensor t+Ateij is defined by: T + A T 6 I J = 2{d^^j + d^x') ( 2 ' 7 ) where are admissible incremental displacement vector and 8 in Equation (2.6) means a variation. The body force and the surface traction at time t + At are given by t+Atff and t+Atf? in Equation (2.6), respectively. t + A t V is the volume of deformed body and t + A t S is the surface of deformed body at time t + At on which surface tractions are prescribed. Utilizing the following relations: / t + A t 0-^+^0*+AtV = j \u00C2\u00AB \" a * ( i ^ + J ^ \ d \u00C2\u00BB \" V Jt+Aty J J Jt+Aty 2 \ O t + A t X j O t + A t X i J t + A t y A V d t + A t x a t+AtrS _ t+At_ t+At Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 27 where dt+AtAj is the projection of a^ + A t5 on the coordinate plane having axis lXj or t + A t X j as normal direction, t+Atrij are the direction cosines of area element dt+AtS with axis lXj or t + A t X j , the Equation (2.6) may be written in the following form: / T + A T ^ ^ ~ ^ + A T V = I t + A t / J W + A 'V + / t+Ataij8uidt+AtAj (2.8) Jt+AtV Qt+^Xj Jt+AtV Jt+&tA. J J \ / In order to solve Equation (2.8), all quantities are transformed into the known computational reference configuration at time t, which is particularly chosen as the material reference domain at time t, i.e., making Pc and P in Figure 2.1 be the same point. This transformation is only due to the change in or motion of the computational reference domain. We have, dSuj _ dSuj dlxk d'+^xj ~ dtxkfr+L'xj [ ' ' t+Atxk = 'xk + ul (2.10) where u% is the mesh point or CRS displacement from time t to t + At. Substituting in Equation (2.9) and reordering, dSui dSui dSui duk d t + A t X j ~ ~\u00C2\u00A5x~ = ~ dtxk 8t+Atxj dividing the above equation by At and let A i \u00E2\u0080\u0094> 0, we obtain, / d5u{ \ A dSui dtvl dtXj J dtxk dtXj which gives, dSui 88 dt+AtXj d*Xj dSui t x otxk dlXj dlXj The volume element dtJrAtV is related to d*V\t through: At The computational reference time rate of the area is [48]: 1 otxk otXj so that: dt+AtAj = c?AJt + ( ^ l d t A j - ^ d i A k ) A t (2.18) Using Equations (2.13), (2.14), (2.17) and (2.18) and considering the relations lx{ = the right hand side of Equation (2.8) may be written as: RHS = [ t+Atf?8uidt+AtV+ [ t+Ataij8uidt+AtAj Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 30 V -^ + f \u00C2\u00A3^d%- - ^ i(r, 0), (f>i(r, 1), c&(0, s) and & (1, s) as shown schematically in Figure 3.1, the mesh coordinate Cj is given by: d(r,s) = (l-s)i(r,0) + si(r,l) + (l-r)(f)i(0,s)+ri(l,s) _ (1 _ r ) ( i _ a)^(0,0) - (1 - r)si{0,1) - rsi{l, 1) - r ( l - * ) & ( l , 0 ) (3.6) 0 < 7\" < 1 0 < 5 < 1 where r, s are normalized coordinates, i = x,y, i.e., (j>x(r, 0) and y(r, 0) etc., represent the x,y coordinates of boundary curves, and cx(r, s), cy(r,s) are nodal x and y coordinates of a point, respectively. The transfinite mapping Equation (3.6) will create a homogeneous mesh when the nodes on the boundaries are distributed homogeneously, and it will allow the boundaries to be represented in a discrete form, e.g., with nodal coordinates on the boundaries. This is a characteristic that is quite convenient and useful in A L E formulation because at each step of simulation, the deformed boundaries are actually described by discrete nodes. The only requirement of the discrete Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 44 (0,D W r , 1 ) (0,1) (1.1) i(0,s) (U\u00C2\u00B0#) ) 4>,(i.s) (0,0) ^(r.O) (1,0) (^i(r.O)) (0,0) *,(r,0) (0,1) Original Region Modified Region Figure 3.1: Transfinite mapping representations in the transfinite method is to have normalized parametric coordinates assigned to the ith point on a curve containing (k + 1) points. The location of internal nodal points are completely determined by the position of boundary nodes. The boundary curves may be described by as many nodes as needed, so that curve fitting error in parametric mapping can be minimized or avoided. In general, we assume that at a given incremental step, at time t, the finite element mesh is distorted. In the next incremental step from time t to t + A i , the mesh will be automatically moved by the program to get a homogeneous one. The new internal nodal coordinates are calculated by Equation (3.6). If the old (distorted) coordinates are c\u00C2\u00B0(r, s), then the mesh velocity is given by: to each node on boundary curves. This process is automated by assigning the coordinate Cj(r,s) - c?(r,a) A i t c* (3.7) (3.8) Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 45 Figure 3.2 gives an example for distorted and automatically modified mesh obtained by applying the above method. It should be noted that, even when the boundaries have discontinuous slopes, for example at A, B and C in Figure 3.2, the method produced a homogeneous mesh. Distorted Mesh at Time t Modified Mesh at Time t+At Figure 3.2: Mesh modification 3.3.2 Consideration of Boundary Motion The above mesh motion scheme will guarantee a homogeneous mesh when boundaries do not move. In a general A L E formulation, however, mesh motion and loading increment occur simultaneously and the boundaries should be allowed to move. In such case, in order to achieve a higher quality mesh, Equation (3.8), has to be modified. First we assume that the boundary velocities are expressed by \u00C2\u00A3;(r, 0),\u00C2\u00A3i(r, l ) , \u00C2\u00A3 ; ( 0 , s ) and s) corresponding to boundaries Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 46 i(r, 0),4>i(r, l),4>i(0,s) and i(l,s) individually, as shown in Figure 3.1. We use the same function as in Equation (3.6) to interpolate the boundary velocities and find the modification term to be used in Equation (3.8); *vf = (l-s)Ci(r,0) + sU(r,l) + (l-r)Ci(0,s) + r^(l,s) - (1 - r)(l - s)$i(0,0) - (1 - r)sZi{0,1) - \u00E2\u0084\u00A2&(1, 1) - r ( l - a ) & ( l , 0 ) (3.9) 0 < r < 1 0 < s < 1 To ensure that the boundary nodes have to be Lagrangian type in the normal direction, the mesh velocities \u00C2\u00A3;(r, 0), \u00C2\u00A3;(r, 1), \u00C2\u00A3;(0, s) and &(1, s) have to be equal to the material velocities tVi(r, 0)/ Vi(r, l),41>;(0, s), and *VJ(1, s); or at least their normal components have to be equal to each other. These material velocities at time t , however, may be unknown, for example as on free surfaces. In order to avoid solving implicit equations and to calculate the normal components, we apply the material velocities at time (t \u00E2\u0080\u0094 At) to approximate the values at time t during the mesh motion. It should be indicated that such assumption is only for mesh motion and does not have any direct effect on the accuracy of solution. Therefore, for Lagrangian boundary points; d(r,0) = t-Atvi(r,0) & ( r , l ) = \u00C2\u00AB - A ' t ; i ( r , l ) 6(0, s) = ' - A t M 0 , *) 6(1,\u00C2\u00AB) = *\" A t f i ( l , *) (3-10) (i = x,y 0 < r < 1 0 < s < l) If higher quality mesh is required, some iterations may be introduced, i.e., using the new material velocities to replace the velocities at time (t\u00E2\u0080\u0094At) in Equation (3.10) until the difference between the material velocities from the two iterations is less than a given tolerance value. In practical applications, this iteration was found to be not necessary since the mesh quality created by Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 47 Equation (3.10) is good enough for most of practical problems. Finally the mesh motion for the nodes in the interior of the domains is given by: t C t c t i t c b C o 1 1 \ vi = vi + vi = ai (3.ll) For the example shown in Figure 3.2, if the distorted mesh at time t has boundaries DE and BE moving at some inhomogeneous speed, using the above mesh motion scheme, we get the mesh at time t + A i as shown in Figure 3.3. So, even when boundaries move at inhomogeneous speeds, which is the usual case in finite element large deformation problems, the scheme may still give homogeneous meshes. Figure 3.3: The mesh motion in A L E Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 48 3.3.3 T h e M o t i o n o f N o d e s on B o u n d a r i e s The boundaries of deformation domain may be classified as one of two categories. The first kind posses no change in boundary shape during the deformation and the shape is known a prior. An example to this type is a boundary with prescribed velocity or displacement conditions, e.g., the tool work-piece interface in indentation or punch forging problems. In the second kind, the deformed boundary shape is unknown in advance and changes during deformation. The boundaries with applied traction forces often lie in this category. Motion of nodes on the first type of boundaries is easy to decide. Because the normal direction of boundaries is known, the mesh velocity in the normal direction should be equal to material velocity, i.e., a>i = 0, bj = 1 in Equation (3.3). In the tangential direction, the node may be moved arbitrarily in order to achieve calculation efficiency. Without losing generality, bj is set to zero and aj may be determined in a way specified by the user, or usually in such a way that the nodes are evenly distributed along the boundaries. On the second type of boundaries, alternative methods have to be employed due to unknown velocity in the direction normal to the boundary. One way of handling this is to transform the second type of boundaries into first type by assuming the normal direction. The disadvantage here is that iterations have to be introduced to make the difference between the assumed normal direction and the calculated one from equilibrium equation less than the preset tolerance. Another method specifies the nodes on the second type of boundaries as Lagrangian points in all directions, i.e., aj \u00E2\u0080\u0094 0 and 6/ = 1 in Equation (3.3). This may result in unexpected or inhomogeneous spacing because of the deformation. In order to have homogeneous mesh pattern, the boundary nodes have to be adjusted in the tangential direction of the boundary. This kind of adjustment can not be done during the incremental step and has to be done after the deformed boundary shape is determined. In this work, the nodes are relocated on the boundaries according to nodal spacing specified by user or are positioned to get an equal distance between \ Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 49 each node. 3 . 3 . 4 I m p l e m e n t a t i o n a n d D i s c u s s i o n The procedures of mesh motion for an incremental step from time t to time t + At may be summarized in the following points: - Get mesh coordinates at time t and the boundary velocities at time t \u00E2\u0080\u0094 At, - From the input data, identify the Lagrangian and Eulerian direction for each boundary node and decide the coefficients \"aj\" and \"6j\" for each D.O.F., - Calculate the value \"aj\" for each internal node according to the boundary coordinates and velocities, using Equations (3.6)-(3.11), - Create the element stiffness matrices, and 'fcjj, - Perform the transformation given by Equation (3.4) to get the equivalent stiffness matrix and load vector, - Assemble the equivalent element stiffness matrices and load vector and solve the assem-bled equations to obtain material velocities, displacements, stresses, etc., - Adjust the position of nodes on the second type boundaries, - Update the material associated properties and mesh coordinates. The application of transfinite mapping method in A L E for mesh motion has many advan-tages. In comparison with other methods discussed above, the transfinite mapping is efficient, simpler to implement and does not introduce curve fitting errors beyond the discretization error that is inherent to the finite element method. The method allows the boundary curves to be represented in discrete form and to have discontinuous slopes. These properties make the Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 50 method directly applicable to domain with irregular polygon boundaries. In transfinite mapping method, the interior node locations are calculated from boundary nodes. The curvatures of the boundary curves and spacing of boundary nodes are accurately reflected in the mesh, as may be shown in Figure 3.3. Furthermore, the boundary curves do not need to be simple functions. Cusps and inflection points provide no special problems. The node motion scheme used in the research for the second type of boundary is simple and more accurate than assuming a normal direction. It is not limited by boundary curve types and it may be combined with transfinite mapping method to handle the mesh motion in the whole deformation domain. Node adjustment in tangential direction on the second kind boundaries is necessary to keep a higher quality mesh. Stiffness matrix processing discussed above is a general method. That may be used with various material models, i.e., rate dependent or independent. 3.4 I N T E G R A T I O N O F S T R E S S R A T E 3.4.1 The Integration Algorithm In large deformation analysis, because of the fact that the configuration of the body is changing continuously, the Cauchy stress can not be integrated simply by adding the stresses and their increments due to deformation directly [47]. The consideration of configuration changing can be incorporated by purely mathematical tensor transformation as given in [58], or alternatively by utilizing the physical definition of Cauchy stress as presented in the following. Assume the material configuration at time t is *C and the corresponding Cauchy stress is Vjj . At time (i + Ai) , the configuration is '+ A'C and the Cauchy stress is t+Ato-ij. As shown in Figure 3.4, the 2nd Piola-Kirchhoff stress at (i + Ai) with respect to the configuration at time Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 51 t, t ^ i j , is related to Cauchy stress l+^crij, at time t + At by [47]: dt+AtXi dtxr *t+At, t *\u00E2\u0080\u00A2 dt+Atz 1 <9'a;n (3.12) where l+AtFi:j = atg^Xi is the deformation gradient tensor or the Jacoban matrix. *3 P ( ' x i ) Configuration at timet ( C Q( T + A TX|) Configuration at time t+At t + A t C Figure 3.4: The material configurations Consider the deformation to be decomposed of two steps; first one is only a rigid body motion and the second one is only straining. Then after the first step, the stress relation may be written as: t+AtJD _ 1 dt+AtXit+Ato(D dt+A% dtxr dlxr (3.13) and since this step is only a rigid body motion, then, t+AtS(I) = t s = t Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 52 where t+Ato-\p and l+AtS^n are Cauchy and 2nd Piolar-Kirchhoff stress at an imaginary intermediate state between t and t + At, respectively, so that Equation (3.13) may be written in the form: i.e., this step only considers the effect of configuration change on Cauchy stress tensor. In the second step, due to straining, the Cauchy stress would have an increment of Atrij = dju'DuAt (3.14) and the stress at time (t + At) would be: t+At 0, gives: \u00E2\u0080\u0094 cr. duk N + CijuDu i ^ \" j t \u00E2\u0080\u009E ^ \" f t t , , /-\u00C2\u00BB t n dlxk dlxk ki (3.17) Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 54 where 'u; is the velocity component in the ith direction of the coordinates. By definition, the Truesdell rate of Cauchy stress may be given by Equation (2.23). Therefore, we may write Equation (3.17) in the following form: The above simple development indicates that the stress integration algorithm derived from the physical definition of the Cauchy stress is actually the same as the integration of Truesdell stress rate equation. 3.4.3 Alternative Bases for Integration Schemes (i) Integration of Truesdell Stress Rate The material configuration at time t is denoted by ' C and at time t + At by t + A t C . A material point in both configuration may be given by the position vectors P(txi,tX2,tx3) and Q(t+Atxi,t+Atx2,t+Atx3), as shown in Figure 3.4, respectively. Tensor quantities in both configurations lC and *+At(7 may be related to each other by proper transformation. Starting from this point, Pinsky et al. [58] derived the implicit integration algorithms for Truesdell stress rate, Jaumann rate and Green-Naghdi rate. The derivation is purely mathematical and according to their development for Truesdell rate, the algorithm is generally implicit and may be expressed as follow: (3.18) * + A t F dtxm \u00C2\u00B0 At dt+AtXi 1 dt+AtXi dt+AtXj + dt+AtXj (3.19) mnkl t+aAt t+aAt d t + a A t X . n where, t+aAt X i = at+AtXi + (1 - afxi Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 55 t+At dt+A X{ F-7. t + a A t \" 1 \u00C2\u00BBJ Qt+cAt Xj t n _ I ( dtyi , when a = 0, it will be explicit and may be expressed as: t + A t _ i &+*xit d^A% i a t + A t ^ A _ a t + A t x j t + A t p A c r m n ^ - \u00E2\u0080\u0094 ^ (3.20) where, A c r m n is defined as in Equation (3.14). Upon comparing Equation (3.20) to Equa-tion (3.15), we find that the only difference is in the second term of the right hand side and if the higher order term of Ui is neglected, the difference will vanish. Numerical experiments show, however, that both equations give almost identical results. In addition, it is necessary to indicate that by assuming that the first step in the deformation is straining and the second is a rigid body motion, as shown in Appendix D, we would get the exact same result as Equation (3.20). (ii) Integration of Jaumann Stress Rate An integration algorithm for Jaumann rate is proposed by Hughes [56]. The purpose of the development in Hughes's work is to keep the \"objectivity\" of the integration. The derivation is based on purely geometrical concepts and the result may be expressed by : t+A^- = r A t ^ > \u00E2\u0084\u00A2 r A t % + . (3.2i) where, ^^R- is the rotation tensor for the increment from time t to t + A i ; ACTJJ is the stress increment given by: AAt) u2 = txxsin^At) + tx^cos^At) - 1) (3.24) Taking derivatives with respect to At and calculating the limit values when At \u00E2\u0080\u0094\u00C2\u00A5 0, gives the velocities at time t as: t t t t V\ = \u00E2\u0080\u0094UJ X2 V2 = UJ X\ From the definition of lDij, the components of the rate of deformation tensor *D at time t is given by: '\"\u00C2\u00BB = | \u00C2\u00A3 = \u00C2\u00B0 <3-25> \u00E2\u0080\u00A2 D a i (p. + \u00C2\u00A3 \u00E2\u0080\u00A2 ) = o 2 \atx2 a*aji/ This verifies that the analytical values of the rate of deformation tensor are identically zero for the case of rigid body rotation, i.e., the tensor 'D is objective with respect to rigid body rotation. Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 59 3.5.2 Forward Difference Algorithm to Calculate *Djj In numerical and finite element analysis, it may not be possible to get the analytical values as indicated above and, in general, a time-discretization procedure is applied to calculate tDij. The conventional method in finite element analysis is to use the average velocity from time t to time t + At as velocity at time t, and employ this to calculate tDij as the rate of deformation tensor at time t. The calculated value of tD{j is also used as the average value of the tensor from time t to time t + At. Mathematically, this is equivalent to a forward difference method. Thus, the velocities are given by; t u\ tXi{cos{ujAt) \u00E2\u0080\u0094 1) \u00E2\u0080\u0094 tx2sin(u>At) V l = A t = A l t u2 tx1sin(u}At) + tx2(cos(u> At) \u00E2\u0080\u0094 1) *2 = At = A l ( 3 - 2 6 ) the tDij components are then obtained by taking derivatives with respect to txi and tx2 as follow; . cos(uAt)-l r>w \u00E2\u0080\u0094 -At t cosjwAt) - 1 \u00C2\u00B0 2 2 ~ A l 4 Z) 1 2 = 0 (3.27) which shows that f D u and tD22 are not equal to zero except when the time step is very small. A finite time increment may cause errors in the calculation of the components of the tensor 'D, especially when the numerical values of these components are of the same order as the right hand side of Equation (3.27). 3.5.3 Forward Difference with a Modification Term In order to remedy the above problem, a modification term is needed [61]. From Equa-tion (3.27), it may be seen that the modification term should be a function of the rigid body Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 60 rotation u;At or, in general, be a function of the components of the spin tensor 'W. In general large deformation , the skew-symmetric spin tensor is defined by; , 3 - 2 8 ) and the modification term for lDij components for general 2D large deformation process may be considered as: *Ji - (WuAty - 1 A i The above modification guarantees the objectivity of lDij in the case of rigid body rotation. With the modification term, numerical calculation of the lDij components are given by: where, Vi - (w 1 2Ai ) 2 - 1 vMHf^-S)^) 2 - 1 A i A i (3.30) and 8ij is the Kronecker delta. For the general rigid body rotation case given above, the spin tensor lWij is calculated in a way similar to Equation (3.27) and the components are given by: *Wn = 0 lW22 = 0 Substituting the above values in Equation (3.30) and (3.29), we have: 'Dn = 0 i2 = 0 Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 61 Therefore, the explicit integration algorithm, Equation (3.16), will satisfy incremental objectivity if is modified as shown above. The *Dij modification presented above may be considered as an alternative way to keep the incremental objectivity of integration algorithms developed by Pinsky [58] and Hughes [56]. The above modification is practically useful when explicit integration algorithms are utilized. More detailed numerical examples on'the effectiveness of the above procedures is given in Chapter 5. 3.6 S T R E S S T R A N S F O R M A T I O N A N D R E T U R N M A P P I N G 3.6.1 General Considerations It may be seen from the above development that in order to keep objectivity of the integration algorithms, the stress at a given time may be first transformed with an appropriate transformation tensor and then it should be added to the stress increment due to pure straining. For integrating the constitutive equation, the algorithm also has to satisfy the material constitutive relations (e.g., Prandtl-Reuss equations) and the incremental plastic consistency. The return mapping algorithm [59,60] is a numerical method used mostly in small deforma-tion problems to satisfy the material constitutive relations and incremental plastic consistency. The method is generally based on two steps. The first one is to calculate an elastic predictor, which obtains the stresses at the end of the increment from the use of the elastic stress-strain relations. The second step is to subsequently map the obtained stresses onto a suitably updated yield surface and restore plastic consistency. This mapping step may be further divided into two steps; one explicit by projecting along the initial plastic flow direction and the other implicit by projecting along the updated plastic flow direction. Although some research has been done on objectivity of stress integration algorithms under finite deformation cases and the algorithms to keep incremental plastic consistency in small deformation cases are well developed, not too much effort has been made to satisfy incremental plastic consistency in finite deformation Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 62 conditions. In the case of large deformation, although the constitutive relation is not the same as for small deformation, the elastic predictor may still be calculated by assuming that the deformation in the current increment is completely elastic. Then the predictor is mapped onto a suitably updated yield surface using the appropriate large strain constitutive relation, thus restoring plastic consistency. For von-Mises material the equivalent stress at time t is defined by: = (3-31) where, ta'i - = tcrij \u00E2\u0080\u0094 | \u00C2\u00A3 ; / c r f c f c is the deviatoric stress and Sij is Kronecker delta. The accumulated equivalent plastic strain is calculated by [64]: \u00C2\u00A3**=Cif0^0^dr (3-32) where; t0 is the initial time, t is the current time of the deformation and TD^is the plastic part of the rate of deformation tensor. For rate-independent materials, the time t is trivial. An important point to be discussed in large deformation analysis is whether the return mapping should be performed before or after the stress transformation. If the return mapping is to be performed before the stress transformation, we should maintain the plastic consistency after the transformation. To investigate this, two different numerical schemes are presented and compared. 3.6.2 Scheme A: Transformation after Return Mapping If at time t, the converged stress is ';Ai; then, calculate Cijki by assuming the increment to be completely elastic; (iii) calculate the intermediate stress state la\p by the use of Equations (3.15), (3.20) or (3.21); (iv) calculate the \"elastic\" stress increments by Acr^ - = t+Atcr[p \u00E2\u0080\u0094 V^-; (v) use 'cr^- and Aaij as the input stress components to regular return mapping algorithm and obtain the mapped stress t + A t c r ; j ; (vi) calculate the deformation history-dependent parameters, such as current yield stress <7 y j e / d and plastic strain epeq. Although this method projects the final stress point on the yield surface, the \"elastic\" stress Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 64 increments include the incremental stress due to proper transformation, i.e. taking into account of the change in configuration. The return mapping procedure has to start with taij instead of t(rlp , because all history dependent parameters such as the current yield stress and the plastic strain are only updated to the state at time t. The above numerical algorithms for keeping plastic incremental consistency, have been combined with the integration algorithms presented in Section 3.4.1. It is found that the combination of integration of Truesdell rate with Scheme A generally gives better results than other cases. Details of this development and calculation are given in Chapter 5. 3.7 UPDATING OF MATERIAL ASSOCIATED PROPERTIES 3.7.1 Basic Updating Scheme In an A L E formulation, the relation between mesh motion and material motion at each point is given by: u\ = di + (3.33) where, ui is material displacement in i direction, ul is mesh displacement in i direction, and no summation on i is observed. Therefore, location of new mesh points and integration points may be easily traced if isoparametric elements are used. For any physical quantity / , the change due to material motion A x / and the change due to mesh motion Axf may be related to each other by (see Appendix A for details): gt+At f Axf = A , / + -\u00C2\u00AB0 (3.34) where t+Atf = /( t + A i x, t + At) is the value of / at material point at time t + At. The material increment A x / may be obtained as in the regular finite element calculations. If a material particle *x at time t moves to t + A t x at time t + At, then: t+A7 = 7 + A x / (3.35) Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 65 Similarly, when mesh point l \ at time t moves to t+Atx at time t + A i , at t + A t x , we have: t+Atfc = tfc + A x f ( 3 > 3 6 ) But at time i , *x = *x, so that * / c = */ = /( 4 x,i) . Therefore, Equations (3.34)-(3.36) may have the form: \u00C2\u00BB u r = ' f + A , / + ^ K - \u00C2\u00AB . ) Equation (3.37) effectively updates any physical quantity / from a material point to a mesh point. The feature that the relation between mesh motion and material motion is known prior is utilized. This procedures effectively eliminates iterations and interpolations or extrapolations used by other researchers [39, 43]. It is important to indicate that the above method may be utilized in updating for Lagrangian and Eulerian method. In updated Lagrangian method, the iteration to detect the relation between new mesh point and old mesh point may still be necessary. However, the updating may be done by Equation (3.37) if the old mesh points are treated as material point in A L E and new mesh points are treated as mesh points in A L E . For Eulerian method, the updating may be carried directly if the imaginary elements are taken as material elements in A L E and mesh points in Eulerian method as mesh points in A L E [7] . 3.7.2 Material Associated Properties at Nodal Points The material associated properties are available at Gaussian integration points by the con-ventional incremental finite element method and the updating scheme presented in the previous section. Nodal values are important, however, for output consideration and for A L E calculation (e.g., Equation (2.20)). In order to obtain these nodal values, local smoothing procedure [62] is applied. Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 66 In this method, assuming the smoothed property / \u00E2\u0080\u009E and the smoothing shape function ha at nodal point a, then: 7 = h j a (3.38) where / is the smoothed property and convention of summation about a is applied. If the unsmoothed property, calculated directly from deformation history is / , a local smoothing equation, by least square method, is applied to minimize the following equation on element level: 4> = J Ja (hafa ~ f) dxdy whereAe is the area of the element and the following condition should be satisfied: d = 0, (3 1,2, ..nn (3.39) (3.40) where, nn is the number of nodes per element. So the smoothed properties may be obtained from the expression: j jf (hah0dxdy) Ja = j jf hp f dxdy (3.41) The above equation may be evaluated using numerical integration procedures. Since the / values at integration points are already known, the values fa can be obtained, if the smoothing function is specified. The values from different elements sharing the same node are averaged and taken as the final property value at the node. For rectangular and parallelogram elements, if ha is a bilinear shape function and a 2 x 2 Gaussian integration rule is adopted in evaluating Equation (3.41), it can be shown that the smoothed nodal properties fi,f2,fz, fi may be expressed as: (3.42) 7i ' ' 2 1 2 1 _ i / I 2 1 2 // 72 1 2 1 + f _ 1 2 1 V3 2 hi < 73 1 _ }/3 2 1 2 ' 2 1 2 fin 1 2 1 _ ^ 2 1 2 2 . Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 67 where, / / , fn, fm and fiv are the unsmoothed values at the Gauss points. For general element shapes and shape functions, solving Equation (3.41) may be costly. To overcome this problem, discrete smoothing is introduced as described in the following. The function in Equation (3.39) is replaced by: = E ( (W.) \K - f \K)2dxdy (3-43) K=l where Ns is the number of sampling points within the element. For 4-node quadrilateral elements, taking the sampling points to be same as the 2 x 2 Gauss points and using bilinear shape functions for ha, the following expression may be obtained, (3.44) ' fx ' < [ ^K=I (h^hi) \K \u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2 ^K=I (^4^4) \K 1 ( J 4 ) K ^K=I yu j \K In the above equation, hi \K,...,h4 \ K are values of shape functions at Gaussian points I, II, III and IV, f \ K (K = I, II,..., IV) is the unsmoothed properties at Gaussian points. Solving for fx, f2, f3 and / 4 , the above equation reduces to Equation (3.42). In this case, it should be noted that f = hafa is an exact least square fit to the unsmoothed stress values fj, fn, fin and fiv for any shape of element. The procedures of updating material associated properties may be summarized as following steps: - Input the properties at material integration points t+Atf, Ui, u\ and t + A t X i = tX{ + ui, - If Ui = ul for all D.O.F. in the element, go to next step; otherwise, calculate the values at mesh integration points using Equation (3.37), then let t+Atf = t+Atf; - Take the values of *+A'/ at different integration points I, II, III and IV, as the values fi, fn, fm and fiv, in Equation (3.42), and get smoothed variables flf f2, f3 and / 4 for nodes 1,2,3 and 4 individually. Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 68 - Add up the values from different elements sharing the same node and divided by the number of elements sharing the node to get the final property values at the mesh node. Chapter 4 TWO-DIMENSIONAL ALE PROGRAM (ALEFE) 4.1 PROGRAM A L E F E 4.1.1 Features of A L E F E Based on the formulation and the numerical algorithms presented in previous chapters, a general 2D finite element program A L E F E is developed with emphasis of application to metal forming problems. The basic features of the developed program may be summarized in the following: - Three formulation options are available: general A L E , updated Lagrangian and Eulerian formulation. The user control parameters are aj and bi in the relation between mesh motion and material motion, *v7 = aj + (1 = 1,..., JV) (4.1) where, no summation on I is observed in the above equation. tvI and tv\ are material and mesh velocities at D.O.F. I respectively, and N is the total number of D.O.F. in the model. If ai \u00E2\u0080\u0094 bj = 0 for all I, the program reduces to Eulerian one; If aj = 0, 6/ = 1 for all I, the program reduces to updated Lagrangian one; otherwise, the program will use a general A L E formulation. - The program incorporates both Jaumann stress rate and Truesdell stress rate in stiffness calculation. The corresponding stress integration algorithms are employed in computing stress. 69 Chapter 4. TWO-DIMENSIONAL ALE PROGRAM (ALEFE) 70 - For general A L E case, the mesh motion for interior nodes is controlled by the nodal motion on boundaries. User only needs to choose the aj and 6j values in Equation (4.1) for boundary nodes. The internal nodes are specified by the program automatically to keep a homogeneous mesh. This is a particularly important features for users and for practical use of the program. It is relatively easy and straight forward to identify the formulation type required, and hence the coefficient a/, bi, for the boundary nodes. - The program handles unsymmetric stiffness matrices and coupled displacements/velocity . boundary conditions. The direct methods for coupled D.O.F. is applied so that no approximation is introduced in this aspect. - Convergence criteria include displacement, residual force and iterative energy. User may choose one, two or all criteria for a given problem. - Rate dependent and rate independent material models are included in the program. For plastic deformation, multi-linear hardening models can be chosen. - Type of elements available is four node quadrilateral element. - The input data format is designed similar to available commercial finite element codes, so that the data generation phase may be directly imported from these programs. - The output data format is designed to be compatible with general graphic simulation and data processing commercial softwares, so that contour, x-y and deformed mesh plots may be easily created from the output data file of A L E F E . - Full and modified Newton-Raphson methods are available. Chapter 4. TWO-DIMENSIONAL ALE PROGRAM (ALEFE) 71 4.1.2 Flow Chart of the Main Program To discuss applications of numerical algorithms in program A L E F E , a general flow chart is shown in Figure (4.1). Two main flow schemes have been tried in the structure of program A L E F E . Both schemes have been implemented in the program, numerically tested and, finally the more efficient one was adopted. These are described briefly in the following. The first scheme for the program flow is summarized in the following steps: 1. Create stiffness matrices from T<7;J, txi; 2. Assemble the matrices and solve for tV{ and *vf; 3. Calculate tV{At, tv?At and use To-ij, tViAt to calculate Aa^; 4. Use 1(Tij and Aa^ as input to stress integration and return mapping algorithm to get the new stress To-ij after the current iteration; 5. Check convergence. If step is not converged, apply the residuals divided by A i as new incremental load and go to step 1; otherwise, leave the iteration loop. The second scheme for the program is summarized in the following steps: 1. Create stiffness matrices from T;Ai, up to the current iteration, i.e., from the beginning of the current incremental step, and the stress at beginning of the current increment V^- are applied to calculate the stress increment Acr;j, so that A. The material's Young's modulus E is taken as 200 GPa. In this example, small deformation but large rotation and displacement occur. The analytical stress solution may be obtained from Hooke's law and Mohr's circle, and it is used as a reference to test the general structure and the correctness of the program. at time t=0 at time t=t Figure 5.1: Element with combined extension and rigid body rotation In the numerical solution, the three different algorithms discussed in Chapter 4 are applied to calculate the 4 Z);J . The stress is integrated using Equation (3.21). Stress results of the above example are shown in Figure 5.2. It is clear that the predictions of the forward difference formula with modification term and the central difference algorithm agree very well Chapter 5. NUMERICAL EXAMPLES 85 with the analytical values, even for large rotation of up to 90 degrees. The conventional forward difference algorithm predicts, however, incorrect values for stress ax. Because the non-zero term in Equation (3.27) is a constant, it brings a constant increment for the normal stresses and in this particular example, the stress increment owing to the rotation is much less than this value, so that the stress ax versus the rotation angle u>At is a straight line. It is found that reducing the incremental step in conventional forward difference algorithm has very little or no effect on the improvement of accuracy of crx. This may be due to the accumulation of the value from non-zero term in Equation (3.27). The above example shows that applying conventional time discretization or forward difference method in calculation of lDij or stress integration may cause significant errors, especially when the boundary conditions are only displacement boundary conditions. The errors can be eliminated by using the central difference method in implicit stress integration. In explicit integration algorithm, however, or when lDij is employed separately, e.g., in visco-plastic deformation, the modification term developed in Chapter 3 may be applied to keep the objectivity of lDij. 5.3 E L A S T I C - P L A S T I C S M A L L D E F L E C T I O N ANALYSIS O F A S Q U A R E E L E M E N T The example is a simple verification problem used by many commercial programs [10]. In this problem, a bar of square cross section is incrementally loaded by tension on two opposite edges, as shown in Figure 5.3. Material properties are taken as Young's modulus E \u00E2\u0080\u0094 18.4 MPa (8/3 ksi), Poisson's ratio v = 0.5, initial yield stress ayieid0 = 55.1 KPa (8.0 psi) and hardening modulus H = 2.6 MPa (8/21 ksi) and the material is assumed to have a bilinear isotropic hardening behavior. The final pressure is increased to 137.8 KPa (20 psi) in 20 steps. Because the geometry and loading are symmetric, only a quarter is modelled with a 4-node er 5. NUMERICAL EXAMPLES o x (Pa) 1 . O O O O E 8 8 . 0 0 0 0 E 7 1 ^ \ --O \u00E2\u0080\u0094 ANALYTICAL SOLUTION \u00E2\u0080\u0094 MODIFIED F O R W A R D \u00E2\u0080\u0094 C E N T R A L D I F F E R E N C E \u00E2\u0080\u0094 FORWARD D IFFERENCE 6 . 0 0 0 0 E 7 \ \ \ \ 4 . 0 0 0 0 E 7 \ - \ \ 2 . 0 0 0 0 E 7 \ \ O . O O O O E O - 1 _\ - 5 . 1 2 0 0 E 2 \ - 1 . O O O O E 9 \u00E2\u0080\u00A2 2 . O O 0 0 E 9 - 3 . 0 0 0 0 E 9 C 1 o 2 0 3 0 -40 S O 6 0 a. Stress distribution o x 7 0 s o g o co At (degree) Xxy(Pa) ANALYTICAL SOLUTION 5.0000E7 - Q MODIFIED F O R W A R D O C E N T R A L D IFFERENCE A FORWARD D IFFERENCE 4.0000E7 -3.0000E7 -2.0000E7 -9.9999E6 --1.4753E2 . . I . . . . i .... i .... i . . . . i . . . . i . , . . I . . . . I . . . 0 10 20 30 40 50 60 b. Stress distribution x Xy 70 80 90 co At (degree) Figure 5.2: The stress results of single element under rotation Chapter 5. NUMERICAL EXAMPLES 87 quadrilateral element. To compare the results, the exact same problem is also simulated by NISA and ANSYS. Figure 5.3: Elastic-plastic small deflection model Figure 5.4 is the curve of tensile load versus displacement in x direction ux. Before a load level of approximate 68.9 KPa (10 psi), results from the developed program, NISA and ANSYS are almost identical. Discrepancy of results starts at a load level of 68.9 KPa (10 psi) and increases with the increasement of load or deformation. This may be attributed to different level of approximation used in each program for handling the effect of finite deformation. Generally, however, the solution from A L E F E is in good agreement with the ones from NISA and ANSYS programs. Chapter 5. NUMERICAL EXAMPLES Figure 5.4: Load versus displacement in x direction Chapter 5. NUMERICAL EXAMPLES 89 5.4 P U R E B E N D I N G O F A R E C T A N G U L A R P L A T E 5.4.1 Problem Description This example is designed to investigate the effect of distorted mesh on the accuracy of simulated results and to show the necessity of mesh motion in A L E . The example is also analysed by NIS A. In this example, a rectangular plane strain plate is subjected to pure bending, as shown in Figure 5.5. The deformation is pure elastic to eliminate effects caused by the differences in handling plasticity related aspects. Young's modulus E is 200 GPa, Poisson's ratio v = 0.3 and concentrated force F is 2000 N/mm. The problem is solved using 5 incremental steps. y F E E o CM 20 (mm) Figure 5.5: Pure bending of a rectangular plate Chapter 5. NUMERICAL EXAMPLES 90 5.4.2 Finite Element Models Finite element models used in the analysis are shown in Figure 5.6. Model-a consists of 100 square quadrilateral elements with 10 elements on each edge. This homogeneous and dense mesh is expected to give more accurate results that will be used as reference [67]. Model-b has only 4 square elements. Although still homogeneous, it is much coarser and is applied as \"standard\" model to compare with model-c and model-d which have distorted meshes. In model-c, although all 4 elements have interior angles less than 180\u00C2\u00B0, so negative Jacobian determinant will not appear, the mesh is severely distorted. Model-d is an extreme case where element 3 has an interior angle larger than 180\u00C2\u00B0. 5.4.3 Results of Lagrangian Formulation The deflections are simulated first by the updated Lagrangian option of the developed program A L E F E and the commercial code NISA. Models a, b and c are run successfully for a maximum value of F equal to 2000 N/mm. The deflection at the lower-right corner of the plate is shown in Figure 5.7. It may be seen that the fine homogeneous mesh, model-a, gives the largest displacements, the distorted mesh, model-c, gives the least displacements; and the coarse homogeneous mesh, model-b, predicts an intermediate displacement values. In the example, the distorted less dense mesh stiffens the solution. This may be generalized to elastic and elastic-plastic problems when displacement methods are applied in finite element analyses, because solution starts from a kinematically admissible condition and a upper bound solution is expected. Model-d did not go through both in A L E F E and NISA because of the negative Jacobian determinant in element 3 with both programs reporting a value of \u00E2\u0080\u00942.6036. The negative Jacobian is due to an interior angle larger than 180\u00C2\u00B0. The only limitation on the shape of the real element is that each interior corner angle must be less than 180\u00C2\u00B0, i.e., the element must be Chapter 5. NUMERICAL EXAMPLES 3 4 1 2 Model -a Model -b 3 4 1 4 1 \ 2 Model-c Model -d Figure 5.6: Finite element models for the pure bending of a plate Chapter 5. NUMERICAL EXAMPLES 92 Figure 5.7: Pure bending of a plate: Lagrangian Formulation Results Chapter 5. NUMERICAL EXAMPLES 93 convex [68]. Practically, however, to achieve higher accuracy, the optimum elements should have interior corner angles of quadrilaterals not far from 90\u00C2\u00B0 [67] and aspect ratios from 1/5 to 5 [68]. 5.4.4 A L E Simulation and Result The same problem with model-c is also simulated by the A L E option of the program A L E F E . The results are compared with the Lagrangian ones from model-a, model-b and model-c in Figure 5.8. In the first incremental step, A L E and updated Lagrangian results coincide. This is obvious since A L E always starts with a updated Lagrangian increment for the problem without contact boundaries. In subsequent steps, the mesh begins to move, in A L E method, according to the scheme specified in Chapter 3 and mesh quality starts to improve, so that the results are closer to those from homogeneous mesh, model-b. After three steps, transition stage due to significant mesh motion seem to end and the response from A L E and from updated Lagrangian method are parallel to each other. The discrepancy is owing to the incremental method which causes some history-dependent effects in the simulation. This example shows that A L E F E and NISA predict identical results when same methods are applied. The mesh distortion has significant effect on the accuracy of the results. Mesh distortion may cause negative Jacobian leading to termination of the simulation, as indicated by model-d. A L E method may eliminate mesh distortion and improve accuracy of simulation. 5.5 P U N C H F O R G I N G P R O C E S S 5.5.1 M o d e l and Updated Lagrangian Results A punch forging process discussed in Chapter 1 is employed to demonstrate the capabilities of the A L E method and to compare with conventional Lagrangian methods. To show possible effects of different meshes on A L E and conventional updated Lagrangian method, homogeneous Chapter 5. NUMERICAL EXAMPLES 94 a. Displacement in x direction 2500 r -0.000 0.010 0.020 0.030 0.040 0.050 0.060 0.070 0.080 0.090 luyl(mm) b. Displacement in y direction Figure 5.8: Pure bending of a plate: A L E Results Chapter 5. NUMERICAL EXAMPLES 95 and inhomogeneous meshes are used in simulation. The material properties, geometry and boundary condition have been described in Chapter 1. The problem is simulated by the developed program A L E F E and commercial codes NISA [10], ANSYS [9] and DEFORM2 [8]. Analysis is performed for up to 60% reduction in the height of the original piece. The homogeneous mesh is shown in Figure 1.4 and is composed of square elements in the undeformed geometry. The results from NISA and ANSYS are obtained using updated Lagrangian method and are shown in Figure 1.5. For proper modeling of this problem, the material point under the punch corner should be free once it moves out of the corner. However, the Lagrangian formulation can not update the boundary condition automatically. This is an important feature in handling most of metal forming problems. This deficiency has the effect of increasing the punch size during the deformation, as shown in Figure 1.5. It is also noted that elements at the punch corner are highly distorted because of the large deformation. Increasing the deformation level causes more severe element distortions and the program eventually terminates the analysis. The deformed shape from DEFORM2 is shown in Figure 1.7. In order to achieve conver-gence for the same deformation level obtained in NISA and ANSYS, a much denser mesh is required, as shown in Figure 1.6. Although remeshing eliminates the punch size increment, it may be seen that the mesh around the punch is highly distorted and more importantly, the predicted deformed boundary of the free surface near the punch corner is unusual and not rational. Another shortcoming is the fluctuations in loads predicted by Lagrangian method, which could not be eliminated by auto remeshing. Figure 1.8 shows a plot for the reduction of the specimen height versus the applied load, from DEFORM2. The result shows noticeable load fluctuation that is pertinent to the updated Lagrangian formulation. Chapter 5. NUMERICAL EXAMPLES 96 5.5.2 A L E Simulation In general, contact boundary condition may be easily handled in A L E by utilizing the feature of arbitrary mesh motion. On the interface of tool and work-piece, two specific directions are considered; perpendicular and parallel to the tool interface. In the direction perpendicular to the tool interface, the nodes will be kept as Lagrangian points, i.e., nodes will be moving with the tool at tool speed. In the other direction, however, nodes will be kept as Eulerian points, i.e., nodes will be fixed parallel to the tool interface. Contact boundary condition may be accurately described, therefore, in any general form of tool workpiece interface. The updating of boundary condition is no longer necessary. The \"punch increasement\", and load fluctuations discussed above, should generally not occur. This will undoubtedly increase the efficiency and accuracy of the simulation of metal forming problems. In the particular case of punch forging process, as shown in Figure 1.4, the nodes on ED are kept as Lagrangian points in Y-direction, i.e., have the same velocity as the punch. The same nodes are treated, however, as Eulerian points in X-direction, i.e., have a zero velocity. This ensures a node remains coincident with each punch corner in the whole deformation process. The nodes on EO are Lagrangian in X-direction and general A L E ones in Y-direction where the speed linearly reduces from v at point E to zero at point O. The nodes on OA are Lagrangian in Y-direction and Eulerian in X-direction. On AB, nodes are treated as Lagrangian in Y-direction, and general A L E in X-direction. Boundaries DC and CB are traction free, so all the nodes are Lagrangian in both directions. The mesh motion scheme described in Chapter 3 is employed to control the internal mesh motion. The simulation is successfully carried out for the whole course of deformation. The results of deformed shape and current yield stress from the A L E method are pre-sented in Figure 5.9. It is obvious that the contact boundary condition for the nodes under the punch is handled automatically and accurately by allowing the degree of freedom in the Chapter 5. NUMERICAL EXAMPLES Figure 5 .9: The current yield stress from A L E Chapter 5. NUMERICAL EXAMPLES 98 vertical direction to be Lagrangian while freezing the degree of freedom in the horizontal di-rection and the punch size increment is eliminated. The deformed shape near the punch corner is acceptable with no element distortion. Both displacements and stresses obtained from the A L E are in agreement with the results reported in [69]. Figure 5.10 shows the load versus reduction curve from the A L E simulation. The load fluctuation caused by updated Lagrangian formulation (Figure 1.8) is completely eliminated in the A L E results. That is believed to be a result of accurately describing the contact boundary condition in the whole deformation pro-cess. Figure 5.11 shows contours of normal stress in punch direction, cry. At the punch corner, the compression stress is maximum. In this area, stress values change dramatically. For the area under the punch, the compression stress reduces from the surface to the center. On the traction-free surface around corner, a small tensile stress area is noticed and values go down from surface to the center in most of the area. 4.0000E7 3.0000E7 If T3 CO O _ l 2.0000E7 1.0000E7 0.1 0.2 0.3 0.4 0.5 0.6 Reduct ion in height Figure 5.10: Load-height reduction result from A L E Chapter 5. NUMERICAL EXAMPLES 99 o y (Pa) 2.6857E9 1.8397E9 9.9365E8 1.4760E8 -6.9844E8 -1.5445E9 -2.3905E9 -3.2366E9 -4.0826E9 -4.9287E9 -5.7747E9 -6.6208E9 -7.4668E9 -8.3129E9 -9.1589E9 0.0 10.0 20.0 30.0 40.0 (mm) Figure 5.11: Contours of cry from A L E Boundary node adjustment in the mesh motion scheme discussed in Chapter 3 is nec-essary. Figure 5.12 shows the deformed shape and current yield stress distribution from A L E simulation without the adjustment. It is clear that although the punch size is not increased in the simulation and the mesh is generally homogeneous, the aspect ratio of the first column of elements outside the punch coiner may not be acceptable. More importantly, most of nodes on the traction-free surface moved to the far right side and no node is located close enough to the punch corner to represent the boundary shape, where the shape changes dramatically. Chapter 5. NUMERICAL EXAMPLES Oyield (Pa) - 1.2799E9 - 1.2114E9 -- 1.1430E9 1.0745E9 1.0060E9 \u00E2\u0080\u0094 9.3761 E8 \u00E2\u0080\u00A2 8.6914E8 8.0068E8 7.3222E8 - 6.6375E8 - 5.9529E8 - 5.2682E8 m 4.5836E8 m 3.8989E8 -3.2143E8 5.0 \u00E2\u0080\u00A2 r , . , i i i i . . . . i . .0 10.0 20.0 30.0 40.0 50.0 (mm) Figure 5.12: The current yield stress without boundary node adjustment Chapter 5. NUMERICAL EXAMPLES 101 5.5.3 Simulations from an Inhomogeneous Mesh An inhomogeneous mesh with same connectivity as the homogeneous mesh is applied as shown in Figure 5.13. Although the mesh is inhomogeneous, elements are still optimally shaped. The simulation is carried out using NISA, ANSYS and A L E F E with A L E option. The displacements of the top right corner, Point C, from all simulations are shown in Table 5.1. Y i i \ 12 mm \ E W ' \ \ \ \ \ D C { mm o i \u00E2\u0080\u0094 1 0 t \ 30 mm B X Figure 5.13: Inhomogeneous mesh for punch forging simulation By comparing the stress and displacement distributions from different programs, it is found that although NISA gives similar displacements from homogeneous and inhomogeneous meshes, stress distributions are much different from the two meshes with stresses from the inhomogeneous mesh 2 to 10 times those obtained from homogeneous mesh. As an example, the current yield stresses are illustrated in Figure 5.14, in which the difference is about 5 times. Chapter 5. NUMERICAL EXAMPLES 102 Table 5.1: Displacements of top right corner from all simulations Program ux mm) Uy mm) Homogeneous Inhomogeneous Homogeneous Inhomogeneous NISA 15.459 14.583 0.003 0.012 ANSYS 18.741 18.652 0.096 0.021 A L E F E 13.929 12.019 0.010 0.040 ANSYS results show closer stress and displacement values and the predicted deformation shape is not improved. Shape and the mesh at the punch corners are severely distorted. The developed A L E F E program gives not only similar displacements, as indicated in Table 5.1, and quite close stress values, but also a realistic deformation shape. The largest difference of stress in A L E exists in rxy where the value from inhomogeneous mesh is 1.6 times that from homogeneous mesh. As an example, the current yield stress contour is shown in Figure 5.15. When compared with the result from homogeneous mesh, shown in Figure 5.9, the difference is less than 15%. The nodes are closer to the punch corner on the traction-free surface, so that the deformed shape around the corner is represented better. The A L E results of the punch forging problem indicates that the \"punch size increase-ment\", load fluctuation and mesh distortion caused by traditional updated Lagrangian method may be eliminated by the A L E scheme. The boundary node adjustment in mesh motion is necessary to keep high quality mesh and to achieve high accuracy. Inhomogeneous meshes and changes in mesh density do not affect the result significantly. Accuracy may be improved if mesh has higher density at locations where deformation is highly inhomogeneous, as at the punch corners. In the traditional updated Lagrangian method, however, such density changes may decrease the calculation accuracy and lead to significantly different results. The discrep-ancy is believed to be caused by the severe mesh distortion which occurrs in dense meshes. Chapter 5. NUMERICAL EXAMPLES 103 a yield (\"1.0E8 Pa) a. The result from homogeneous mesh Figure 5.14: Comparison of current yield stresses from NISA Chapter 5. NUMERICAL EXAMPLES Figure 5.15: The A L E result of current yield stress from inhomogeneous mesh Chapter 5. NUMERICAL EXAMPLES 105 5.6 T E N S I L E N E C K I N G ANALYSIS OF A N E L A S T I C - P L A S T I C B A R The initiation of tensile necking of an elastic-plastic bar in plane strain is simulated by updated Lagrangian method, which is degenerated from the A L E formulation. The objective is to check the correctness of the method and to compare the integration algorithms discussed in Chapter 3. The different integration algorithms used in this example are: integration algorithm of Truesdell rate discussed in Chapter 3, Equation (3.15), integration algorithm of Truesdell rate developed by Pinsky [58], Equation (3.20), and integration of Jaumann rate developed by Hughes [56], Equation (3.21), combined with transformation A and B discussed in the chapter. The bar has a length to width ratio of 3 and the central one-sixth of the length is thinned slightly so that the necking analysis may be formulated as a standard deformation problem instead of an eigenvalue bifurcation problem. Only one quarter of the bar is modelled due to symmetry and the finite element model is shown in Figure 5.16. There are 256 quadrilateral 4-node elements and 297 nodes in the model. The material has a Young's modulus to yield stress ratio of 7.5 x 102 , Poisson's ratio of 0.3 and the ratio of hardening modulus to yield stress is constant at 1.25 with a yield stress value of 250 MP a. The length L is taken as 30.0 mm and W0 = 10.0 mm in calculation. The amplitude S of the neck is defined by the width reduction at center section normalised by the initial length L [3]. The simulation is carried out until S is about 0.07. The same problem is also solved using the commercial finite element code ANSYS [9]. Results are compared with the results reported by McMeeking and Rice [3]. The necking parameter S is plotted against the engineering strain 7 = A.L/L. In Figure 5.17, Jaumann stress rate results are labelled as Jaumann-A and Jaumann-B. These correspond to the use of Jaumann rate, Equation (3.21), with transformation-A and transformation-B, respectively. Similarly, Truesdell-A and Truesdell-B are from the integration of Truesdell rate, Equation (3.15), with transformation-A and B, respectively. Chapter 5. NUMERICAL EXAMPLES Figure 5.16: The finite element model for tensile necking er5. NUMERICAL EXAMPLES Figure 5.17: The necking S and engineering straining 7 Chapter 5. NUMERICAL EXAMPLES 108 It can be seen from Figure 5.17 that for integration of Jaumann rate, the transformation algorithm A and B lead to same results. For Truesdell rate, transformation-A and B give quite different results. This indicates that Truesdell rate is more sensitive to the location of stress transformation. It should be mentioned here that the results from Truesdell A scheme have almost same slopes and very close values to those reported by McMeeking [3]. The results from ANSYS show very much different values of necking. The ANSYS run was terminated at 7 = 0.57 because of divergence, whereas all the runs of the developed A L E program have been successful for values up to 7 = 1.0. Figure 5.18 shows results from different schemes developed by Pinsky et al. [58] for the integration of Truesdell rate, Equation (3.20) in this thesis, (curves: Truesdell-A, P and Truesdell-B, P). Figure 5.18 also shows corresponding results developed in this research, Equation (3.15), ( curves: Truesdell-A and Truesdell-B). The results indicate that schemes developed in this research give same results as those developed by Pinsky et al. [58]. This is expected since the difference between Equation (3.15) and Equation (3.20) is only in the second term on right hand side and is caused only by higher order terms of A i . Figure 5.19 shows the contours of the equivalent von-Mises stress and the ratio of the equivalent stress to current yield stress at 7 = 0.75. Although no stress distribution is reported in [3] , the shape of unloading zone is discussed. From Figure 5.19, the zone where the ratio of the equivalent stress to current yield stress is less than 1.0 is experiencing unloading. The unloading zone shown here has the same shape and location as that reported in [3]. The example indicates that the transformation algorithm A and B may generally be combined with any integration algorithm. The integration of Jaumann rate is less sensitive to the location of stress transformation. The combination of Truesdell rate with Transformation-A generally gives better results than other cases. Algorithms for integration of Truesdell rate developed in this research gives same results as those developed by Pinsky [58]. Similar conclusions to those reported above have been obtained through simulation of the punch forging Chapter 5. NUMERICAL EXAMPLES 109 Chapter 5. NUMERICAL EXAMPLES Figure 5.19: Contours of equivalent stress (Truesdell rate, transformation Chapter 5. NUMERICAL EXAMPLES 111 problem by A L E method. 5.7 SHEET METAL EXTRUSION 5.7.1 Extrusion through a Curved Die Plane strain sheet metal extrusion is another typical metal forming problem in which large strains are expected and boundary condition updating as well as re-meshing are required in traditional updated Lagrangian methods. The A L E method is applied to simulate two extrusion cases to show the advantage of A L E method in such deformation problems where the deformation extent is independent of the displacement of the tool after certain steps, i.e., when piston moves a distance equal to the length of the die. The first simulation (curved die) has been taken as a bench mark problem by many researchers [7, 70] whereas the second simulation has some limited experimental results available in the literature. An aluminium billet in an initially stress free state is forced into the die while sliding between two smooth rigid plates. A rigid smooth piston pressing against the rear face of the billet and moving with prescribed velocity provides the driving force. The die produces 25% thickness reduction over a distance 1.2a, where a is the half thickness of the original sheet and taken as 25.4 mm (1 in) as shown in Figure 5.20. Material properties of the Aluminium billet are: Young's modulus E = 6.89 x 104 MPa (1 x 107 psi); hardening modulus H = 1.14 x IO3 MPa (1.65 x 105 psi), initial yield stress tr^wo = 3.95 x 102 MPa (5.7 x 104 psi), and Poisson's ratio v \u00E2\u0080\u0094 0.3. The example is the same as the one discussed in Chapter 1, except that the die is shaped in the form of 5th order polynomial curve with zero curvature and slope at the ends. The same process was simulated with updated Lagrangian method by Lee [12] and by Yamada [70], where the billet was assumed to have been initially machined to fit the die. Figure 5.20 shows the geometry and the finite element mesh used in A L E simulation. In the A L E simulation, the mesh and the corresponding nodes confined to the die area are taken Chapter 5. NUMERICAL EXAMPLES Figure 5.20: Original geometry and mesh in extrusion through a curved die Chapter 5. NUMERICAL EXAMPLES 113 as Eulerian during the whole course of deformation. The nodes on traction free boundaries are designed to be completely Lagrangian, whereas those on symmetric axes, contacting with piston, or sliding on plate are taken as Lagrangian in the normal direction and the general A L E in tangential direction. All other nodes are assumed to be A L E ones and are controlled by the mesh motion scheme detailed in Chapter 3. -eq 0.502172 0.468694 0.435216 0.401738 0.368259 0.334781 0.301303 0.267825 0.234347 0.200869 0.167391 0.133913 0.100434 0.0669563 0.0334781 x/a Figure 5.21: Contours of plastic strain after piston moved 2a Figure 5.21 shows the plastic strain distribution and the deformed shape after a piston movement of 2a. As shown from the figure, a piston displacement of 2a corresponds to approximately 2.0 times the width of the die. The figure shows the shape of plastic zone and how it is expanding. The distribution of longitudinal stress ax at different lateral positions, labelled \"lateral levels\", is also shown in Figure 5.22. A lateral level corresponds to the initial distance from the central line of the sheet. Similar results from reference [12] are shown in Chapter 5. NUMERICAL EXAMPLES 114 Figure 5.22 as well. The values of shear stress rxy obtained from the two methods are illustrated in Figure 5.23. Generally, good agreement is observed. The results of reference [12] were obtained from a sequentially updated Lagrangian scheme with special introduction of initial stress stiffness terms as well as dilatation and rotation terms in the stress updating scheme. Special attention was given when applying the boundary condition at the die face and after the material leaves the die face. Numerical and computational difficulties are reported in [12] in the process of correctly applying and updating the boundary conditions. The treatment described is essentially unique for the given application. 5.7.2 Extrusion through a Straight Die Plane strain metal extrusion is investigated experimentally and analyzed with slip line field method by Farmer and Oxley [71]. Measuring the die separation force is not practical and would greatly complicate the design of the apparatus. Therefore, only the extrusion force is measured and the friction on the die interface is minimized by effective lubrication. The half die-angle is 30\u00C2\u00B0 and the thickness reduction ratio is 50%. The work-piece has cross-section 2.54 mm (0.1 in) thick x 10.16 mm (0.4 in) wide, the thickness being reduced to 1.27 mm (0.05 in) during extrusion. In order to reduce the peak load and to reduce the time taken to reach steady-state condition, the work-piece is pre-formed to the angle of die. The material is 5052-H34 tempered aluminium. A stress /strain curve for this material obtained from compression test is applied to derive the initial yield stress of \u00E2\u0080\u0094 o 100 \u00E2\u0080\u0094-f\u00E2\u0080\u0094 , i , 200 t --1-300 Step No. . . f . - a \u00E2\u0080\u00A2 CASE1.A - - O - - CASE1 ,B \u00E2\u0080\u00A2 CASE2 A CASE3 V CASE4 O CASE5 ' Parameters lor different cases are shown in Table:5.2 ' ' \u00E2\u0080\u00A2 ' \u00E2\u0080\u00A2 1 1 ' 400 500 600 Figure 5.31: Locations of neutral points when velocity criterion is applied The correct and final location of the neutral point should keep the external forces self-balanced. That is a unique feature of the rolling process and generally makes the simulation of Chapter 5. NUMERICAL EXAMPLES 129 (MPa) 400 0.00 (MPa) 30 0.02 0.04 0.06 0.08 Arc of contact p a. The normal stress 0.10 0.12 (Rad.) 0.00 0.02 0.04 0.06 0.08 Arc of contact (3 b. The friction stress 0.10 0.12 (Rad.) Figure 5.32: Rolling pressure and friction stress when velocity criterion is applied Chapter 5. NUMERICAL EXAMPLES 130 rolling more complicated. The above methods of deciding neutral points introduced resultant forces on either the free edge HA or at the neutral point itself. These forces are not present in the real system and so the methods are unacceptable. A more rational definition of the neutral point is thought to be given by the following force balance equation in the rolling direction [77]: / (TRCOSB - pRsinB) d8 - {TRCOS8 + pRsinB) d8 = 0 (5.5) J/3N JO where, 8 is the angle from the vertical line to the line connecting the contact point and the roll center as shown in Figure 5.29, 8N is the value of 8 at the neutral point N, and 8G is the value of 8 at point G. In this proposed simulation of the neutral point, the velocity boundary condition may be either applied at the neutral point itself or on the external boundary HA as discussed above. Figure 5.33 shows how the location of the neutral point moves from the initial guess to the final steady state location when this simulation is used. It is observed that the neutral points moved from the initial guesses to steady state locations for all reduction cases, and the difference in the initial guess has no impact on the final result. To show this, two different initial guesses are considered for Case 1 and the results are shown in Figure 5.33 as Case la and Case lb. It is observed that the final neutral point has the same position after the 2nd step. Applying velocity boundary condition on the stress-free boundary H A gives the same results for neutral point in all 5 cases. Figure 5.34 shows the impact on rolling pressure when applying the velocity boundary conditions on the stress-free boundary HA (curve ALE-BOU) or on the neutral point itself (curve ALE-NEU). It is obvious that the two methods give almost identical results and the \"jumps\" caused by the previous methods are eliminated. It is also important to investigate the reaction forces in this new simulation. The reaction force at HA, when the velocity boundary condition is applied at the same boundary, is 2.31 N/mm, which is quite small and may be ignored if compared to the rolling force of 875.11 N/mm. The above results suggest that this Chapter 5. NUMERICAL EXAMPLES method is more reliable and more accurate than the ones previously discussed. 131 \u00E2\u0080\u0094l 10 100 200 Step No. 300 400 500 Figure 5.33: Locations of neutral points when force condition is applied 5.8.5 Simulation of Rolling Pressures As far as the rolling pressure distribution patterns is concerned, the classical slab method gives a distribution with maximum pressure, known as \"friction hill\", at the neutral point, on both sides of which the pressure decreases monotonically, as shown in Figure 5.36. However, some circumstances were observed experimentally [76] in which the pressure distribution possess double peaks. In the present computed results, double-peak pressure distribution curves, as well as \"friction hill\" distribution are obtained depending on the parameters involved. Figure 5.35 and Figure 5.36 show the pressure variations for two typical cases, Case 1 and Case 5, respectively. For small roll ratios, R/h0, two-peak roll pressures are obtained with maximums near entrance and exit and with pressure drop in the middle of contact arc. For cases of large R/h0, Figure 5.36 shows a friction-hill distribution. It may be observed that the computed results are, generally, Chapter 5. NUMERICAL EXAMPLES 350 0.04 0.06 0.08 Arc of contact p (Rad.) a. The normal pressure 0.12 20 15 _ 10 (0 0_ 2 to 0 c o o \Y -5 -10 -15 -20 i 7 ! 0.00 0.02 0.04 0.06 0.08 Arc of contact P (Rad.) b. The friction stress 0.10 0.12 Figure 5 .34: Rolling pressure and friction stress when force criterion is applied Chapter 5. NUMERICAL EXAMPLES 133 in good qualitative agreement with the experimental results. The discrepancies at the entrance and exit points may be attributed to inaccurate description of friction, the roll flattening and the spring back of work-piece after exit. The experimental results show significant increase in contact angle at both entrance and exit, while the contact angle in the numerical simulation is smaller, obtained only from pure geometric relations and has the same value as in a classical slab method. 350 i~ -0.02 0.00 0.02 0.04 0.06 0.08 Arc of contact p (Rad.) 0.10 0.12 Figure 5.35: Comparison of rolling pressure, small R/ho 5.8.6 Simulation of Rolling Force and Torque The roll separating force and the roll torque under different deformation reduction but same R/hQ values form Case 1 to Case 4, are simulated by applying the velocity boundary condition at neutral point and entrance. Figure 5.37 and Figure 5.38 show a comparison of the numerically obtained results and experimental ones from reference [76]. Once again, the Chapter 5. NUMERICAL EXAMPLES 134 Figure 5.36: Comparison of rolling pressure, large R/h0 two boundary conditions (ALE-BOU and ALE-NEU) gave almost identical results as in the roll pressure distribution. The simulated results for rolling force are smaller than the measured values for some reductions and larger for others whereas the simulated torque values generally overestimate the experimental values. The difference between the computed and measured values may be attributed to the assumption of rigid roll in the finite element analysis, to the uncertainty in friction modeling and measurements and to the strain hardening modeling and measurements. Figure 5.39 is the distribution of equivalent plastic strain for Case 1. It shows the spread of plastic zones within the strip. Prior to entry, a large plastic zone is observed and yielding is indicated to take place before the point of contact. At the center, yielding occurs very near to the point of contact. Continued elasto-plastic flow is observed after yielding. Figure 5.40 shows contours of current yield stress for Case 5 at different steps. It may be observed that the maximum values are moving horizontally from entrance to the exit. The figure also shows how Chapter 5. NUMERICAL EXAMPLES 135 2500 r-2000 , 1000 500 I ...^^Or. , O EXPERIMENT \u00E2\u0080\u00A2 - \u00E2\u0080\u00A2 - - ALE-BOU \u00E2\u0080\u00A2\u00E2\u0080\u0094* ALE-NEU I _L 0.15 0.20 0.25 0.30 Reduction in height 0.35 0.40 Figure 5.37: Comparison of rolling force 1.4000E4 1.2000E4 h 1.0000E4 8.0000E3 0 o r -6.0000E3 4.0000E3 2.0000E3 0.0000E0 0.10 0.15 0.20 0.25 Reduction in height 0.30 0.35 0.40 Figure 5.38: Comparison of rolling torque Chapter 5. NUMERICAL EXAMPLES 136 the plastic zone is developed. In the simulation, a steady state condition was reached after 300 incremental steps. x (mm) Figure 5.39: Equivalent plastic strain (Case 1: R/ho = 12.5, Reduction = 14.17%) 5.8.7 Concluding Remarks One of the most complicated deformation processes in metal forming: the rolling process, is successfully simulated by the Eulerian method reduced from the general A L E formulation developed in this research. Program A L E F E is directly applied to simulate the process without special developments or enhancements. Various methods for applying correct boundary con-ditions and locating neutral points are investigated and proved to be helpful in simulating such kind of problems efficiently and accurately. The method used in handling friction forces in this Chapter 5. NUMERICAL EXAMPLES 137 a. Step 120 b. Step 210 c. Step 300 'yield 82.0175 79.903 77.7885 75.674 73.5595 71.445 69.3305 67.216 65.1015 62.987 60.8725 58.758 56.6435 54.529 52.4145 x (mm) \u00C2\u00B0yield (MPa) 90 4475 87.771 B5.0945 82.418 79.7415 77.065 74.3885 71.712 69.0355 66.359 63.6825 61.006 58.3295 55.653 52.9765 yield 10 x (mm) (MPa) 95.0366 92.0541 89.0717 86.0892 83.1068 80.1244 77.1419 74.1595 71.1771 68.1946 65.2122 62.2297 SO 2473 56.2649 53.2824 x (mm) Figure 5.40: Current yield stress at different steps (Case 5: R/h0 Reduction = 22.83%) 39.0, Chapter 5. NUMERICAL EXAMPLES 138 work is simple and proved to be effective in conjunction with the Eulerian simulation. The computed results showed both double-peak as well as \"friction hill\" distributions for rolling pressure. In general, the results are in good agreement with experimental ones. 5 . 9 C O M P R E S S I O N B E T W E E N W E D G E - S H A P E D D I E S 5 .9.1 Deformation Conditions Plane strain compression of pre-shaped material between wedge-shaped dies is simulated by A L E method. The same process is investigated experimentally and analysed with slip line method by Chitkara and Johnson [78]. The dies are considered to be perfectly rough, so sticking friction conditions are applied in the simulation. Another distinct characteristic of the process is that the dies are wedge shaped which brings out two consequences. A more general coupled displacement / velocity boundary condition on the die-metal interfaces is needed because the dies themselves move in the process and the deformation is no longer symmetric about the vertical axis. The general shape and size of specimen and die are indicated in Figure 5.41. In order to get the curves of mean vertical die pressure and end velocities with respect to strip thickness, different specimen dimensions are employed as shown in Table 5.3. The unique tool angle r/ = 3\u00C2\u00B0 is used because most of the comprehensive experimental results are available in literature for this angle [78]. The specimen ABCDEFGP is compressed between two inclined dies BC and F G with angle of inclination n with the horizontal direction. The specimen is pre-shaped such that surface BC and F G are to be matched with the tool face. The dotted lines indicate the approximate shape after deformation. Velocities vt and vs are at large and small end of the specimen, respectively, while each tool is advancing with a velocity v. Experiments were conducted on 300-ton Denson hydraulic testing machine and loads of less than 3.45 MPa (500 lb/in2) were measured independently by means of two calibrated Chapter 5. NUMERICAL EXAMPLES y Initial shape Deformed shape Figure 5.41: General shape of the specimen and die Table 5.3: Specimen dimensions in compression between wedge-shaped dies No m Ho degree mm in mm in 1 3\u00C2\u00B0 25.40 1.000 22.74 0.895 2 3\u00C2\u00B0 22.23 0.875 19.56 0.770 3 3\u00C2\u00B0 19.05 0.750 16.39 0.645 4 3\u00C2\u00B0 15.88 0.625 13.21 0.520 5 3\u00C2\u00B0 12.70 0.500 10.04 0.395 6 3\u00C2\u00B0 9.53 0.375 6.86 0.270 7 3\u00C2\u00B0 6.35 0.250 3.69 0.145 8 3\u00C2\u00B0 4.45 0.175 1.78 0.070 Chapter 5. NUMERICAL EXAMPLES 140 pressure gauges, so that the die load may be recorded. Dial gauges were used for measuring downward displacement of the die and lateral movements of the specimen. Therefore, the deformed specimen thickness and velocities at two ends can be measured. The rough dies were prepared by machining a square grid on the surface and using no lubrication in compression [78]. The material of specimen chosen in A L E simulation is half hard aluminium. The stress-strain curve is obtained by experiment [78]. Using bilinear approximation of the constitutive curve, we use initial yield stress eryieid0 = 276 MPa and hardening modulus H = 25 MPa. Other material properties used in the simulation are: Young's modulus E = 69 GPa and Poisson's ratio v = 0.33 [72]. 5.9.2 Finite Element Model Due to symmetry about the horizontal axis only upper half of the specimen is discretized. The finite element model for Case 1 is shown in Figure 5.42. A total of 240,4-node quadrilateral elements are used in the model. The material nodes on horizontal axis are restricted to move vertically, however they are free to move horizontally. The material points on die-metal interface F G are restrained to move horizontally to achieve the sticking friction condition, but they are allowed to move in vertical direction along with the die movement. In A L E simulation, the mesh nodes on edge F G are treated as Lagrangian points in both horizontal and vertical direction. The mesh nodes on GB and FC are taken as Eulerian points, i.e., fixed in horizontal direction, and general A L E points in vertical direction where the velocity changes linearly from die velocity v at F and G to zero at C and B. The mesh points on AB and CD are kept as Lagrangian points in vertical direction but as general A L E point in horizontal direction where velocity increases linearly from zero at B and C to the material velocities of points A and D at A and D, respectively. The mesh points on BC are Lagrangian points in vertical direction and Eulerian points in horizontal direction. All other boundary nodes Chapter 5. NUMERICAL EXAMPLES 141 y (mm) M M I T T T T T T T T T T T T I ^ E 11. I .1 . Ii 1 .1 .1 .1 i l . 1 . 1 . 1 .1 I I. I.I .1 ii I.I .1 I I. Ill 1 I l. l. l il J ,1 ,1 ,1 i l , I, 1 ,1 l I, ,1 .1. ll 71 -30 -20 -10 0 10 20 30 40 50 60 x (mm) B e D Figure 5.42: The finite element model for Case 1 are on traction-free boundaries, and, therefore, they are Lagrangian points in both horizontal and vertical directions. The motion of internal nodes is specified by the mesh motion scheme discussion in Chapter 3. The compression is performed for up to 50% reduction in initial height in 500 steps. 5.9.3 Slip Line Field Solution The slip line field solution [78] assumes that centred fans radiate from corners of the dies and that the starting slip lines are at 45\u00C2\u00B0 to the plane of symmetry of the specimen. The usual technique of approximating small arcs by their chords is applied to build up the slip-line field. The total normal and net tangential forces on the dies are evaluated. Then, vertical die pressures are obtained by resolving these forces and summing algebraically [78]. It should also be noted that the average yield stress is used in slip line analysis. The velocities of the specimen at the large and small end are achieved from the hodograph constructed according to the slip line 15 10 5 n Chapter 5. NUMERICAL EXAMPLES 142 field. 5.9.4 Comparison of Results The results predicted by A L E method are compared with the experimental and slip line field analysis results published by Chitkara and Johnson [78]. Figure 5.43 shows the result of Case 1 giving the variation of load on the die with respect to compressive strain In Hc/hc. The load-strain plot of A L E result matches closely the experimental results when strain is less than 0.30. An obvious discrepancy starts to appear after this level and increases as the strain increases. The discrepancies at large reduction levels may be due to a reduction in friction at the interfaces as the compression proceeds in experiment, which was reported in [78]. Besides, the experiments used a finite width 25.4 mm (1 in) of the specimen, compared with 25.4 mm (1 in) of projected length of die. Therefore, the plane strain conditions might not be achieved experimentally. Non-dimensionalized mean vertical pressures pv/2k as a function of the current strip thickness of the specimen at large side hi is presented in Figure 5.44. The result from A L E matches the experimental values much better than the slip line results. Figure 5.45 shows results of end velocities of specimen at different thickness hi. The values from A L E and slip line theory agree with the experimental results very well. 5.9.5 Development of Deformed Shapes and Plastic Zones Figure 5.46 shows the deformed shapes and contours of equivalent plastic strains at steps when reduction in height is 15%, 30% and 45% for Case 1. It is observed that even when reduction is 45%, the mesh is still homogeneous everywhere in the deformation domain. The punch or die size increasement which occurs in updated Lagrangian method does not appear here. From the equivalent plastic strain distributions at different steps, it is may be seen that Chapter 5. NUMERICAL EXAMPLES Figure 5.43: The variation of vertical load on the die with compressive strain Chapter 5. NUMERICAL EXAMPLES 144 Figure 5.44: The variation of mean vertical die pressure with current strip thickness Chapter 5. NUMERICAL EXAMPLES 145 Chapter 5. NUMERICAL EXAMPLES 146 the plastic zones are initially formed at the die corners and the strain values are highest through the whole deformation. The central parts of the interfaces have lower values of plastic strain, which means that the deformation is small and is incurred by the non-slip friction condition. This may be the reason why in slip line fields, a rigid body block is usually assumed there. It is interesting to indicate that plastic zones are already formed outside of the dies, but in slip line field constructed by Chitkara and Johnson [78], the plastic deformation begins only inside the dies, which may contribute to the inaccuracy of slip line results. Plastic zones are mainly located between and around the two dies. For most area outside the dies, plastic deformation is almost zero. er 5. NUMERICAL EXAMPLES L P eq 0.376818 0.351698 0.326576 0.301465 0276334 0251212 0.226091 0.20097 0.175849 0160727 0125606 0.100485 0.0753637 0.0502425 0.0251212 -30 -20 -10 0 10 20 30 40 50 60 a. The plastic zones at 15 % reduction in height b. The plastic zones at 30 % reduction in height p ^eq B 1.48481 1.38582 1.28684 1.18785 c. The plastic zones at 45 % reduction in height Figure 5.46: Equivalent plastic strain at different steps Chapter 6 C O N C L U S I O N S A N D F U T U R E W O R K 6.1 C O N C L U S I O N S It is shown that when the conventional Lagrangian methods are applied to simulate finite strain deformation processes, many problems arise: extensive mesh distortion or entanglement, load fluctuations and incorrect description of contact boundary conditions with sharp edges or corners. Available remeshing and rezoning techniques may overcome some mesh distortion and boundary description problems but they are generally not robust and they create additional problems in convergence and consume extensive CPU time. In this research, an A L E formu-lation is shown to have the potential of eliminating the above problems and of providing a promising tool by combining the merits of both Lagrangian and Eulerian schemes. A critical discussion of A L E formulation in solid mechanics is given in this work. It is shown that many authors have concentrated on particular problems that led to a special form or a restricted formulation. A general A L E formulation is proposed in this thesis and it is developed from the principal of virtual work after discretization for a linear increment. The presented formulation is applicable to various material constitutive models and to general loading. A mesh motion scheme is developed based on the transfinite mapping method and is extended to consider the effect of general boundary motion. The scheme is simple and efficient one that does not cause curve fitting errors. It allows the boundary curves to be represented in discrete form and to have discontinuous slopes. Therefore, it is directly applicable to domains with irregular polygons as boundaries. This scheme is shown to be computationally more 148 Chapter 6. CONCLUSIONS AND FUTURE WORK 149 efficient and is straight forward to implement in an A L E formulation. In the thesis, motion of nodes on boundaries, the classification of boundary types and the corresponding detailed motion schemes are presented and shown to be general and efficient in obtaining a high quality mesh. The effect of distorted mesh on the accuracy of the results and the necessity for an efficient A L E mesh motion scheme are demonstrated through numerical examples. It is shown that the developed mesh motion scheme maintains a homogeneous mesh under very large deformation levels. The lDij components calculated from the conventional forward difference algorithm may not vanish under rigid body rotation when the time increment is finite. This may lead to excessive error accumulation in practice. A modified forward difference algorithm or a central difference one provides an alternative for calculating the components of 'D^- in large deformation and large rotation problems. This ensures the objectivity of the 'D tensor with respect to rigid body rotations and maintains objectivity in explicit stress integration algorithm when applied separately. The modification term for the components of the 'D tensor is presented in this thesis and is also shown to maintain objectivity of the tensor with respect to rigid body motions. An alternative stress integration algorithm developed in this thesis is shown to be equivalent to an integration of Truesdell rate equation. The development is based on physical meaning of 2nd Piola-Kirchhoff stress and finite deformation. The difference between the developed algorithm and the one based on purely mathematical transformation is only in the higher order terms of time increment At, and it is demonstrated that both algorithms produce identical numerical results. Two numerical implementation schemes for stress integration coupled with a return mapping method are presented to keep incremental plastic consistency. The proposed schemes may generally be combined with any integration algorithm. It is shown that the integration of Jaumann rate is less sensitive to the location of stress transformation than that of Truesdell rate. The combination of Truesdell rate with transformation A presented in this work generally gives better results than other cases. Chapter 6. CONCLUSIONS AND FUTURE WORK 150 In the updating of material associated properties, the relation between material derivative and referential time derivative is applied to create the relation between increment due to mesh motion and increment due to material motion. Using the developed relation, the material associated properties may be updated directly from material points to mesh points without the need for time consuming interpolations and the creation of a hypothetical mesh as was done by other authors. Combined with local smoothing by the least square method, the updating scheme provides an efficient and accurate routine to update material associated properties. Practical application of the developed A L E formulation is presented in the metal forming area. Simulation of metal extrusion and punch forging processes is performed. It is shown that the load fluctuation and punch size increasement, normally associated with updated Lagrangian formulation, are eliminated. A L E results indicate that the mesh motion algorithm maintains homogeneous mesh throughout the deformation process. It is also shown that description of changing contact boundary conditions is straight forward and easy to perform in A L E formulation. A necking bifurcation of a bar in tension is simulated by an updated Lagrangian method that is reduced from the A L E formulation. Through this example, it is shown that the combination of Truesdell rate with transformation after return mapping gives consistently better results than use of Jaumann rate. One of the most important and complicated metal forming processes is rolling. The steady state condition of rolling of sheet metals is simulated with the Eulerian option, i.e. the particular case of A L E method when the mesh does not move. Application of correct boundary conditions and simulation of neutral points are investigated. The simulation is performed with the general 2D finite element program A L E F E developed in this research without the need to append special routines for the application. No particular interface elements or friction layers are needed. For different ratio of initial thickness of work piece and roll radius cases, the program predicts both double and single peaks of friction-hill of pressure distribution, as indicated by experiments. Rolling pressure, force and torques are in reasonable agreement with the experimental results. Chapter 6. CONCLUSIONS AND FUTURE WORK 1 5 1 Compression between wedge-shaped dies is also simulated by the A L E method. It is noted that although the deformation process is not symmetrical, the A L E method still maintains a homogeneous mesh during deformation and eliminates boundary condition problems caused by updated Lagrangian method. Simulated results are in good agreement with experimental ones. 6.2 F U T U R E W O R K The A L E formulation presented in this research may be extended to simulate dynamic problems. Explicit time integration schemes should be investigated and applied in this area to achieve higher calculation efficiency. Some effort should be concentrated on mesh motion schemes with dynamic effects. The application of A L E in fracture mechanics is also a promising area since it may facilitate the simulation of crack propagation in a continuous manner. However, some details, e.g., tracing material particles and updating material-associated properties at crack tip, are not well developed. The combination of transfinite method with other methods, such as isoparametric interpo-lation method, is another area worth to explore. Such combination makes it possible to take the advantage of both techniques. Therefore, meshes with complicated topologies or excessively distorted meshes may be moved easily to obtain better quality . Other type of elements, like triangular, 8-node quadrilateral element should be incorporated in the program A L E F E , so that user can choose different element types according to the feature of a specific problem. In addition, curved boundaries may be described more accurately. In metal forming, deformation often occurs at elevated temperature, so that the deformation rate may have significant effect. Large deformation normally produces heat and temperature rise. Therefore, rate dependent material model should be applied in the simulation. Fur-thermore, thermal effects should be coupled with the deformation analysis to achieve higher Chapter 6. CONCLUSIONS AND FUTURE WORK accuracy. Bibliography [1] Hill, R., \"Some Basic Principles in the Mechanics of Solids without a Natural Time\", J . Mech. Phy. Solids, Vol.7, pp. 209-225, 1959. [2] Hibbitt, H . D., Marcal, P. V. , Rice, J . R., \"A Finite Element Formulation for Problems of Large Strain and Large Displacement\", Int. J . Solids Struct., Vol.6, pp. 1069-1086, 1970. [3] McMeeking, R. M . , Rice, J . R., \"Finite Element Formulations for Problems of Large Elastic-plastic Deformation\", Int. J . Solids Struct., Vol.11, pp. 601-616, 1975. [4] Maniatty, A. M . , Dawson, P. R., Weber, G. G. , \"An Eulerian Elasto-viscoplastic Formulation for Steady-state Forming Process\", Int. J . Mech. Sci., Vol. 33, pp. 361-377, 1991. [5] Gadala, M . S., Dokainish, M . A. , Oravas, G. A E . , \"Formulation Methods of Geo-metric and Material Nonlinearity Problems\", Int. J . Nume Meth. Engng., Vol.20, pp. 887-914, 1984. [6] Gadala, M . S., Oravas, G. A E . , Dokainish, M . A. , \"A Consistent Eulerian Formu-lation of Large Deformation Problems in Static and Dynamics\", Int. J . Non-linear Mechanics, Vol.18, pp. 21-35, 1983. [7] Abo-Elkhier, M . , Oravas, G. A E . , Dokainish, M . A. , \"A Consistent Eulerian Formu-lation for Large Deformation Analysis with Reference to Metal-extrusion Process\", Int. J . Non-linear Mechanics, Vol.23, pp. 37-52, 1988. [8] Scientific Forming Technologies Corporation, Design Environment for Forming Users Manual (Version 4.1), Scientific Forming Technologies Corporation, Columbus, USA, 1995. [9] Swanson Analysis Systems Inc., Engineering Analysis System User's Manual, Swan-son Analysis Systems Inc., Houston, USA, 1992. [10] Engineering Mechanics Research Corporation, Users Manual for NISA Systems Analysis, Engineering Mechanics Research Corporation, Michigan, USA, 1992. [11] Lee, C. H . , Kobayashi, S., \"Elastoplastic Analysis of Plane-strain and Axisymmetric Flat Punch Indentation by the Finite Element Method\", Int. J . Mech. Sci., Vol.12, pp. 349-370, 1970. 153 Bibliography 154 [12] Lee, E . H . , Mallett, R. L . , \"Stress and Deformation of Metal Extrusion Process\", Comp. Mech. Appl. Mech Engrg., Vol.10, pp. 339-353, 1977. [13] Voyiadjis, C. Z., Foroozesh, M.,\"A Finite Strain Total Lagrangian Finite Element Solution for Metal Extrusion Problems\", Comp. Appl. Mech. Engrg., Vol.86, pp. 337-370, 1991. [14] Lush, A. M . , \"Numerical Grid Generation Used for Remeshing Finite Element Anal-yses of Metal Forming\", In: Numerical Grid Generation in Computational Fluid Me-chanics' 88 ( Sengupta, S., Hauser, J . , Eiseman, P. R., Thompson, J . F . eds), 2nd Conference on Grid Generation in Computational Fluid Dynamics, pp. 1029-1038, Swansea, U K , 1988. [15] Cheng, J . , Kikuchi, N. , \"A Mesh Re-zoning Technique for Finite Element Simula-tions of Metal Forming Process\", Int. J . Nume. Meth. in Engrg., Vol.23, pp. 219-288, 1986. [16] Gelten, I. C. J . M . , Konter, I. A. W. A. , \"Application of Mesh-rezoning in the Updated Lagrangian Method to Metal Forming Analyses, In: Numerical Methods in Industrial Forming Processes (Pittman, J . F . T. eds), pp. 511-521, Pineridge Press, Swansea, U K , 1982. [17] Wang, J . , Gadala, M . S.,\" Remarks on the Finite Element Solution of Large Strain and Forming Problems\", Proceedings of Canadian Congress of Applied Mechanics-C A N C A M 95, pp. 262-263, Victoria, B C , Canada, 1995. [18] Kanto, Y . , Tagawa, G. , \"A Dynamic of Contact Buckling Analysis by the Penalty Finite Element Method\", Int. J . for Nume. Methods in Engrg., Vol.29, pp. 755-774, 1990. [19] Nour-Omid, B. , Wriggers, P., \"A Two-level Iteration Method for Solution of Contact Problems\", Comp. Meth. Appl. Mech. Engrg., Vol.54, pp. 131-144, 1986. [20] Padovan, J . , Moscarello, R., \"Pantographing Self-adaptive G A P Elements\", Com-puters & Structures, Vol.20, pp. 745-758, 1985. [21] Abo-Elkhier, M . , On the Numerical Solution of Nonlinear Problems in Continuum Mechanics, Ph. D. thesis, McMaster University, Hamilton, Canada, 1985. [22] Derbalian, K . A. , Lee, E . H . , Mallet, R. L . , McMeeking, R. M , \"Finite Element Metal Forming Analysis with Spatially Fixed Mesh\", In: Application of Numerical Methods for Forming Process. (Armen, H. eds.) A S M E - A M D Vol.28, pp. 39-47, 1978. Bibliography 155 [23] Oden, J . T. , Braucrdi, H . J . , \" On the Calculation of Consistent Stress Distributions in Finite Element Approximations\", Int. J . Nume. Mech. Engrg., Vol.3, pp. 317-325, 1971. [24] Ghosh, S., Kikuchi, N. ,\" An Arbitrary Lagrangian-Eulerian Finite Element Method for Large Deformation Analysis of Elastic-viscoplastic Solids\", Comp. Meth. Appl. Mech. Engrg., Vol.86, pp. 127-188, 1991. [25] Cough, R., Sharp, R., Otero, I., Tipon, R., McCallen. R., \"Application of A L E Techniques to Metal Forming Simulations\", Advanced Computational Methods for Materials Modeling, American Society of Mechanical Engineers, Applied Mechanics Division, A M D Vol.180, pp. 133-140, 1993. [26] Wang, J . , Gadala, M . S., \"Aspects in the Formulation and Application of Finite Elements to Metal Forming Problems\", Proceedings of the International Conference, Modelling and Simulation in Metallurgical Engineering and Material Science, pp. 577-582, Beijing, P. R. China, 1996. [27] Noh, W. F . , \"A Time-dependent Two-space-dimensional Coupled Eulerian-Lagrangian Code\", In: Methods in Computational Physics (Alderb, B. et al eds), Vol.3, Academic Press, New York, USA, 1964. [28] Wang, J . , Gadala, M.S., \"Formulation and Survey of A L E Method in Nonlinear Solid Mechanics\", Finite Elements in Analysis and Design, Vol.24, pp. 253-269, 1997. [29] Donea, J . , Fasoli-stella, D. P., Giuliani, S., \" Lagrangian and Eulerian Finite Element Techniques for Transient Fluid-structure Interaction Problems\", Trans. SMiRT-4, Paper B l /2 , pp. 115-19, San Francisco, USA, 1977. [30] Belytschko, T. , Kennedy J . M . , \"Computer Methods for Subassembly Simulation\", Nuclear Engrg. Des., Vol.47, pp. 17-38, 1978. [31] Lynch, P. R., O'Neill, K . , \"Continuously Deforming Finite Elements for the Solution of Parabolic Problems with and without Phase Change\", Int. J . Nume. Meth. Engrg., Vol.17, pp. 81-96, 1981. [32] Haber, R. B. , Abel, J . F . \"Contact-slip Analysis Using Mixed Displacements\", J . Engrg. Mech. Vol.109, pp. 441-429, 1983. [33] Huetink, J . , \"Analysis of Metal Forming Processes Based on a Combined Eulerian-Lagrangian Finite Formulation\", In: Numerical Methods in Industrial Forming Pro-cesses ( Pittman, J . F . T. eds), pp. 501-509, Pineridge, Swansea, U K , 1982. Bibliography 156 [34] Huetink, J . , Vreede, P. T. , Lugt, J . V . D., \"Progress in Mixed Eulerian-Lagrangian Finite Element Simulation of Forming Process\", Int. J . Nume. Meth Engng., Vol.30, pp. 1441-1457, 1990. [35] Benson, D. J . , \"An Efficient, Accurate, Simple A L E Method for Nonlinear Finite Element Programs\", Comp. Meth. Appl. Mech. Engrg., Vol. 72, pp. 305-350, 1989. [36] Haber, R. B. , \"A Mixed Eulerian-Lagrangian Displacement Model for Large-deformation Analysis in Solid Mechanics\", Comp. Meth. Appl. Mech. Engng., Vol.43, pp. 277-292, 1984. [37] Yamada, T . , Kikuchi, F . , \"An Arbitrary Lagrangian-Eulerian Finite Element Method for Incompressible Hyper elasticity\", Comp. Meth. Appl. Mech. Engrg., Vol.102, pp.149-177, 1993. [38] Schreurs, P. J . G. , Veldpaus, F . E . , Brekelmans, W. A. M . , \"An Arbitrary Eule-rian Lagrangian Finite Element Model for the Simulation of Geometric Nonlinear Hyperelastic and Elasto-plastic Deformation Processes\", In: Numerical Methods in Industrial Forming Processes (Pittman, J . F . T. eds), pp. 491-500, Pineridge, Swansea, 1982. [39] Schreurs, P. J . G. , Veldpaus, F. E . , Brekelmans, W. A. M . , \"Simulation of Forming Processes, Using the Arbitrary Eulerian-Lagrangian Formulation\", Comp. Meth. Appl. Mech. Engng., Vol.58, pp. 19-36, 1986. [40] Hughes, T. J . R., Liu, W. K. , Zimmerman, T. K.: \"Lagrangian-Eulerian Finite Element Formulation for Incompressible Viscous Flows\", Comp. Meth. Appl. Mech. Engrg., Vol.29, pp. 329-349, 1981. [41] Liu, W. K . , Chang, H . , Chen, J . , Belytschko, T. , \"Arbitrary Lagrangian-Eulerian Petrov-Galerkin Finite Elements for Nonlinear Continua\", Comp. Meth. Appl. Mech. Engrg., Vol.68, pp. 259-310, 1988. [42] Liu, W. K . , Chen, J . S. Belytschko, T. , Zhang, Y . F . ,\" Adaptive A L E Finite El-ements with Particular Reference to External Work Rate on Frictional Interface\", Comp. Mech. Appl. Mech. Engrg., Vol.93, pp. 189-216, 1991. [43] Ghosh, S., Kikuchi, N. , \"Finite Element Formulation for the Simulation of Hot Sheet Metal Forming Process\", Int. J . Engrg. Sci., Vol.26, pp. 143-161, 1988. [44] Ghosh, S., \"Finite Element Simulation of Some Extrusion Processes Using the Arbi-trary Lagrangian-Eulerian Description\", J . Mat. Shaping Technol., Vol.8, pp. 53-64, 1990. Bibliography 157 [45] Huerta, A. , Casadei, F ., \"New A L E Applications in Non-linear Fast-transient Solid Dynamic\", Engineering Computations, Vol.11, pp.317-345, 1994. [46] Bathe, K . J . , Ramm, E . , Wilson, E . L . , \"Finite Element Formulation for Large deformation Dynamic Analysis\", Int. J . Nume. Meth. Engrg., Vol.9, pp. 353-386, 1975. [47] Bathe, K . , Finite Element in Engineering Analysis, Prentice-Hall Inc., New Jersey, USA, 1982. [48] TruesdeU, C. , Toupin, R., \"The Classical Field Theories\", In: Encyclopedia of Physics, Vol.III/1, pp. 552-610, Spriner, Berlin-New York, 1960. [49] Truesdell, C , \"The Simplest Rate Theory of Pure Elasticity\", Communications in Pure and Applied Mathematics, Vol.8, pp. 123-132, 1955. [50] Dabke, P.D., Haque,L, Srikrishna, M . and Jackson, J .E . , \"A Knowledge-Based Sys-tem for Generation and Control of Finite-Element Meshes in Forging Simulation\", J . Mat. Eng. and Perf., Vol.1, pp. 415-428, 1992. [51] Lush, A. , \"Numerical Grid Generation Used for Remeshing Finite Element Anal-yses of Metal Forming\", In: Numerical Grid Generation in Computational Fluid Mechanics' 88 (Sengupta, S., et al, eds.), pp.1029-1038, Swansea, U K , 1988. [52] Sparis, P. D., \"A Method for Generating Boundary-orthogonal Curvilinear Coordi-nate Systems Using the Biharmonic Equation\", Journal of Computational Physics, Vol.61, pp. 445-462, 1985. [53] Harber, R., Shephard, M . S., Abel, J . E . , Gallagher, R. H . , Greenberg, D. P.,\"A General Two Dimensional, Graphic Finite Element Preprocessor Utilizing Discrete Transfinite Mapping\", Int. J . Nume. Meth. Eng. Vol.17, pp. 1015-1044, 1981. [54] Gadala, M.S., Wang, J . , \"A Practical Procedure for Mesh Motion in Arbitrary Lagrangian-Eulerian Method\", Submitted for publication in Engineering with Com-puters July, 1997. [55] Gadala, M . S., Wang, J . , \"Numerical Treatment of Mesh Distortion Problems in F E A \" , Submitted for publication in Int. J . of Computers and Structures, Sept., 1997. [56] Hughes, T. J . R., Winget, J . , \"Finite Rotation Effects in Numerical Integration of Rate Constitutive Equations Arising in Large-deformation Analysis\", International Journal for Numerical Methods in Engineering, Vol.15, pp. 1862-1867, 1980. Bibliography 158 [57] Hibbitt, Karlsson and Sorensen, Inc., ABAQUS Theory Manual, Hibbit, Karlsson & Sorensen, Inc., Pawtucket, USA, 1992. [58] Pinsky, P. M . , Ortiz M . , Pister, K . S., \"Numerical Integration of Rate Constitutive Equations in Finite Deformation Analysis\", Computer Methods in Applied Mechan-ics and Engineering, Vol.40, pp. 137-158, 1983. [59] Ortiz, M . , Popov, E . P., \"Accuracy and Stability of Integration Algorithms for Elastoplastic Constitutive Relation\", International Journal for Numerical Methods in Engineering, Vol.21, pp. 1561-1576, 1985. [60] Ortiz, M . , Simo, J . C , \"An Analysis of a New Class of Integration Algorithms for Elastoplastic Constitutive Relations\", International Journal for Numerical Methods in Engineering, Vol.23, pp. 353-366, 1985. [61] Wang, J.,Gadala, M.S., \"Computational Implementation of Stress Integration in F E Analysis of Elasto-plastic Large Deformation Problems\", Submitted for publication in Engineering Computations, March, 1997. [62] Hinton, E . , Campbell, J . S., \"Local and Global Smoothing of Discontinuous Finite Element Functions Using a Least Square Method\", Int. J . of Nume. Mech. Eng., Vol.8, pp. 461-480, 1974. [63] Oden, J . T. , Braucnli, H. J.,\"On the Calculation of Consistent Stress Distributions in Finite Element Application\", Int. Jour. Nume. Meth. Engrg., Vol.3, pp. 317-325, 1971. [64] Malvern, L. E . , Introduction to the Mechanics of a Continuous Medium, Prentice-Hall Inc, New Jersey, USA, 1969. [65] Bathe, K . J . , Cimento, A. P., \"Some Practical Procedures for the Solution of Non-linear Finite Element Equations\", Computer Methods in Applied Mechanics and Engineering, Vol.22, pp.59-85, 1980. [66] Owen, D. J . R., Hinton, E . , Finite Element in Plasticity, Theory and Practice, Pineridge Press Limited, Swansea, U K , 1980 [67] Cook, R. D., Malkus, D. S., Plesha, M . E . , Concepts and Applications of Finite Element Analysis, John Wiley & Sons, New York, USA, 1989. [68] Burnett, D. S., Finite Element Analysis, from Concepts to Applications, Addison-wesley Publishing Company, Menlo Park, California, USA, 1987. [69] Huerta, A. , Casadei, F . , \"New A L E Applications in Non-linear Fast-transient Solid Dynamic\", Engineering, Computations, Vol.11, pp. 317-345, 1982. Bibliography 159 Yamada, Y . , Hirakawa, FL, \"Large Deformation and Instability Analysis in Metal Forming Processing\", Application of Numerical Methods to forming Processes, A S M E , AMD-Vol.28, (Edited by H . Armen), pp.27-38, 1978. Farmer, L. E . , Oxley, P. L. B., \"A Slip Field for Plane Strain Extrusion of a Strain-hardening Material\", J . Mech. Phys. Solids, Vol.19, pp.369-388, 1971. Lyman, T. , Metals Handbook, Vol.1, Properties and Selection of Metals, American Society for Metals, Metals Park, Ohio, USA, 1961. Li , G . , Kobayashi, S., \"Rigid-plastic Finite Element Analysis of Plain Strain Rolling\", J . of Eng. for Industry, Vol.104, pp.55-64,1982. Hwang, S.M., Joun, M.S., \"Analysis of Hot Strip Rolling by a Penalty Rigid-viscoplastic Finite Element Method\", Int. J . Mech. Sci., Vol.34, pp.971-985, 1992. Liu, C , Hartley, P., Sturgess, C. E . N. , Rowe, G. W., \"Finite Element Modelling of Deformation and Spread in Slab Rolling\", Int. J . Mech. Sci., Vol.29, pp.271-283, 1987. Al-Salehi, F . A. R. , Fibrank, T. C , Lancaster, P. R., \"An Experimental Determi-nation of the Roll Pressure Distributions in Cold Rolling\", Int. J . Mech. Sci. Vol.15, pp.693-710, 1973. Wusatowski, Z., Fundamentals of Rolling, Pergamon Press, New York, USA, 1969. Chitkara, N. R., Johnson, W., \" Plane Strain Compression of pre-shaped material between Wedge-shaped Dies\", Int. J . Mech. Sci., Vol.4, ppl51-164, 1974. Appendix A Material and Referential Time Rate A physical quantity'/ = /('x) at Material Reference System (MRS) point P ( ' : c i , ' x2,' x3) at time t becomes t + A t f = / ( ' + A 'x) at time t+At as the point P moves to Q ( t + A t x u t + A t x 2 , t + A t x3 The MRS change of ' / may be expressed by: A x / = / ( ' + A ' x , i + A i ) - / ( ' x , i ) (A.l) Similarly, the change of ' / at the corresponding Computational Reference System (CRS) point P^Xif X2,1 X3) w n e n t n e P o i n t moves to R(t+Atx\,t+At X 2 , t + A t X3) is denoted by: Ax/ = /rA'x,* + A0-/(fX,*) (A.2) To derive the relation between A x / and Axf, we assume that the MRS and CRS points coincide at point P at time t, as shown in Figure (A.l), so that: and, at time t + At, f(tX,t) = f(tx,t) = tf (A.3) t + A ' x =' x + A x (A.4) t + A tX =*X + &X (A.5) 160 Appendix A. Material and Referential Time Rate 161 Figure A . l : The description of relation between mesh motion and material motion Appendix A. Material and Referential Time Rate 162 Assuming that the MRS point located at R at time t + At comes from a material point P*(tx\,t a;*./ X 3 ) with displacement Ax*, we have: = * x * + Ax* = t + A V (A.6) f(t+AtX,t + At) = f(t+Atx,t + At) (A.7) Substitute Equation (A.3) and (A.7) into Equation (A.2): Axf = fC+^^t + A^-fC+^t + A^ + fC+^t + A^-fC^t) = A x / + ( / (t + A t x * , t + At)- / ( t + A t x , t + At)) = A x / + V / (t + A f x , t + At)\u00C2\u00BB ( i + A t x * - t + A t x) + h.o.t. where V is gradient symbol and h.o.t. is higher order term of ( *+ A t x * \u00E2\u0080\u0094 t + A t x ) . But, t+A* x* _t+At x _ 4 + A * ^ _*+At x = 'x + A x - ' x - A x = A x - A x Therefore, A x / = Axf + Vf{t+Atx,t + At)*(AX-Ax) + h.o.t. Where h.o.t. is higher order term of (A% \u00E2\u0080\u0094 Ax) . If h.o.t. is ignored, we have: A * / = A */ + | ^ K - ^ ) (A-8) Dividing the above equation by At and taking the limit when At \u00E2\u0080\u0094> 0, 7A = f f + ^Cvl - *vk) (A.9) Appendix B Updated Lagrangian Formulation An updated-Lagrangian formulation given by Rice [3] is based on the following equation: / v ( v ^ - I v \u00E2\u0080\u009E , ( , W - | | ^ ) ) ^ = jtvtfiB-8vidtV + f tfS-SvutS (B.10) where 'Aj = \ (\u00C2\u00A7r^ + ft^) is the deformation rate at time t and ' r ^ is Jaumann rate of Kirchoff stress which is defined by: where 'cr^ is Cauchy stress and ' J is the ratio of volume in the current state to the volume in reference state. The time rate can be written as: ' Ttj ~ TJt<7ij + TJt\u00C2\u00B0'ij using the Jaumann rate for both 'r^- and V^-, we get: t - J t T t _ I t Tt-J T *3 ran + ' J4 < ( B . l l ) Since ' r ^ is symmetric, the right hand side of Equation (B.10) may be expressed by: f (t j d S v i 1 (dSvk dSvt\ t ldlvk d'vA d8vkt dW\j U { T i j &Xi 2 ^ 0*3; + a\u00C2\u00ABaj J a i j {d^ + d*xk) d*Xi ^frxj f dSvi (t-J (&v, , d*vk\ *ajk ( 5 ^ f c \ \ = lv\u00C2\u00A5x~AT--^{o^k + o^J+^{d^k-\u00C2\u00A5x-Jrv (B-12) 163 Appendix B. Updated Lagrangian Formulation 164 Substituting * J = [48], and choosing the current state as a reference so that * J = 1. The Equation (B.l 1) will take the form: V J - *\u00C2\u00AB\u00E2\u0080\u00A2 9 t V k 4 - rrJ Substituting the above equation in Equation (B.12), we get: LHS = ( p ( ^ + - ^ - ' ^ ( p . + p L ) + ' ^ ( p . . p L ) ) < t V hv otxj \ J otxk 2 ya'xfc otXj) 2 \otxk dtxi J J So the updated Lagrangian formulation Equation (B.10) takes the following form: dSvj L .j _tvik( #S dW\ lfjk[ d%_ _ dW\ t dW\ j v d'xj \ a i j 2 [d'xk + d'xj + 2 [frxk d'xij + (Tii^xk) = JtvSv\f?*V + jtg8v\f?*S (B.13) Equation (B.13) is same as Equation (2.25) in Chapter 2 except the load correction terms introduced in Equation (2.25). Appendix C Mapping Between Material and Mesh Domain It is important to ensure that the mesh motion scheme proposed Equation (3.3) guarantees one-to-one mapping between material points and mesh points, i.e., TJC = ^ 0. This may be verified as follow: At time t \u00E2\u0080\u0094 At, assume '~A ta:; =*-A t ^ i , i.e., material points and mesh points coincide and from t \u00E2\u0080\u0094 At to t the material velocity is t _ A t f i , then, tx t-At, lxi + At'-^vi = l-atXi + At1 Xi= Xi + At Vi = Xi + At (a; + b{i) v{i) j From which we obtain a relation between txi and tXi as: H = t X i - A t ( a i + b{i)t-Atv[i))+Att-Atvi = *Xi - Ate* + (l - 6 ( i )) At'-^vu and the determinant of the Jacobian transformation becomes: dlXi dlXj 8^ + (l - 6 ( 0) At d'Xj No summation on (i). Since * AtVi is independent of we have d d t ^ x ) tjc I u 0 and: (C.14) (C.15) (C.16) (C.17) 165 Appendix D Stress Increment Procedure We consider decomposing the general deformation into two steps; straining and rigid body motion. In the first step, due to straining, the Cauchy stress would have an increment of Ao-ij = dju'DuAt (D.18) where, Ciju is the element constitutive tensor, At is the time increment and *Dki is the rate of deformation tensor. The stress after the first step would be: t+At*l? = V \u00C2\u00AB + Aery = \+**S\u00E2\u0084\u00A2 (D.19) where, t+Ats[p is the 2nd Piola-Kirchhoff stress at an imaginary intermediate state with respect to the configuration at time t and t+Ato-\p is the corresponding Cauchy stress. The 2nd Piola-Kirchhoff stress is related to Cauchy stress by: t+*t 1 d f + A t a : i t + A t dt+AtXj t+At-g 3 xm d xn and since the second step is only a rigid body motion, then, l+AtSmn = ?\"S\u00E2\u0084\u00A2 (D.21) Substituting Equation (D.21) and (D.19) into Equation (D.20), we get: 1 d t + A t x i t d t + A t x j 1 dt+AtXi. dt+At t+At\u00E2\u0080\u009E _ \u00E2\u0080\u00A2*- w ~lt\u00E2\u0080\u009E *\u00C2\u00B1 0 i ' \" m 991 1 0 <+*'F d*xm \" m n d*xn ' t+AtY d*xm 77171 d*xn Equation (D.22) is same as the one developed by Pinsky [58] and both are consistent with Truesdell stress rate constitutive equation. 166 "@en . "Thesis/Dissertation"@en . "1998-05"@en . "10.14288/1.0080928"@en . "eng"@en . "Mechanical Engineering"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "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."@en . "Graduate"@en . "Arbitrary Lagrangian-Eulerian method and its application in solid mechanics"@en . "Text"@en . "http://hdl.handle.net/2429/8616"@en .