<http://www.w3.org/2009/08/skos-reference/skos.html#note>
"ARBITRARY LAGRANGIAN-EULERIAN APPLICATION  IN SOLID  M E T H O D A N D ITS  MECHANICS  By Jin Wang B. A . Sc., M . A . Sc. University of Science and Technology Beijing  A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T O F T H E REQUIREMENTS FOR T H E D E G R E E OF D O C T O R O F PHILOSOPHY  in T H E F A C U L T Y OF G R A D U A T E STUDIES D E P A R T M E N T O F MECHANICAL ENGINEERING  We accept this thesis as conforming to the required standard  T H E UNIVERSITY O F BRITISH COLUMBIA  April 1998 \u00A9 Jin Wang, 1998  \u000CIn presenting  this  degree at the  thesis  in  partial fulfilment  of  University of  British Columbia,  I agree  freely available for reference copying  of  department  this or  publication of  and study.  thesis for scholarly by  this  his  or  her  Department The University of British Columbia Vancouver, Canada  DE-6 (2/88)  Ajrfd  l6&gt;  requirements that the  I further agree  purposes  representatives.  may be It  thesis for financial gain shall not  permission.  Date  the  that  advanced  Library shall make it  by the  understood be  an  permission for extensive  granted  is  for  that  allowed without  head  of  my  copying  or  my written  \u000CAbstract  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  \u000Ctreatment for rate of deformation tensor  Dij is presented to maintain incremental objectivity  t  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, A L E F E , 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 A L E F E . Various benchmark and practical problems are simulated by the developed program. Practical 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  \u000CTable of Contents  Abstract  ii  List of Tables  ix  List of Figures  x  Nomenclature  xiii  Acknowledgement  xix  1 INTRODUCTION  1  1.1  NONLINEAR FINITE E L E M E N T METHODS  1  1.2  L A G R A N G I A N FORMULATIONS  2  1.2.1  Sample Examples  2  1.2.2  Limitation of Lagrangian Type Formulation  12  1.3  E U L E R I A N FORMULATION  14  1.4  ARBITRARY LAGRANGIAN-EULERIAN M E T H O D  15  1.5  OBJECTIVES A N D SCOPE  16  2 A R B I T R A R Y L A G R A N G I A N 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  \u000C2.3 3  DISCUSSION A N D CHARACTERISTICS OF PROPOSED FORMULATION  31  NUMERICAL IMPLEMENTATION OF A L E  35  3.1  INTRODUCTION A N D 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  M E S H MOTION S C H E M E 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  3.5  3.6  3.7  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  INCREMENTAL OBJECTIVITY A N D N U M E R I C A L T R E A T M E N T OF * D  U  56  3.5.1  General Considerations  56  3.5.2  Forward Difference Algorithm to Calculate D y  59  3.5.3  Forward Difference with a Modification Term  59  t  STRESS TRANSFORMATION A N D 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  UPDATING OF MATERIAL ASSOCIATED PROPERTIES  64  v  \u000C4  Basic Updating Scheme  64  3.7.2  Material Associated Properties at Nodal Points  65  TWO-DIMENSIONAL 4.1  4.2  4.3  4.4 5  3.7.1  A L EP R O G R A M (ALEFE)  69  PROGRAM ALEFE  69  4.1.1  Features of A L E F E  69  4.1.2  Flow Chart of the Main Program  71  MAIN SUBROUTINES  74  4.2.1  Functions of Subroutines  74  4.2.2  Flow Chart of Mesh Motion  76  C O N V E R G E N C E CRITERIA  76  4.3.1  General Discussion  76  4.3.2  Displacement Criterion  78  4.3.3  Force Criterion  4.3.4  Energy Criterion  V  78 79  C O U P L E D DISPLACEMENT CONDITION  80  NUMERICAL EXAMPLES  83  5.1  OVERVIEW  83  5.2  ROTATION OF A N E L A S T I C A L L Y D E F O R M E D E L E M E N T  83  5.3  ELASTIC-PLASTIC S M A L L DEFLECTION ANALYSIS OF A SQUARE  5.4  ELEMENT  85  PURE BENDING OF A R E C T A N G U L A R 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  \u000C5.5  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 B A R  105  5.7  SHEET M E T A L EXTRUSION  Ill  5.7.1  Extrusion through a Curved Die  Ill  5.7.2  Extrusion through a Straight Die  114  5.7.3  Discussion  121  5.8  5.9  6  P U N C H FORGING PROCESS  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  COMPRESSION B E T W E E N 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  CONCLUSIONS 6.1  AND FUTURE WORK  CONCLUSIONS  148 148  vii  \u000C6.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  viii  \u000CList 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  \u000CList of Figures  1.1  The meshes before and after deformation  3  1.2  Not updated boundary condition  4  1.3  Distribution of stress a 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 P  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  l  x  \u000C5.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 a from A L E  99  y  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 a from A L E and reference [12]  115  5.23 Comparison of r  116  x  xy  from A L E and reference [12]  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  \u000C5.29 The definition of geometry in rolling 5.30 Finite element model (Case 1: R/h  0  123 = 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/h  133  5.36 Comparison of rolling pressure, large R/h  134  5.37 Comparison of rolling force  135  5.38 Comparison of rolling torque  135  0  0  5.39 Equivalent plastic strain (Case 1: R/h  0  = 12.5, Reduction = 14.17%)  5.40 Current yield stress at different steps (Case 5: R/h  Q  136  = 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  161  The description of relation between mesh motion and material motion  xii  \u000CNomenclature  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  d  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  p  f  the value of / at computational referential point  f?  components of body forces  ff  components of surface tractions  /-  unsmoothed quantity at node a  c  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  Xlll  \u000Chu  values of shape function for displacements in i direction at D.O.F. I  h  smoothing shape function at node a  H  hardening modulus  i  index for directions / coordinate axes (i =  /  index for D.O.F.  I  unit matrix  J  ratio of material volume to the volume in referential state  a  J  l,2,3ora:,y,z)  Jacobian matrix  c  k,k  current and initial shear yield stresses, respectively  K , ku  stiffness matrix related to 'v  K , kjj  stiffness matrix related to v  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  n  number nodes per element  rii  components of the unit vector normal to boundary surface  N  total number of D.O.F in the model  0  c  n  4  c  N  number of sampling points in an element  p  normal pressure  pi  external nodal force at D.O.F. /  p^  external nodal force in k  p  equivalent load vector  q  distributed tensile force  qi  equivalent nodal force at D.O.F. /  s  th  iteration at D.O.F. /  xiv  \u000Cqj^  equivalent nodal force in k  v  normalised parametric coordinate 0 &lt; r &lt; 1  TI  residual nodal force at D.O.F. /  th  residual nodal force in k  th  iteration at D.O.F. /  at D.O.F. /  R , Rij  rotation tensor  s  normalised parametric coordinate 0 &lt; s &lt; 1  S  surface of deformed body on which tractions are prescribed  Su  2  nd  Piola Kirchhoff stress  Sip  2  nd  Piola Kirchhoff stress at an intermediate state  t  time  td,tf,t  displacement, force and energy tolerances, respectively  T  transformation matrix  e  U  total radial displacement  U  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\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\u00B0  computational reference velocity vector  v1  mesh velocity in direction i  r  x  Vi( )  material velocity of node a at direction i in global coordinates  i(a)  material velocity of node a at direction i in local coordinates  a  v  xv  \u000Cv\\  mesh velocity at D.O.F. J  v'  material particle velocity vector in local coordinate system  l  V  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 '  displacement increment in k  A /  change of physical quantity due to material motion  A&lt;Tij  increment of Cauchy stress  Af  change of physical quantity due to computational reference system motion  6  angle for coupled displacements  n  tool angle for wedge shaped dies  5  th  x  x  e  iteration at D.O.F. /  equivalent plastic strain,  p eq  &amp;(r, s)  boundary velocity of a geometric domain  crij  components of Cauchy stress tensor  v\\ 5 0\"2, o-  the principal Cauchy stresses  3  o'yieidi o-yi ido current and initial yield stresses, respectively e  cr'l  Jaumann rate of Cauchy stresses  xvi  \u000C&lt;rj-J  Truesdell rate of Cauchy stresses  a[-  deviatoric Cauchy stress tensor  o~\\p  Cauchy stresses at an intermediate state I  o-W  equivalent stresses at an intermediate state /  r  friction stress  r'  Jaumann rate of Kirchoff stress  H  friction coefficient  v  Poisson's ratio  4&gt;  objective function in the least square method  4&gt;i{r,s)  boundary coordinate of a geometric domain  X  position vector of referential coordinate  Xi  mesh point coordinate in axis i  u&gt;  angular velocity  V/  gradient of /  J  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  \u000Cmaterial time rate of a quantity A  computational reference time rate of a quantity  T  transformation of a matrix  xviii  \u000CAcknowledgement  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 U B C - U G F 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.  xix  \u000CChapter 1  INTRODUCTION  1.1  N O N L I N E A R FINITE 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  \u000CChapter 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 F E 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  LAGRANGIAN  1.2.1  FORMULATIONS  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 10 MPa, Poisson's ratio v = 0.33, 4  hardening modulus H = 138 MPa with initial yield stress o- i o = 90 MPa. Axi-symmetric yie  d  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 d = 0.15w. The simulation was carried out using both NISA and ANSYS programs. The mesh p  topology before and after deformation is presented in Figure 1.1. The radial displacement U , T  the second principle stress cr , and the equivalent plastic strain e  p  2  eq  are used for comparison of the  \u000CChapter 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\u20224-  ~~1  \\.  Punch Size Increasement  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 10 MP 4  v = 0.3, strain hardening i f = 1.1 x 10 MPa and initial yield stress o- i 3  yie  a, Poisson's ratio = 4.0 x 10 MP 2  d0  a.  The metal extrusion problem is modelled by plane strain 8-node isoparametric elements and  \u000CChapter 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 U and x  stress a , both in the extrusion direction, are examined. The distributions of U and a x  x  x  from  NISA and ANSYS are very similar and the values of U are close to each other with a difference x  of less than 0.4%. The values of a are, however, quite different; and show a difference of up x  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 A N S Y S , are significantly different. As an example, distribution of a on the outer x  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.  \u000CChapter 1.  INTRODUCTION  5  Ox(MPa)  3.5E2  0.0E0  -3.5E2  -7.0E2 -  -10.5E2 -  -14.0E2  , i, . , , , -2.0 -1.5  ^  .  .  i  .  -1.0  .  \u2014\u2022\u2014  Ox(MPa) N o G a p Element  \u2014A\u2014  Ox(MPa) G a p Element  , .i. , . . -0.5  i  0  .  .  .  .  i .... i 0.5 1.0  Figure 1.3: Distribution of stress &lt;r on outer surface x  \u2122a.  \u000CChapter 1.  INTRODUCTION  6  In simulation of the punch indentation, once the material point under the punch corner moves out of the corner, it should be allowed to move up. However, the Lagrangian formulation can not update the boundary condition automatically. This is an important feature needed in handling most of metal forming problems. Also, it may not be feasible to manually update boundary conditions after each incremental step. This deficiency has the effect of increasing the punch size during the deformation, as shown in Figure 1.1. For metal extrusion, the situation is worse. Figure 1.2 shows a typical deformed mesh (superimposed on the original one) obtained from such programs. The reduction in thickness depicted by the program does not agree with the specified value because of the absence of geometric boundary updating in the program. Thus the predicted deformation is unacceptable.  (iii) Plane Strain Forging  ,  Plane strain forging is also simulated by NISA and ANSYS. The material model is chosen to be an elasto-plastic one with Young's modulus E = 200.0 GPa, Poisson's ratio v = 0.30, initial yield stress &lt;T i  yie d0  = 250.0 MPa and plastic modulus i f = 1.0 GPa.  deformed by a rigid frictionless tool with a prescribed vertical velocity.  The body is  Only a quarter of  the domain is studied because of symmetry, and plane strain condition is assumed. The finite element model used in NISA and ANSYS is shown in Figure 1.4. The analysis is performed for up to 60% reduction in the height of the original piece to show the effect of large deformation. Figure 1.5 shows the results from ANSYS [9] and NISA [10].  Because the boundary  conditions are not updated automatically, the tool or punch width will be actually increasing with the deformation process as in the indentation example. It is also noted that elements at the punch corner are highly distorted because of large deformation. Increasing the deformation level causes more severe element distortions and the program terminates the analysis. The use of more advanced capabilities such as the point-to-surface gap element in ANSYS may eliminate some of the above problems such as the punch size one.  These elements have their own  \u000CChapter 1.  INTRODUCTION  Y  i  V  * 12 mm  X \\  \\\\ \\ \\ \\ \\ \\  D  C  10 mm  E  7  0  A  30 mm  B  X  Figure 1.4: The finite element model problems, however, such as convergence difficulties, user definition of stiffness parameters and violation of local equilibrium at the punch corner. For this particular example, contact elements were used on the tool work-piece interface and several trial cases with different reduction ratios, different tolerance limits and different stiffness parameters were performed. All cases tried ceased to give successful convergence beyond a reduction of 30%. The same process is simulated by DEFORM2 [8] program, which is one of the most widely used metal forming programs because of its capability of automatic remeshing and updating contact boundary conditions. The finite element model used in DEFORM2 is shown in Figure 1.6. A much denser mesh was required to achieve convergence for the same deformation level obtained in NISA and ANSYS.  Figure 1.7 shows the deformed shape from DEFORM2.  Although the punch size increment is eliminated by remeshing, it may be seen that the mesh around the punch is highly distorted and the punch corner is \"cut off\". More importantly,  \u000CChapter 1.  INTRODUCTION  8  Punch Size Increasement  ANSYS S.2 HUG Z9 1996 12:32:13 DISPLACEMENT STEP=1 SUB =750 TIME=1 RSYS=0 DMX =.018891 \u00ABDSCA=1 2U =1 DIST-.02SB17 XF -.0Z1379 YF =.005005  60* deformation, ANSYS  ANSYS  EMRCNlSA/TllSPLAY  TIME: .1IXKHIL+01 60 % DEFORMATION, WITHOUT T1MEAMP. CAUCI1Y STRESS  \\UGnXMUA6:$0 ROTJ x ROTJ  -0  NISA  Figure 1.5: The deformed shapes by updated Lagrangian method  \u000CChapter 1.  INTRODUCTION  \u2022 12 mm  i  \u2022+  30 mm \u2022&gt;  Figure 1.6: F i n i t e element model used i n D E F O R M 2 simulation  \u000CChapter 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 loadfluctuationthat is pertinent to the updated Lagrangian formulation. Similarfluctuationis 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 \u2014j..r.r.|..T.r.rT..rr.y.T..rr.\u00A5.T.T.r^  -8.7  2.9  14.5  26.0  37.6  Figure 1.7: Deformed mesh from D E F O R M 2  (mm)  \u000CChapter 1.  INTRODUCTION  2.200  1.760  1.320  11  j  A  I ,.  j ;  / f  /  T  i  j \\  . j . .  1  r  !  A/  / '  o  T3  0.440  -\u2022  \\  4 , 1 ' '  :  0.0  0.12  0.00 0.24  0.36  0.48  0.60  Reduction  Figure 1.8: Load-height reduction result from D E F O R M 2  \u000CChapter 1.  1.2.2  INTRODUCTION  12  Limitation of Lagrangian Type Formulation  Most of the available commercial F E 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  \u000CChapter 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 F E 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  \u000CChapter 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  EULERIAN  FORMULATION  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 materialassociated properties (e.g., strain, stress) of the material points momentarily occupying the  \u000CChapter 1.  INTRODUCTION  integration points at the beginning of each step.  15  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  ARBITRARY LAGRANGIAN-EULERIAN  METHOD  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  \u000CChapter 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 loadfluctuations,may describe precisely any contact boundary conditions and make boundary condition updating no longer necessary after each incremental step.  1.5  OBJECTIVES A N D SCOPE A L E method has the potential to eliminate problems caused by Lagrangian or Eulerian  methods in simulation of generalfinitedeformation 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 historydependent 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  \u000Cer 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 consistency 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 2Dfiniteelement program and to simulate bench mark problems and metal forming processes to show the features of A L E method.  \u000CChapter 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 (ALE) F O R M U L A T I O N  2.1  SURVEY OF ALE METHODS 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 twodimensional hydrodynamics problems with moving fluid boundaries. Later, the A L E method was introduced into the finite element method [29] fromfluidmechanics for modelling the solidfluid 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 afiniteelement 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  \u000CChapter 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  \u000CChapter 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  \u000CChapter 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:  '/\u2022 = 7 + A  (2.1)  (J X{  where xi is material coordinate, vi and v^ are material and mesh velocities individually. t  t  t  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 reformulated 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 deformation 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  \u000CChapter 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. d*vi dt  dt  \u2022x  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 introduced 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  \u000CChapter 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 A L E 3 D 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.  \u000CChapter 2.  2.2  ARBITRARY  LAGRANGIAN  EULERIAN  (ALE) FORMULATION  24  CONSISTENT A L E FORMULATION  2.2.1  Geometric and Kinematic Description  B o u n d a r y of C o m p u t a t i o n a l D o m a i n  t  F  X3,  X3  l  Figure 2.1: Schematic diagram for domains and mapping in A L E description  Assuming a material particle P^xif  x * x ) in a material reference system (MRS) at time 2  3  t has velocity '^(t = 1,2,3). The particle moves to Q(  x,  t+At  Similarly, a grid point P\u00B0( xi, t  t  (CRS) at time t has velocity  X2? X3)  m  x,  t+At  t+At  1  2  x ) at time t + At. 3  the corresponding computational reference system  = 1,2,3) and moves to Q { Xi, c  t+At  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  \u000CChapter 2. ARBITRARY  LAGRANGIAN  EULERIAN  (ALE) FORMULATION  25  the material domain and computational domain at any time t, i.e., for each:  'xi  (i = 1 , 2 , 3 )  = x ( xu X2, X3) t  t  t  i  (2.3)  we have only one:  (i = 1 , 2 , 3 )  *Xi = Xi{ xi, x , x ) t  t  t  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 rii is the component of the unit vector normal to the boundary. The physical interpretation t  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 xi t  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 Xit  Assume that the material time derivative and the computational reference time derivative are denoted by a superscripted dot \"\u2022\" 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.  \u000CChapter 2. ARBITRARY  2.2.2  LAGRANGIAN  EULERIAN  (ALE) FORMULATION  26  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  Jt+Aty  = [  f?8u d V+  t+At  [  t+At  i  Jt+Aty  J  J  \"\"ffSuJ+^S  Jt+At  1  J  (2.6)  i  '  v  S  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,  aij  t+At  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^^  =  d^x')  + j  where  (  2  '  7  )  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 f?  in Equation (2.6), respectively.  t+At  t + A t  ff  t+At  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:  / Jt+Aty  t + A t  0-^+^0*+ V At  J  =  j  \u00AB \"  t+Aty  t+AtrS  _  t+At_ t+At  * ( i ^  a  Jt+Aty  J  2 \\O A  V  d  t + A t  t + A t  Xj  a x  +  J ^ \\ d \u00BB \" V O  t + A t  XiJ  and  \u000CChapter 2. ARBITRARY  where d Aj t+At  LAGRANGIAN  is the projection of a^  as normal direction, rij Xj,  /  (ALE) FORMULATION  27  5 on the coordinate plane having axis Xj or l  are the direction cosines of area element d S  t+At  t + A t  +At  EULERIAN  t + A t  Xj  with axis Xj or  t+At  l  the Equation (2.6) may be written in the following form: T  +  A  T  Jt+At  V  ^ ^ ~ ^  +  A  T  t+At  = I  V  Qt+^Xj  /JW 'V+ + A  a 8u d A  t+At  /  Jt+At  Jt+&amp;t .  V  t+At  ij  i  j  J  A  (2.8) 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 P and P in Figure 2.1 be the same point. This transformation is only c  due to the change in or motion of the computational reference domain. We have, dSuj  _  dSuj  d'+^xj  ~  dtxkfr+L'xj  x  dx l  k [  = 'xk + ul  t+At  k  ' '  (2.10)  where u% is the mesh point or CRS displacement from time t to t + At.  Substituting in  Equation (2.9) and reordering,  d  dSui  dSui  Xj  ~ ~\u00A5x~  t + A t  dSui du  k  ~ dx  =  8+ x  t  t  k  At  j  dividing the above equation by At and let A i \u2014 &gt; 0, we obtain, / d5u{ \\  dSui d vl  A  t  d Xj J  dx  t  d Xj  t  t  k  which gives, dSui d Xj t+At  88 d*  Xj  dSui d Xj  ox t  l  The volume element d V tJrAt  t x  d Xj l  k  is related to d*V\\ through:  &lt;f+^V = a*v\\  t  t  + (cfV^At  (2.12)  \u000CChapter 2. ARBITRARY  LAGRANGIAN  EULERIAN  (ALE) FORMULATION  28  The computational reference time rate of volume is [48]:  {cfvy and it follows that: &lt;? V +At  =&lt;TV  *X  +  ^-AtcTV  (2.13)  O Xk t  The stress component in Equation (2.8) may be expressed by: t+At  o~ij =  (Tij  t  t  +  (^(Tijj  At  aij  =  l  ^ N  +  ^  (  J  ^ r  -  (2  14)  where, (*o-jj) = {^o-- + ('i;\u00A3 \u2014' f f c ) ! ^ ) denotes the computational reference time rate of A  a j.  t  i  Since the computation reference domain is chosen to be the material reference domain at  time t, i.e., x  \u2014 x%, the following relations may be obtained: l  t  i  dSui d'xj  dSu; i  &lt;TV  dSui  d*XA  = (TV  = (TV _t  _  X  l  Substituting Equations (2.11), (2.13), (2.14) and the above equations into the left hand side of Equation (2.8), the following linearized form of the internal virtual work may be obtained: dSui  LHS  J+Aty  d'+^xj  Jt+AtV  ' dSu;  dSui d vi  d Xj  d Xk d Xj  t  N  \u2014-\u2014=\u2022 Ai t  v  JV  +  d xi t  l  t  ^AtcfV) pp] O Xj t  +  j  AtcfV OX I t  i  \u000CChapter 2. ARBITRARY  hv  LAGRANGIAN  o Xj  EULERIAN  hv d Xj  t  \\  t  (ALE) FORMULATION  + V.,  d*x  ' '  k  l  \\  29  AtcfV  All terms with order higher than linear term of At are ignored. Considering the right hand side of Equation (2.8), we have:  7f = 7f  t+A  +* ff At  (2.16)  A  t  x  From Equation (2.1), t rBA  ti  Where ff' l  =  t rB- . (tc  ti  +{ p v  is the material time rate of ff  -  at fB  Ji p)-K-  t \u201E, \\  u  v  t  at time t. Substituting into Equation (2.16) and  l  considering *Xi =' Xi  7f ='/f+(7/- + (  t+A  ,  Similarly, the area element d Aj  t  f \u00A3 - f * ) ^ ) A t  (2.17)  t  may be described as:  t+At  d* Aj  =  +M  d*Aj  +d A &gt;At t  /  The computational reference time rate of the area is [48]:  ox  1  o Xj  t  t  k  so that: d A  =  t+At  j  c?A  +(^ld A -^d A )At t  Jt  (2.18)  i  j  k  Using Equations (2.13), (2.14), (2.17) and (2.18) and considering the relations x l  {  right hand side of Equation (2.8) may be written as: RHS  =  [  f?8u d V+  t+At  t+At  i  [  a 8u d A  t+At  t+At  ij  i  j  =  the  \u000CChapter 2. ARBITRARY  LAGRANGIAN  EULERIAN (ALE) FORMULATION  V^- + f \u00A3^d%- - ^&lt;i%) \\d'x  K  d'xj  k  I  30  M ,  where all higher order terms of At are neglected. Since d*Aj =' rijcPS and V^-TIJ =' / j ' 5  the traction (or load) rate with the displacements and gradients held fixed, corresponding to the A ^ T ? in [2], we have: RHS  =  +  J  Su ffd V t  t  i  + J  SufajfAj  t iuiCvl-'^^pi-'njAUtS  (2.19)  O Xji  JS t  Finally, using LHS Equation (2.15) and RHS Equation (2.19), and applying the relation:  hv  O Xj t  Jty  JA t  i  \u000CChapter 2. ARBITRARY  LAGRANGIAN  EULERIAN  (ALE) FORMULATION  31  We obtain the final A L E formulation in the form:  = l SulfPd V + l Sulffd S t  v  (2.20)  t  ts  The above formulation is a complete and general one in the sense that it may be applied to various types of finite deformation problems with generalized types of loadings and boundary conditions. The formulation is also suitable for implementing various types of rate dependent and rate-independent material constitutive relations. Another unique feature in the formulation, is the existence of generalized load correction terms that facilitates handling of deformationdependent loads, e.g., follower loads. These features are discussed in the following subsection.  2.3  DISCUSSION A N D CHARACTERISTICS OF PROPOSED  FORMU-  LATION The above formulation has special features and may differ from similar ones in the literature in the following main aspects: - No specific assumptions were introduced after the principle of virtual work.  Also,  velocities, instead of displacements, are used as primary variables. This makes it more straight forward to implement rate-dependent material laws. - Consistent load correction term to handle deformation-dependent loads are introduced. - The integral volume and surface have clear physical meanings which are the material domain after the last incremental step.  This is arrived at by choosing the arbitrary  reference domain coincident with the material domain after each incremental step.  \u000CChapter 2. ARBITRARY  LAGRANGIAN  EULERIAN  (ALE) FORMULATION  32  The formulation represents a strict A L E method because the mesh motion and deformation can occur independently at the same time, i.e., both material velocity and mesh velocity are included in Equation (2.20). Different objective stress rates can be applied to the formulation. If Jaumann stress rate [47],  is introduced, the formulation (2.20) takes the following particular form:  +  lv\u00A5x-[- \u00A5x~  + \u00A5x-  aik  aij  k  -I,*v hs  +  {  V  k  k  V  9x  \\  ox  ~  Vx-  Vk)  i  v  k  dx /  l  t  k  k  a'zfe  t  k  ox ) t  k  = JfV 8u\\f?-*V+ [ Snlftc^S Js l  (2.22)  t  ty  If Truesdell stress rate [49],  \u00AB  a  ~  ~ \u00A5x~  ~ \u00A5x~  akj  aik  k  \u00A5x~  +  aij  k  ( 2 k  -  2 3 )  is applied, Equation (2.20) takes the form: dSui ('t_T rp t  Iv \u00A5x~ {  r ds ( Ui  d Vi d i l  ,t_  - + \u00A5x~  a  akj  + k  t  d'v]  b^vi  ,t ( .  tv  \u00A5x~  aik  k  d'vi  = / &amp;u\\f?\u00A3V + I SulffJS J'V  J*S  -  t  dv l  t  k  \u00A5x- )  aij  k  t  *  V  .d^A  (2.24)  \u000CChapter 2.  ARBITRARY  LAGRANGIAN  EULERIAN  (ALE) FORMULATION  33  - The formulation reduces to updated Lagrangian one when the material velocity is chosen equal to mesh velocity, i.e., when V{ = vf. Equation (2.22) takes the form: f hv  dSui,(  PXJ \\ ' \\  f 6u\\f *-fV tY  t  &lt;9^fc\\  t(  2  %  hv  =  rik (d'vj  -J  t  [frxk  d^kj  +  frxj  +  J's  \\  f  cr  jfe  / &amp;VJ  d^fe\\  \\a*z  Pxi)  2  fc  d'zfc  t  d'vk \\  j.  ^d'xk)  o Xk) t  f Wfi'*S  (2.25)  ts  which is the same as the updated Lagrangian formulation given by McMeeking [3], as verified in Appendix B, except for the load correction terms:  and  J*s  \\  d xk l  o Xk) t  - On the other hand, if mesh velocity is chosen to be equal to zero i.e., v? = 0, the formulation (2.20) reduces to a general Eulerian formulation, as shown in the following equation:  - The problem of mesh distortion is easier to handle in A L E method. The A L E method can not only handle the mesh distortion but also improve the calculation accuracy. Because the motion of the mesh is arbitrary, the new mesh is always chosen to give an optimal results. - In the application of Equation (2.20), various domains in the structure may be specified as Lagrangian, Eulerian or A L E ones. A practical example for the importance of this point  \u000CChapter 2. ARBITRARY  LAGRANGIAN  EULERIAN  (ALE) FORMULATION  34  may be seen when we consider the interface between the tool and work-piece in a metal forming problem. In this case, 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 same 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 size increasement\", and load fluctuations discussed in Chapter 1, should generally not occur. This will undoubtedly increase the efficiency and accuracy of the simulation of metal forming problems.  \u000CChapter 3 NUMERICAL IMPLEMENTATION OF A L E  3.1  INTRODUCTION A N D BRIEF SURVEY  3.1.1  Mesh Motion and Stiffness Matrix Processing  In general, the stiffness matrix in A L E formulation may be rectangular and nonsymmetric. The rectangular nature of the matrix results from the fact that there are two unknowns for each degree of freedom, one related to the material and the another related to the mesh, while only one equation may be derived from the equilibrium for each degree of freedom. There are generally two ways to get a solvable set of equations in an A L E formulation. One is to specify the mesh displacements or velocities before solving the linearized equilibrium equations. In each incremental step, the mesh motion is decided by certain algorithms and data from previous increments. It is generally independent of the deformation in the current step. This method is straight forward and has been applied by many researchers [24, 39,44]. The method, however, is not a general one and may not be capable of considering the effect of the current material motion. Another method is to set up supplementary constraint equations, i.e., relations between material displacements or velocities and mesh displacements or velocities. In this method, the mesh motion is coupled with the material motion in the current step by the constraint equations. This method is a more general one and has the potential of producing higher quality mesh, but it will generally double the number of unknowns and increase the calculation time significantly. The mesh motion algorithm is an important and critical aspect in A L E formulation and it affects the computation efficiency and accuracy. Schreurs [39] introduced a mesh motion  35  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  OF ALE  36  technique for an A L E formulation that requires the solution of a simultaneous set of equations. Other researchers applied regularfiniteelement mesh generation techniques to move the computational mesh in A L E . For example, Benson [35] used relaxation stencil method derived from Laplace's equation and Ghosh [44] employed local elliptic and algebraic mesh generation methods. Generating mesh by solving differential equations may create a high quality mesh, especially when applying Poisson's equation [50] or biharmonic equation [51, 52], but it is time consuming and may introduce other types of errors in the process. Other techniques utilize algebraic interpolation and introduce a set of algebraic equations. Algebraic interpolation methods are generally simpler to apply, but they introduce curvefittingerrors in the description of regions whose boundary may not be exactly described by polynomials of the same order as those appearing in the interpolation functions [53]. In this chapter, a new method is introduced to handle the supplementary constraint equations that are produced in a mesh motion algorithm [54]. In this method, the equations are processed and incorporated in the stiffness matrix on an element rather than global level. The element equations are then reduced before assembly so that the number of unknowns on the global level are kept same as updated Lagrangian method. This reduces the need for larger computer space and C P U time significantly.  A mesh motion algorithm [55] is developed based on  transfinite mapping method [53]. This technique provides a homogeneous mesh and matches the boundary of a given region at an infinite number of points, so that no curvefittingerrors may be introduced. The new mesh location is directly determined by an explicit formula without solving any equation, which also expedites the calculation speed. Another distinct advantage of the transfinite mapping technique is that it allows the discrete representation of boundary curves and discontinuity of slope of boundary curves. This makes the method convenient and efficient since in A L E , the region boundaries would be normally expressed as discrete curves formed by nodes. Unlike the regular transfinite mapping scheme, the new procedure may take account of the motion of the boundary and may also adjust the nodes on boundaries to obtain a  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  OF ALE  37  higher quality mesh.  3.1.2  Stress Integration  In finite deformation analysis, due to the continuous change in material configuration, the constitutive equation is usually expressed as a relation between some objective rates of stress and the rate of deformation tensor. The stress integration algorithms have to keep incremental objectivity, i.e., it has to be invariant with respect to superimposed rigid body motions within a given increment. For small deformation problems, the standard time-discretization procedures may be applied to rate constitutive equations without causing too much error, which makes the integration simple and straight forward. For finite deformation analysis, however, the standard time-discretization only achieves objectivity in the limit of vanishingly small time steps [3]. This may lead to excessive error accumulation in practice when a finite time increment is used. An algorithm for integrating rate constitutive equation where Jaumann rate is used is presented in [56] and later employed in ABAQUS [57]. One of the weak-points of the method is the requirement of additional computation of rotation tensor \u00A3  + A t  R from t to t + At.  A  more general discussion about the requirements that should be satisfied by an algorithm for the integration of rate constitutive equation is given by Pinsky [58]. A theoretical derivation of a general implicit integration algorithm is also presented in [58] by pure mathematical mapping between deformed and undeformed configuration. The consistency and incremental objectivity of the integration algorithms are clearly verified. Details of the numerical implementations are not, however, discussed. Another aspect of the integration of rate constitutive equations is to satisfy the material plastic constitutive relations (e.g. Prandtl-Reuss equations) and the incremental plastic consistency, i.e., to keep the state of stress within the elastic domain or on the yield surface. In small deformation analysis, this may be achieved by a general return mapping method [59, 60].  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  OF ALE  38  For large deformation analysis, the return mapping methods should, however, be applied consistently with the integration algorithm. Numerical investigation of this point and consistent procedures for applying these algorithms are discussed in this chapter. The presented algorithms are combined consistently with the general return mapping method to satisfy the above mentioned requirements. The incremental objectivity of stress integration algorithms needs the rate of deformation tensor D j to be zero when only rigid body motion occurs [58]. It is generally achieved in t  i  implicit integration algorithms by using a central difference integration scheme, i.e., a = 0.5 [56, 58], where a is the coefficient applied in the implicit integration schemes. In an explicit integration algorithm, a numerical treatment based on forward difference method (a = 0) is necessary to keep the incremental objectivity. In this chapter, two explicit numerical integration algorithm for stress rate equation are derived based on the physical definition of stress instead of the purely mathematical tensor transformation presented in [58].  It is verified that the  developed algorithms are equivalent to integration of Truesdell stress rate equation.  One  algorithm is exactly the same as the equivalent algorithm in [58], and another one is slightly different. Practically, however, it is shown that the two algorithms give almost exactly the same results [61].  3.1.3  Updating Material Associated Properties  In A L E simulation, the material associated properties, such as strains, stresses have to be updated after each incremental step because the mesh motion is independent of the material motion. Such updating is necessary in Eulerian method and updated Lagrangian method if re-meshing is applied. One method of achieving this in Eulerian formulation is to create an imaginary finite element mesh through the material points at time t + At and use this mesh to calculate properties at integration or nodal points by interpolation [7]. This process is very  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  OF ALE  39  time consuming and requires different interpolation logic every step and certain simplifications may be applied to the procedure [21, 22]. In updated Lagrangian formulation with remeshing capabilities, the updating of material associated properties is done through three steps [16]. First, properties at the old mesh nodal points are obtained from the old mesh integration points by extrapolation. This step is not consistent when integration points undergo material nonlinearity. A weighted nodal averaging may then be obtained. Normally, a direct relation between old and new meshes is not known. Therefore, in the second step, the old mesh element to which a new mesh node belongs is found. This is done by checking iteratively for the normalized local coordinates of the new mesh node in old mesh elements. If the absolute magnitudes of all the normalized coordinates are less than or equal to 1, then the old mesh element in which a new mesh node belongs may be deciphered. After that, the values of variables at the new mesh nodes may be defined by simple interpolation within the old mesh element. In the third step, the values at the integration points of the new mesh elements are determined by interpolation from the new nodal values. In A L E formulation, researchers use an updating algorithm called \"follower point \" or \"pseudo-material element\" [39, 43]. In this method, all integration point variables are extrapolated to the nodal points of the \"pseudo material elements\" and a weighted averaging may be used. After that, iterations are used to find the \"pseudo material element\" to which mesh points will be related. The new nodal and integration point values are then calculated by interpolation within the \"pseudo material element\". The main disadvantage of this method is the use of one interpolation plus one extrapolation schemes. This may greatly degrade the accuracy of simulation and the calculated parameters. This procedures is basically similar to the one used after re-meshing in updated Lagrangian method. Interpolation is necessary for various methods of updating. It is usually done by interpolation functions similar to those used for displacement or velocity in the finite element discretization [62]. A slightly different scheme is to create a continuous field for all variables by a method  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  OF ALE  40  called \"local smoothing with triangulation\" [16], then, interpolate the variables linearly within one triangle. A different procedure for interpolation is developed by Derbalian [22] using conjugate approximation method [63].  By constructing a set of conjugate shape functions  (bi-orthogonal) to those used to represent the displacement or velocity field, a consistent approximation for the stress field may be obtained which is continuous across inter-element boundaries and which involves less mean error. Because this method introduces a large system of equations, it is usually applicable to Eulerian formulation since the large system of equation need to be computed only once. For other cases, the interpolated elements are different in each step and the large system of equations must be solved iteratively during each incremental step. In this research, an updating scheme based on continuum mechanics equations is used to update material associated properties. The relation between mesh motion and material motion is employed to avoid iterations required in previous schemes. The procedure is described in Appendix A. The method is more accurate, simpler to implement and eliminates the necessity of interpolations and extrapolations mentioned above.  3.2  STIFFNESS M A T R I X PROCESSING Using the standard finite element procedure to discretize the A L E Equation (2.20), we get, 4  K ' v +* K  C  V = V  (3.1)  where, *v is the vector of material velocities at nodal points, *v is the vector of mesh velocities c  at nodal points, ' K is the stiffness matrix related to *v, ' K i s the stiffness matrix related to *v , c  c  'p is the external loading rate vector. In Equation (3.1), the total number of unknowns is twice the total number of D.O.F. and the equation number is only half of the number of unknowns. In order to solve the equations by conventional methods, supplementary equations equal to number of total D.O.F. are needed. These equations are supplied by the relations between v 4  and 'v , i.e., an explicit mesh motion scheme in A L E . c  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  OF ALE  41  The basic idea is to avoid setting up the supplementary constraint equations to solve Equation (3.1) on a global level. Instead, we set up the constraint on the element level and reduce the equations by condensing out the mesh motion variables before assembly. This will reduce the global D.O.E to be only the material velocities. On element level, Equation (3.1) may be written as: %/vj  + %/vj  =%  {I = 1,2,  where, n is the number of D.O.E in the element.  n)  (3.2)  We use a transfinite mapping for mesh  velocities, *Vj in the following general form: *v j = aj + c  6 J)*T;(J) (  (J  = 1,2,  n)  (3.3)  where, no summation on \"J\" is observed. It should be noted that the mesh motion scheme of Equation (3.3) guarantees one-to-one mapping between material domain and computational reference domain, as verified in Appendix C and that the choice of the coefficients aj and bj will identify the type of formulation as follow: If aj = bj = 0, the formulation reduces to Eulerian one, If aj \u2014 0 and bj = 1, the formulation reduces to Lagrangian, and Otherwise,  the formulation will be a general A L E one.  The above three cases may be mixed in a given problem so that certain regions of the model may have different formulation. It is also important to realize that at each node, different D.O.Es are allowed to have different formulation type. This property makes handling contact boundary condition much easier. The only apparent constraint here is that the boundary points should be always kept as Lagrangian type in the direction normal to the boundary so that the mesh and material points will always represent the same boundary. Substituting Equation (3.3) into Equation (3.2), we get: (%J  + 'kpj^fvj  = % - %jaj  (3.4)  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  where there is no summation on (J) in k ^b^ y t  42  In matrix form, the equation may be written  c  I  OF ALE  J  as: 'K'v =  (3.5)  where *K is equivalent stiffness matrix and ' p is equivalent load rate vector. In Equation (3.4), 7  the only unknowns are the material velocities vj. i  Thus, by introducing the supplementary  constraint equations on the element level and modifying the elemental stiffness matrices, the mesh velocities may be condensed out of the element equilibrium equations, so that the number of unknowns will be finally the same as in the traditional finite element formulation. In Equation (3.4), the equivalent stiffness matrix is ku l  load is ' / / \u2014 tk^jdj.  +  and equivalent element  The conventional finite element assembly and elimination method may  now be applied directly to solve for material velocities vj. This procedure is more efficient and t  the usual finite element solution routines may be used with minor changes. The only limitation to the procedure is that the mesh velocities may be coupled with the material velocities only at the same D . O . E . In practice, such limitation is trivial, because more complicated relations between material and mesh velocities may not guarantee a better mesh quality, but will generally increase the computation time significantly. If for some particular reason, the mesh velocity *vj has to be coupled with the D.O.E in other elements, the above method is still applicable. Only difference is that processing would create \"fictitious\" nodes whose D . O . E is coupled by the mesh velocity equations and which are not connected with the element physically.  3.3  M E S H M O T I O N S C H E M E IN A L E  3.3.1  Transfinite M a p p i n g M e t h o d  To complete the mesh motion procedure and finalize the supplementary constraint equations, it is essential to decide the values of the coefficient aj and 6/ in Equation (3.3). As discussed  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  OF ALE  43  above, on boundaries of a deformation domain, the nodes have to be Lagrangian type in the normal direction in order to keep a one-to-one mapping between the material and computational reference domain, i.e., aj = 0,6/ = 1. These boundary points may move, however, in the tangential direction with a velocity different from the mesh velocity. All the internal nodes are assumed as general A L E points moving at speed vj = a/, which is specified by transfinite l  mapping method, and the coefficients hi are set to zero for these general points. The transfinite mapping method is employed to decide the internal nodal speed [55]. Originally, transfinite mapping is applied to create a mesh on a geometric domain when the boundaries are specified [53].  Assuming a region with four boundary curves specified as  4&gt;i(r, 0), (f&gt;i(r, 1), c&amp;(0, s) and &amp; (1, s) as shown schematically in Figure 3.1, the mesh coordinate Cj is given by:  d(r,s)  =  (l-s)&lt;f&gt; (r,0)  _  (1 _ ) ( i _ )^(0,0) - (1 - r)s&lt;f&gt;i{0,1) - rs&lt;j&gt;i{l, 1)  -  r(l-*)&amp;(l,0)  i  r  + s&lt;f&gt; (r,l) +  (l-r)(f) (0,s)+r&lt;f&gt; (l,s)  i  i  i  a  (3.6)  0 &lt; 7\" &lt; 1  0 &lt; 5 &lt; 1  where r, s are normalized coordinates, i = x,y, i.e., (j&gt; (r, 0) and &lt;j&gt; (r, 0) etc., represent the x  x,y coordinates of boundary curves, and c (r, s), c (r,s) x  y  y  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  \u000CChapter 3. NUMERICAL  (0,D  W  IMPLEMENTATION  OF ALE  44  (0,1)  r,1)  (1.1)  &lt;t&gt;i(0,s)  4&gt;,(i.s)  (U\u00B0#) )  (0,0)  ^(r.O)  *,(r,0)  (0,0)  (1,0)  (0,1)  (^i(r.O)) Modified Region  Original Region  Figure 3.1: Transfinite mapping representations in the transfinite method is to have normalized parametric coordinates assigned to each node on boundary curves. This process is automated by assigning the coordinate 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\u00B0(r, s), then the mesh velocity is given by: Cj(r,s) - c?(r,a)  Ai  t c*  (3.7) (3.8)  \u000CChapter 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 \u00A3;(r, 0),\u00A3i(r, l ) , \u00A3 ; ( 0 , s ) and  s) corresponding to boundaries  \u000CChapter 3. NUMERICAL  &lt;/&gt;i(r, 0),4&gt;i(r, l),4&gt;i(0,s)  IMPLEMENTATION  OF ALE  46  and &lt;/&gt;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) - \u2122&amp;(1, 1)  -  r(l-a)&amp;(l,0)  (3.9)  0&lt; r&lt; 1  0&lt; s &lt;1  To ensure that the boundary nodes have to be Lagrangian type in the normal direction, the mesh velocities \u00A3;(r, 0), \u00A3;(r, 1), \u00A3;(0, s) and &amp;(1, s) have to be equal to the material velocities Vi(r, 0)/ Vi(r, l), 1&gt;;(0, s), and *VJ(1, s); or at least their normal components have to be equal  t  4  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 \u2014 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) = - v (r,0) t  At  i  6(0, s) = ' - M 0 , *) A t  (i = x,y  0&lt; r &lt; 1  &amp;(r,l) = \u00AB- 't; (r,l) A  i  6(1,\u00AB) = * \" f i ( l , *) At  (3-10)  0 &lt; s &lt; 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\u2014At) 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  \u000CChapter 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  i = i  v  v  c  t  i  t  + i v  c  Co  b  = i a  1 1  \\  (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  \u000CChapter 3. NUMERICAL  3.3.3  IMPLEMENTATION  OF ALE  48  T h e M o t i o n of Nodes on Boundaries  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&gt;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 \u2014 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  \\  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  OF ALE  49  each node.  3.3.4  Implementation and Discussion  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 \u2014 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 assembled 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 advantages. 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  \u000CChapter 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 OF STRESS  3.4.1  RATE  T h e Integration A l g o r i t h m  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 V j j . At time (i + Ai), the configuration is '+ 'C and the Cauchy stress is o-ij. A  t+At  As shown in  Figure 3.4, the 2nd Piola-Kirchhoff stress at (i + A i ) with respect to the configuration at time  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  OF ALE  51  t, t ^ i j , is related to Cauchy stress l+^crij, at time t + At by [47]: d Xi t+At  d z t+At  *t+At, t *\u2022  dx t  1  r  where l F +At  = g^ at  i:j  &lt;9'a;  (3.12)  n  is the deformation gradient tensor or the Jacoban matrix.  Xi  Q(  P('xi)  Configuration at timet C  T+AT  X|)  Configuration at time t+At C  (  t + A t  *3 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: AtJD  t+  _  1  d Xi Ato(D t+At  t+  r  and since this step is only a rigid body motion, then, S  =  t  s  %  t+A  l  r  t+At (I)  d  dx  dx t  =  t  (3.13)  \u000CChapter 3. NUMERICAL  where  o-\\p  t+At  and  l S^ +At  n  IMPLEMENTATION  OF ALE  52  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: ij  t+At&lt;T  =  t+At  o~lp + Ao~ij  i.e.,  + Ao-ij  (3.15)  or, H-A*  .%3  d  t + A t  x  d  t+At  i t  Xj  + djkfDuAt  (3.16)  where, Cijki is the element constitutive tensor, At is the time increment, Wki is the rate of deformation tensor,  - |  +J^).  In decomposion of the deformation, if the first step and the second step are switched, i.e., if the first step is only straining and the second one is only rigid body motion, we get another integration scheme as shown in Appendix D.  \u000CChapter 3. NUMERICAL  3.4.2  IMPLEMENTATION  OF ALE  53  Consistency of the Integration Algorithm  Integration algorithms have to be consistent with the constitutive equations. It may be shown that the integration algorithm developed in Equation (3.16) is equivalent to using Truesdell stress rate in the rate constitutive equation. We consider the determinant of the deformation gradient tensor as: t+At p  \\F  \\F = \\F + \\F  + d  du  k  dx l  k  where, |f,F| = 1 and u =  x  \u2014 x is the displacements from time t to t + At, so that;  t+At  k  l  k  k  1 t+Atp  du  k  d u d xi  d  t  c  x  t  dx  du\u201E  m  d  m  t  l  k  x \u201E  with the higher orders of Ui being ignored. Substituting the above equation into Equation (3.16), gives: du  \\ d { xi + Ui)  UJ)  d ('xj +  l  k  o-ij = 1  t  +  Ciju'DuAt  Rearranging the above equation and neglecting all the higher order terms of uf, we obtain: dx  d Xj  t  *+At\u201E..  _  dui  t  i  t  0~rr,n,~. L  x  m  d Xi  t  k  dx dx t  t  n  rr  t  n  *  0x  dUj_  dx l  r  n  t  mn  \u2014  l  0x  m  d Xj du  t  d Xi  t  0~mn ~7^1 ~f~  ft x  ox  l  d xj  t  \"777  A \u2022  +  Cij i D iAt t  k  k  k  dui  4 -  X . \\ X. t  t  mn^nj  O X, t  I i m ^ranr\\+ u  ox t  n  r  t  duk_  r  dx l  k  or, &lt;r  t+At  - - n-t  At  _  I  \u2014 T T  At  d u*  7^ Vd&lt;z  t_  I  (?mi + m j  m  t_ I  t  du  ^ \"ftt,  , /-\u00BB  o  du  P. d x, m  \u2014 cr.  N k  +  CijuDu  t  which, upon taking the limit when At \u2014 &gt; 0, gives: i  dx l  ^ \" j t \u201E  dx l  k  k  t  n  ki  (3.17)  \u000CChapter 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:  (3.18) 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 A material point in both configuration may be given by the position vectors and Q( xi, x , x ), t+At  t+At  C.  P( xi, X2, x ) t  t  t  3  as shown in Figure 3.4, respectively. Tensor quantities in both  t+At  2  t + A t  3  configurations C and *+ (7 may be related to each other by proper transformation. Starting l  At  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: 1  *  + A t  d  t+At  d Xi t+At  F  dx  m  At  Xj  \u00B0  t  d Xi t+At  +  d  t+At  mnkl t+aAt t+aAt  where, t+aAt  = a Xi t+At  X i  + (1 -  afxi  d  Xj  t + a A t  X. n  (3.19)  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  t+At t+aAt\"  d  t+A  F1  55  X{  Qt+cAt  \u00BBJ 7.  OF ALE  Xj  t n_I(  dtyi  ,  when a = 0, it will be explicit and may be expressed as: t + A t  _  i  &amp;+*x  d^ %  i  A  it  a  t + A t  ^ _  Acr  t+Atp  where, A c r  m n  a  A  m n  t + A t  xj  ^-\u2014^  (3.20)  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  where, ^^R-  ^- = r  A  t  ^ &gt; \u2122 r  A  t  % + .  is the rotation tensor for the increment from time t to t + A i ;  ACTJJ is the stress increment given by: A&lt;T;J = Cijkil+ox^D At kl  The l+o.s&amp;tDki is calculated by central difference method: *  n  t+O.BAt^W  _ l _ ( _ d ^ Ut+O.SAt^ 2  -  d%_\\  +  at+0.5At J x  (3.2i)  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  OF ALE  56  and t+0.5At  If displacement  Ui =  X{  t+At  \u2014 Xi l  = 0.b x t+At  X i  {  + 0.5*3*  and the gradient of displacement with respect to  t + 0 5 t  Xi  is:  dui Qt+0.5At  R can be expressed in matrix form as [56]:  where, I is unit matrix and G has components of Gij. Upon comparing Equation (3.21) and Equation (3.15), it is important to note that in Equation (3.21) the stress at time t is transformed by a rotation tensor l R j +At  {  tion (3.15), it is transformed by deformation gradient tensor, ^^F-.  whereas in Equa-  For a pure rigid body  motion, the Truesdell rate is the same as Jaumann rate, the deformation gradient tensor reduces to the rotation tensor and determinant of deformation gradient tensor is equal to unity, so that Equation (3.15) and Equation (3.21) lead to the same results. The objectivity of the algorithms have already been verified in [56] and [58].  3.5  I N C R E M E N T A L OBJECTIVITY A N D NUMERICAL T R E A T M E N T OF *Dy  3.5.1  General Considerations  The rate of deformation tensor is not only used in integration of rate type constitutive equations as in Equations (3.15, 3.20 and 3.21), but it is also used in other applications, such as visco-elasticity and creep involving large strains and rotations. The incremental objectivity of Dij l  means that under finite incremental step, its components will be independent of any  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  OF ALE  57  rigid body rotation. For large deformation analysis, this is shown to be achieved when the integration scheme is used with a value a = 0.5 ( central difference method) in implicit stress integration algorithms [58]. In many applications, however, Wij is employed separately, not with stress integration, or is used with an explicit stress integration, e.g., Equation (3.16). In this section it is shown that using the conventional forward difference method in such analyses may only achieve objectivity when the time steps are very small and it may lead to excessive error accumulation in practice. A forward difference method with a correction term is presented to keep the incremental objectivity in situations where rigid body rotation is presented.  It may be instructive to show first that analytically, the components of the rate of deformation tensor do identically vanish in the case of rigid body rotation. This may be  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  OF ALE  58  achieved by considering a general two-dimensional rigid body rotation with an angular velocity en. If the rotation center O is chosen as the origin of a fixed Cartesian coordinate system, a material particle located at P ( x , x ) t  t  at time t moves to anew position P  t  1  2  t + A t  (  t + A t  xi,  t + A t  x2)  at time t + At, as shown in Figure 3.5. The new position may be expressed in the following relations:  t + A t  Xi  x  t+At  2  = txiCos^At)  \u2014  = txisin^At)  + x cos(u;At)  x sin(uiAt)  t  2  t  (3.23)  2  so that, the displacements are given by: ui \u2014 *a;i(cos(a;Ai) \u2014 1) \u2014 x sin(u&gt;At) t  2  u  2  + tx^cos^At)  = txxsin^At)  -  1)  (3.24)  Taking derivatives with respect to At and calculating the limit values when At \u2014\u00A5 0, gives the velocities at time t as: t  t  t  V\\ = \u2014UJ X  t  V2 = UJ X\\  2  From the definition of Dij, the components of the rate of deformation tensor *D at time t is l  given by:  '\"\u00BB = | \u00A3 = \u00B0 \u2022  D  a  i  (p. + \u00A3 \u2022 ) o  &lt;- &gt; 3 25  =  2 \\a x t  2  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.  \u000CChapter 3. NUMERICAL  3.5.2  IMPLEMENTATION  OF ALE  59  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 D j. t  i  The conventional method infiniteelement analysis is to use the average velocity from time t to time t + At as velocity at time t, and employ this to calculate D j as the rate of deformation t  i  tensor at time t. The calculated value of D{j is also used as the average value of the tensor t  from time t to time t + At.  Mathematically, this is equivalent to a forward difference method.  Thus, the velocities are given by; Xi{cos{ujAt) \u2014 1) \u2014 x sin(u&gt;At)  u\\  t V  l  At u  =  t  t  2  *  =  At  Al  =  x sin(u}At) + x (cos(u&gt; At) \u2014 1)  t  2  2  t  t  1  2  Al  =  ( 3  the Dij components are then obtained by taking derivatives with respect to x t  t  and x t  i  2  -  2 6 )  as  follow; .  r&gt;w \u2014 -  \u00B0  which shows that D u and D f  t  22  At  cosjwAt) - 1  t  4  cos(uAt)-l  2 2  Z)  12  ~ =  Al 0  (3.27)  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 Equation (3.27), it may be seen that the modification term should be a function of the rigid body  \u000CChapter 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 Dij components for general 2D large deformation process may l  be considered as: *Ji - (WuAty Ai  -1  The above modification guarantees the objectivity of  Dij in the case of rigid body rotation.  l  With the modification term, numerical calculation of the Dij components are given by: l  where,  Vi - (w  12  Ai  vMHf^-S)^) 2  -1  )2  Ai  Ai  1  (3.30)  and 8ij is the Kronecker delta. For the general rigid body rotation case given above, the spin tensor Wij is calculated in a l  way similar to Equation (3.27) and the components are given by: *Wn  = 0  W  =0  l  22  Substituting the above values in Equation (3.30) and (3.29), we have:  'Dn = 0 &lt;D = 0 22  4  /J&gt;i2 = 0  \u000CChapter 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  STRESS T R A N S F O R M A T I O N  3.6.1  AND RETURN  MAPPING  G e n e r a l 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 befirsttransformed 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 deformation 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  \u000CChapter 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:  = where, a' - = crij \u2014 | \u00A3 ; / c r t  t  i  fcfc  (3-31)  is the deviatoric stress and Sij is Kronecker delta. The accumulated  equivalent plastic strain is calculated by [64]:  **=Cif ^ ^  \u00A3  0  0  -  dr  (3 32)  where; t is the initial time, t is the current time of the deformation and D^is the plastic part T  0  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 '&lt;jjj and material velocity is Vi, with other quantities as t  shown in Figure 3.4, the procedure may be summarized as follow: (i) input data for integration: o~ j, v , At, X{, t  t  i  (ii) update material coordinate:  t+At  t  i  ai; = X{ + viAt; then, calculate Cijki by assuming the  increment to be completely elastic;  t  t  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  OF ALE  63  (iii) calculate the elastic stress increments A&lt;Tij by the use of Equation (3.14); (iv) use stress quantities 'cr^ and AUJJ as input to the return mapping algorithm and calculate an intermediate stress state ' c r j p that satisfies the incremental plastic consistency; then, calculate deformation history-dependent parameters, such as the current yield surface  cr  y i e W  and plastic  strain e ; p  eq  (v) calculate the stress increments by Acr^- =  cr\\f  - aij, use Equations (3.15), (3.20)  t+At  l  or (3.21) to update stresses and obtain a second intermediate stress state corresponding equivalent stress  t + A t  1  _  7  to the yield surface by using the relation  t+At-( ) H  Vyield u  cr|J ^ and its  a ^ by Equation (3.31);  (vi) scale down the equivalent stress point VJJ ' t+At_  t+At  eq  In the procedure above, it may be seen that step (vi) is necessary to keep the final stress point  t + A t  3.6.3  aij  plastically consistent, i.e., on the yield surface.  Scheme B: Transformation Before Return Mapping  The second procedure may be summarized as follow: (i) input data for stress integration: V^-, Vi, A i , Xi; l  l  (ii) update material coordinates: a;j = Xi + 't&gt;;Ai; then, calculate Cijki by assuming the t+At  l  increment to be completely elastic; (iii) calculate the intermediate stress state a\\p by the use of Equations (3.15), (3.20) or (3.21); l  (iv) calculate the \"elastic\" stress increments by Acr^- = (v) use  'cr^-  cr[p \u2014 V^-;  t+At  and Aaij as the input stress components to regular return mapping algorithm and  obtain the mapped stress  t + A t  cr;j;  (vi) calculate the deformation history-dependent parameters, such as current yield stress &lt;7 j y  e/d  and plastic strain e . p  eq  Although this method projects the final stress point on the yield surface, the \"elastic\" stress  \u000CChapter 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 aij instead of t  (rlp , because all history dependent parameters such as the current yield stress and the plastic  t  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 mesh motion A f x  x  / and the change due to  may be related to each other by (see Appendix A for details): gt+At f  xf  A  where f t+At  = /(  increment A  x  t+Ai  = A,/ +  --\u00AB0  (3.34)  x , t + At) is the value of / at material point at time t + At. The material  / 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:  7=7+ A  t+A  x  /  (3.35)  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  Similarly, when mesh point \\ at time l  OF ALE  t moves to x  t+At  t+At  fc  =  c  tf  +  A  x  at time  t+  65  A i , at  f  t + A t  x , we have: (  3  &gt;  3  6  )  But at time i , *x = *x, so that * / = * / = / ( x , i ) . Therefore, Equations (3.34)-(3.36) may c  4  have the form:  \u00BB r u  = ' f + A , / + ^ K - \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 conventional incrementalfiniteelement 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.  \u000CChapter 3. NUMERICAL  IMPLEMENTATION  OF ALE  66  In this method, assuming the smoothed property / \u201E and the smoothing shape function h  a  at nodal point a, then: 7 = hj  (3.38)  a  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&gt; = J J  (h f  a  a  a  ~ f)  (3.39)  dxdy  whereA is the area of the element and the following condition should be satisfied: e  d&lt;f&gt; = 0,  (3  1,2, ..n  (3.40)  n  where, n is the number of nodes per element. So the smoothed properties may be obtained n  from the expression: j jf (h h dxdy) a  J  0  a  = 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 f  can be obtained, if the smoothing  a  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 h is a bilinear shape function and a 2 x 2 a  Gaussian integration rule is adopted in evaluating Equation (3.41), it can be shown that the smoothed nodal properties fi,f ,fz, 2  7i ' 7 7 2  3  fi may be expressed as: 1 2  ' 2  1 2  1+  f 1 2  1 _ }/3 2  1 2  1 _  1 _  i/I  _ 1 2 ' 2  ^ 2  1 2  2  1 2  1  //  V3 2  hi &lt;  1 2  fin 2  .  (3.42)  \u000CChapter 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: &lt;/&gt; = E ( ( W . ) \\K - f \\ )  y  (3-43)  2dxd  K  K=l  where N is the number of sampling points within the element. s  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 h , the following expression may be obtained, a  ' fx ' (3.44)  &lt;  (h^hi)  [ ^K=I  \\K  \u2022\u2022\u2022 ^K=I (^4^4) \\K 1 ( J 4 )  In the above equation, hi \\ ,...,h4 K  \\  K  K ^K=I  yu j \\K  are values of shape functions at Gaussian points  I, II, III and IV, f \\ (K = I, II,..., IV) is the unsmoothed properties at Gaussian points. K  Solving for f , f , f and / , the above equation reduces to Equation (3.42). In this case, it x  2  3  4  should be noted that f = h f a  a  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  f,  t+At  Ui,  u\\ and  t + A t  Xi  = X{ t  + 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 f t+At  =  f;  t+At  - Take the values of *+ '/ at different integration points I, II, III and IV, as the values A  fi, fn, fm and fiv, in Equation (3.42), and get smoothed variables f  lf  for nodes 1,2,3 and 4 individually.  f , f and / 2  3  4  \u000CChapter 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.  \u000CChapter 4 TWO-DIMENSIONAL A L E P R O G R A M (ALEFE)  4.1  PROGRAM ALEFE  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, *v = aj +  (1 = 1 , . . . , JV)  7  where, no summation on I is observed in the above equation. v t  (4.1) and v\\ are material t  I  and mesh velocities at D.O.F. I respectively, and N is the total number of D.O.F. in the model. If ai \u2014 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  \u000CChapter 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.  \u000CChapter 4.  4.1.2  TWO-DIMENSIONAL  ALE PROGRAM  (ALEFE)  71  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 &lt;7;J, x ; T  t  i  2. Assemble the matrices and solve for V{ and *vf; t  3. Calculate V{At, t  and use  v?At  t  o-ij, ViAt  T  t  to calculate  Aa^;  4. Use (Tij and Aa^ as input to stress integration and return mapping algorithm to 1  get the new stress  after the current iteration;  o-ij  T  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  &lt;Jij,  T  x;  t  i  2. Assemble the matrices and solve for V{ and v\\; t  t  3. Calculate accumulated displacements in the current load increment step tt\u00BB = u{ + V{At  t  4. Use  and u\\ \u2014 u\\ +  &lt;Tij  t  v At  t  c  i  and use 'er^, U{ to calculate  Aa^;  and Aaij as input to stress integration and return mapping algorithm to get  the new stress  &lt;Tij  T  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, end the iteration loop.  \u000Chapter 4.  TWO-DIMENSIONAL  ALE PROGRAM  (  Start  (ALEFE)  )  Input d a t a , initialize t, tag =0, Uj=0,  Uj=0,  a,,  bj, t  x  j  A p p l y l o a d increment, l e t T g ^ t a j j  T C r e a t stiffness matrices f r o m % j ] x Q. O O  _l  P r o c e s s the m a t r i c e s , a s s e m b l y a n d s o l v e v f a n d v j  c o  Uj=Uj+VjAt, U;=  to  uf+\\ZfAt, u s e t a j j , u * to calculate A o n  T  E .a  U s e integration algorithm a n d return m a p p i n g to get CL O o \u2014I  o.  the n e w s t r e s s Toy after the current iteration  cr LU  I  C a l c u l a t e the residual f o r c e , internal f o r c e , c h e c k equilibrium  55  J J p d a t e L a g r a n g i a n m e t h o d Uj=_uj_2_ \"|~~Yes~ U p d a t e r e f e r e n c e configuration a n d s t r e s s  U p d a t e material a s s o c i a t e d properties to m e s h points  Yes O u t p u t results  find) D e c i d e hte n e w m e s h s c h e m e a , a n d bj  A p p l y the n e w residual a n d l o a d incremental a s n e w load i n c r e m e n t a l , let ui=0 a n d  ^=0  Figure 4.1: Flow chart of the main program  \u000CChapter 4.  TWO-DIMENSIONAL  ALE PROGRAM  (ALEFE)  73  The difference between the two flow schemes is mainly in steps 3 and 4. In first scheme, the displacements of the current iteration 'v;A\u00A3 and the updated stress variables to calculate the stress increment  A&lt;T;J,  &lt;Tij  T  are taken as basic  so that ACT;.,- is only the increment in current  iteration and is caused by * ^ A i . The return mapping and stress integration starts from the stress state after last iteration,  &lt;j{j.  T  This state may not be generally a converged stress state and  may not reflect a proper deformation history. In the second scheme, however, the accumulated displacement, U{ =  + 4t&gt;;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&lt;T;J  is the \"accumulated\" increment from the  beginning of the current step. The return mapping and stress integration starts from the stress state after the last incremental step  'cr^-.  This stress state is always a converged state and reflects  a proper deformation history. Another disadvantage of the first method is that the current yield surface or yield stress may be overestimated, which may lead to wrong results. To explain this, we consider an iteration load that takes the current yield surface at a point to a higher level. If this iteration load is not right, but is then corrected during the increment, the yield surface at the point will still be registered at the higher value obtained from the iteration. This may, in general, not reflect a true deformation history. In the second scheme, the integration and return mapping are always performed from the last converged or accepted stress state. Thus the final converged results are not affected by errors or overshoots that may occur in a specific iteration. Numerical experiments show that the first scheme is less robust from the convergence point of view and gives less accurate results upon comparing to experiments or published data than the second scheme. Consequently, the second scheme is employed in program A L E F E .  \u000CChapter 4.  4.2  TWO-DIMENSIONAL  ALE PROGRAM  (ALEFE)  74  M A I N SUBROUTINES  4.2.1  Functions of Subroutines  A total of twenty-three subroutines are compiled with the main program to perform the A L E simulation. The function of each subroutine is described as:  M A I N Input analysis data such as total number of incremental steps, tolerances, maximum iteration number; specify the boundary nodes as either Lagrangian or general A L E nodes. Output results according to user's requirement. IND ATA Input initial data from the text file created by other commercialfiniteelement codes, e.g., NISA. The data include finite element model, material properties and boundary conditions.  INCRL Evaluates the incremental load vector according to the total load, total number of incremental steps and residual forces from previous iteration. G A U S S Q Sets up the integration point position in natural coordinate system and calculates weighting factors for numerical integration.  SHADRI Calculates shape function values and their derivatives with respect to natural coordinates at a particular point.  J A C O B D Calculates Cartesian coordinates of a Gaussian points, Jacobian determinant and the Cartesian derivatives of shape functions at Gaussian points. DORDEP  Evaluates the constitutive matrix Cijki for plane strain problems.  STIFMA Evaluates A L E stiffness matrices for each element.  F R O N T Frontal solver for unsymmetric stiffness matrix.  \u000CChapter 4.  TWO-DIMENSIONAL  ALE PROGRAM  (ALEFE)  75  S R S 1 N Evaluates stresses, plastic strains and current yield stress at Gaussian points and calculates residual forces at nodal points. E L E P L 2 Calculates the stress increment due to elastic deformation and the part due to plastic deformation. I N V 2 Evaluates stress invariants and the equivalent stresses from given stress components. P O P V 2 Uses return mapping method to update the stress increment due to plastic deformation. S T R N 2 Calculates the rate of deformation tensor. S T R U P D Updates the material associated properties from material points to mesh points. A L E D E C Assigns mesh motion for boundary nodes according to the user's specification. M E S M O T Calculates mesh motion for all interior nodes. B O U A D H Adjusts boundary nodes after each incremental step to keep homogeneous distribution of boundary nodes. C O N V E R Calculates convergence parameters and checks convergence. D E C O N F Stores the output data in a format compatible with other commercial graphic processing software for post processing. W R I D A Write data files for further processing. RE ADA  Restores data from direct access files.  CHADIM  Converts stress, rate of deformation from one-dimension array to two dimension  array. MMATX  Calculates the product of two matrices.  \u000CChapter 4.  4.2.2  TWO-DIMENSIONAL  ALE PROGRAM  (ALEFE)  76  Flow C h a r t of M e s h M o t i o n  In A L E F E , mesh motion is completed by three subroutines, A L E D E C , M E S M O T and B O U A D H which, individually, calculates the aj and 6/ values for boundary and interior nodes and adjusts boundary nodes after each incremental step. User only needs to specify each D.O.F. on boundaries as Lagrangian, Eulerian or general A L E . Users may specify particular areas as Lagrangian, Eulerian or A L E ones. Such specifications will be maintained during the whole simulation. Figure 4.2 shows a flow chart for the mesh motion scheme.  4.3  CONVERGENCE  4.3.1  CRITERIA  G e n e r a l Discussion  In an incremental solution strategy based on iterative methods, realistic convergence criteria are necessary for the termination of the iteration. If the convergence tolerance is too loose, inaccurate results may be obtained, whereas if the tolerance is too tight, much computational effort is spent to obtain needless accuracy [65]. ANS YS [9] uses displacements and residual forces as convergence criteria. The maximum value, absolute value and Euclidean norm options are included. In residual force check, the total applied external force instead of the increments of external forces are employed as the base. In NISA [10], displacements, residual forces and iterative energy are checked in each iteration. In the displacement criterion, a load step will be assumed converged if the ratio of the maximum absolute iterative displacement to the maximum absolute displacement at first iteration in the current step is less than or equal to the preset tolerance. As in ANSYS, Euclidean norm of residual forces is checked, but the base is the incremental external forces. In the energy criterion, a load step will be assumed converged if the ratio of the iterative energy to the energy at first iteration is less than or equal to the preset tolerance. In [66] residual force criterion similar to ANSYS and energy criterion (ratio of the sum of absolute values of work done by the  \u000Cer 4.  TWO-DIMENSIONAL  ALE PROGRAM  (ALEFE)  (Start)  _r  Input a i , b i for all boundary nodes and some particular areas if needed Increment step n=1, t=At u, =0, uf =0 Load vector Iteration k=1 Calculate stiffness matrices and solve for v i and reaction force rate 1  I  u ^ u i + 'viAt,  T  =uf+arr 'Vji) b(i)  New load vector  Use u,to update stress, strain, yield surface and calculate residual forces .,onvergec[_ Yes  No  k=k+1  t  For all nodes other than specified by users, set ai =bi =0  Adjust nodes on boundaries by subroutine BOUADH and update u f  +  :  Use U| and u | to update material associated properties to mesh points  X  t=t+At, t = t , + u f X|  x  \u2014f  Yes  According to users's specification, decide the values ai and bi for all boundary nodes, using subroutine ALEDEC Use 'xi and ui /At as input to calculate ai and b| for all interior points, call subroutine MESMOT v n=n+1  Figure 4.2: The flow chart for mesh control  \u000CChapter 4. TWO-DIMENSIONAL  ALE PROGRAM  (ALEFE)  78  iterative residual forces to the sum of absolute values of work done by the accumulated external forces in the first iteration) are used. In program A L E F E , displacements, residual forces and energy criteria are employed for checking convergence. Brief description of these criteria is described in the following.  4.3.2  Displacement Criterion  In this criterion, it is assumed that the ratio of the norm of the current iterative displacement to the norm of the displacement in first iteration be within a preset tolerance. This may be written as: (4.2)  where N is total number of D.O.F. in the model, Au\\ ' = v At  is the displacement in k  t  th  I  iteration, vj is the material velocity in the k t  th  iteration at global D.O.F., I, t is displacement d  tolerance. Although this criteria is effective in some analyses, numerical experience showed that other criteria should be, in general, incorporated with this one [65].  4.3.3  Force Criterion  Another criterion is based on the out of balance forces. Substituting Equation (2.7) and the shape function approximation into Equation (2.6), we get the equivalent nodal force at D.O.F. T. t+At  - o-ij t+At  \u2022+Aty  (huj + h ) d*  V  +  jIti  2  and the external force at D.O.F. I: t+At  Pi = /  Jt+Aty  ffh d V  t+At  t+At  iI  + [  Jt+Atg  fihu&lt;f S  t+At  +At  \u000CChapter 4. TWO-DIMENSIONAL  ALE PROGRAM  (ALEFE)  79  Therefore, the residual force at D.O.F. / can be expressed as: t+At  n =  t+At\u201E  t+At\u201E  Pi -  qi  where hu is the value of shape function in i direction at D.O.F. I.  If the variables in the  expressions for internal equivalent nodal force and external nodal force are calculated at the current iteration k, we get the residual force at the k  iteration in the incremental step expressed  th  as: t+At (k)  _ +Atp(k)  r  _ t+Atq{k)  t  ^  ^  A force convergence criterion is chosen as the ratio of the norm of the iterative residual force to the load increment, i.e.,  ,  V  f  =  l  (  1  ]  VS^*^-^)  2  &lt; t,  (4.4)  where'qi is the converged equivalent nodal force at D.O.F. / at time t, and tf is force tolerance. In order that the above criterion does not become too restrictive for small load increments, the maximum load increment during the course of deformation up to the current increment is used. Equation (4.4) becomes:  ^-J  &lt;  t  f  max{\\/sf ('+ pJ - tqi) ) At  0)  (4.5)  2  =1  It should be noted that one of the disadvantages in using a force criterion is the possible existence of inconsistent units in the force vector, e.g., forces and moments in beam elements.  4.3.4  Energy Criterion  In order to provide some indication of when both the displacements and forces are near to equilibrium values, the iterative energy, i.e., the work done by the iterative residual forces, is  \u000CChapter 4.  TWO-DIMENSIONAL  ALE PROGRAM  (ALEFE)  80  compared to the initial internal energy increment, so that:  sf (A4 =1  XN (t+At f) =1  where  t + A t  p  fc)t+At  4 ) fc_1)  _t  &lt; t  e  (  qi)Au  (4.6)  f  p j ^ is the accumulated external force at D.O.F. I at t + A i . 0  It is generally accepted that a combination of any two criteria provide an effective conver gence criteria for most engineering problems [65].  4.4  COUPLED DISPLACEMENT CONDITION  Figure 4.3: The relation between global and local velocities  The direct transformation method is applied to impose coupled displacement conditions. As shown in Figure 4.3, if a nodal point a is confined to move along a surface AB, we have  \u000CChapter 4.  TWO-DIMENSIONAL  ALE PROGRAM  (ALEFE)  81  the following relation between local and global displacement at node a: cos9  Vx(a) y(&lt;*) )  J  sinO  \u2014sinQ cosQ  (4.7)  V  In the local coordinates x', y', the boundary condition will be Vy(a)  0 and the conventional  elimination algorithm may be employed. It is first necessary to transform the stiffness matrix for the pertinent node using the transformation tensor: cos8  sinO  \u2014sinO cosO  = T (\u00AB)  (4.8)  In matrix form, the following equations are utilized: (4.9) where *v is the global velocity vector, *v' is local velocity vector and T is transformation matrix. Substituting Equation (4.9) into the final finite element equations (Equation (3.5)), we have:  'KTV  (4.10)  The above equation may be solved for 'v' by the conventional assembly and elimination procedure. It should be noted that the stiffness matrix resulted in Equation (4.10), i.e., *KT, is an unsymmetric one. To keep the symmetry characteristic, Equation (4.10) should be premultiplied by T '.. This treatment is not necessary in general A L E formulation since *K is T  originally unsymmetric. Figure 4.4 shows a flow chart for the frontal solution scheme.  \u000Cer 4. TWO-DIMENSIONAL ALE PROGRAM (ALEFE)  (Start) Detect the last apperance of each node Initialize the arrays to store global stiffness matrix, modified equations and loads Element L=1 For each D.O.F., find a|and b i I Get'K, ' K and T if any node coupled, calculate'K then 'KT c  I For each D.O.F. find the destination in wave front I Assemble the element stiffness matrix *KT 1  Check each node to find the ones which can be eliminated, extract the equations and right hand sides Modify the right hand side and the active stiffness matrix  Back-substitution, starting from last variable M Get the equation and right hand side for variable M i Calculate velocity v or displacement u or reaction force l  M  M  Transfer to get the corresponding global displacement or velocity  Figure 4.4: The flow chart for frontal solution  \u000CChapter 5  NUMERICAL  5.1  EXAMPLES  OVERVIEW Numerical examples are employed to demonstrate the correctness of the developed for-  mulation, numerical algorithm and its implementation and to highlight features of A L E in comparison with other Lagrangian and Eulerian formulation. All examples are plane strain deformation problems and may be divided into two categories.  The first category is sam-  ple problems which include rotation of an elastically deformed element, elastic-plastic small deflection analysis of a square element and pure bending of rectangular plate. This kind of problem is invoked to verify the correctness and accuracy of the program and the numerical algorithms. The second category of problems is more practical in nature and includes punch forging process, tensile necking of a elastic-plastic bar, sheet metal extrusion through curved and straight die, rolling and compression between wedge-shaped dies. These problems are generally more complicated and are applied to investigate the advantages of A L E in simulation of real industrial problems. All examples are simulated by the developed A L E F E program and one or more available commercial codes. Results are compared with available numerical and/or experimental ones from literature.  5.2  ROTATION OF A N ELASTICALLY DEFORMED  ELEMENT  A combined extension and rotation of single element is designed to show the necessity of the special numerical treatment of rate of deformation tensor 'Dij. The two-dimensional plane  83  \u000CChapter 5. NUMERICAL  EXAMPLES  84  strain element shown in Figure 5.1 at time t = 0, is subjected to displacement in the x direction, u  = 4.5 x 10~ fe and displacement in the y direction, u = \u20141.125 x 10~ 6, where b is the 4  x  4  y  dimension of the square element. Choosing a value of Poisson's ratio v = 0.2, means that the stress o-y will be equal to zero at time t = 0. A rigid body rotation of 90 degrees is imposed on the element as shown in the figure through a constant angular velocity u&gt;. 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 Z);J. The stress is integrated using Equation (3.21). Stress results 4  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  \u000CChapter 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 a . x  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 a  x  versus the rotation angle u&gt;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 cr . This may be due to the accumulation of the value x  from non-zero term in Equation (3.27). The above example shows that applying conventional time discretization or forward difference method in calculation of Dij or stress integration may cause significant errors, l  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 Dij is employed separately, e.g., in viscol  plastic deformation, the modification term developed in Chapter 3 may be applied to keep the objectivity of Dij. l  5.3  ELASTIC-PLASTIC 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 ELEMENT 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 \u2014 18.4 MPa  (8/3 ksi), Poisson's ratio v = 0.5, initial yield stress a i yie  d0  = 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. Thefinalpressure 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  \u000Cer 5.  NUMERICAL  EXAMPLES  o (Pa) x  1 .OOOOE8  \u2014  --O 1 ^  8.0000E7  ANALYTICAL SOLUTION  \u2014  MODIFIED F O R W A R D  \u2014 \u2014  CENTRAL DIFFERENCE FORWARD DIFFERENCE  \\ \\ \\  6.0000E7  \\ \\ \\  4.0000E7  -  \\ \\ \\  2.0000E7  \\ O.OOOOEO -  1  _\\  -5.1200E2  \\  -1 . O O O O E 9  \u20222.OO00E9  -3.0000E9  C  1 o  20  30  -40  SO  60  a. Stress distribution o  so  7 0  go  co At (degree) x  Xxy(Pa) 5.0000E7 -  Q  ANALYTICAL SOLUTION MODIFIED F O R W A R D  O A  CENTRAL DIFFERENCE FORWARD DIFFERENCE  4.0000E7 -  3.0000E7 -  2.0000E7 -  9.9999E6 -  -1.4753E2  .  0  . I . .  10  . .  i  20  .... i .... i . 30 40  . . .  i . . . .  50  i . ,  60  b. Stress distribution x y  .  .  I  70  .  .  .  .  I  80  .  . .  90  co At (degree)  X  Figure 5.2: The stress results of single element under rotation  \u000CChapter 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 u . Before a x  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.  \u000CChapter 5. NUMERICAL  EXAMPLES  Figure 5.4: Load versus displacement in x direction  \u000CChapter 5. NUMERICAL  5.4  EXAMPLES  PURE BENDING OF A RECTANGULAR  5.4.1  89  PLATE  P r o b l e m 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  \u000CChapter 5. NUMERICAL  5.4.2  EXAMPLES  90  Finite Element M o d e l s  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\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\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 \u20142.6036. The negative Jacobian is due to an interior angle larger than 180\u00B0. The only limitation on the shape of the real element is that each interior corner angle must be less than 180\u00B0, i.e., the element must be  \u000CChapter 5. NUMERICAL  EXAMPLES  3  4  1  2  Model-a  Model-b  4 3  4  1  1  Model-c  \\  2  Model-d  Figure 5.6: Finite element models for the pure bending of a plate  \u000CChapter 5. NUMERICAL  EXAMPLES  Figure 5.7: Pure bending of a plate: Lagrangian Formulation Results  92  \u000CChapter 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\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 FORGING PROCESS  5.5.1  M o d e l and U p d a t e d 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  \u000CChapter 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  lu l(mm) y  b. Displacement in y direction  Figure 5.8: Pure bending of a plate: A L E Results  0.080  0.090  \u000CChapter 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 convergence 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 loadfluctuationthat is pertinent to the updated Lagrangian formulation.  \u000CChapter 5. NUMERICAL  5.5.2  EXAMPLES  96  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 E D 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 E O 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 A B , nodes are treated as Lagrangian in Y-direction, and general A L E in X-direction. Boundaries D C 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 presented 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  \u000CChapter 5. NUMERICAL  EXAMPLES  Figure 5.9: The current yield stress from A L E  \u000CChapter 5. NUMERICAL  EXAMPLES  98  vertical direction to be Lagrangian while freezing the degree of freedom in the horizontal direction 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 process. 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  0.1  0.2  0.3  0.4  0.5  1.0000E7 R e d u c t i o n in height  Figure 5.10: Load-height reduction result from A L E  0.6  \u000CChapter 5. NUMERICAL  99  EXAMPLES  o (Pa) y  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 cr from A L E y  Boundary node adjustment in the mesh motion scheme discussed in Chapter 3 is necessary. 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.  \u000CChapter 5. NUMERICAL  EXAMPLES  Oyield  (Pa)  -  1.2799E9 1.2114E9 -- 1.1430E9 1.0745E9 1.0060E9 9.3761 E8 \u2014 8.6914E8 8.0068E8 7.3222E8 - 6.6375E8 5.9529E8 - 5.2682E8 m 4.5836E8 3.8989E8 3.2143E8  \u2022  m  -  \u2022  5.0  r , .0  .  ,  10.0  i  i  i  i . . . 20.0  . i . 30.0  40.0  50.0  (mm)  Figure 5.12: The current yield stress without boundary node adjustment  \u000CChapter 5. NUMERICAL  5.5.3  EXAMPLES  101  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  ii  \\ 12 mm  E  W  '\\  \\ \\ \\\\  D  C  mm  {  \\  o i\u2014  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.  \u000CChapter 5. NUMERICAL  EXAMPLES  102  Table 5.1: Displacements of top right corner from all simulations Program  u  x  Homogeneous NISA ANSYS ALEFE  mm) Inhomogeneous  15.459 18.741 13.929  14.583 18.652 12.019  mm) Homogeneous Inhomogeneous 0.003 0.012 0.096 0.021 0.010 0.040 Uy  A N S Y S 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 r  xy  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 increasement\", 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 discrepancy is believed to be caused by the severe mesh distortion which occurrs in dense meshes.  \u000CChapter 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  \u000CChapter 5. NUMERICAL  EXAMPLES  Figure 5.15: The A L E result of current yield stress from inhomogeneous mesh  \u000CChapter 5. NUMERICAL  5.6  EXAMPLES  105  T E N S I L E N E C K I N G ANALYSIS O F 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 10 , Poisson's ratio of 0.3 and the ratio of hardening modulus to yield 2  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 W = 10.0 mm in calculation. 0  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 commercialfiniteelement 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.  \u000CChapter 5. NUMERICAL  EXAMPLES  Figure 5.16: The finite element model for tensile necking  \u000Cer5.  NUMERICAL  EXAMPLES  Figure 5.17: The necking S and engineering straining 7  \u000CChapter 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 TransformationA 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  \u000CChapter 5.  NUMERICAL  EXAMPLES  109  \u000CChapter 5. NUMERICAL  EXAMPLES  Figure 5.19: Contours of equivalent stress (Truesdell rate, transformation  \u000CChapter 5. NUMERICAL  EXAMPLES  111  problem by A L E method.  5.7 5.7.1  SHEET M E T A L EXTRUSION 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 10 MPa (1 x 10 psi); hardening modulus H = 1.14 x IO 4  7  3  MPa (1.65 x 10 psi), initial yield stress tr^wo = 3.95 x 10 MPa (5.7 x 10 psi), and 5  2  4  Poisson's ratio v \u2014 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  \u000CChapter 5. NUMERICAL  EXAMPLES  Figure 5.20: Original geometry and mesh in extrusion through a curved die  \u000CChapter 5. NUMERICAL  113  EXAMPLES  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 a at different lateral positions, x  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  \u000CChapter 5. NUMERICAL  EXAMPLES  114  Figure 5.22 as well. The values of shear stress r  xy  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\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 &lt;r MPa.  yieW0  = 267.7 MPa and strain hardening modulus H = 85.9  Other material properties are: Young's modulus E = 6.9 x 10 MPa and Poisson's 4  ratio v = 0.32 [72]. Because the geometry and boundary condition are symmetric, only half of the billet is considered. The finite element model is shown in Figure 5.24. The same problem is simulated by NISA and ANSYS. In the A L E simulation, all the nodes in patch B C F G in Figure 5.24  \u000Cer 5. NUMERICAL  EXAMPLES  a. Result from A L E  Figure 5.22: Comparison of o~ from A L E and reference [12] x  \u000CChapter 5. NUMERICAL  EXAMPLES  a. Result from A L E  T xy  (ksi)  \" 7  5  \"  iz -A I REDUCTION'  b. Result from reference  Figure 5.23: Comparison of r  xy  from A L E and reference [12]  \u000CChapter 5. NUMERICAL  EXAMPLES  117  are kept as Eulerian points, so that the contact boundary conditions of the die work-piece are satisfied throughout the deformation history. The nodes on boundaries A B , C D and H G are kept as Lagrangian points in Y direction and general A L E points in X direction. The nodes on A H , F E and E D are treated as Lagrangian points to keep one-to-one mapping between mesh domain and material domain. Once again, the internal points are handled by the mesh motion scheme discussed in Chapter 3. The piston motion is maintained until some distance after the steady deformation is achieved.  y ( mm) 2.0  './//////////  A/ / X ,  1  G  1.0 F  E  0.5  |  1.27 (mm)  1.5  c)  f  J 1  L  J 2  1  \u2022  _L_ 3  m^f  ii\u2014L  \u2022  &lt;[  E3  x (mm)  5  C  D  Figure 5.24: The finite element model for straight die extrusion  Figure 5.25 shows contours of current yield stress (plastic zones) and the deformed mesh from A L E after a piston motion of 1.2 mm. It may be seen that the mesh is homogeneous in the whole domain with the highest deformation area just before the exit around the die corner. This is different from the curved die, that was considered in the previous section, where the largest deformation is located just after the exit. This may be attributed to the zero curvature  \u000CChapter 5. NUMERICAL  118  EXAMPLES  and slope at the exit.  Oyieid (MPa)  y (mm) 2.0 r 1.5 f1.0  371.863 364.916 357.969 351.022 344.076 337.129 330.182 323.235 316.288 309.341 302.394 295.448 288.501 281.554 274.607  0.5 0.0  x (mm)  Figure 5.25: The current yield stress (plastic zone) from A L E  With the same mesh used above, both NISA and ANSYS analyses were terminated automatically by the programs after a piston motion of only about 0.18 mm because the mesh was entangled around the exit of the die. Therefore, a coarser mesh is used for NISA and ANSYS simulations, as shown by the gray lines in Figure 5.26. With the new model, NISA results were successfully obtained for a piston motion of 0.54 mm. The final distorted mesh at the end of the run is shown by the black lines in Figure 5.26. It may be noted that most of the distortion is due to the in-ability to update the boundary conditions. For the ANSYS run, the program did not converge at 0.19 mm of piston displacement and the simulation was terminated due to excessive mesh distortion. Results of the extrusion forces from A L E , ANSYS and NISA are compared with the  \u000CChapter 5. NUMERICAL  EXAMPLES  Initial m e s h Distorted m e s h  Figure 5.26: The initial and deformed mesh from NISA  \u000CChapter 5. NUMERICAL  EXAMPLES  120  experimental one from [71], as shown in Figure 5.27. It is observed that the extrusion force from A L E is in close agreement with the experimental one. The A L E and experimental one show a steady force after the piston moves about 0.2 mm.  The extrusion forces from NISA  and A N S Y S clearly did not show this behavior. An important reason for this is believed to be that the contact boundary condition is not updated automatically, so that the reduction in height increased with the piston movement, which is not practical. Figure 5.28 shows a comparison of the pressure distribution on the piston face from A L E with the one from slip line method [71]. It is observed that the pressure from A L E is in good agreement with the result from slip line theory except at points close to the surface of the work-piece (i.e., y = 1.27 mm). The reason of this difference is not clear, but at y = 1.27mm, the result from slip line theory give a tension stress, which is not reasonable.  1500  i--  ALEFE Experiment NISA ANSYS |  1000  0 0.0  //  0.2  0.4  0.6  0.8  1.0  1.2  Displacement of piston (mm)  Figure 5.27: Comparison of extrusion force  1.4  \u000CChapter 5. NUMERICAL  0.0  EXAMPLES  0.2  0.4  121  0.6  0.8  1.0  1.2  1.4  Location in thickness direction (mm)  Figure 5.28: Pressure distribution on piston face 5.7.3  Discussion  The metal extrusion examples indicates that A L E formulation is able to handle changes in the kinematic boundary conditions easily and keep mesh homogeneous in the whole deformation process. In traditional updated Lagrangian method, this may only be achieved by particular boundary updating scheme and a specially developed program, as in [12]. The effect of die shape on the plastic zones is clearly reflected in the results. The stress, extrusion force and pressure distribution generally agrees with other numerical results from the literature [12], as well as experiment and slip line [71] results.  \u000CChapter 5. NUMERICAL  5.8  EXAMPLES  122  STRIP ROLLING  5.8.1  Particular Difficulties in Rolling Process  Rolling is one of the oldest industrial metal working processes. In view of the tremendous amount and wide variety of rolled products manufactured every year, it may be considered as one of the most important. For more than half a century, numerous investigations, both analytical and experimental, have been carried out on rolling [73]. Slab analysis, slip-line field and upperand lower-bound methods have been widely used in theoretical analyses of rolling. However, owing to the complexities of deformation involved in rolling, various degrees of simplification and idealization are necessary. These methods have provided useful but limited information in rolling. The increase in quality requirements and the concern for maintaining and improving the competitiveness of rolled products have led to demands for more detailed information and understanding of the process.  Under more realistic conditions, accurate determinations of  detailed metal flow in rolling processes became possible only when the finite element method was introduced into the analysis. Numerical simulation of rolling is challenging. Besides boundary condition updating and mesh distortion due to large deformation, friction on the roll work-piece contact boundary is more important and complicated than other forming processes. Because of the change in relative speed between the roll and the work-piece, the friction force along the contact surface (arc G F in Figure 5.29) changes its direction at the neutral point N, as shown in Figure 5.29. The location of neutral point changes during the non-steady state deformation and becomes steady at somewhere between the exit section and the halfway point of the contact surface. The location of the neutral point depends mainly on the friction condition and is not known a prior. A difficulty arises when attempting to implement Coulomb friction model in which the tangential frictional stress distribution is dependent on the normal pressure distribution. Some researchers developed particular formulations to handle the friction [74]; others employ  \u000CChapter 5. NUMERICAL  EXAMPLES  123  friction layer or interface elements [42, 75]. In this research, we do not introduce the particular formulation to handle rolling simulation. Instead, an approximate method based on the general program used above is applied. Rolling process is simulated by Eulerian method reduced from the general A L E formulation.  y  Direction of rolling  Figure 5.29: The definition of geometry in rolling  \u000CChapters.  NUMERICAL  EXAMPLES  124  Table 5.2: Geometric dimensions of rolling process, (Length unit: mm) Case No. 1 2 3 4 5  5.8.2  Initial Height ho 6.274 6.274 6.274 6.274 2.057  Final Height hi 5.385 4.902  R/ho 12.5 12.5 12.5 12.5 39.0  4.445 4.115 1.588  Reduction 14.17 % 21.86 % 29.40 % 34.41 % 22.83 %  Deformation Conditions in Rolling  A series of cases of strip rolling of aluminium with different dimensions have been numerically analyzed with the developed program A L E F E . For comparison, the material properties and geometry are chosen to be the same as in reference [76]. The material is elasto-plastic with Young's modulus E = 6.894 x 10 MPa, Poisson's ratio v = 0.30 and initial yield stress 4  fyieido  = 50.3 MPa. The material has isotropic strain hardening with the following relation: /  O-yield =  where, c r \u00A3  yieW  eq  p  \u00A3  (TyieldO I 1 +  X0.26  j  (5.1)  is current yield stress, 1 S  equivalent plastic strain.  In the simulation, a piece-wise linear model for the plastic hardening modulus H is used to approximate Equation (5.1). The geometric dimensions of simulated rolling process are shown in Table 5.2. Two values of the ratio of the roll diameter to the initial height of work-piece, Case 1 and Case 5, are employed to show the different types of the normal stress distribution on the contact surface. A series of reduction under same ratio R/h  0  are. applied to examine the  change of rolling force and torque with reductions. The friction condition on the interface is determined by experiment [76]. It can be described by Coulomb law: r  =  pp  (5.2)  \u000CChapter 5. NUMERICAL  EXAMPLES  125  where, T is friction stress; p is normal pressure; and p is equivalent friction coefficient, (p = 0.1 for the studied cases). The roller is assumed rigid and the elastic deformation of work-piece before getting into and after leaving the roll gap is ignored. The deformation is simulated under plane-strain condition. Because of symmetry, only one-half of the model is analyzed. The finite element model for case-1 is shown in Figure 5.30. The model consists of 480 four-node quadrilateral isoparametric elements.  Figure 5.30: Finite element model (Case 1: R/ho = 12.5, Reduction 14.17%)  5.8.3  N u m e r i c a l Treatments i n the Simulation  The steady state strip rolling process is simulated by Eulerian formulation Equation (2.26) which is degenerated from the general A L E formulation Equation (2.20). Referring to Figure 5.29, it is assumed that there is a spatially fixed boundary, A B C D E F G H , and the material crosses this boundary and goes out of the rolls. New material being continuously convected into the fixed region has prescribed initial state, usually no deformation and stresses. The simulation is carried on until the steady state is achieved, in which all the variables, i.e., velocities, stresses  \u000CChapter 5.  NUMERICAL  EXAMPLES  126  and nodal forces etc., do not change as the new material enters the rollers. The Coulomb friction is calculated from the normal pressure distribution in the previous incremental step. At the first incremental step, the friction stress is assumed to be homogeneous on the whole contact surface and given by: T = 0.02fc  (5.3)  o  where, k is the initial shear yield strength, k = o- i / \\/3. 0  0  yie  d0  Such treatment should not impact the final results, since the simulation does not finish until differences of all variables between two consecutive iterations, including normal pressure, are less than the preset tolerances. This treatment, however, makes the handling of the friction effect easier and more straight forward. Specification of boundary condition on contact surfaces is now discussed. A l l the nodes on the contact surface are sliding except the neutral points which may spread over a small area. A tangential stress equal to \"pp\" is applied to sliding nodes. The direction is dependent on whether the node is in forward sliding area or backward sliding area. The velocity boundary condition in the neutral area is that the neutral point has a circumferential speed equal to that of the roller, i.e., there is no relative slip between the roller and the neutral point. However, a prescribed velocity boundary condition may cause significant reaction force in tangential direction which is not rational. An alternative method is to specifiy a prescribed horizontal velocity at the boundary H A , as shown in Figure 5.29. The velocity is calculated by: v  ho  = \u2014r\u2014  (5.4)  where, h is the height of work-piece at neutral point, n  v  rh  is the horizontal component of roll speed.  Equation (5.4) assumes homogeneous v t  x  distribution on the cross section passing through the  neutral point N M and a volume consistency between N M and HA. This assumption, is trivial  \u000CChapter 5. NUMERICAL  EXAMPLES  127  and should not impact the simulation results. Another requirement to be considered here is the total reaction force caused on edge HA. This force should be much smaller than the reaction forces elsewhere if no other external traction is applied in the rolling. In simulation, the above two ways of applying the velocity boundary condition are employed and compared.  5.8.4  Simulation of Neutral Point Location  The existence of a neutral point is the unique feature of rolling. Its location directly affects the results of simulation, so it is important to correctly detect or simulate it. However, finding the neutral point is not easy. The neutral point has the same speed as the roller, so it is natural to think that the neutral point should be the point on the contact surface that has a velocity equal or close to the tangential roller velocity. To find such a point, an additional iteration loop needs to be nested outside the equilibrium iteration in each incremental step. At the beginning of a new step, we apply the boundary conditions according to the results from the last converged step, then find the new location of the neutral point from the calculated nodal velocities. If the difference between the new and old neutral points is not within the specified tolerance, the boundary conditions have to be updated according to the new neutral point. The process continues until convergence occurs. Results of the neutral point simulation by the above procedure are shown in Figure 5.31. Thefigureshows the x-location of the neutral point as a function of the incremental step. It is observed that the neutral points for all five-reduction cases did not change from the initial guess throughout the deformation process and that different initial guess lead to different results. Two different initial guesses are made for Case 1 and shown in Figure 5.31, labelled as Case la and Case lb. The distribution of rolling pressure and friction stress on the contact surface with respect to angle 8 in Figure 5.29 is also examined. When velocity boundary condition is applied at neutral point, both rolling pressure and friction stresses show a sudden \"jump\". The results from Case  \u000CChapter 5. NUMERICAL  EXAMPLES  128  1 are shown in Figure 5.32 (curve ALE-NEU). The \"jump\" occurs at the neutral point. It is believed that this is due to the reaction resulting from the application of the boundary condition. On the other hand, if the velocity boundary condition is applied on the external boundary H A , the neutral point fluctuates and may not in general converge. Figure 5.32 (curve ALE-BOU), shows the distribution of rolling pressure and friction stress for Case 1 from such simulation. We see that the \"jumps\" are eliminated, but the distribution patterns are quite different. Another important point to be considered here is that the reaction force on H A is quite high (108  N/mm)  and may not, in general, be ignored if compared to the total rolling force (845 N/mm).  In  addition to these points, we see that the above simulations of neutral points may not maintain equilibrium of external forces.  14  12 I  |  y-  l  y  \u00A3  '-a  8  \u2014f \"\"  --1-  t  y  10  ..f.-a \u2014f-  O)  c  -  5  \u00A3 cb  4 h  ^ o  ^o-  \u2022 CASE1.A - - O - - CASE1 ,B \u2022 CASE2  --&lt;!&gt; \u2014 o  CASE3  A  CASE5  are shown in Table:5.2  , i , 200  CASE4  O  ' Parameters lor different cases  \u2014-f\u2014  100  V  300  '  ' \u2022 ' \u2022 400  1  1  '  500  600  Step No.  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 selfbalanced. That is a unique feature of the rolling process and generally makes the simulation of  \u000CChapter 5. NUMERICAL  EXAMPLES  129  (MPa) 400  0.00  0.02  0.04  0.06 Arc of contact  0.08  0.10  0.12 (Rad.)  p  a. The normal stress (MPa) 30  0.00  0.02  0.04  0.06 Arc of contact  0.08  0.10  0.12  (Rad.)  (3  b. The friction stress  Figure 5.32: Rolling pressure and friction stress when velocity criterion is applied  \u000CChapter 5. NUMERICAL  EXAMPLES  130  rolling more complicated. The above methods of deciding neutral points introduced resultant forces on either the free edge H A 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  J/3  N  - pRsinB)  d8 -  {TRCOS8 + pRsinB)  d8 = 0  (5.5)  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 H A 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 H A (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  \u000CChapter 5. NUMERICAL  EXAMPLES  131  method is more reliable and more accurate than the ones previously discussed.  \u2014l  10 100  200  300  400  500  Step No.  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/h , two-peak roll pressures are obtained with maximums near entrance 0  and exit and with pressure drop in the middle of contact arc. For cases of large R/h , Figure 5.36 0  shows a friction-hill distribution. It may be observed that the computed results are, generally,  \u000CChapter 5. NUMERICAL  EXAMPLES  350  0.04  0.06  0.12  0.08  Arc of contact  p  (Rad.)  a. The normal pressure 20  i  15  7  !  _  (0 0_  10  2 to c  0  o o  \\Y  -5  -10  -15  -20 0.00  0.02  0.04  0.06 A r c of contact  0.08 P  0.10  0.12  (Rad.)  b. The friction stress  Figure 5 . 3 4 : Rolling pressure and friction stress when force criterion is applied  \u000CChapter 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  0.10  0.12  Arc of contact p (Rad.)  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/h  Q  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  \u000CChapter 5. NUMERICAL  EXAMPLES  134  Figure 5.36: Comparison of rolling pressure, large R/h  0  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  \u000CChapter 5. NUMERICAL  EXAMPLES  135  2500 r-  2000  , 1000  ...^^Or.  I  ,  500  O EXPERIMENT \u2022 - \u2022 - - ALE-BOU \u2022\u2014* ALE-NEU I  0.15  0.20  _L 0.30  0.25  0.35  0.40  Reduction in height  Figure 5.37: Comparison of rolling force  1.4000E4  1.2000E4 h  1.0000E4  8.0000E3 0  6.0000E3  o r-  4.0000E3  2.0000E3  0.0000E0 0.10  0.15  0.20  0.25  0.30  0.35  Reduction in height  Figure 5.38: Comparison of rolling torque  0.40  \u000CChapter 5. NUMERICAL  136  EXAMPLES  the plastic zone is developed. In the simulation, a steady state condition was reached after 300 incremental steps.  x  Figure 5.39: Equivalent plastic strain (Case 1: R/ho = 12.5, Reduction  5.8.7  (mm)  = 14.17%)  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 conditions 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  \u000CChapter 5. NUMERICAL  EXAMPLES  137  'yield &lt; &gt; MPa  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)  a. Step 120 \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  x (mm)  b. Step 210  yield (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  10  x (mm)  c. Step 300  Figure 5.40: Current yield Reduction = 22.83%)  stress at different  steps (Case 5:  R/h  0  39.0,  \u000CChapter 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  COMPRESSION  5.9.1  BETWEEN WEDGE-SHAPED  DIES  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\u00B0 is used because most of the comprehensive experimental results are available in literature for this angle [78]. The specimen A B C D E F G P is compressed between two inclined dies B C and F G with angle of inclination n with the horizontal direction. The specimen is pre-shaped such that surface B C and F G are to be matched with the tool face. The dotted lines indicate the approximate shape after deformation. Velocities v and v are at large and small end of the t  s  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/in ) 2  were measured independently by means of two calibrated  \u000CChapter 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 1 2 3 4 5 6 7 8  m in  degree  mm  3\u00B0 3\u00B0 3\u00B0  25.40 22.23 19.05  1.000 0.875 0.750  3\u00B0 3\u00B0 3\u00B0 3\u00B0 3\u00B0  15.88 12.70 9.53 6.35 4.45  0.625 0.500 0.375 0.250 0.175  Ho mm  in  22.74 19.56 16.39 13.21 10.04 6.86 3.69 1.78  0.895 0.770 0.645 0.520 0.395 0.270 0.145 0.070  \u000CChapter 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 stressstrain curve is obtained by experiment [78]. Using bilinear approximation of the constitutive curve, we use initial yield stress er i  yie d0  = 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 F C 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 A B and C D 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 B C are Lagrangian points in vertical direction and Eulerian points in horizontal direction. All other boundary nodes  \u000CChapter 5. NUMERICAL  EXAMPLES  141  y (mm) 15 M M  10  ITTTTTTTTTTTTI  ^  E  5  n  11. I .1 . Ii 1 .1 . 1. 1 il.1.1.1 . 1 I I. I.I . 1 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  B  10  20  30  40  50  e  60  x (mm)  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\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  \u000CChapter 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  H /h . c  c  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 p /2k v  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  \u000CChapter 5. NUMERICAL  EXAMPLES  Figure 5.43: The variation of vertical load on the die with compressive strain  \u000CChapter 5. NUMERICAL  EXAMPLES  144  Figure 5.44: The variation of mean vertical die pressure with current strip thickness  \u000CChapter 5.  NUMERICAL  EXAMPLES  145  \u000CChapter 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.  \u000Cer 5. NUMERICAL  EXAMPLES  P eq 0.376818  L  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  B  p  ^eq 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  \u000CChapter 6  CONCLUSIONS A N D F U T U R E W O R K  6.1  CONCLUSIONS 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 formulation 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  \u000CChapter 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  Dij  l  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 A t , 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.  \u000CChapter 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.  \u000CChapter 6. CONCLUSIONS AND FUTURE WORK  151  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 interpolation 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  \u000CChapter 6.  accuracy.  CONCLUSIONS  AND FUTURE  WORK  \u000CBibliography  [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 Geometric 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 Formulation 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 Formulation 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, U S A , 1995. [9] Swanson Analysis Systems Inc., Engineering Analysis System User's Manual, Swanson Analysis Systems Inc., Houston, U S A , 1992. [10] Engineering Mechanics Research Corporation, Users Manual for NISA Systems Analysis, Engineering Mechanics Research Corporation, Michigan, U S A , 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  \u000CBibliography  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 Analyses of Metal Forming\", In: Numerical Grid Generation in Computational Fluid Mechanics' 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 Simulations 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 MechanicsC 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\", Computers &amp; 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.  \u000CBibliography  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 EulerianLagrangian Code\", In: Methods in Computational Physics (Alderb, B . et al eds), Vol.3, Academic Press, New York, U S A , 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, U S A , 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 EulerianLagrangian Finite Formulation\", In: Numerical Methods in Industrial Forming Processes ( Pittman, J . F . T . eds), pp. 501-509, Pineridge, Swansea, U K , 1982.  \u000CBibliography  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 Largedeformation 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 Eulerian 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 E l 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 Arbitrary Lagrangian-Eulerian Description\", J . Mat. Shaping Technol., Vol.8, pp. 53-64, 1990.  \u000CBibliography  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, U S A , 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 System 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 Analyses 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 Coordinate 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 Computers 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.  \u000CBibliography  158  [57] Hibbitt, Karlsson and Sorensen, Inc., A B A Q U S Theory Manual, Hibbit, Karlsson &amp; Sorensen, Inc., Pawtucket, U S A , 1992. [58] Pinsky, P. M . , Ortiz M . , Pister, K . S., \"Numerical Integration of Rate Constitutive Equations in Finite Deformation Analysis\", Computer Methods in Applied Mechanics 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, PrenticeHall Inc, New Jersey, U S A , 1969. [65] Bathe, K . J . , Cimento, A . P., \"Some Practical Procedures for the Solution of Nonlinear 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 &amp; Sons, New York, U S A , 1989. [68] Burnett, D . S., Finite Element Analysis, from Concepts to Applications, Addisonwesley Publishing Company, Menlo Park, California, U S A , 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.  \u000CBibliography  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 Strainhardening 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, U S A , 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 Rigidviscoplastic 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 Determination 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, U S A , 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.  \u000CAppendix A Material and Referential Time Rate  A physical quantity'/ = /('x) at Material Reference System (MRS) point P ( ' : c i , ' x ,' x ) 2  at time t becomes  t + A t  f  = / ( ' ' x ) at time t+At as the point P moves to Q ( +A  t + A t  x  t + A t u  3  x ,  t+At  2  x  The MRS change of ' / may be expressed by: A / x  =  /('  + A  'x,i + Ai)-/('x,i)  (A.l)  Similarly, the change o f ' / at the corresponding Computational Reference System (CRS) point P^Xif  X2, X3) 1  w  n  e  n  t  n  e  P  o i n t  moves to R( x\\, t+At  t+At  X2,  t + A t  X3)  is denoted by:  A / = /r'x,* + A0-/( X,*) A  f  x  To derive the relation between A  / and A f,  x  x  (A.2)  we assume that the MRS and CRS points  coincide at point P at time t, as shown in Figure (A.l), so that:  and, f( X,t)  = f( x,t)  t  t  = f t  (A.3)  at time t + At, t + A  t  +  A  ' x =' x + A x  (A.4)  =*X + &amp;X  (A.5)  t X  160  3  \u000CAppendix A. Material and Referential Time Rate  161  Figure A . l : The description of relation between mesh motion and material motion  \u000CAppendix 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*( x\\, a;*./ X 3 ) with displacement A x * , we have: t  t  = * * + Ax* =  t  +  A  x  f( X,t  V  (A.6)  + At) = f( x,t  t+At  + At)  t+At  (A.7)  Substitute Equation (A.3) and (A.7) into Equation (A.2): Af x  =  fC+^^t  =  A / + (/(  =  A /+V/(  + A^ + fC+^t  + A^-fC+^t t + A t  x  x * , t + At)- / (  t + A f  x  x , t + At)\u00BB (  t + A t  i + A t  x  x  _  = =  4  A^-fC^t)  x , t + At))  x* -  t + A t  x) + h.o.t.  where V is gradient symbol and h.o.t. is higher order term of ( * + t+A* * _ t + A t  +  + * ^ _*+At  At x  * \u2014  t + A t  x ) . But,  A  x  'x + A x - ' x - A x Ax-Ax  Therefore, A / x  =  A f + Vf{ x,t t+At  x  + At)*(A -Ax) X  + h.o.t.  Where h.o.t. is higher order term of (A% \u2014 A x ) . If h.o.t. is ignored, we have:  A  */ = */ + | ^ K - ^ ) A  ( -) A 8  Dividing the above equation by At and taking the limit when At \u2014 &gt; 0,  7  A  = f + ^Cvl f  - *v ) k  (A.9)  \u000CAppendix B  U p d a t e d Lagrangian Formulation  An updated-Lagrangian formulation given by Rice [3] is based on the following equation:  / (v^-Iv\u201E,(, v  =  j f -8v d V t  tv  B  i  +f  t  i  -||^))^  W  tfS-SvutS  (B.10)  where ' A j = \\ (\u00A7r^ + f t ^ ) 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: ' tj T  T  J \u00B0'ij  +  ~ J &lt;7ij t  T  t  using the Jaumann rate for both 'r^- and V^-, we get: t-J T *3  t  T  t  _  I  ran  t  +'J  Tt-J 4  &lt;  (B.ll)  Since ' r ^ is symmetric, the right hand side of Equation (B.10) may be expressed by:  f  (t j  U { =  f  i  dSv  T i j  d S v  i  &amp;Xi  1  ( k  dSv \\  dSv  2 ^ 0*3;  t  +  ld v l  t  a\u00ABaj J  k  a  i  (&amp;v, , d*v \\  (t-J  lv\u00A5x~A --^{o^  k  T  +  k  j  {d^  *a  jk  +  d'vA  d8v  d*x )  d*Xi  k  ( 5 ^  o^J ^{d^ -\u00A5x-Jr  v  +  k  163  kt  f  c  \\ \\  dW\\j ^frxj  -  (B 12)  \u000CAppendix  B. Updated Lagrangian  Substituting * J =  Formulation  164  [48], and choosing the current state as a reference so that * J = 1.  The Equation (B.l 1) will take the form:  V  -  J  *\u00AB\u2022  9  t  V  4 - rr  k  J  Substituting the above equation in Equation (B.12), we get: LHS  =  ( p ( ^ hv o x \\ t  +  J  - ^ - ' ^ ( p . ox 2 ya'xfc  j  p L )  +  t  k  o Xj) t  +  ' ^ ( p . . p L ) ) 2 \\o x dx J J t  &lt;  t  V  t  k  i  So the updated Lagrangian formulation Equation (B.10) takes the following form: dSvj L .j _ vik( t  v d'xj \\  =  a  J Sv\\f?*V tv  i  j  2 [d'xk  + j 8v\\f?*S tg  dW\\  #S +  d'xj  +  fjk[  l  2  d%_ _ dW\\ [frxk  d'xij  +  t  dW\\ j ^x )  (Tii  k  (B.13)  Equation (B.13) is same as Equation (2.25) in Chapter 2 except the load correction terms introduced in Equation (2.25).  \u000CAppendix C  M a p p i n g Between M a t e r i a l and M e s h D o m a i n  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., J = T  ^ 0. This may  C  be verified as follow: At time t \u2014 At, assume '~ a:; =*At  from t \u2014 At to t the material velocity is t  f i , then,  t _ A t  l t-At,  xi + At'-^vi  x  Xi=  ^ i , i.e., material points and mesh points coincide and  At  Xi + At  Vi  =  = - Xi l  + At  at  1  Xi + At (a; + b  (C.14)  v j  {i)  {i)  From which we obtain a relation between x and Xi as: t  t  i  H  = =  t X i  -At(a  *  Xi  i  +  b - v )+At - v t  {i)  - Ate* + ( l - 6 ) (i)  At  t  At  [i)  i  (C.15)  At'-^vu  and the determinant of the Jacobian transformation becomes: d Xi l  d Xj l  8^ + ( l - 6 ) A t (0  No summation on  (i).  Since * Vi At  is independent of tjc  I  (C.16)  d'Xj we have  d d t  ^  x )  0 and:  u  (C.17)  165  \u000CAppendix 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:  *l? = V \u00AB + Aery = \\+**S\u2122  t+At  where,  t s[p +At  (D.19)  is the 2nd Piola-Kirchhoff stress at an imaginary intermediate state with respect  o-\\p is the corresponding Cauchy stress. The 2nd Piola-  t+At  to the configuration at time t and  Kirchhoff stress is related to Cauchy stress by: 1  t+*t  d  a:  f + A t  d Xj t+At  i t + A t  3x  t+At-g  dx  m  n  and since the second step is only a rigid body motion, then,  l S +At  mn  = ?\"S\u2122  (D.21)  Substituting Equation (D.21) and (D.19) into Equation (D.20), we get: t+At\u201E  _ 1 0  1  \u2022*-  &lt;+*'F  d  w  t + A t  x  d  i t  ~lt\u201E  d*x  m  \"  t + A t  *\u00B1  m  n  x  j  0 i  d*x  1  ' t+ Y At  n  d  .  t+At  '  d*x  \" 77171  m  d  t+At  Xi  m 991  d*x  n  Equation (D.22) is same as the one developed by Pinsky [58] and both are consistent with Truesdell stress rate constitutive equation. 166  "@en
<https://open.library.ubc.ca/cIRcle/collections/ubctheses/831/items/1.0080928>