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 © Jin Wang, 1998 In 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> 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 Abstract The applicability and accuracy of existing finite element formulation methods for finite strain deformation and metal forming problems are investigated. It is shown that the existing formulation methods, both Lagrangian and Eulerian type, are not suitable nor efficient for large deformation problems especially when boundary conditions change during deformation, as is the case in most metal forming problems. This creates a need for a more general and efficient type of formulation. An Arbitrary Lagrangian-Eulerian (ALE) method is presented for the general application in solid mechanics and large deformation problems. A consistent A L E formulation is developed from the virtual work equation transformed to an arbitrary computational reference configuration. The formulation presents a general approach to A L E method in solid mechanics applications. It includes load correction terms and it is suitable for both rate-dependent as well as rate-independent material constitutive laws. The proposed formulation reduces to both updated Lagrangian and Eulerian formulations as special cases. The formulation is presented in a form that makes the programming an extension to existing Lagrangian and Eulerian type programs. An efficient mesh motion scheme for A L E formulation is developed with a procedure for handling boundary motion within the scheme, which can ensure homogeneous mesh results. A practical and more efficient numerical method is presented to handle supplementary constraint equations on element level rather than on the global level. Different numerical algorithms for the integration of the rate type constitutive equation are investigated and coupled with the return mapping algorithm to provide plastic incremental consistency. A numerical procedure for stress integration is developed based on the physical meaning of stress. Jaumann and Truesdell rates are taken as the objective stress rates in the constitutive equation. An alternative numerical ii treatment for rate of deformation tensor 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 Table 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 2.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 4 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 5.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 6.2 FUTURE WORK 151 Bibliography 153 Appendices 160 A Material and Referential Time Rate 160 B Updated Lagrangian Formulation 163 C Mapping Between Material and Mesh Domain 165 D Stress Increment Procedure 166 viii List o f Tables 5.1 Displacements of top right corner from all simulations 102 5.2 Geometric dimensions of rolling process, (Length unit: mm) 124 5.3 Specimen dimensions in compression between wedge-shaped dies 139 ix List of Figures 1.1 The meshes before and after deformation 3 1.2 Not updated boundary condition 4 1.3 Distribution of stress 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 5.3 Elastic-plastic small deflection model 87 5.4 Load versus displacement in x direction 88 5.5 Pure bending of a rectangular plate 89 5.6 Finite element models for the pure bending of a plate 91 5.7 Pure bending of a plate: Lagrangian Formulation Results 92 5.8 Pure bending of a plate: A L E Results 94 5.9 The current yield stress from A L E 97 5.10 Load-height reduction result from A L E 98 5.11 Contours of 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 5.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 Nomenclature General Variables: coefficients in relation of material motion and mesh motion area of an element b dimension of a square element Ci coordinates of a point P old coordinates of a point P Cijki material constitutive tensor 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 hu 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 qj^ equivalent nodal force in k v normalised parametric coordinate 0 < r < 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 < s < 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°j 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° 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 v\ 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<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 &(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 <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> objective function in the least square method 4>i{r,s) boundary coordinate of a geometric domain X position vector of referential coordinate Xi mesh point coordinate in axis i u> angular velocity V/ gradient of / 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 material time rate of a quantity A computational reference time rate of a quantity T transformation of a matrix xviii Acknowledgement I would like to express my deep appreciation to my supervisor, Dr. Mohamed S. Gadala, for his invaluable help, judicious advice and understanding throughout my graduate studies. His patient guidance and continuous encouragement effectively helped me to proceed with this work. I express my gratitude and appreciation to Dr. Henry Vaughan for his invaluable suggestions and advice. I would also like to thank Dr. Stanley G. Hutton and Dr. Reza Vaziri for their help and valuable discussions and suggestions at various stages of this work. I acknowledge the financial support of Natural Science and Engineering Research Council of Canada, of 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 Chapter 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 Chapter 1. INTRODUCTION 2 not easy to simulate non-steady static or dynamic behavior. Trials have been made to relate the moving material points (in terms of F E Gauss points) to the fixed spatial mesh [6, 7]. Much work is still required, however, to refine and establish this approach. This chapter investigates the applicability and accuracy of existing formulation methods to the finite strain deformation problems. The objective is to examine some of the above aspects through the use of available 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 Chapter 1. INTRODUCTION 3 results. The results obtained from the two programs showed close agreement in displacement and strains but revealed noticeable discrepancies when comparing stress values. Also, there is noticeable disagreement between the results from both programs and similar ones reported in the literature [11]. w/2 i•4- ~~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 Chapter 1. INTRODUCTION 4 Figure 1.2: Not updated boundary condition 2-node gap elements. The nodal and deformed shape after the piston moves 0.3a are shown in Figure 1.2 by grey and black lines, respectively. The distributions of displacement 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. Chapter 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 . —•— Ox(MPa) N o G a p Element —A— Ox(MPa) G a p Element , .i. , . . -0.5 i 0 . . . . i .... i 0.5 1.0 Figure 1.3: Distribution of stress <r on outer surface x ™a. Chapter 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 <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 Chapter 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, Chapter 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 «DSCA=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 Chapter 1. INTRODUCTION • 12 mm i •+ 30 mm •> Figure 1.6: F i n i t e element model used i n D E F O R M 2 simulation Chapter 1. INTRODUCTION 10 the predicted deformed boundary of the free surface near the punch corner is unusual and not rational. Another shortcoming is the fluctuations in loads predicted by Lagrangian method. Figure 1.8 shows a plot for the reduction of the specimen height versus the applied load, from DEFORM2. The result shows noticeable 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 —j..r.r.|..T.r.rT..rr.y.T..rr.¥.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) Chapter 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 -• \ 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 Chapter 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 Chapter 1. INTRODUCTION 13 higher quality finite element meshes, but it is time consuming and, more significantly, may still produce poorly shaped meshes near boundaries in some cases [17]. The problems of contact and friction boundary conditions pose another serious concern in the solution of metal forming applications. Contact condition, which is defined from the geometrical compatibility on the contact surface or the impenetrability condition, has been introduced by the direct method, elements with special mechanical properties, the Lagrangian multiplier method and the penalty function method [18]. Since the contact area is a prior unknown, the boundary conditions of the contact problems are determined as part of the solution. Although the contact surfaces obtained by the Lagrangian multiplier method satisfy the contact condition in the integral sense, additional unknowns are required and the total number of unknowns in the system equations increases. In addition, the associated tangent matrix may be indefinite and has zero diagonal entries that pose some difficulties in the solution steps [19]. For the penalty function method, the solution results satisfy the contact conditions only approximately. The accuracy of the approximate solution depends strongly on the penalty parameter, so the correct choice for these parameters is the essence of the algorithm. When the penalty parameter is chosen to be too large, it leads to numerical problems in the form of loss of accuracy in the solution. On the other hand, a too small choice, results in unacceptable penetration of one body into the other. The node to node or node to surface gap or contact elements have been recently introduced in commercial 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 Chapter 1. INTRODUCTION 14 stiffness characterizations, i.e., poorly conditioned Jacobians. They are also awkward to apply in situations requiring friction effects and they require significant amount of equilibrium iterations. Furthermore, parameters of the contact element have to be chosen by the user, and the convergence speed and the accuracy of the solution are quite dependent on the choice of these parameters. Sometimes, especially when friction effect is incorporated, convergence may be difficult to achieve, as indicated by the experience from solving the plane strain forging problem. 1.3 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 Chapter 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 Chapter 1. INTRODUCTION 16 of the reference domain is represented by a set of grid points, which may be interpreted as the movement of a finite element mesh. Therefore, in an A L E formulation, the finite element mesh need not adhere to the material or be fixed in space but may be moved arbitrarily relative to the material. A proper A L E formulation should reduce to Lagrangian formulation if we choose to use the same motion for the computational and material meshes. On the other hand, if we choose to fix the computational mesh, an A L E formulation should reduce to Eulerian formulation. Combining the merits of both Lagrangian and Eulerian formulation, A L E is more suitable to handle mesh distortion and entanglement and special boundary condition changes in metal forming problems. If the nodes on tool-workpiece interface are specified as Eulerian points in tangential direction, it may eliminate 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 er 1. INTRODUCTION 17 for differences and their applicability to large deformation metal forming problems. Develop a complete A L E formulation that may be applied to general finite deformation cases and degenerate to either updated Lagrangian or Eulerian formulation as extreme cases [28]. Emphasise on the generality of the developed formulation to use various material laws, stress rates and on its special application to metal forming problems. Investigate an efficient motion scheme for the A L E method and a numerical algorithm to process the supplementary equations arising from mesh motion. Existing schemes tend to introduce higher degree of complexity in the implementation and require more CPU processing time. Study the effect of various stress integration methods on objectivity and plastic 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. Chapter 2 A R B I T R A R Y L A G R A N G I A N E U L E R I A N (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 Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 19 mesh motion do not occur at the same time and the discretized equilibrium equations do not include simultaneously mesh velocities and material velocities as unknowns. In this regard, Huetink formulation may be considered as an updated Lagrangian formulation (ULF) coupled with remeshing after each incremental step. Benson [35] described a method called simple A L E . It seems that the method was developed mainly for fluid mechanics. The method was introduced to reduce the computer time of large scale analyses involving finite deformation through the use of explicit finite element schemes. The time step size of the explicit program is, therefore, limited by element mesh dimensions. In order to modify the mesh, the A L E method was implemented and the equations were derived by substituting the relationship between the material time derivative and reference configuration time derivative into the governing equations for a continuum in a Lagrangian coordinate system. An operator split was advocated to decouple the equations so that for every increment, two steps were needed. First, a Lagrangian step was performed, in which the mesh moved with the material during the step. The solution was then mapped from the Lagrangian mesh to the reference one in order to complete the A L E step. This method may be also considered as an updated Lagrangian formulation coupled with remeshing after increment. An improved version of A L E formulation was presented by Haber [36]. In Haber's work two displacement variables were considered to be primary unknowns, i.e., the Lagrangian and Eulerian type of displacements. This allowed some distinction between the material and the computational motion. The new technique was applied to large-deformation frictional contact and fracture mechanics. It was demonstrated that the mixed Eulerian-Lagrangian description (i.e., ALE) may be used to vary the crack length in a continuous way and may easily handle contact boundary condition. The paper does not, however, clearly indicate which displacement stands for mesh motion and there is no discussion about how to move the computational finite element mesh. Furthermore, one of the most demanding aspects in A L E , the updating of Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 20 material-associated properties was not addressed. It is also assumed that the formulation is limited to linear elastic material and the load is independent of displacement, therefore, the final formulation may not be directly applicable to non-linear material and deformation-dependent loads. Similar to the above work, Yamada introduced two types of displacements called the material and spatial incremental displacements and derived an A L E formulation for plane deformation of hyperelastic material [37]. For such kind of materials, the stress can be uniquely determined from the strain energy density and it is independent of strain histories. Thus, no state quantity of the particles needs to be introduced except the deformation gradient tensor. The application of formulation to only path-independent hyperelastic material models eliminated the need for rigorously addressing the problem of calculating and updating material associated properties, and makes the formulation not applicable to the general solid mechanics deformation problems. An A L E formulation specifically derived for solid mechanics was described by Schreurs [38, 39]. In [39], the difference between the CRS (Computational Reference System) derivative and the MRS (Material Reference System) derivative of a physical quantity was clearly stated. Starting from the equilibrium equation in material domain, the weak form was set up by the principle of weighted residuals. The final discretized equilibrium equations were obtained by transforming integration over material domain to a reference domain. However, the physical meaning of the reference domain was not clearly addressed. A computational mesh moving method was created by assuming the grid point as a material point of a linear isotropic fictitious body. Assuming one shape of a stress-free element being optimal, linear set of algebraic equations were set up to have the deformed mesh elements recover to the optimal shape. By solving the equations, the grid nodal displacements were determined. This method is time consuming and does not guarantee an optimal mesh because the real material is generally elastic-plastic. Follower points or cells are used to represent the material associated properties and deformation history. These properties at new integration points decided by CRS points Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 21 were updated by extrapolation and interpolations from the follower points. Hughes [40] elaborated on some basic concepts related to A L E . An important one is the relationship between the material time derivative and referential time derivative * / A of a physical quantity * / . These should satisfy the following relation: '/• = 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 Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 22 transport theorem widely used in fluid mechanics was applied to an arbitrary moving control volume to obtain the field equations of mechanics with respect to an arbitrarily moving grid point. The weak form of the differential equations were obtained by using appropriate weighting functions and integrating over the current grid volume. In Ghosh's work it was assumed that the motion of the material with respect to the grid is quasi-static, i.e., | = 0. In the opinion of the author, this assumption is not warranted. In quasi-static problems, it is the material point acceleration, as expressed by Equation (2.2), that may be neglected instead of the motion of the material with respect to the grid points, because the term contributes only partially to the material point acceleration. d*vi dt dt •x 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 Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 23 to run a pure Lagrangian mode, no further action is taken and the code proceeds to the next time step. If a pure Eulerian calculation is desired, the nodes are placed back in their original positions. The user has options available to tailor the evolution of the mesh in order to maximize either efficiency or accuracy. The mesh updating scheme implemented in 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. Chapter 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°( 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 Chapter 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 "•" and a cap " A " , respectively. The relationship between the two derivatives is given by Equation (2.1) (See Appendix A for details). Equation (2.1) is important since it makes it possible to track the material deformation history when A L E is used. Chapter 2. ARBITRARY 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 « " 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 » " V O t + A t XiJ and Chapter 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+&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 ~ ~¥x~ t + A t dSui du k ~ dx = 8+ x t t k At j dividing the above equation by At and let A i — > 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: <f+^V = a*v\ t t + (cfV^At (2.12) Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION 28 The computational reference time rate of volume is [48]: {cfvy and it follows that: <? V +At =<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;£ —' 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 — x%, the following relations may be obtained: l t i dSui d'xj dSu; i <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 —-—=• Ai t v JV + d xi t l t ^AtcfV) pp] O Xj t + j AtcfV OX I t i Chapter 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 „, \ 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 £ - f * ) ^ ) A t (2.17) t may be described as: t+At d* Aj = +M d*Aj +d A >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 Chapter 2. ARBITRARY LAGRANGIAN EULERIAN (ALE) FORMULATION V^- + f £^d%- - ^<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 Chapter 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. Chapter 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¥x-[- ¥x~ + ¥x- 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], « a ~ ~ ¥x~ ~ ¥x~ akj aik k ¥x~ + aij k ( 2 k - 2 3 ) is applied, Equation (2.20) takes the form: dSui ('t_T rp t Iv ¥x~ { r ds ( Ui d Vi d i l ,t_ - + ¥x~ a akj + k t d'v] b^vi ,t ( . tv ¥x~ aik k d'vi = / &u\f?£V + I SulffJS J'V J*S - t dv l t k ¥x- ) aij k t * V .d^A (2.24) Chapter 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 <9^fc\ t( 2 % hv = rik (d'vj -J t [frxk d^kj + frxj + J's \ f cr jfe / &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 Chapter 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. Chapter 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 Chapter 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 Chapter 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 £ + 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]. Chapter 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 Chapter 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 Chapter 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 Chapter 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 — 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) Chapter 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 ' / / — 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 Chapter 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>i(r, 0), (f>i(r, 1), c&(0, s) and & (1, s) as shown schematically in Figure 3.1, the mesh coordinate Cj is given by: d(r,s) = (l-s)<f> (r,0) _ (1 _ ) ( i _ )^(0,0) - (1 - r)s<f>i{0,1) - rs<j>i{l, 1) - r(l-*)&(l,0) i r + s<f> (r,l) + (l-r)(f) (0,s)+r<f> (l,s) i i i a (3.6) 0 < 7" < 1 0 < 5 < 1 where r, s are normalized coordinates, i = x,y, i.e., (j> (r, 0) and <j> (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 Chapter 3. NUMERICAL (0,D W IMPLEMENTATION OF ALE 44 (0,1) r,1) (1.1) <t>i(0,s) 4>,(i.s) (U°#) ) (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°(r, s), then the mesh velocity is given by: Cj(r,s) - c?(r,a) Ai t c* (3.7) (3.8) Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 45 Figure 3.2 gives an example for distorted and automatically modified mesh obtained by applying the above method. It should be noted that, even when the boundaries have discontinuous slopes, for example at A, B and C in Figure 3.2, the method produced a homogeneous mesh. Distorted Mesh at Time t Modified Mesh at Time t+At Figure 3.2: Mesh modification 3.3.2 Consideration of Boundary Motion The above mesh motion scheme will guarantee a homogeneous mesh when boundaries do not move. In a general A L E formulation, however, mesh motion and loading increment occur simultaneously and the boundaries should be allowed to move. In such case, in order to achieve a higher quality mesh, Equation (3.8), has to be modified. First we assume that the boundary velocities are expressed by £;(r, 0),£i(r, l ) , £ ; ( 0 , s ) and s) corresponding to boundaries Chapter 3. NUMERICAL </>i(r, 0),4>i(r, l),4>i(0,s) IMPLEMENTATION OF ALE 46 and </>i(l,s) individually, as shown in Figure 3.1. We use the same function as in Equation (3.6) to interpolate the boundary velocities and find the modification term to be used in Equation (3.8); *vf = (l-s)Ci(r,0) + sU(r,l) + (l-r)Ci(0,s) + r^(l,s) - (1 - r)(l - s)$i(0,0) - (1 - r)sZi{0,1) - ™&(1, 1) - r(l-a)&(l,0) (3.9) 0< r< 1 0< s <1 To ensure that the boundary nodes have to be Lagrangian type in the normal direction, the mesh velocities £;(r, 0), £;(r, 1), £;(0, s) and &(1, s) have to be equal to the material velocities Vi(r, 0)/ Vi(r, l), 1>;(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 — 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< r < 1 &(r,l) = «- 't; (r,l) A i 6(1,«) = * " f i ( l , *) At (3-10) 0 < s < l) If higher quality mesh is required, some iterations may be introduced, i.e., using the new material velocities to replace the velocities at time (t—At) in Equation (3.10) until the difference between the material velocities from the two iterations is less than a given tolerance value. In practical applications, this iteration was found to be not necessary since the mesh quality created by Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 47 Equation (3.10) is good enough for most of practical problems. Finally the mesh motion for the nodes in the interior of the domains is given by: t C t 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 Chapter 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>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 — 0 and 6/ = 1 in Equation (3.3). This may result in unexpected or inhomogeneous spacing because of the deformation. In order to have homogeneous mesh pattern, the boundary nodes have to be adjusted in the tangential direction of the boundary. This kind of adjustment can not be done during the incremental step and has to be done after the deformed boundary shape is determined. In this work, the nodes are relocated on the boundaries according to nodal spacing specified by user or are positioned to get an equal distance between \ Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 49 each node. 3.3.4 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 — 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 Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 50 method directly applicable to domain with irregular polygon boundaries. In transfinite mapping method, the interior node locations are calculated from boundary nodes. The curvatures of the boundary curves and spacing of boundary nodes are accurately reflected in the mesh, as may be shown in Figure 3.3. Furthermore, the boundary curves do not need to be simple functions. Cusps and inflection points provide no special problems. The node motion scheme used in the research for the second type of boundary is simple and more accurate than assuming a normal direction. It is not limited by boundary curve types and it may be combined with transfinite mapping method to handle the mesh motion in the whole deformation domain. Node adjustment in tangential direction on the second kind boundaries is necessary to keep a higher quality mesh. Stiffness matrix processing discussed above is a general method. That may be used with various material models, i.e., rate dependent or independent. 3.4 I N T E G R A T I O N 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 Chapter 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 *• dx t 1 r where l F +At = g^ at i:j <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) Chapter 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<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. Chapter 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 — 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„ m d m t l k x „ 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„.. _ 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 — l 0x m d Xj du t d Xi t 0~mn ~7^1 ~f~ ft x ox l d xj t "777 A • + 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, <r t+At - - n-t At _ I — T T At d u* 7^ Vd<z t_ I (?mi + m j m t_ I t du ^ "ftt, , /-» o du P. d x, m — cr. N k + CijuDu t which, upon taking the limit when At — > 0, gives: i dx l ^ " j t „ dx l k k t n ki (3.17) Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 54 where 'u; is the velocity component in the ith direction of the coordinates. By definition, the Truesdell rate of Cauchy stress may be given by Equation (2.23). Therefore, we may write Equation (3.17) in the following form: (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 ° 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) Chapter 3. NUMERICAL IMPLEMENTATION t+At t+aAt" d t+A F1 55 X{ Qt+cAt »J 7. OF ALE Xj t n_I( dtyi , when a = 0, it will be explicit and may be expressed as: t + A t _ i &+*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 ^-—^ (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 ^ > ™ 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<T;J = Cijkil+ox^D At kl The l+o.s&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) Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 56 and t+0.5At If displacement Ui = X{ t+At — 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 Chapter 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 Chapter 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) — = txisin^At) + x cos(u;At) x sin(uiAt) t 2 t (3.23) 2 so that, the displacements are given by: ui — *a;i(cos(a;Ai) — 1) — x sin(u>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 —¥ 0, gives the velocities at time t as: t t t V\ = —UJ 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: '"» = | £ = ° • D a i (p. + £ • ) o <- > 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. Chapter 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) — 1) — x sin(u>At) u\ t V l At u = t t 2 * = At Al = x sin(u}At) + x (cos(u> At) — 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>w — - ° 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 Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 60 rotation u;At or, in general, be a function of the components of the spin tensor 'W. In general large deformation , the skew-symmetric spin tensor is defined by; , 3 - 2 8 ) and the modification term for 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 <D = 0 22 4 /J>i2 = 0 Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 61 Therefore, the explicit integration algorithm, Equation (3.16), will satisfy incremental objectivity if is modified as shown above. The *Dij modification presented above may be considered as an alternative way to keep the incremental objectivity of integration algorithms developed by Pinsky [58] and Hughes [56]. The above modification is practically useful when explicit integration algorithms are utilized. More detailed numerical examples on'the effectiveness of the above procedures is given in Chapter 5. 3.6 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 Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 62 conditions. In the case of large deformation, although the constitutive relation is not the same as for small deformation, the elastic predictor may still be calculated by assuming that the deformation in the current increment is completely elastic. Then the predictor is mapped onto a suitably updated yield surface using the appropriate large strain constitutive relation, thus restoring plastic consistency. For von-Mises material the equivalent stress at time t is defined by: = where, a' - = crij — | £ ; / 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 ^ ^ £ 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 '<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 Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 63 (iii) calculate the elastic stress increments A<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>;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 — 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 <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 Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 64 increments include the incremental stress due to proper transformation, i.e. taking into account of the change in configuration. The return mapping procedure has to start with 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,/ + --«0 (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) Chapter 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 > 3 6 ) But at time i , *x = *x, so that * / = * / = / ( x , i ) . Therefore, Equations (3.34)-(3.36) may c 4 have the form: » r u = ' f + A , / + ^ K - « . ) 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. Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 66 In this method, assuming the smoothed property / „ 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> = 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<f> = 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 < 1 2 fin 2 . (3.42) Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 67 where, / / , fn, fm and fiv are the unsmoothed values at the Gauss points. For general element shapes and shape functions, solving Equation (3.41) may be costly. To overcome this problem, discrete smoothing is introduced as described in the following. The function in Equation (3.39) is replaced by: </> = E ( ( W . ) \K - f \ ) 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) < (h^hi) [ ^K=I \K ••• ^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 Chapter 3. NUMERICAL IMPLEMENTATION OF ALE 68 - Add up the values from different elements sharing the same node and divided by the number of elements sharing the node to get the final property values at the mesh node. Chapter 4 TWO-DIMENSIONAL 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 — bj = 0 for all I, the program reduces to Eulerian one; If aj = 0, 6/ = 1 for all I, the program reduces to updated Lagrangian one; otherwise, the program will use a general A L E formulation. - The program incorporates both Jaumann stress rate and Truesdell stress rate in stiffness calculation. The corresponding stress integration algorithms are employed in computing stress. 69 Chapter 4. TWO-DIMENSIONAL ALE PROGRAM (ALEFE) 70 - For general A L E case, the mesh motion for interior nodes is controlled by the nodal motion on boundaries. User only needs to choose the aj and 6j values in Equation (4.1) for boundary nodes. The internal nodes are specified by the program automatically to keep a homogeneous mesh. This is a particularly important features for users and for practical use of the program. It is relatively easy and straight forward to identify the formulation type required, and hence the coefficient a/, bi, for the boundary nodes. - The program handles unsymmetric stiffness matrices and coupled displacements/velocity . boundary conditions. The direct methods for coupled D.O.F. is applied so that no approximation is introduced in this aspect. - Convergence criteria include displacement, residual force and iterative energy. User may choose one, two or all criteria for a given problem. - Rate dependent and rate independent material models are included in the program. For plastic deformation, multi-linear hardening models can be chosen. - Type of elements available is four node quadrilateral element. - The input data format is designed similar to available commercial finite element codes, so that the data generation phase may be directly imported from these programs. - The output data format is designed to be compatible with general graphic simulation and data processing commercial softwares, so that contour, x-y and deformed mesh plots may be easily created from the output data file of A L E F E . - Full and modified Newton-Raphson methods are available. Chapter 4. 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 <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 <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» = u{ + V{At t 4. Use and u\ — u\ + <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 <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. hapter 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 —I 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 Chapter 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£ and the updated stress variables to calculate the stress increment A<T;J, <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, <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>;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<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 . Chapter 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. Chapter 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. Chapter 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 er 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 —f 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 Chapter 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 •+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<f S t+At +At Chapter 4. TWO-DIMENSIONAL ALE PROGRAM (ALEFE) 79 Therefore, the residual force at D.O.F. / can be expressed as: t+At n = t+At„ t+At„ 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 < 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 < 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 Chapter 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 < 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 Chapter 4. TWO-DIMENSIONAL ALE PROGRAM (ALEFE) 81 the following relation between local and global displacement at node a: cos9 Vx(a) y(<*) ) J sinO —sinQ 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 —sinO cosO = T («) (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. er 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 Chapter 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 Chapter 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 = —1.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>. 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 Chapter 5. NUMERICAL EXAMPLES 85 with the analytical values, even for large rotation of up to 90 degrees. The conventional forward difference algorithm predicts, however, incorrect values for stress 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>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 — 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 er 5. NUMERICAL EXAMPLES o (Pa) x 1 .OOOOE8 — --O 1 ^ 8.0000E7 ANALYTICAL SOLUTION — MODIFIED F O R W A R D — — CENTRAL DIFFERENCE FORWARD DIFFERENCE \ \ \ 6.0000E7 \ \ \ 4.0000E7 - \ \ \ 2.0000E7 \ O.OOOOEO - 1 _\ -5.1200E2 \ -1 . O O O O E 9 •2.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 Chapter 5. NUMERICAL EXAMPLES 87 quadrilateral element. To compare the results, the exact same problem is also simulated by NISA and ANSYS. Figure 5.3: Elastic-plastic small deflection model Figure 5.4 is the curve of tensile load versus displacement in x direction 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. Chapter 5. NUMERICAL EXAMPLES Figure 5.4: Load versus displacement in x direction Chapter 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 Chapter 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°, 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°. 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 —2.6036. The negative Jacobian is due to an interior angle larger than 180°. The only limitation on the shape of the real element is that each interior corner angle must be less than 180°, i.e., the element must be Chapter 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 Chapter 5. NUMERICAL EXAMPLES Figure 5.7: Pure bending of a plate: Lagrangian Formulation Results 92 Chapter 5. NUMERICAL EXAMPLES 93 convex [68]. Practically, however, to achieve higher accuracy, the optimum elements should have interior corner angles of quadrilaterals not far from 90° [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 Chapter 5. NUMERICAL EXAMPLES 94 a. Displacement in x direction 2500 r - 0.000 0.010 0.020 0.030 0.040 0.050 0.060 0.070 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 Chapter 5. NUMERICAL EXAMPLES 95 and inhomogeneous meshes are used in simulation. The material properties, geometry and boundary condition have been described in Chapter 1. The problem is simulated by the developed program A L E F E and commercial codes NISA [10], ANSYS [9] and DEFORM2 [8]. Analysis is performed for up to 60% reduction in the height of the original piece. The homogeneous mesh is shown in Figure 1.4 and is composed of square elements in the undeformed geometry. The results from NISA and ANSYS are obtained using updated Lagrangian method and are shown in Figure 1.5. For proper modeling of this problem, the material point under the punch corner should be free once it moves out of the corner. However, the Lagrangian formulation can not update the boundary condition automatically. This is an important feature in handling most of metal forming problems. This deficiency has the effect of increasing the punch size during the deformation, as shown in Figure 1.5. It is also noted that elements at the punch corner are highly distorted because of the large deformation. Increasing the deformation level causes more severe element distortions and the program eventually terminates the analysis. The deformed shape from DEFORM2 is shown in Figure 1.7. In order to achieve 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. Chapter 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 Chapter 5. NUMERICAL EXAMPLES Figure 5.9: The current yield stress from A L E Chapter 5. NUMERICAL EXAMPLES 98 vertical direction to be Lagrangian while freezing the degree of freedom in the horizontal 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 Chapter 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. Chapter 5. NUMERICAL EXAMPLES Oyield (Pa) - 1.2799E9 1.2114E9 -- 1.1430E9 1.0745E9 1.0060E9 9.3761 E8 — 8.6914E8 8.0068E8 7.3222E8 - 6.6375E8 5.9529E8 - 5.2682E8 m 4.5836E8 3.8989E8 3.2143E8 • m - • 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 Chapter 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— 1 0 t\ 30 mm B X Figure 5.13: Inhomogeneous mesh for punch forging simulation By comparing the stress and displacement distributions from different programs, it is found that although NISA gives similar displacements from homogeneous and inhomogeneous meshes, stress distributions are much different from the two meshes with stresses from the inhomogeneous mesh 2 to 10 times those obtained from homogeneous mesh. As an example, the current yield stresses are illustrated in Figure 5.14, in which the difference is about 5 times. Chapter 5. NUMERICAL EXAMPLES 102 Table 5.1: Displacements of top right corner from all simulations Program 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. Chapter 5. NUMERICAL EXAMPLES 103 a yield ("1.0E8 Pa) a. The result from homogeneous mesh Figure 5.14: Comparison of current yield stresses from NISA Chapter 5. NUMERICAL EXAMPLES Figure 5.15: The A L E result of current yield stress from inhomogeneous mesh Chapter 5. NUMERICAL 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. Chapter 5. NUMERICAL EXAMPLES Figure 5.16: The finite element model for tensile necking er5. NUMERICAL EXAMPLES Figure 5.17: The necking S and engineering straining 7 Chapter 5. NUMERICAL EXAMPLES 108 It can be seen from Figure 5.17 that for integration of Jaumann rate, the transformation algorithm A and B lead to same results. For Truesdell rate, transformation-A and B give quite different results. This indicates that Truesdell rate is more sensitive to the location of stress transformation. It should be mentioned here that the results from Truesdell A scheme have almost same slopes and very close values to those reported by McMeeking [3]. The results from ANSYS show very much different values of necking. The ANSYS run was terminated at 7 = 0.57 because of divergence, whereas all the runs of the developed A L E program have been successful for values up to 7 = 1.0. Figure 5.18 shows results from different schemes developed by Pinsky et al. [58] for the integration of Truesdell rate, Equation (3.20) in this thesis, (curves: Truesdell-A, P and Truesdell-B, P). Figure 5.18 also shows corresponding results developed in this research, Equation (3.15), ( curves: Truesdell-A and Truesdell-B). The results indicate that schemes developed in this research give same results as those developed by Pinsky et al. [58]. This is expected since the difference between Equation (3.15) and Equation (3.20) is only in the second term on right hand side and is caused only by higher order terms of A i . Figure 5.19 shows the contours of the equivalent von-Mises stress and the ratio of the equivalent stress to current yield stress at 7 = 0.75. Although no stress distribution is reported in [3] , the shape of unloading zone is discussed. From Figure 5.19, the zone where the ratio of the equivalent stress to current yield stress is less than 1.0 is experiencing unloading. The unloading zone shown here has the same shape and location as that reported in [3]. The example indicates that the transformation algorithm A and B may generally be combined with any integration algorithm. The integration of Jaumann rate is less sensitive to the location of stress transformation. The combination of Truesdell rate with 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 Chapter 5. NUMERICAL EXAMPLES 109 Chapter 5. NUMERICAL EXAMPLES Figure 5.19: Contours of equivalent stress (Truesdell rate, transformation Chapter 5. NUMERICAL EXAMPLES 111 problem by A L E method. 5.7 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 — 0.3. The example is the same as the one discussed in Chapter 1, except that the die is shaped in the form of 5th order polynomial curve with zero curvature and slope at the ends. The same process was simulated with updated Lagrangian method by Lee [12] and by Yamada [70], where the billet was assumed to have been initially machined to fit the die. Figure 5.20 shows the geometry and the finite element mesh used in A L E simulation. In the A L E simulation, the mesh and the corresponding nodes confined to the die area are taken Chapter 5. NUMERICAL EXAMPLES Figure 5.20: Original geometry and mesh in extrusion through a curved die Chapter 5. NUMERICAL 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 Chapter 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° 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 <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 er 5. NUMERICAL EXAMPLES a. Result from A L E Figure 5.22: Comparison of o~ from A L E and reference [12] x Chapter 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] Chapter 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 • _L_ 3 m^f ii—L • <[ 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 Chapter 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 Chapter 5. NUMERICAL EXAMPLES Initial m e s h Distorted m e s h Figure 5.26: The initial and deformed mesh from NISA Chapter 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 Chapter 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. Chapter 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 Chapter 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 Chapters. 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 £ yieW eq p £ (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) Chapter 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 Chapter 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 = —r— (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 Chapter 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 Chapter 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 £ '-a 8 —f "" --1- t y 10 ..f.-a —f- O) c - 5 £ cb 4 h ^ o ^o- • CASE1.A - - O - - CASE1 ,B • CASE2 --<!> — o CASE3 A CASE5 are shown in Table:5.2 , i , 200 CASE4 O ' Parameters lor different cases —-f— 100 V 300 ' ' • ' • 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 Chapter 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 Chapter 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 Chapter 5. NUMERICAL EXAMPLES 131 method is more reliable and more accurate than the ones previously discussed. —l 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, Chapter 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 Chapter 5. NUMERICAL EXAMPLES 133 in good qualitative agreement with the experimental results. The discrepancies at the entrance and exit points may be attributed to inaccurate description of friction, the roll flattening and the spring back of work-piece after exit. The experimental results show significant increase in contact angle at both entrance and exit, while the contact angle in the numerical simulation is smaller, obtained only from pure geometric relations and has the same value as in a classical slab method. 350 i~ -0.02 0.00 0.02 0.04 0.06 0.08 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 Chapter 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 Chapter 5. NUMERICAL EXAMPLES 135 2500 r- 2000 , 1000 ...^^Or. I , 500 O EXPERIMENT • - • - - ALE-BOU •—* 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 Chapter 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 Chapter 5. NUMERICAL EXAMPLES 137 'yield < > 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 °yield (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, Chapter 5. NUMERICAL EXAMPLES 138 work is simple and proved to be effective in conjunction with the Eulerian simulation. The computed results showed both double-peak as well as "friction hill" distributions for rolling pressure. In general, the results are in good agreement with experimental ones. 5.9 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° 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 Chapter 5. NUMERICAL EXAMPLES y Initial shape Deformed shape Figure 5.41: General shape of the specimen and die Table 5.3: Specimen dimensions in compression between wedge-shaped dies No 1 2 3 4 5 6 7 8 m in degree mm 3° 3° 3° 25.40 22.23 19.05 1.000 0.875 0.750 3° 3° 3° 3° 3° 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 Chapter 5. NUMERICAL EXAMPLES 140 pressure gauges, so that the die load may be recorded. Dial gauges were used for measuring downward displacement of the die and lateral movements of the specimen. Therefore, the deformed specimen thickness and velocities at two ends can be measured. The rough dies were prepared by machining a square grid on the surface and using no lubrication in compression [78]. The material of specimen chosen in A L E simulation is half hard aluminium. The 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 Chapter 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° 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 Chapter 5. NUMERICAL EXAMPLES 142 field. 5.9.4 Comparison of Results The results predicted by A L E method are compared with the experimental and slip line field analysis results published by Chitkara and Johnson [78]. Figure 5.43 shows the result of Case 1 giving the variation of load on the die with respect to compressive strain In 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 Chapter 5. NUMERICAL EXAMPLES Figure 5.43: The variation of vertical load on the die with compressive strain Chapter 5. NUMERICAL EXAMPLES 144 Figure 5.44: The variation of mean vertical die pressure with current strip thickness Chapter 5. NUMERICAL EXAMPLES 145 Chapter 5. NUMERICAL EXAMPLES 146 the plastic zones are initially formed at the die corners and the strain values are highest through the whole deformation. The central parts of the interfaces have lower values of plastic strain, which means that the deformation is small and is incurred by the non-slip friction condition. This may be the reason why in slip line fields, a rigid body block is usually assumed there. It is interesting to indicate that plastic zones are already formed outside of the dies, but in slip line field constructed by Chitkara and Johnson [78], the plastic deformation begins only inside the dies, which may contribute to the inaccuracy of slip line results. Plastic zones are mainly located between and around the two dies. For most area outside the dies, plastic deformation is almost zero. er 5. NUMERICAL EXAMPLES 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 Chapter 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 Chapter 6. CONCLUSIONS AND FUTURE WORK 149 efficient and is straight forward to implement in an A L E formulation. In the thesis, motion of nodes on boundaries, the classification of boundary types and the corresponding detailed motion schemes are presented and shown to be general and efficient in obtaining a high quality mesh. The effect of distorted mesh on the accuracy of the results and the necessity for an efficient A L E mesh motion scheme are demonstrated through numerical examples. It is shown that the developed mesh motion scheme maintains a homogeneous mesh under very large deformation levels. The 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. Chapter 6. CONCLUSIONS AND FUTURE WORK 150 In the updating of material associated properties, the relation between material derivative and referential time derivative is applied to create the relation between increment due to mesh motion and increment due to material motion. Using the developed relation, the material associated properties may be updated directly from material points to mesh points without the need for time consuming interpolations and the creation of a hypothetical mesh as was done by other authors. Combined with local smoothing by the least square method, the updating scheme provides an efficient and accurate routine to update material associated properties. Practical application of the developed A L E formulation is presented in the metal forming area. Simulation of metal extrusion and punch forging processes is performed. It is shown that the load fluctuation and punch size increasement, normally associated with updated Lagrangian formulation, are eliminated. A L E results indicate that the mesh motion algorithm maintains homogeneous mesh throughout the deformation process. It is also shown that description of changing contact boundary conditions is straight forward and easy to perform in A L E formulation. A necking bifurcation of a bar in tension is simulated by an updated Lagrangian method that is reduced from the A L E formulation. Through this example, it is shown that the combination of Truesdell rate with transformation after return mapping gives consistently better results than use of Jaumann rate. One of the most important and complicated metal forming processes is rolling. The steady state condition of rolling of sheet metals is simulated with the Eulerian option, i.e. the particular case of A L E method when the mesh does not move. Application of correct boundary conditions and simulation of neutral points are investigated. The simulation is performed with the general 2D finite element program A L E F E developed in this research without the need to append special routines for the application. No particular interface elements or friction layers are needed. For different ratio of initial thickness of work piece and roll radius cases, the program predicts both double and single peaks of friction-hill of pressure distribution, as indicated by experiments. Rolling pressure, force and torques are in reasonable agreement with the experimental results. Chapter 6. CONCLUSIONS AND FUTURE WORK 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 Chapter 6. accuracy. CONCLUSIONS AND FUTURE WORK Bibliography [1] Hill, R., "Some Basic Principles in the Mechanics of Solids without a Natural Time", J . Mech. Phy. Solids, Vol.7, pp. 209-225, 1959. [2] Hibbitt, H . D . , Marcal, P. V . , Rice, J . R., "A Finite Element Formulation for Problems of Large Strain and Large Displacement", Int. J . Solids Struct., Vol.6, pp. 1069-1086, 1970. [3] McMeeking, R. M . , Rice, J . R., "Finite Element Formulations for Problems of Large Elastic-plastic Deformation", Int. J . Solids Struct., Vol.11, pp. 601-616, 1975. [4] Maniatty, A . M . , Dawson, P. R., Weber, G . G . , "An Eulerian Elasto-viscoplastic Formulation for Steady-state Forming Process", Int. J . Mech. Sci., Vol. 33, pp. 361-377, 1991. [5] Gadala, M . S., Dokainish, M . A . , Oravas, G . A E . , "Formulation Methods of 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 Bibliography 154 [12] Lee, E . H . , Mallett, R. L . , "Stress and Deformation of Metal Extrusion Process", Comp. Mech. Appl. Mech Engrg., Vol.10, pp. 339-353, 1977. [13] Voyiadjis, C . Z., Foroozesh, M.,"A Finite Strain Total Lagrangian Finite Element Solution for Metal Extrusion Problems", Comp. Appl. Mech. Engrg., Vol.86, pp. 337-370, 1991. [14] Lush, A . M . , "Numerical Grid Generation Used for Remeshing Finite Element 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 & Structures, Vol.20, pp. 745-758, 1985. [21] Abo-Elkhier, M . , On the Numerical Solution of Nonlinear Problems in Continuum Mechanics, Ph. D. thesis, McMaster University, Hamilton, Canada, 1985. [22] Derbalian, K . A . , Lee, E . H . , Mallet, R. L . , McMeeking, R. M , "Finite Element Metal Forming Analysis with Spatially Fixed Mesh", In: Application of Numerical Methods for Forming Process. (Armen, H . eds.) A S M E - A M D Vol.28, pp. 39-47, 1978. Bibliography 155 [23] Oden, J . T . , Braucrdi, H . J . , " On the Calculation of Consistent Stress Distributions in Finite Element Approximations", Int. J . Nume. Mech. Engrg., Vol.3, pp. 317-325, 1971. [24] Ghosh, S., Kikuchi, N.," An Arbitrary Lagrangian-Eulerian Finite Element Method for Large Deformation Analysis of Elastic-viscoplastic Solids", Comp. Meth. Appl. Mech. Engrg., Vol.86, pp. 127-188, 1991. [25] Cough, R., Sharp, R., Otero, I., Tipon, R., McCallen. R., "Application of A L E Techniques to Metal Forming Simulations", Advanced Computational Methods for Materials Modeling, American Society of Mechanical Engineers, Applied Mechanics Division, A M D Vol.180, pp. 133-140, 1993. [26] Wang, J . , Gadala, M . S., "Aspects in the Formulation and Application of Finite Elements to Metal Forming Problems", Proceedings of the International Conference, Modelling and Simulation in Metallurgical Engineering and Material Science, pp. 577-582, Beijing, P. R. China, 1996. [27] Noh, W . F . , "A Time-dependent Two-space-dimensional Coupled 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. Bibliography 156 [34] Huetink, J . , Vreede, P. T . , Lugt, J . V . D . , "Progress in Mixed Eulerian-Lagrangian Finite Element Simulation of Forming Process", Int. J . Nume. Meth Engng., Vol.30, pp. 1441-1457, 1990. [35] Benson, D. J . , "An Efficient, Accurate, Simple A L E Method for Nonlinear Finite Element Programs", Comp. Meth. Appl. Mech. Engrg., Vol. 72, pp. 305-350, 1989. [36] Haber, R. B . , "A Mixed Eulerian-Lagrangian Displacement Model for 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. Bibliography 157 [45] Huerta, A . , Casadei, F ., "New A L E Applications in Non-linear Fast-transient Solid Dynamic", Engineering Computations, Vol.11, pp.317-345, 1994. [46] Bathe, K . J . , Ramm, E . , Wilson, E . L . , "Finite Element Formulation for Large deformation Dynamic Analysis", Int. J . Nume. Meth. Engrg., Vol.9, pp. 353-386, 1975. [47] Bathe, K . , Finite Element in Engineering Analysis, Prentice-Hall Inc., New Jersey, 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. Bibliography 158 [57] Hibbitt, Karlsson and Sorensen, Inc., A B A Q U S Theory Manual, Hibbit, Karlsson & 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 & 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. Bibliography 159 Yamada, Y . , Hirakawa, FL, "Large Deformation and Instability Analysis in Metal Forming Processing", Application of Numerical Methods to forming Processes, A S M E , AMD-Vol.28, (Edited by H . Armen), pp.27-38, 1978. Farmer, L . E . , Oxley, P. L . B., "A Slip Field for Plane Strain Extrusion of a 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. Appendix 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 + &X (A.5) t X 160 3 Appendix A. Material and Referential Time Rate 161 Figure A . l : The description of relation between mesh motion and material motion Appendix A. Material and Referential Time Rate 162 Assuming that the MRS point located at R at time t + At comes from a material point P*( 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)» ( 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 * — 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% — 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 — > 0, 7 A = f + ^Cvl f - *v ) k (A.9) Appendix 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„,(, v = j f -8v d V t tv B i +f t i -||^))^ W tfS-SvutS (B.10) where ' A j = \ (§r^ + 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 °'ij + ~ J <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 < (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 &Xi 1 ( k dSv \ dSv 2 ^ 0*3; t + ld v l t a«aj J k a i (&v, , d*v \ (t-J lv¥x~A --^{o^ k T + k j {d^ *a jk + d'vA d8v d*x ) d*Xi k ( 5 ^ o^J ^{d^ -¥x-Jr v + k 163 kt f c \ \ dW\j ^frxj - (B 12) Appendix 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 *«• 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 < 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). Appendix 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 — At, assume '~ a:; =*At from t — 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 Appendix D Stress Increment Procedure We consider decomposing the general deformation into two steps; straining and rigid body motion. In the first step, due to straining, the Cauchy stress would have an increment of Ao-ij = dju'DuAt (D.18) where, Ciju is the element constitutive tensor, At is the time increment and *Dki is the rate of deformation tensor. The stress after the first step would be: *l? = V « + Aery = \+**S™ 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™ (D.21) Substituting Equation (D.21) and (D.19) into Equation (D.20), we get: t+At„ _ 1 0 1 •*- <+*'F d w t + A t x d i t ~lt„ d*x m " t + A t *± 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
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Arbitrary Lagrangian-Eulerian method and its application...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Arbitrary Lagrangian-Eulerian method and its application in solid mechanics Wang, Jin 1998
pdf
Page Metadata
Item Metadata
Title | Arbitrary Lagrangian-Eulerian method and its application in solid mechanics |
Creator |
Wang, Jin |
Date | 1998 |
Date Issued | 2009-06-02T19:42:10Z |
Description | 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 ALE formulation is developed from the virtual work equation transformed to an arbitrary computational reference configuration. The formulation presents a general approach to ALE 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 ALE 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 treatment for rate of deformation tensor [sup t]D[sub ij] is presented to maintain incremental objectivity of the tensor. It is shown by numerical examples that the use of Truesdell stress rate with a developed numerical integration procedure gives consistently more accurate results than other procedures presented. An algorithm for updating material associated properties is presented and applied in simulation of various metal forming problems. A 2D finite element program, ALEFE, based on the presented formulation is developed and tested. The program may reduce to an updated Lagrangian or Eulerian methods as special cases. The mesh motion for the whole domain is controlled by the motion of the boundary nodes. The program can handle unsymmetric stiffness matrices and coupled displacements/velocity boundary conditions. The input data is designed to be similar to available commercial finite element codes, so that the data generation phase may be directly imported from these programs. The output data format is designed to be compatible with general graphic simulation and data processing commercial softwares, so that contour, x-y and deformed mesh plots may be easily created from the output data of ALEFE. Various benchmark and practical problems are simulated by the developed program. 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. |
Extent | 8790981 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
Date Available | 2009-06-02 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080928 |
URI | http://hdl.handle.net/2429/8616 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1998-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- ubc_1998-272672.pdf [ 8.38MB ]
- Metadata
- JSON: 1.0080928.json
- JSON-LD: 1.0080928+ld.json
- RDF/XML (Pretty): 1.0080928.xml
- RDF/JSON: 1.0080928+rdf.json
- Turtle: 1.0080928+rdf-turtle.txt
- N-Triples: 1.0080928+rdf-ntriples.txt
- Original Record: 1.0080928 +original-record.json
- Full Text
- 1.0080928.txt
- Citation
- 1.0080928.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
India | 46 | 4 |
China | 41 | 4 |
United States | 30 | 1 |
Germany | 19 | 4 |
Canada | 12 | 0 |
United Kingdom | 10 | 2 |
Russia | 8 | 0 |
Brazil | 6 | 1 |
France | 6 | 1 |
Italy | 5 | 1 |
Poland | 4 | 0 |
South Africa | 4 | 0 |
Netherlands | 4 | 2 |
City | Views | Downloads |
---|---|---|
Unknown | 56 | 12 |
Shenzhen | 23 | 3 |
Bangalore | 13 | 2 |
Beijing | 9 | 1 |
Kharagpur | 9 | 0 |
Montreal | 8 | 0 |
Ashburn | 7 | 0 |
Chennai | 4 | 0 |
Shanghai | 4 | 0 |
Atlanta | 4 | 0 |
Pretoria | 4 | 0 |
Ranchi | 4 | 0 |
Mountain View | 4 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080928/manifest