.23.“ :3: “.7. Em, A. i936. . v .4 . .rry... :li..:t...:..( 1.1.5.5.. 13...... r... 5. AL... I. a nah“. "u. . .1.st An... I.“ .3. .‘.\.1. 4 . 1:13.): .31.... 3.3.. .. x as... 7 till. : L. .5! :a. 111...! I 9., u! vi: 3.7.1:: v éiiw 1:: . .v . :17... .9 E33. .3. .1 .12: 4:1. 3.. ‘ a....H.r.€. uwbyuwun. v i... .x. . v \ . . .. pity.) at...» D xx: txs-z 3 5:... s9... (xxuv 1 a! n 2 . 929.4... :vl‘iw yr. f if :1, 2:. : .Vir’\.9 1 1.5.3,}.- r :rv 3.. 3 . . . :7. 3 .m . 3:: .11.! v. c.) . ~ 7. may? .2. w I33»... .5‘5 I .1 iii . V . $6.5; _....):: . 25a... . ..............¢.r».. ._. v x . .. _ i .v haul. urn.“ u’ 3' L. ‘ .. h 4 .C a :31; A \IEit. a y :31 251. .Cv-S. {1 $3.}: :2... . .r v.1: a, $2.92: .31.. z t .13.. .6 15%;: . Ln‘l~. {.1 in 1.. I: . IX... .703: )9 WOxI. 3...; lb‘fibii 53.5.. I. 9. [3| it E x .1... z \ 5.5:»..- .- 3! Q; 3. n it}... I 3.... 1...... \n. : .arxv szuJi L» .. . \~ Pfivzv.’ 239‘». 1'0} . 3 n 1. vi $.10 9.1-}. i “8635 .LIBRARY Michigan State University This is to certify that the thesis entitled Arbitrary Lagrangian Eulerian (ALE) Description for Large Deformation Metal Forming Problems presented by Muhammad Mukarrum Raheel has been accepted towards fulfillment of the requirements for MasteLaLSciencdegree in Mechanical. Engineering (Kw I! V . aJor professor Date 20 MIN 2.0m 0-7639 MS U is an Affirmative Action/Equal Opportunity Institution PLACE IN RETURN Box to remove this checkout from your record. To AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE iii“ 0. S," 9 11007 0'6 2’: i3 9.? 6/01 cJCIFiC/DateDue.p65-p.15 ARBITRARY LAGRANGIAN EULERIAN (ALE) DESCRIPTION FOR LARGE DEFORMATION METAL FORNIIN G PROBLEMS By Muhammad Mukarrum Raheel A THESIS Submitted to Michigan State University In partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Mechanical Engineering 2001 ABSTRACT ARBITRARY LAGRANGIAN EULERIAN (ALE) DESCRIPTION FOR LARGE- DEFORMATION, METAL-FORMING PROBLEMS By Muhammad Mukarrum Raheel Arbitrary Lagrangian Eulerian (ALE) description has been used to simulate large- deforrnation, metal-forming problems. Certain mesh governing algorithms are implemented in ALE formulations to control mesh convection. A punch problem is simulated and advantages of ALE simulations over pure Lagrangian simulations are shown. Traditionally, Lagrangian and Eulerian kinematic descriptions are used for finite-element analysis of large-deformation problems. In the Lagrangian description, the grid points adhere to the material points throughout the course of deformation. Lagrangian description has certain limitations in large-deformation problems. In the Eulerian description the reference system is kept fixed in space, and every spatial point is identified by an invariable set of three independent co—ordinates. Because of excessive movement of interfaces and boundaries, it is arduous to take material associated boundary conditions into account. This research pursues Arbitrary Lagrangian Eulerian (ALE) description which combines the advantages of both Lagrangian and Eulerian descriptions. ALE introduces an additional (reference) mapping, which brings additional flexibility in the formulation so that mesh neither adheres to the material nor the space; rather the mesh is in relative motion with the material and it rearranges itself during the deformation process. Dedicated to the Late Dr. Dinesh Balagangadhar iii ACKNOWLEDGEMENT I am grateful to my thesis advisor the late Dr. Dinesh Balagangadhar for his valuable academic and non-academic assistance, continuous guidance, and motivation. He helped me a great deal in all the formulations and simulations discussed in the thesis. His support and encouragement was key to the completion of this work. My interaction with him gave me, not only a high level of academic knowledge, but also provided me a chance to work with such a wonderful person. I would like to thank him for his technical and financial help. My committee members comprising of Dr. Diaz, Dr. Benard, and Dr. Pourbograt deserve special thanks for suggesting and encouraging me after the untimely demise of Dr. Dinesh. Moreover, it is my duty to mention the contribution of all my instructors in the National University of Sciences and Technology (NUS'I‘), Pakistan who taught and encouraged me to reach this point. iv TABLE OF CONTENTS LIST OF TABLES .................................................................................. vii LIST OF FIGURES ................................................................................ viii CHAPTER 1 .......................................................................................... 1 INTRODUCTION .................................................................................... 1 Numerical simulations ...................................................................... 2 Scope of the thesis ........................................................................... 3 Outline of thesis .............................................................................. 4 CHAPTER 2 .......................................................................................... 6 KINEMATIC DESCRIPTIONS ................................................................... 6 Lagrangian and Eulerian descriptions .................................................... 6 Arbitrary Lagrangian Eulerian descriptions ............................................. 9 CHAPTER 3 ......................................................................................... 11 CONTINUUM MECHANICS FOR ALE KINEMATIC DESCRIPTION ................ 11 Co-ordinate systems ....................................................................... 11 Displacement, velocity and acceleration fields ........................................ 13 Transformation equations ................................................................. 14 Conservation of mass and momentum in reference frame ............................ 14 CHAPTER 4 .......................................................................................... 19 KINEMATIC FORMULATION FOR LARGE DEFORMATION PROBLEMS. . . . 19 Constitutive relations for large deformation Hyperelasticity ........................ 2O Variational description for large deformation .......................................... 21 Lagrangian formulation ................................................................... 24 Arbitrary Lagrangian Eulerian formulation ............................................. 26 Numerical implementation ............................................................... 36 CHAPTERS .......................................................................................... 39 PUNCH INDENTATION PROBLEM-A NUMERICAL EXAMPLE ...................... 39 Boundary conditions ....................................................................... 41 Results and discussion ..................................................................... 42 Energy modes ............................................................................... 50 Applications ................................................................................. 54 CHAPTER6 .......................................................................................... 57 CONCLUSIONS AND RECOMMENDATIONS ............................................. 57 APPENDD( A ....................................................................................... 60 MATRIX DEFINITIONS .......................................................................... 60 APPENDIX B ....................................................................................... 63 IMPORTANT RESULTS USED IN FINITE ELEMENT DISCRETIZATION ........... 63 REFERENCES ...................................................................................... 67 vi LIST OF TABLES Table 5.1. Nodesets of the work piece ......................................................... 40 Table 5.2. Boundary conditions for Lagrangian based simulation ......................... 41 Table 5.3. Boundary conditions for ALE based simulation .................................. 42 vii LIST OF FIGURES Figure 3.1. Reference, Lagrangian and Eulerian configurations .............................. 12 Figure 4.1. Mappings for Lagrangian and ALE formulations ................................. 19 Figure 4.2. Nine node finite element ............................................................. 36 Figure 4.3. Flow chart of FEADSA ............................................................... 37 Figure 5.1. Discretized finite element domain ................................................... 40 Figure 5.2. Initial, reference, and deformed configuration for distortion energy constraint ..................................................................................... 44 Figure 5.3. Deformed configuration for element geometry constraint ....................... 45 Figure 5.4. Comparison of Lagrangian and ALE based simulation at various time steps .............................................................................................. 47 Figure 5.5. Reference frame configuration at various time steps ............................. 50 Figure 5.6. Zero energy modes .................................................................... 52 Figure 5.7. Non-zero energy modes .............................................................. 53 Figure 5.8. Lagrangian based simulation of extrusion process ............................... 55 Figure 5.9. Eulerian simulation of mould filling process ...................................... 55 viii Chapter 1 INTRODUCTION The finite element method is a widely used technique to simulate large deformation metal forming processes like forging, extrusion, drawing, rolling, upsetting etc. In this technique, the domain defining the continuum is discretized into a number of geometric elements on which the solution is approximated by simple basic functions. The material properties and constitutive relationships are expressed over these finite elements in terms of unknown values at element nodes. A set of equation results after assembling all the elements of the domain. The equations are then solved to obtain the approximate behavior of the continuum. The ability of the finite element method to solve problems with complex geometries makes it popular. The domain over which the finite element problem is formulated can vary and this leads to a variety of kinematic descriptions. Initially, Lagrangian kinematic description for finite element meshes dominated the field of large deformation analysis. In a Lagrangian description, a set of material points are identified with the same set of grid points throughout the course of deformation, thus allowing no material motion relative to the convected mesh [1]. Another popular description which has been extensively used over the years to simulate various manufacturing processes is the Eulerian description. In the Eulerian description, there is a material motion relative to the stationary mesh. We keep the reference system fixed in space and each point in space is unambiguously identified by an invariable set of three independent co-ordinates. In other words, in this description the material associated with a given reference volume changes during a deformation, which implies a material motion relative to the stationary mesh. The particles can translate across element boundaries and the particle associated with a node can change at each stage of deformation [2]. A more recently developed kinematic description is the Arbitrary Lagrangian Eulerian (ALE) description, which combines the above mentioned descriptions. As opposed to Lagrangian and Eulerian description, the ALE description introduces an additional (reference) mapping, which brings flexibility in the formulation. In particular, the ALE finite element mesh neither adheres to the material nor the space; rather the mesh is in relative motion with the material and it rearranges itself during the deformation process [3]. This brings more flexibility than the pure Lagrangian and Eulerian formulations. One of the far reaching results obtained by ALE is that it has effectively reduced mesh distortion in the interior of the deforming body. Thus, we can monitor the grid motion at sharp edges and comers of the work piece using ALE, resulting in computationally efficient and accurate results. 1.1 Numerical simulations The significant development in computational capacity has made it possible for us to analyze any physical process through the use of numerical simulations. Traditionally, engineers tend to use experimental techniques more than any other analysis tool; but they are time consuming and expensive. The acceptance of computational methods in industry is now increasing due to improved awareness, availability of appropriate software and reduced computational cost. The simulations in contrast to traditional analysis tools are more extensive and fast, they give more insight within a relatively short computer time. They allow much more inexpensive parameter studies than experimental procedures used previously. A numerical simulation usually consists of three steps: the first being the description of the effects that occur during the process (for example material behavior and frictional behavior etc) in a mathematical model. One must know the effects that occur in the process in order to construct a mathematical model that describes the real process accurately. The second step is to solve the mathematical model with a numerical method, which requires a good understanding of numerical mathematics. One should know the limitations of the numerical method, too. The last step is the interpretation of the results, and here it is important to have a good understanding of both numerical mathematics and the mechanics of the process. For example, effects that are due to numerical algorithm should not be interpreted as physical effects of the process and vice versa. Although research in computational methods has had numerous successes, considerable work still needs to be done. For example, there is a need to develop models which better represent material behavior and friction phenomenon. Work also needs to be done to improve existing computational tools. For instance, in the finite element method new computational procedures for reducing computational expense have to be developed. As such, issues regarding implicit versus explicit time integration schemes, mesh distortion and adaptivity must be researched. In this thesis, emphasis is on reducing the mesh distortion so that we get computational efficiency. 1.2 Scope of the thesis In this thesis the Arbitrary Lagrangian Eulerian description is investigated with the intention of making improvements. The merits and demerits of Lagrangian and Eulerian descriptions are discussed. Based on various shortcomings of these descriptions, ALE is implemented to solve large deformation metal forming processes. A comparison of pure Lagrangian based simulation and ALE based simulation is shown in a punching problem, with the results showing the fact that certain limitations of Lagrangian formulations can be overcome by using the ALE formulations. Two mesh governing algorithms are suggested, which can control the convection of mesh during the large deformation process. The focus is to improve the geometry of distorted mesh. When we use the Lagrangian description, the elements start losing their symmetry and thus affect the precision of the results. After formulating the punch problem using ALE description, the element formulation was coded in an existing FEM code written in FORTRAN 77 and the boundary value problem was solved. Appropriate material properties and integration schemes based on Newton Raphson method are used to get computational accuracies. 1.3 Outline of thesis The outline of thesis is as follows. In Chapter 2, various advantages and disadvantages of Lagrangian and Eulerian description are analyzed. It is emphasized, how ALE circumvents the limitations of the two descriptions and results in an improved formulation. Chapter 3 discusses the co-ordinate systems used in Lagrangian, Eulerian and ALE descriptions. The transformation equations are developed, and it is shown how the conservation of mass and momentum equations are expressed in reference frame. In Chapter 4, two different ALE mesh governing algorithms are presented along with iterative equations and finite element discretization. In Chapter 5, the two proposed ALE mesh governing algorithms are used to simulate a punch problem. The extent to which each mesh governing algorithm is successful in controlling convection of mesh is qualitatively analyzed and discussed. A comparison of the Lagrangian and ALE formulations is shown in the punch problem. The Lagrangian mesh suffers from excessive mesh distortion while ALE mesh rearranges itself during the course of deformation and results in a better geometry of final mesh. In Chapter 6, various conclusions are drawn about ALE formulations and future research directions are suggested. Appendices that are pertinent to the thesis are provided at the end. Appendix A discusses the matrix definitions and Appendix B gives some of the mathematical results that were used in finite element formulations. Chapter 2 KINEMATIC DESCRIPTIONS 2.1 Lagrangian and Eulerian descriptions In the last two decades, tremendous progress has been made in the area of numerical methods like the finite element and boundary element methods to analyze complicated large deformation problems involving a wide variety of non-linear materials. Significant contributions in this area are made by Hibbit, Marcal, and Rice [4]; McMeeking and Rice [5]; and Kikuchi and Cheng [6]. One of the most challenging aspects of finite element methods for non-linear problems is the determination of the proper kinematic description for the specific problem at hand. Most of these computational methods use a pure Lagrangian (total or updated) kinematic description for the finite element mesh. The Lagrangian description has enjoyed popularity because of the following advantages [7] 1. Pure Lagrangian approach has the advantage of having to satisfy less complex governing equations because of the absence of convection terms in the formulation. 2. When the Lagrangian description is used, the material points coincide with finite element mesh and quadrature points throughout the deformation; thus it enables to accurately define the material properties, boundary conditions, and stress and strain rates. Although, these merits of Lagrangian description seem attractive initially, however serious limitations are encountered when metal forming operations are simulated with Lagrangian approach. A few of these limitations are discussed by Liu [7], and Ghosh and Kikuchi [3]: 1. Lagrangian description neglects to take care of the convection effects, which is a serious limitation in using it with large deformation problems. Lack of control over mesh movement results in distorted and entangled mesh with large changes in element dimensions. Sometimes the entanglement and element distortion becomes so bad that even negative volumes are observed. Because of elements moving with deforming body, excessively distorted elements occur, which causes numerical problems and adversely affects the accuracies of solution. 2. Another drawback of using Lagrangian description comes when the boundary condition has to be specified on a material point, which might move itself. Situations in which we have contact boundaries, sharp edges or comers, complex shapes or abrupt surface discontinuities; the pure Lagrangian description is incapable of accurately representing the contact boundary condition. Moreover, using the Lagrangian description effectively alters the dimension and shape of the work piece, which is highly undesirable. 3. Moreover, Balagangadhar [8] pointed out that Lagrangian methods are computationally inefficient since they require large meshes, transient analyses, complicated contact algorithms and transient adaptive meshing strategies. Despite the popularity of the Lagrangian description, there are certain problems in solid mechanics that are described most naturally by an Eulerian description. For example, in modeling an extrusion process, it may be more convenient to monitor the flow of material through a fixed region of space in the vicinity of a forming die, rather than to describe the motion of a fixed set of material particles as they pass the die. In effect, the current configuration is discretized and the initial configuration is treated as unknown [2]. Because of the idea of relative motion of material to convected mesh in Eulerian description, we get some advantages which make it preferable in some situations, e. g., The Eulerian description is preferred when it is convenient to model a fixed region in space for situations that may involve large flows, large distortions, and mixing of materials. It has enjoyed great popularity to solve fluid flow problems [7]. The Eulerian description is preferred when we have to obtain the internal deformation accurately, a common situation in most fluid flow problems. Moreover, Eulerian methods have been widely used for free surface flow simulations. The pure Eulerian approach with its characteristic of keeping the mesh stationary has some difficulties as well, few of them are indicated in [7] and [9], e. g.: 1. Numerical difficulties due to convective effects arise because of the relative motion between the flow of material and fixed mesh. In the Eulerian description, the boundary of the deforming body does not coincide everywhere and always with an element side, thus making it arduous to take material associated boundary conditions into account. Therefore, special accommodation is needed in Eulerian description because material interfaces and boundaries may move through the mesh. Thus, in such situations where boundaries or interfaces move substantially, the Eulerian description is not really suited. 3. Another issue in Eulerian formulations comes in the treatment of free surfaces. The governing equations are expressed on the deformed configuration; however the free surface locations are not known a priori. To determine their location, time consuming iterative updating or successive recalculation are used [10]. 2.2 Arbitrary Lagrangian Eulerian descriptions In view of the above discussion of Lagrangian and Eulerian descriptions, the shortcomings of each one of these call for a kinematic description that should combine the advantages of both the above approaches into a single description. Thus, in recent years an alternative description called Arbitrary Lagrangian Eulerian (ALE) has been suggested to overcome these shortcomings. In ALE, the finite element mesh is neither constrained to move with the associated material nor is required to remain stationary in space. Thus in ALE, the finite element mesh need not adhere to the material but may be in general motion relative to the material. We get a totally flexible mesh, which can move with pointwise arbitrary velocity while reserving the potential to represent a Lagrangian or Eulerian description as limiting cases at points where such descriptions are desired [1]. ALE formulations can thus be effectively used on contact problems in solid mechanics (e. g. punching) to achieve computationally efficient and accurate results. This description has been employed by Haber [2], Liu [7] and by Ghosh and Kikuchi [3] to execute large deformation analysis of elastic-plastic solids. One of the far reaching results obtained by ALE is that it circumvents major difficulties associated with Lagrangian and Eulerian descriptions and it reduces the mesh distortion in the interior of the deforming body. Therefore, we can monitor the grid motion at sharp edges and comers of the work piece in an effective manner using ALE. The potential of the ALE description remains to be explored and is the aim of this thesis. A large volume of literature exists for the finite element solution of large deformation problems involving contact. Among the various methods for incorporating the constraint conditions are the introduction of contact forces through equilibrium conditions by Chan and Tuba [11], through the use of Lagrange multiplier techniques by Hughes [12], through mixed elements by Tseng and Olsen [13], and through penalty functions by Kikuchi and Oden [14]. In this thesis, two new approaches are introduced to govern the convection of mesh. One of the approaches is to introduce a constraint in the formulation that refers to the minimization of distortion energy. The other approach puts a constraint on the geometry of the deforming element so that elements preserve their geometry. The two approaches for ALE formulations are qualitatively compared to each other by implementing the proposed formulations in a punch problem. Moreover, another comparison is shown to highlight the considerable improvements brought by ALE formulations compared to the pure Lagrangian approach. 10 Chapter 3 CONTINUUM MECHANICS FOR ALE KINEMATIC DESCRIPTION In this chapter, continuum mechanics for the ALE description in general has been developed. Large displacements and deformations can occur in the simulations of forming processes. In simulating any solid mechanics problem with FEM, we have to formulate a mathematical model, using a set of independent variables and co-ordinates within a reference frame in order to identify material points of the body or points in space. 3.1 Co-ordinate systems In the analysis of processes with large deformations, the motion and deformation can be specified with respect to several co-ordinate systems. Depending on the choice of the co- ordinate system used in the formulation we can have three descriptions. 1. Lagrangian description where state variables are functions of material co- ordinates. 2. Eulerian description where state variables are function of spatial co-ordinates. 3. Referential or ALE description where state variables are a function of reference co-ordinates. In Figure 3.1 the relations between the three domains in the referential description are sketched. Since it represents the referential description, the referential co-ordinate x is the independent variable. As a result the mappings onto the spatial domain and the material ll domain are functions of time, which means that a referential point corresponds at each time to a different spatial point and a different material point. The mappings (I) and ‘1' are one-to-one mappings between three domains [15]. This means that every point in one of the domains corresponds to one point in the other two domains. Note that the Lagrangian and Eulerian descriptions can be seen as special cases of the referential description. The Lagrangian and Eulerian descriptions are obtained when 2! equals X and x respectively. Thus, in the Lagrangian description the mapping ‘I’is the identity function and in the Eulerian description the mapping (1) equals the identity function. X2 Xi Eulerian Co-ordinate (¢°‘P") Reference Co-ordinate X1 Lagrangian Co-ordinate Figure 3.1. Reference, Lagrangian and Eulerian configurations 12 3.2 Displacement, velocity and acceleration fields In this section, a general case is considered and displacement, velocity and acceleration fields are expressed in three co-ordinate systems described above. In the Lagrangian description the material particles X are marked with the initial position that will also be referred to as X = x(X ,0). The current position is expressed as a function of the initial position X and time t as x(X ,t) . The displacement field can be written as x=(,(,t) (3.1) X = m 1,0 (3.2) The Lagrangian description is especially attractive for describing physical variables that are associated to material points. The velocity field in three co-ordinate systems can be given as v(X,t) = E’gi) ad) ‘2’" X t t (33) mm = ( < . >. >> a: In the Eulerian configuration v(x,t) = v(X ,t) (3.4) v(x,t) = v(‘P(;{,tl),t) or (3.5) 00:, t) = v(‘I’(' (x,t),t),t) Similarly in the referential configuration mm) =v(x,t) (3.6) i"(.1’.t)= 9((¢(2’,t).t) (3.7) Defining the acceleration field in the three configurations in the same way 13 E)v(‘l"l (X ,t), 1)) a(X,I) = at (3.8) 5(x,t) = a(X,t) (3.9) 5(x,t) = a(‘I’("(x,t),t),t) (3.10) Eta/,1) = fi(x,t) (3.11) aunt) = fi((¢(z.t),t) (3.12) 3.3 Transformation equations The deformation gradient is defined as the transformation of the initial differential line element dX to the deformed differential line element dx. dx = FdX (3.13) where, F is given as F = I + Vu (3,14) The relative change in volume between the undeforrned and deformed state is equal to J. J = det F (3.15) For a scalar, vector and tensor fields ( a, a , and A) respectively, we have the transformation equations V Xa= 17ng [a (3.16) an=VlaFJ (3.17) VXA=VZAFJ (3.18) 3.4 Conservation of mass and momentum in reference frame The conservation of mass in referential form is of particular interest to solve any continuum problem. A detailed derivation of conservation laws in referential form can be 14 The conservation of mass in referential form is of particular interest to solve any continuum problem. A detailed derivation of conservation laws in referential form can be found in [16]. Consider an arbitrary volume (21 fixed in the referential domain, Rx bounded by a surface F1 and a continuous medium with a density [)(zJ). We can write the volume Q! with respect to fix and 9x representation of co-ordinates. The total mass in Q: at time t is given as M: [MD]: jdex= jp°d§2x (3.19) 91 (2x ox where p°(X .t.) = Jp(x.t) (3.20) fi(z.t)=jp(x.t) (3.21) 7 = det[§-’5-] (3.22) 81]- 3x,- J _ 111215;] (3.23) 1 Similarly, the principle of conservation of momentum can also be derived [16]. The total rate of change of momentum of the medium occupying at time t in the referential domain is given as [ fi(;(,t)v(x,t)dQZ= [1351+ [ figdflz (3.24) x o,z 1",. oz 8 a: where? is the force per unit area acting on the surface I"! and g is the body force per unit mass acting in (21. The force on the deformed spatial surface per unit of referential area 15 ? may be written as a function of the first Piola-Kirchoff stress tensor 7“ and the outward unit normal ii to the referential surface as It should be noted that the first Piola-Kirchoff stress tensor T is defined here in the referential sense, i.e., it is defined with respect to the fixed referential domain. Moreover, T is related to the Cauchy stress tensor 0 and to the first Piola-Kirchoff stress tensor in its classical sense T° (i.e., defined with respect to the material surface at to) by the fact that all of them give the same force on the deformed surface, de , but use different exterior unit normals and unit surfaces [16], namely (aims, = (no)de = (n°.T°)de (3.26) 01‘ , ax, TI'J' = Ja—x""0'kj (3.27) x x 32’. '1}!- = 151:0” (3.28) where n and 71" are the exterior unit normals to the deformed surface de and to the material surface de at time to respectively. Substituting equation (3.25) into equation (3.24) and using the divergence theorem to transform the surface integral to volume integral, one obtains A i1| 1 man]: I flwg, d522, (3.29) 31x 92 a, 81,- 16 The right hand side of above equation is transformed using Reynolds transport theorem and the divergence theorem into A " , Btu/w. 8T.- J 9.2/4. +__J_'€_'. (192:! _1_‘+figi dgz (330) Q, at 2' an 9: 32'; where ax]. (0-: 3.31 x and a) is defined as the particle velocity with respect to the referential co-ordinates. Equation (3.30) is reduced to ~ . aw.“ . 313.. ap_v,l + ’pv' = fl+figi inRZ (3.32) at Z 311- 32’,- After noticing that ins arbitrarily chosen, we can simplify the above equation by using continuity, thus the final form of the equilibrium equation in the referential form may be written in R2, as “av. . av. 3T}, . 1 ,_'=_ . 3.33 pat +ijazj [alj+pg.] ( ) Z If the Lagrangian description is used, the preceding equation is transformed using 1 = X (3.34) w = O (3.35) . ax. ] =J =d t -—' 3.36 e [axj] ( ) l7 ,5 =P° (3.37) The corresponding momentum equation for equation (3.33) in Rx in the Lagrangian description is given as . av. 3T3- . . p a“ = [$21—+ p gil 1n Rx (3.38) X i where we have used the fact T“ = T° (3.39) If the Eulerian description is taken 1 = x (3.40) a) = v (3.41) i =1 (3.42) 15 = .0 (3-43) The corresponding momentum equation for equation (3.32) is given as ._|' ..__'_= _ . R 3.44 where we have used the fact H. u q (3.45) For more details the reader is referred to Appendix A of [16]. Based on the basic concepts developed in this chapter for ALE description, large deformation formulations are developed in the next chapter. 18 Chapter 4 KINEMATIC FORMULATION FOR LARGE DEFORMATION PROBLEMS In this chapter, a rigorous development of pure Lagrangian and Arbitrary Lagrangian Eulerian kinematic formulations for large deformation problems is presented. We seek to formulate the constitutive relations and the finite element matrices for large deformation problems. Figure 4.1 shows the reference configuration used in the ALE kinematic formulations. It can be seen that the ALE formulation introduces a reference configuration, which consists of a set of grid points in arbitrary motion in space.1 /u\ x=X+u M (a). Lagrangian kinematic formulations x = X + u /g\ /u\« W (b). ALE kinematic formulations X=r+g x=r+g+u Figure 4.1. Mappings for Lagrangian and ALE formulations IIt should be noted that a different notation is chosen in figure (4.1) compared to notation in figure (3.1). All formulations in the subsequent chapters are based on figure (4.1). 19 The relationship between current and initial co-ordinates for any central point in pure Lagrangian description is given as x = X + u (4.1) The same relationship is given in the ALE description as X=r+g (4.2) x=X+u x=r+g+u (4.3) 4.1 Constitutive relations for large deformation Hyperelasticity In this section, various constitutive relations to be used in the formulations are presented. A hyperelastic material model is used in all formulations. The strain energy density function e for large deformation hyperelasticity is given as e023) = 2(5) = grog) wt, (156+ gage?» (4.4) where it and p are the Lame’s constants. 1} (E) = Tr(E) ((152) = nag?) g(i3(lf’)) is assumed to be zero in the formulation.1 The second Piola-Kirchoff stress tensor is symmetric and it can be computed using 3e _—_ 4' I g(i3 (F )) can be chosen based on problem at hand, e.g. g(i3(E)) = K (log(det E ))2 with K been an arbitrary constant. 20 or s = 2Tr(1§)1 + 2211;; (4.6) The first Piola-Kirchoff stress tensor is unsymmetrical, and it is given as p = F— (4.7) E = g (4.8) or If = I + Vu (4.9) The Cauchy stress tensor is given as a = 1 FT (4.10) - J -2 2 where J is the Jacobian given as J = det(E) (4-11) The Strain tensor can be given in terms of the deformation gradient as E= (ETE-I) (4.12) Nit- 01' g =%(Vu+VTu+VTuVu) (4.13) 4.2 Variational description for large deformation In order to construct finite element approximations for the solution of large deformation problems it is necessary to write the formulation in a variational form. The integral forms can be written either in the initial configuration or in the current configuration. The 21 simplest approach is to start from an initial configuration because integrals here are all expressed over domains, which do not change during the deformation process. The potential energy can be written in the initial configuration as [I = Ie(E)dQ—W (4.14) (I where e(l§) is the strain energy density function for large deformation hyperelasticity given in equation (4.4). Let V be the set of all kinematically admissible solutions v ={ue H1(Q):u=00n 1“,} where 1““ is the position of the boundary where displacements are prescribed. If u is a kinematically admissible equilibrium solution, then u is a stationary point of the potential energy. This implies that at u 5H=0V5ueV (4.15) The potential of the external work is given as W: [146,210+ jugidr, (4.16) o r, where t,- denotes specified tractions in the initial configuration and F, is the traction boundary surface in the initial configuration. Thus, we can write equation (4. 14) as [1: je(g)do— Iuibidfl— juitidr, (4.17) n n r, Taking the variation of equation (4.17), we obtain an e jaEijsijdo- Iduibidfl— j6uitidI‘, :0 v (Sue v (4.18) o o r, 22 where 6a,- is a variation of the initial configuration displacement (i.e. a virtual displacement), which is arbitrary except at the kinematic boundary condition locations 1‘“ where it vanishes. We note that by using equation (4.5) and constructing the variation of (SE. U , the first term in the integrand of equation (4.18) can be expressed in alternate forms as 51505:} = 6F,“.ij.5‘,-j (4.19) where symmetry of Sij has been used. The variation of deformation gradient may be expressed directly in terms of the current configuration displacement as 35u 6F“ E axk or, (4.20) i 5F“ 5 Vé'u (4.21) Using the above results, after integration by parts using Green’s theorem, the variational equation (4.18) can be written as 1J 5n = — [511k [(ijsij), + 6kg), 1219+ I5uk[ijS--ni - 6,), mm = o v 614 e V (4.22) a r, This gives the equations of equilibrium in initial configuration as (ijSl-j),i + (51.1)171 = 0 (4.23) PkiJ +171. = 0 (4-24) The initial configuration traction boundary condition is obtained as SUijni -(6ki)t,- = 0 V u e V or, (4.25) 23 4.3 Lagrangian formulation The Lagrangian formulations are obtained by considering the equilibrium equations derived above based on variational principles. The weak form and finite element formulations for pure Lagrangian description are derived by considering (5 H = 0 V 5a 6 V Equation (4.22) can be written as jau-(Divg+b)d§2— Ida-(En—IMI“, =0 v 5ue v o r, From vector algebra we know the relation Div(ETu) = 13 : Vu + u.DivE Using the above relation in equation (4.28), we obtain IDiv(ETJu)dQ— jgzvaudo+ jaubdo— [au-(gn—tmr, :0 v (Sue v n a a r, The divergence theorem IDiv v (19 = Jundl“ Q 1‘ yields ngaumr— jgzvaudm jaubdo— [661312.113 + J‘ 611.1211“, :0 v 5ue v r n a r, r, 01' jau.13ndr— j 13 : V5udfl+ jaubdo— jau.13ndr,+ j 611.1211“, = 0 v (Sue v P Q Q rt I‘, Canceling out terms and using the fact that 611 = 0 on 1““ we obtain [5:V6udo— [aubdm I5u1d1‘,=0 v aue v o n r, (4.27) (4.28) (4.29) (4.30) (4.31) (4.32) (4.33) (4.34) Assuming there are no body forces and tractions, we get rid of second and third terms. 24 Thus, the equilibrium equation in its weak form is written as 6H=IV6qudQ=OV6ueV or a 8e (4.35) an: [FTvauz—do=0 v due v {2 3E 4.3.1 Iterative equations for Lagrangian description Let u be the equilibrium solution and u°be a neighboring solution then, defining r 5 51-104, (5a) then r(u,5u) = O V 5u e V (4.36) we can write 0 ar 0 0 r(u,5u) = r(u ,5u)+[—).(u —u )+0("u-u )= 0 V due V (4.37) Bu As r is a linear function of 5a , ignoring higher order terms, there is a linear operator L and a bilinear operator a such that L(§u)+a(5u,Au) =0 V 5146 V (4.38) where Au 5 (u —u°) , indeed 2 FTV Au +VT A F a(§u,Au)= IVT(Au)V§u:—a-f-d§2+ IETV5u: a 82 (" ( ) ( u) ")dQ (4.39) n 315 o 81;: 2 2 Invoking the symmetry ofgki, the major and minor symmetries of SE:- , and rearranging the terms in the above equation yields 8e aze a (611,211.) = JVT(Au)V5u :ggdo+ JETvau :-a—E—2 ETV(Au)do (4.40) 25 2 where 83?: = 11(1 ® I ) + 22111 is a 4th order tensor. Upon discretization by finite element method we can express a(§u,Au) = (5u)T [K]Au T (4.41) L(6u) = (614) R where 5u and Au are vectors containing finite element degrees of freedom and K is the tangent stiffness matrix evaluated at u, and R is the residual vector at uo. The equation (4.38) results in the (incremental) finite element equations {R} +[K]{Au} = 0 (4.42) 4.4 Arbitrary Lagrangian Eulerian formulation In addition to the equilibrium equation, ALE formulations require equations to govern the convection of mesh. Various constraints can be imposed to achieve control over the mesh distortion. ALE augments the equilibrium equation by a functional, which helps to control the convection of mesh. This thesis presents two mesh governing algorithms to be used with ALE formulations and discuss about the ability of each algorithm in controlling mesh convection. Therefore, for ALE formulations 1'[ = [419219 ._W + HAW, (4.43) o anwnm = 0 v a. e v (4.44) Thus, the difference in Lagrangian and Eulerian formulation comes from the augmented mesh governing functional HMesh as indicated in equation (4.43). The choice of the mesh governing algorithm depends on problem at hand. 26 The deformation gradients and Jacobians for the ALE mappings shown in figure (4. 1) are given as dx F =— 4.45 -1 dr ( ) dX = — 4.46 ~8 dl' ( ) J , = det(F,) (4.47) 18 = det(Fg) (4.48) Moreover, the total deformation gradient can be expressed as E. = 553 (4.49) The residual corresponding to equilibrium equation (4.35) can be expressed in ALE in the following way. Starting from an: IVdu:EdQ=OV 5ueV 0 Using, E = ES and If = Ff; in above equation we get an: IVJu:(F,Eg"§)dQ (4.50) (2 Now from the transformation equations, we know that V§u = 161.58" d9 = J 8 d9, we get the weak form of equilibrium equation in ALE description given as 611 = j (V,6u5;) : (155848)]ng =0 v 514,55: 6 v (4.51) a, 27 In the following section we discuss the variational formulation and finite element discretization for the two mesh governing algorithms. 4.4.1 Minimize distortion measure In addition to the equilibrium equation, we introduce a functional developed by minimizing the distortion energy measure given as 1 l-IMesh = ISL/t —1)2er (4.52) 9r The variation of the above equation yields the weak form to be used in the finite element formulations. The details of this can be found in Appendix B. 511M“, = j (J, —1)J,(E,’T: 613mg, = 0 v 5u,§g e v (4.53) a. By noting that in the ALE formulations 15 = I+V,g +V,u (4.54) E8 = I +Vrg (4.55) 515, = V,6g + V,§u (4.56) 5E, = 17,58 - (4.57) Equation (4.53) can be rewritten as 511M“, = j (J, —1)J,(E,‘T :V,5g)dQ, + a. (4.58) I (J, -1)J,(If,'T :V,5u)dfl, = o v 514,6g e v 0, Replacing the values from equations (4.51) and (4.58) in equation (4.44) we get 28 an: [ (V,6u5;‘):(151§g"‘s)Jgdo,+ j(J,-1)J,(13‘T.V,6g)do, Q, Q,- (4.59) + j (J,—1)J,(1,E‘,’T :V,6u)d£2, = 0 v 5u,5g e V 4.4.1.1 Iterative equations for ALE formulation Let u , g be the equilibrium solution and 14° , g° be a neighboring solution then, defining r = 5H(u, Eu, g, 6g) (4.60) with 6g = 0 on 1"“ (4.61) we can write, 0 o a o a , r(u,§u,g,6g)=r(u ,6u,g ,5g)+(5£].(u—u )+[3;-].(g-g ) +0(||u-u°ll)+0(llg-s° (4.62) )=0V 5u,5geV As r is a linear function of 6a and 5g, there is a linear operator L and a bilinear operator a such that L(5u)+L(5g)+a(§u,Au,5g)+a(5g,Ag,§u) =0 V 5u,6g e V (4.63) Indeed (6a, Au, 6g): 0] (V, Jqu"') (95— 1:1.{15 }Fg' ‘5)1 do -T +QJ;(V,5uF-l)1(FF-laS—zhAuDJ d!) + Q[(V 6a): ((.9 BF, uu{A })(J —1)JdQ (4.64) + I(V,6u:E-T)a((1‘aul)1’){Au}dQ + n[(V 6g): (31; :t{Au})(J —1)JdQ n . -r a((J,-1)J,) + ertVrfig-F. ) a“ {Au}do, and 29 F-l a; {Ag}):(F,Fg’1S)Jng, a a(5u,6g,Ag) = J(V,5u a -1 31" gg{A }F g“ + j (V,5qu"): (— a. + a] (V,§qu“) : (1‘ng-l a—{Ag} Fg‘lS)Jgd:, (4.65) -r + AW, 5uFP—l) (FFg-IS)-5§-{Ag}d§2 + “[(V 6a): (a 81“, g{Ag})(J _1)J dQ -r +j:(V,u6:E‘T)a((1’a_l)1'){Ag}dQ,+n_[(V 5g): (8:: {Ag})(] _1)J dQ a 8 _ 3((1 —1)J) + (V,5 IF,T) ‘ ' A do, o‘[ g 38 { 8} Upon discretization by finite element method we can express a(5u,Au,6g) =(6u)T[Kl]{Au}+(6g)T[K2]{Au} a(5u,Ag.6g)=(5u)T[K3]{Ag}+(58)T[K4]{A8} (4.66) L<6u> = (6a)T R1 L(58) = (63? R2 where 611 ,6g ,Au and Ag are vectors containing finite element degrees of freedom and K1 , K 2 , K 3 and K 4 are the tangent stiffness matrices evaluated at uo and go , and R1 and R2 are the residual vectors at uo and go. The equation (4.66) results in the (incremental) finite element equations {R1} +[K1]{A"} +[K2]{A“} +{R2}+[K3]{A8} +[K4]{Ag} =0 4.4.1.2 Finite element discretization In order to get the element stiffness matrix we need to manipulate the above equations in such a way so that they can be easily implemented in a FORTRAN code. The transformed equations are given below. For more details refer to Appendix B. 30 4.4.1.2.1 Residual terms a. j (J, —1)J,(F," :V,6g)dQ, = {6g‘1’ j [B]T{F," }(J, —1)J,do, Q, Qr j (J, —1)J,(F," :V,§u)d£2, = {611‘ }’ j [B]T{F,'T )(J, —1)J,dQ, 0, 4.4.1.2.2 Stiffness terms Q[(V, 6a F,“): (— —u‘{Au}(F,' 1S))J, do: 16a 1’ [131 ’11oF,"1’1Io1Fg“s1’11B1J,da.1Au1 0r ‘3—5-{Au})J, do: [ (V,§uF,“) : (F,F,' 12. 16u‘1’j1B1’110F,”118F,"0111C11F’01110113118412.1411} 8F," [1V6u11 9r {MIX-l "DerQr : —16u1’:]181’18"oF 11T11B11J— 111412 1Au1 n[(V,6g:1«‘,") 3((11— 010”“er = Bu 16a 1’ [1B] T{F.‘T}{1’,"T}T[Bl(21,-1)J,dfl,{Au} 98F," [1V 6g1za1 12. 14u111J 411.49, = —16g 1’ [1B]’1F;’08“11T11B11J—1>J 41214111 3((11— ‘1)11) Bu 16g 1’ [131’18"118"1’131121.—1>J.dn.14u1 0r QJ'(V, 6g :F,") {Au}dQ, = 31 j(V,§uF,“) : (F,F,“S)J,d12, = {611‘ 1’ j [B]T[I o F,"]’1F, o S]{F,’1}J,dQ, 12. (4.67) (4.68) (4.69) (4.70) (4.71) (4.72) (4.73) (4.74) (4.75) I(V,§u Q1' 3178-] . -l ,8 1Ag11.1F.Fg S1Jgdfl.= —16u‘1’ j1B1’1F,F,“S 011111111,“1 0 F,"’11B1J,d£2.14g1 0r j (V,5uF,") : (£{Ag}F,"S))J,dQ, = a, 8g 16::‘1’ ] 181’11oF,"1’1I01845101814412.1411} Q —l I r ~16u‘1’ j 1B1’1I o F," 1’1F, 0511510 F,"11B1J,d0,14g1 a, F_1_BS j (V,5uF,"):(F,, {——Ag})J, (112,: 1614‘ 1’ [IBM an,"1’ 1F.F,"1 0111611F’ o 111P11BngdQ.1Ag1 j (V auF,.") (F,F,%—“S) 8 g{Ag}dQ,=- 9r 16a‘1’ [181’1IoF:"1’1F: os11F,“11F," 1’131Jgda,14g1 -T [1V 6a): (LE, 1Ag111J. 411.49, = —16u 1’ ] 131’1F." OE“][T1181(J, 411.412.1411} 9. I(V,5u : FIT) 9r 8((11—1)Jr) Q : ag {Agld r 16u‘1’ [181’1F."1 1F." 1’1B1121. — 111.419.1411} 8F," [1V,6g:14g111J 411.419. = 9r —16g 1’8 ] 181’1F." oF.“117"11B1121, —11J.dn.14g1 nr [(17,811 , E") 810. —11J,1 38 {138}er = 16g‘1’ j1B]’1F."1 1F."1’181121. — 111.42.14.11 32 (4.76) (4.77) (4.78) (4.79) (4.80) (4.81) (4.82) (4.83) (4.84) 4.4.2 Constraint on element geometry The second approach used to control the convection of mesh during large deformation problems was to preserve the element geometry during the course of deformation. For this constraint 611,,,,,, = j {6g}.{VJ,}d12, = 0 V§g e V or (4.85) a. 5W... = [ 6g,J,,do, = 0 V§g e V 1:, HM. = j (6g,J,),do, — [ 6g,,J,do = 0 V5g e V o, a, 611,,,,,, = j 6g,J,n,dr, - j (V.5g)J,dQ, = 0 V5g e V r, r The first term on left hand side of above equation is 0, so we get 611m, = j (V.5g)J,dQ, = 0 vag e V (4.86) 9r Combining the equilibrium and mesh governing equations, we get 61] = [ (V,5uF,“) : (13F,"8)J,d12, + j (V,.6g)J,do, =0 vau,6g e V (4.87) o, 12, Indeed a(§u,Au,6g) = I (V,6uF,“) : (%{Au}F,'IS)J,dQ, + a, u 4. ,as ( 88) _ % Bu AudQ au{ } , j (V,6uF,“1):(F,F,‘ {Au})J,dQ, + j (V,6g.1) Q, 9r 33 Fg-11A 'FF‘ISJdQ a, 311.1,,1, ,+ 8 a(§g,Ag,5u) = Q_[(V,6u j (V,5uF,7‘) :(——'16g}F,,7‘S)J :10, + a] (V, auF,7‘): }S)J,dQ,+ 0. ag (4.89) 91w, 5a F,7:‘) (33— E{Ag}F,,’lS)J :10, + j (V auF,7‘ :F,F,aJ7‘S) —1Ag}do, + 88 j1V.6gJ1%J—'1Ag1 a, 8 The above equations are transformed in the following way for easy implementation in FORTRAN code. For more details refer to Appendix B. 4.4.2.1 Residual terms [ (V,§g)J,dQ, =15g‘1’ [ [B]’11}J,do, (4.90) a. r [(V,6uF,7‘):(F,F,7‘S)J,d12,=16u'1’ j [B]T[IOF,‘T]T[F, 0511F,7‘}J,d12, (4.91) 0, ar 4.4.2.2 Stiffness terms n[(V, 5a F,:7‘) (9,,— 1:"{,Au}(F ‘S)J, dQ= (4.92) {511 }’ [[81’110F,7’]’[101F,7‘S}’][B]J,d12,1Au} 12. j(V,6uF,7‘):(F,F,‘—1Au})J,aSdo: r (4.93) 16u‘1’ ]1B1’IIoF,7’1’ 1F.F,7‘ 6111C11F’6111011811,40.1Au1 -T j (V 5a): (3:; {Au })(J, —1)J, dQ,= (4.94) —16u 1’ [181’187’ oF.“11T11B11J.~11J.d0,14u1 0. 34 I (Vr58 :FFT) 12. a“ " 1)") 1Au}d12, = Bu 16u‘1’ I131’187’1187’1’1B]12J. —11J,dn,14u1 9r —T j (17.581131; 12. {Au})(J1-I)erflr : u —16g‘1’ [181’187’ o F.7‘11T11311J, — 11141114111 0 f j (V158 = Fir) 0. a((Jr _1)Jr) {Au}dQ, : au 16::‘1’ I 181’18"11F7’1’1B1121, 411.412.1411 9r aF7‘ _, j (V,6u 8 {Ag}) : (F,F, S)J,d£2, = o, 38 _W. 1’ I181’1F.F,7‘s 0111T11Fg7‘ o F,7’11B]J,dfl,14g1 Gr 1 (V,5uF,7‘) : (%E’—{Ag}F,"S)J,dQ, = 12. 16u‘1’ 1131’11oF,7’1’1Io1F,7‘S1’11B1J,d!2,1Ag1 9r -1 j (V,6uF7‘) : (F, 8F“ 1Ag}S)J,d12, = o, g 38 —16u‘1’ j1B1’1I o F,,7’1’1F, oS11Fg7‘ o F,7’11B1J,d£2,14g1 fir _ - BS j (V,5uF,1):(F,F,IE{A8})Jngr = r 16u‘1’ [181’11 o F,,7’1’18Fg7‘ 0111C11F’ 0111P1181Jgdn,14g1 0, -l -l 318 [ (V,6uF, ):(F,F, S)—1Ag}do, = o, 88 (My [131’11 o F,,7’1’1F, 6511Fg7‘11F,7’1’1B]J,da,14g1 Dr 31, (V1581) o"; 36’ 14:149. =16g‘1’ 1131’1111Fg‘ 1’18140.1Ag1 9r 35 (4.95) (4.96) (4.97) (4.98) (4.99) (4.100) (4.101) (4.102) (4.103) 4.5 Numerical implementation A wide variety of numerical details have to be considered in the construction of a global solution strategy for solving large deformation problems using the formulations discussed above. The physical grid point domain is discretized into isoparametric nine-node elements with variables interpolated by the Lagrangian shape functions. It is to be remembered that the elements are essentially made up of grid points rather than of material points. A schematic of the nine node element is shown in figure (4.2). The node numbering scheme is also shown in the figure. Each node has four degree of freedom, namely uI , u2 , g1 and g2 . The degrees of freedom for node 1 are shown for illustration. “241 l l 1, '5’ 2> “1 1 g 1 Figure 4.2. Nine node finite element A FORTRAN 77 computer code named FEADSA has been developed for the implementation of the aforementioned algorithm. The mesh generation, shape functions, derivatives of shape functions, deformation gradients, residuals and stiffness terms are coded up as an element in FEADSA and then a frontal solver is used to do the computations. A flow chart describing the overall structure of the code is shown in figure (4.3). 36 [ Initialize ; 1 Yes I Assemble K, M, R 1 Evaluate Pseudo-loads ] Converged [ Solve for sensitivities ] No Yes Figure 4.3. Flow chart of FEADSA The above figure depicts how the code can be used for a transient, nonlinear problem to perform both the analysis and the design sensitivity analysis. For steady problems, the transient part of the code (the outer loop on the left hand side of figure (4.3)) will not be active. For linear problems, Newton Raphson iterations (the inner loop .on the left hand side of figure (4.3)) will not be invoked during the solution procedure. If the design sensitivity analysis is not required, the evaluation of the pseudo-loads and the solution of sensitivities (the loop on the right hand side of figure (4.3)) will be eliminated. 37 An element “LRGDEF2D” was written in this code based on the ALE formulations discussed in previous section. An input file is generated and boundary conditions are specified to solve the Boundary value problem. The element is capable of doing the calculations based on Lagrangian description and ALE description at the same time. If a boundary condition is specified at all nodes imposing g1 and g2 to be zero, the formulation becomes a pure Lagrangian description without any control over mesh convection. This facilitates the comparison between the Lagrangian and ALE simulations. The element “LRGDEFZD” is capable of solving any large deformation metal forming problem using ALE description. In the next chapter a simple academic test case is performed. A comparison of Lagrangian based simulation and ALE simulation is shown in a simple 2D unsteady punching problem. 38 Chapter 5 PUNCH INDENTATION PROBLEM-A NUMERICAL EXAMPLE In the previous chapter the Lagrangian and ALE methods were described in detail. The numerical simulation of large deformation problems can be performed using the relatively simpler Lagrangian method. However, generally in the case of large deformation, distortion of mesh occurs. As a result the calculation becomes inaccurate or it may even crash. The ALE method can be used to reduce mesh distortion. Several large deformation problems including upsetting, punching, extrusion, ring rolling, etc, can be simulated using ALE based formulation discussed in chapter 4. In this chapter the applicability of ALE method in simulation of a punch indentation problem is shown. The punch problem has been widely used by many authors to validate their results of ALE formulations before. Some of the noteworthy contributions in this area are Haber [2], and Brekelmans, Veldpaus, Schreurs [9]. The simulation of punching problem is carried out using both the Lagrangian and ALE descriptions and results are compared to validate the concept of ALE formulations. The figure (5.1) shows the work piece represented by a discretized domain. The nine node finite elements with node numbers can be seen in figure (5.1). 39 '7.l-~..68§-l$92-i?91—-I7t'-*;72.1—-:?33—T74l—;755—261—577 I I I I I I I I I I I-—57—-158——59I—-160r—-161‘v—-1 62L—63—64—<651—-66 ILILIIIILI 4 1—46—{47—481—49—(501—51:——52——~53—-154~—55 I I I___3|7 3'9 4'gf—41—142-«43r—44 7*1 TV 111’“ 1. 1 I 2 31—24—425—261—27—128I—29F—30 —131}—-I32,I—33 0 ‘1—77 21—113 —.7—;4‘7—: 511—1161—1777—1—181—31 4177771101—1— 111 0 2 4 6 8 1 0 x Figure 5.1. Discretized finite element domain The Table (5.1) below shows different nodes of the dicretized domain comprising the various nodesets. This facilitates the enforcement of boundary conditions on the work piece. Table 5.1. Nodesets of the work piece NODESET NODES LEFI‘ 1,12,23,34,45,56,67 RIGHT l l,22,33,44,55,66,77 TOP 67,68,69,70,7 l ,72,73,74,75,76,77 BOTTOM l ,2,3,4,5,6,7,8,9, 10,1 1 PUNCH 67,68,69 40 5.1 Boundary conditions The boundary conditions that were enforced on the domain while performing the Lagrangian simulation are shown in Table (5.2). It must be noted that in Table (5.2), all nodes have 3rd and 4th degrees of freedom equal to 0. In other words, the effect of ALE mesh governing equations is cancelled out. Thus, there is no reference frame to be used in the simulation; making it pure Lagrangian based simulation. The displacement boundary condition enforced on the punch nodes make the domain move in the same way as the work piece moves during the forging process. The punch displacement is carried out in 20 time steps. A displacement of —2.0 means a 33% reduction in height of the original work piece. Table 5.2. Boundary conditions for Lagrangian based simulation NODESET DOF 1 DOF 2 DOF 3 DOF 4 LEFT 0.0 FREE 0.0 0.0 RIGHT FREE FREE 0.0 0.0 TOP FREE FREE 0.0 0.0 BOTTOM FREE 0.0 0.0 0.0 PUNCH 0.0 -2.0 0.0 0.0 41 The Table (5.3) gives the boundary conditions to be associated with the ALE based simulation. It must be noted that the 3rd and 4‘h degrees of freedom are able to move in this case and this enables mesh rearrangement during the course of deformation. Table 5.3. Boundary conditions for ALE based simulation NODESET DOF l DOF 2 DOF 3 DOF4 LEFT 0.0 FREE 0.0 FREE RIGHT FREE FREE 0.0 FREE TOP FREE FREE FREE 0.0 BOTTOM FREE 0.0 FREE 0.0 PUNCH 0.0 -2.0 0.0 0.0 5.2 Results and discussion The punch indentation problem is solved using Lagrangian and ALE based formulations. When ALE formulation is used, both mesh governing algorithms discussed in chapter 4 are implemented to see which one behaves better in controlling convection of mesh. In the next section, first a comparison of the two mesh governing algorithms used in ALE formulations is presented and then a comparison of Largrangian and ALE simulation is shown to see mesh distortion in both kinematic descriptions. 5.2.1 Comparison of mesh governing algorithms Figure (5.2) is the result obtained by using the distortion energy algorithm in ALE formulation. It shows the initial configuration, reference configuration, and deformed 42 configuration of the work piece at 20th time step. The mesh governing algorithm based on minimization of distortion energy is successful in controlling the mesh movement to some extent, however it can be seen from the deformed configuration of figure (5.2) that there can be certain improvement made in the quality of mesh beneath the punch. The reference frame moves during the course of deformation and helps the mesh to adjust and avoid excessive distortion. However, the rearrangement process is still not enough. There is certainly more room for rearrangement. This approach may have serious limitations because when plasticity is added to the problem, the constraint will struggle to maintain the mesh quality. Instead we might end up in getting shear distortion and entanglement. This motivated the use of another algorithm based on putting a constraint on the geometry of the element. 43 Figure 5.2. Initial, reference and deformed configuration for distortion energy constraint The problem was also solved using the element geometry constraint discussed in chapter 4. Figure (5.3) is the result of the implementation of this constraint in ALE formulation. It only shows the deformed configuration at 20th time step. A detailed discussion of mesh rearrangement and reference frame motion is discussed in the next section. However, it can be easily seen that the mesh distortion is much less than what we had with the distortion energy constraint. The elements preserve their geometry even after significant 33 % reduction in height of the work piece. 44 \l TYl~17l~111l1111T1Tf:T'1.YIT “II \\ A I / // #P/ I O I. )— Figure 5.3. Deformed configuration for element geometry constraint The reference frame shows lot more rearrangement in this case, so the final mesh quality is improved. This constraint gives us much better results as compared to the minimization of distortion energy constraint, so ALE formulation based on this constraint will be used while making comparison with Lagrangian description. 5.2.2 Comparison of ALE and Lagrangian simulations The next important comparison is between the two kinematic descriptions discussed in the thesis. In this section, ALE simulation of the punch problem based on element geometry constraint is compared to the results of simulation of exactly same problem simulated with Lagrangian description. Figure (5.4) shows a comparison of the two simulations at 15', 6th, 11th, 16'“, and 20th time step. 45 ALE Lagrangian Time Step = 1 Lagrangian Time Step = 6 Lagrangian Time Step = 11 46 Lagrangian Time Step = 16 Lagrangian Time Step = 20 Figure 5.4. Comparison of Lagrangian and ALE based simulation at various time steps Clearly, the results of the Lagrangian description differ strongly from the results of the ALE calculations in the above simulation. It can be easily seen that the Lagrangian description captures the displacement for first 10 time steps almost comparable to ALE description. However, when the punch really gets deep, it starts losing precision and severe mesh distortion can be seen at 16'h time step. When it gets to 17th time step, which is about 25% reduction in height, we can see that the Lagrangian mesh crashes. Severe entanglement can be seen in the final time step. 47 In the Lagrangian method the elements become severely distorted in the region with large strains. This is a general characteristic of the Lagrangian method, i.e., the elements become distorted or stretched in regions with large strains. Especially in these regions, one would like to describe the solution more accurately. In the ALE calculations the elements are made to preserve their geometry in the high strain region. The effect of geometry can be clearly seen when the Lagrangian based simulation crashes in 17th time step, whereas the constraint of element geometry helps the mesh to convect smoothly during the course of motion. Even at the end of 20th time step, the elements are still in good shape in the case of the ALE formulation; whereas they can be seen totally collapsed for the Lagrangian simulation. This mesh quality improvement in ALE formulations is attributed to the introduction of the reference frame. The mesh rearrangement during the course of deformation is shown in figure (5.5). The reference configuration is shown at corresponding time steps. It can be seen that the rearrangement in case of element geometry constraint is lot more effective than the one seen for distortion energy constraint. The comparison of reference configuration for the element geometry constraint and distortion energy constraint at 20th time step reveals that there is more rearrangement in case of element geometry constraint. This helps to reduce mesh distortion in large deformation analysis. 48 ilII'FIT‘rF‘II'III' 11W III11117 1' I 11111 I. l L 1 LL 1 LL 1 O J N W # O 1111 1111-1111! O .b lr'1'111rrrY l-lI'TlI-l 111117.! 1 4 6 Time Step = l 10 '1 I l I 1 1 L 1 l J 1 I I'II‘I'ITIIIIII '11:. 4 Time Step = 6 10 1111 il.i J l L 4 4 Time Step = l 1 49 6 1O o 2 4 6 a 10 Time Step =18 CD I 0') 11'11‘1'1111111 I'TTII‘IT1I« 1 I —-——- A l A l l l I I 1 l 1 L i 1 Time Step = 20 10 @1- Figure 5.5. Reference frame configuration at various time steps Energy modes The above formulation is successful in terms of controlling mesh movement. However, before doing actual forming problems; it is useful to do the modal analysis of the problem. The eigenvalues and eigenvectors of the global stiffness matrix were plotted to 50 study the zero and non-zero energy modes. Figure (5.6) gives the zero energy modes and figure (5.7) gives non-zero energy modes of the system. 51 Figure 5.6. Zero energy modes 52 Figure 5.7. Non-zero energy modes 53 The study of zero and non-zero energy modes gives an insight into the study of element deformation. It is essential to make a few changes in the element formulation to remove zero energy modes. The reason being that when those modes come into play, they can cause undesired behavior, which affects the convection of the mesh. It remains one of the future explorations to come up with a even better constraint condition that would remove zero energy modes and give a better mesh quality. 5.4 Applications In the above section, we simulated an example of academic interest. It was a simple case illustrating the use of Arbitrary Lagrangian Eulerian description in numerical simulations and its advantages over pure Lagrangian based simulations. There are various industrial manufacturing processes that can be modeled and simulated using ALE methodology. The analysis and optimization of various forming processes like extrusion, forging, and rolling are facilitated by the use of ALE based numerical simulations. For instance, the simulation of the extrusion process with a sharp angle at the die exit can be simulated with ALE simulations. The elements at the die exit are locally distorted, so it becomes necessary to use ALE to capture the behavior close to the sharp edge. This situation is shown in figure (5.8) taken from Stoker [15], where an extrusion process is simulated based on Lagrangian formulations. The elements at the die exit are locally distorted because of the sharp angle at the die exit. Such a situation can be effectively handled by ALE formulations. 54 a M1 WSW: '1 11mm” ‘ mama AMWIM Fig 5.8.Lagrangian based simulation of extrusion process Another problem is that it is difficult to obtain an accurate description of the free surface during various manufacturing processes. The ALE method is found very useful to capture free surface motion during processes like injection molding of a polymer. The free surface motion in the Eulerian simulation of mould-filling process is shown in figure (5.9) taken from Stoker [15], where the shaded portion is the material. (a) “II (C) Fig 5.9.Eulerian simulation of mould filling process Significant improvement is made in this area of research by Stoker [15]. The proposed algorithms in this thesis can be used to model the extrusion and mould filling process described above. 55 The ring rolling process has been investigated using ALE by Liu and Hu [17], there is certainly a lot of work to be done in simulations of plane strain ring rolling process based on ALE formulations. Typically, complex contact conditions, and variation in frictional conditions on an interface for ring rolling process can be analyzed using ALE description. 56 Chapter 6 CONCLUSIONS AND RECOMIVIENDATIONS This thesis elucidates the power of the Arbitrary Lagrangian Eulerian finite element method to solve large deformation problems. The algorithms described in this thesis may find direct application in metal forming analysis and fracture mechanics problems where many compromises are often made in accordance with the limitations of the conventional method. The overall algorithm proved to be numerically stable and reasonable in the sense that the results obtained are in very good accordance with physical experience. The ALE approach is implemented for the analysis of nonlinear solid mechanics problems, and it has been shown that this approach is useful for problems where severe mesh distortion is anticipated. The limitations of the Lagrangian description, when it comes to large deformation solid mechanics problems are presented. The Lagrangian description lacks control over the mesh movement resulting in distorted (sometimes entangled) meshes with large changes in element dimensions, adversely affecting the accuracy of solutions. Moreover, problems involving certain contact boundaries, especially those with sharp edges or comers, may not be represented precisely in this description. This is due to the fact that the boundary condition has to be specified on a material point that might itself move. Situations of this kind are frequently encountered in the numerical simulation of metal forming processes, e. g., extrusion, drawing, etc., where the punch or die faces may have acute edges or abrupt surface discontinuities. The suggested algorithms based on ALE formulations are successful in circumventing these problems. A punch problem is simulated based on Lagrangian and ALE formulations 57 with the result suggesting that an ALE based simulation is computationally more stable and gives better results. The Lagrangian description is unable to sustain large deformation and ends up in severe entanglement. A major objective in metal forming is the determination of the initial shape of the work piece (preform) and of the process parameters (e.g. the die shape) that lead to a final product with desired geometry and material properties. The solutions to these inverse problems are usually obtained by trial and error methods, using the results of analysis for each set of preforrns and process parameters. The ALE based finite element formulation presented in the thesis can be very useful for preform design problems in metal forming. For example, an ALE formulation is better suited to be used in design of preform in open die forging of a work piece, where we want to achieve a final product without barreling. Such preform design problems have been solved by Zabaras and Badrinarayanan [18], Grandhi and Gao [19], Zabaras and Srikanth [20]. The preform design and the die design problems can be formulated under a rigorous mathematical basis by posing them as optimization problems. The objective function for these optimization problems can be defined as an appropriate measure of the error between the desired final state and the numerically calculated state for a given set of design variables. A sensitivity analysis can be performed to evaluate the gradients of the objective function. The present work can be extended to develop a finite element method for large deformation analysis of elastic-viscoplastic (time dependent plastic) material behavior. The ALE finite element model for large deformation analysis of rate sensitive materials can be very useful for industrial applications. The response of materials subjected to loading at elevated temperatures (~0.5 melting point) is best represented by rate 58 dependent plastic models exhibiting the characteristic of combined creep and plasticity. ALE based simulations can be made to model various sheet metal forming operations, where grain orientation due to pre-rolling gives rise to structure anisotropy. The presented models can be modified in future to include temperature dependence, rate dependent yield strength and plastic anisotropy. Moreover, it will be interesting to develop the conservation laws, constitutive equations, and the equation of state for path- dependent materials using ALE finite element method. 59 APPENDIX A MATRIX DEFINITIONS Any tensor is represented in a vector form with ‘ (A11 Arz A:A=< vec1111 421' 1422. For any two tensors, A and B, with components A=[All A12:IB=|:BH BIZ] A21 A22 321 322 The tensor multiplication AB AB = [41311 + Arszr 411312 + 1412322 ] 421311 +42sz1 A21312 +A2232 can be expressed as -All 0 A12 0 - 0 A 0 A, AB = “ 2 B { } ,2, 11,, 0 { 1 _ 0 A21 0 A22_ PBII 321 0 O 7 B B 0 0 {AB}: 12 22 {A} 0 0 BH 32, _ 0 0 3,2 322- OI' {AB}=[A01]{B}=[IOB’]{A} 60 (A.1) (A2) (A3) (A4) (A5) (A6) where the Kronecker product 0 is defined as [AoB]=[ 4.8 4:8 42.8 42.8 The Kronecker product 0 has the property {BXA} = {BoA’} {X} Similarly we can also define and [A1813]C=(B.C)A [A813], =43, Two matrices [B] and [B]T are defined such that (W) = 1811161 {Vru} = [B]T {ue} Moreover the double dot product of two tensors A and B is defined as We can write the double dot product of two tensors A and B in terms of two vectors A and B. Let A:B=A,.,B,,. IandB=< 61 I (A.7) (A8) (A9) (A.10) (A.11) (A.12) (A.13) (A.14) 1311‘ A'B-{A A A A )13‘21 (A15 ' _ 11 21 12 22 B ' ) 21 1322. 01' which is the same result as in equation (A. 1 3). Thus the 2“" order tensors in the formulations are considered as vectors and the double dot product is taken in the same way as shown above. Thus, j (A):(B)dQ, = [11’qu (A.17) 0. a. The Transpose matrix [T] is defined to obtain the transpose of any vector. ME IF” F T F If F=I "1 and F =1") { 1 F,, I I F,, We. 1%: then 1F,,‘ ’1 o o 0‘1F,,‘ F 0 0 1 0 F , 21,= , 12, (A.18) E2 0 l 0 0 F2| IF22. -0 0 0 1- IF22. '1 0 o 0‘ 0 0 1 0 , . where [T] = rs called the Transpose matrix. 0 1 0 O _0 0 0 l, 62 APPENDIX B IMPORTANT RESULTS USED IN FINITE ELEMENT DISCRETIZATION In Appendix B, few of the terms used in finite element formulations are derived. These equations are based on matrix definitions given in Appendix A. E=%[FTF—I] (8.1) S=[C][E] 113.21 F=I+V,u (8.3) F, =I+V,g (B.4) F} = I+V,g+V,u (B.5) 3111184311151.) 18.61 Bu 8F, Tg'IASI-IBIIASI (3’7) @41414311431 (13.81 38 317—1 -1 —T $=—[F OF ] (B.9) aF-T =—[F" oF7‘][T] (B.10) BF 887‘ -1 -T 1Au}=—[F, OF, ][T][B]{Au} (8.11) Bu 8F,“ -1 —T ag 1Ag}=-[F, 0F, ][T][B]1Ag} (8.12) 63 aF7‘ 32:14:11 =—[F.“ 06"][TllBl1Ag1 (13-13) 31714 —T —r a, {MP-[E 0F. ][T][B]{Au} (8.14) aE-T -T —1 ag 1Ag)=—[F, OF, ][T][B]1Ag} (8.15) 31:84 —T -1 ag 1Ag}=—[F, 0F, ][T][B]{Ag} (8.16) $1M} = [I o F,7’][B]{Au} (8.17) where Q =[1 O F,,-T] (B.18) g—g-{Ag}=[[lOF,"T]—[FOF,'T]][B]{Ag} (B.19) where P=[[IoF,7’]—[F@F,7’]] (8.20) 33141=1c11flc>q12113114u1 (3.2:) BS 7 $1Ag1=ICl[F GIIIPllBlMg} (13-22) 3‘18 -T a]: -T. $148}=J,(F, IBIIASI) (1324) Based on above relationships, the finite element discretization of the residual and stiffness terms can be obtained. To illustrate the general procedure, one residual term is transformed into finite element form below. Consider L(6u) : [ (V,6uF,7‘) : (F,F,7‘S)J,do, or n. (8.25) L(6u) : j (V,6uF,7‘)’(F,F,7‘S)J,do, see (A.17) 0. Since V ,6,uF7‘ =I(V, 6u)F,=7‘ -,[1@F7’][8]16u} see (A. 8) and (A. 12) (8.26) F,F,7‘s : [17, os]{F,7‘} see (A8) (8.27) Combining the above two equations and substituting in equation (B.25), we get the finite element discretization for the given residual. L(6u)= {611”} I [B] [IoF,7’] [F os]{ ,7‘}J, an, (8.28) where we have used the identity 11411811611’ =1C1’1B1’141’ Similarly, an example of discretization of a stiffness term is given below. Consider one of the stiffness terms of the above residual (6:1, Au): 1 (V 6a F,7‘): (9,51 1Au}(F,7 ‘5))J, do or (8.29) a(6u, Au): nj1V,6uF,’7‘) (63— u'1Au}(F,7 ‘S))J,do Since V ,,=6uF7‘ -,[1@F7’][B]{6u} (8.30) 3F _1 -l T —u—'a{Au}F, s: [101F, S} ][B]1Au} see (A8) and (B6) (8.31) 65 Combining the above two equations and substituting in equation (B.29), we get the finite element discretization for the given residual. a(6u,Au) = {6u°}’ j [B]T [IOF,'T]T [I o1F,7‘S}’][B] 1,612, {Au} (8.32) Qr 66 REFERENCES [1] Ghosh, S. (1990), “ Finite Deformation Analysis Involving Contact of Deformable Bodies Using the Finite Element Method with Arbitrary Lagrangian-Eulerian Description”, Winter ASME conference, pp. 273-280. [2] Haber, R. B. (1984), “A Mixed Eulerian-Lagrangian Displacement Model for Large Deformation Analysis in Solid Mechanics,” Comp. Meth. Appl. Mech. Engng., V01. 43, pp 277-292. [3] Ghosh, S. and Kikuchi, N. (1990), “An Arbitrary Lagrangian Eulerian Finite Element Method for Large Deformation Analysis of Elastic-Viscoplastic Solids,” Comp. Meth. Aappl. Mech. Engng., Vol. 86, pp 127-188 [4] Hibbit, H. D., Marcal, RV. and Rice, J. R. (1970), “A Finite Element Formulation for Problems of Large Strain and Large Displacement,” Int. J. Solids and Struc, Vol. 6, pp 1069-1086. [5] McMeeking, R. M. and Rice, J. R. (1975), “Finite Element Formulation for Problems of Large Elastic Plastic Deformation,” Int. J. Solids and Struc., Vol. 11, pp 601-616. 67 [6] Cheng, J. H. and Kikuchi, N. (1985), “An Analysis of Metal Forming Processes Using Large Deformation Elastic-Plastic Formulations,” Comp. Meth. Appl. Mech. Engng., Vol. 49, PP 71-108. [7] Liu, W. K., Belytschko, T . and Chang, H. (1986), “An Arbitrary Lagrangian- Eulerian Finite Element Method for Path Dependent Materials,” Comp. Meth. Aappl. Mech. Engng., Vol. 58, pp 227-245. [8] Balagangadhar, D. (1998), “A Reference Frame Formulation for the Analysis and Design Steady Manufacturing Processes”, PhD thesis, University of Illinois at Urbana- Champaign. [9] Brekelmans, W. A. M., Veldpaus, F. E. and Schreurs, P.J.G. (1986), “ Simulation of Forming Processes, Using The Eulerian-Lagrangian Formulation,” Comp. Meth. Aappl. Mech. Engng., Vol. 58, pp 19-36. [10] Dawson RR, and Thompson E.G. (1978), “ Finite Element Analysis of Steady- state Elasto-viscoplastic Flow by the Initial Stress—rate Method,” Int. J. Num. Meth. Engng., Vol. 12, pp 47-57. [11] Chan, SK. and Tuba, 1. S. (1971), “A Finite Element Method for Contact Problems of Solid Bodies, Part 1, Theory and Validation,: Int. J. Mech. Sci., Vol. 13, pp 615-625. 68 [12] Hughes, T.J.R., Taylor, R.L., Sackman, J. L., Cumier, A . and Kanokhukulchai, W.(1976), “ A Finite Element Method for A Class of Contact-Impact Problems,” Comp. Meth. Appl. Mech. Engng., Vol. 8 , pp 249-276 [13] Tseng, J . and Olsen, M. D. (1981), “The Mixed Finite Element Method Applied to Two Dimensional Elastic Contact Problems,” Int. J. Num. Meth. Engng., Vol. 17, pp 991-1014. [14] Kikuchi, N. and Oden, J.T. (1988), “ Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, “ SIAM, Philadelphia. [15] Stoker, C., (1998) “ Developments of the Arbitrary Lagrangian-Eulerian Method in non-linear Solid Mechanics.” [16] Liu, W.K. and Huerta, A. (1988) “ Viscous Flow With Large Free Surface Motion,” Comp. Meth. Aappl. Mech. Engng., Vol. 69, pp 277-324. [17] Liu, W. K and Hu, Y.K. (1992), “ ALE Finite Element Formulations For Ring Rolling Analysis,” Int. J. Num. Meth. Engng., Vol. 33, pp 1217-1236. [18] Zabaras, N. and Badrinarayanan, S. (1996), “ A Sensitivity Analysis For The Optimal Design Of Metal-Forming Processes,” Comp. Meth. Aappl. Mech. Engng., Vol. 129, pp 319-348. 69 [19] Grandhi, R. V. and Gao,Z.Y. (1999), “Sensitivity Analysis And Shape Optimization For Preform Design in Therrno-Mechanical Coupled Analysis ,” Int. J. Num. Meth. Engng., Vol. 45, pp 1349-1373. [20] Zabaras, N. and Srikanth, A. (2000), “ Shape Optimization and Preform Design In Metal Forming Processes,” Comp. Meth. Aappl. Mech. Engng., Vol. 190, pp 1859-1901. [21] Benson, D. (1989), “ An Efficient, Accurate, Simple ALE Method For Nonlinear Finite Element Method,” Comp. Meth. Aappl. Mech. Engng., Vol. 72, pp 305-350. [22] Balagangadhar, D. and Tortorelli, D. A. (2000) “Design of Large Deformation Steady Elastoplastic Manufacturing Processes. Parts 1 and 2,” Int. J. Num. Meth. Engng., Vol. 49, pp 950-950. [23] Gurtin, M.E., “An Introduction to Continuum Mechanics,” London: Academic Press, 1981. [24] Ghosh, S. and Kikuchi, N. (1988), “ Finite Element Formulation for the Simulation of Hot Sheet Metal Forming Processes,” Int. J. Engng. Sci., Vol. 26, pp. 143- 161. 70 111111111111111111111111111