-J ., .. v'.' . v 3i .c‘ .‘ u : t _.. - ‘Itzl ' ‘ ; ‘n;.vvgl, i ‘35»; f'n. ..‘1 “L. «<- -... 9;...‘3. —» A ' - x.- , J ‘ «I £:;2:: - " ‘ . .--:¢1£.3_ 3 In a p u . ’.~ 1 ' ,". 5 ~ Ix ~! ‘1 ‘W . u, ‘i‘t'z'. 1,“ We?!" ’ 526% fits} A {:5 Mi. lg. . . xr‘h‘z‘ ".‘ $.u.[ '2"‘ ll "a " 1" gt: .L "i: F" 1‘ .. “j ,b , grfirfiggq * ‘ . .tH ' fa": ~ - mac"? ‘3‘ ‘a : 4* 3:14?le ' ,35- 3 - P e’:v;'.;§;i§:m;$§§~m g m {91* 3 M i’ ' ’f‘-'2--u§= i .J'zfii ‘ : mow mu 4 , M 9m. :4 v. = - tap-'9 Tm iia4=3p¥913722i WM 1;; was; 45’? ._ “3:5: ‘ " W, fig! “LA \" afl-ZJ ugh: ‘, aw it?" " . ; 1&1? :m ‘ ’ :it Ii! 3; 1' u- ...- . '. ...”.w «.3 DIE} 1'3"?” wt: 3: 1.25:: .4, -~ :« p u‘ “‘ -.—4 w gr 4-..- ‘41,.- .~ .‘1 .1 ~ - . i ‘ am ‘0- ... . ,, V ,... A ‘ MW w.—— »-«vnfirm... " I- .. -..-m l 1 "m: “‘. ... w ....) ‘ ‘ :5 f. .1! _ a; 131139;?3'flifié t§§:3!53;fi'§9fi T 151; Ll‘xil‘QSQJEQé‘E‘ 215$} ‘ k z 1'. 2 1 . é 2 n 0 I waif .; A "' vi-‘éfimw '4: J0: If I? ml.“ u—- .u‘ : . .ds's: ......u.. =. {f 4 . 4......" ..‘..- u.‘ .¢,. f ‘2‘71' :,{;!i§:;'L\'g?;‘ §‘:§{~i‘ , “fixing ‘f h "“fl‘fit‘gb‘. 1'. ' ' 7' “f0? Y" . ‘ M x ' " I .13 1' -‘ téj :' ' ,‘V " Nah " $313. \ ......» , 4 ‘ w W '€:......- » ‘ u. ’3 '1 I L! .. 5 .:qu~ .- . - ....n- ’ a a; . . u v - ...... ...-.-- < < ' - . w .... .... ... ~ . M n - J24..un\~. ..‘ v .-s. g m - ... ...... 1‘ ..N. I ._ “1“2. I: i";:2 " '*,’ .f'ii" ‘ "$.22 ... _ m .., ' _. ’3‘: ‘ . .. "NM“ or ....- . ; " r ; ‘ Sig-1 ... @315? ' Aria ." ‘4‘ fir. 3?? ,5-.. a.“ ___ .. _ _ LC ‘1 .Mp— -m -.. 5M .- A a? .v _ +m-u— .: -.- . -......~ -..~."4>- — -' ”7:5; J: .": 1 . ...... - . I...“ " -..4. v.4..- .mn—n .4,» - ...- v ‘1 H9" ‘ W" .\ I? i x 1 Jr. .‘2' - m .. w: , ‘ wen-:7 .3; u >-.. '2..— - ,3 m ...p. "‘r2 vv~o—- ... .... - .« - -5\-O ’f‘vfifl‘: r-v .— .——. ' ~ , ' ,.: ' A.‘ . ' 0-140 ‘ ,_ 3353‘ w _'. . ‘ - > '- --'-:'2 l- .... _ _ . ' ,. .. n:_- U, . — on» > ' '. ‘ , . < fl.- ‘ 1" ._ . ...w... .. . -A r V.ov-. '3‘ -‘ - 'v-::;:=1-— -’ , ‘ - - 'l" w ' . '0.» 0‘10 .fl." .. ,4. ,..,.. .33"' ~g::;;:::*::. . v . u ... K A .2-.. . “ , .."" . ... ,5‘==5E53:.V .“ _-m.~ 5‘; 1.- ... ”-‘Zffl “‘L.“ THESIS 1K3 (F 3 C) '\._ tr: lelllllllllll' This is to certify that the dissertation entitled Solid Finite Elements for Sheet Metal ' Forming Simulation presented by Lorenzo Marco Smith has been accepted towards fulfillment of the requirements for Ph.D. degree in Engineering Mechanics Mam Major professor Gamay, f gouty Major professdr MM; my MS U is an Affirmative Action/Equal Opportunity Institution 0-12771 LIBRARY Michigan State University PLACE IN REFURN Box to remove this checkout from your record. To AVOID FINES return on or before date due. MAY BE RECAHED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 11/00 WWW." SOLID FINITE ELEMENTS FOR SHEET METAL FORMING SIMULATION By Lorenzo M. Smith A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Materials Science and Mechanics 1999 ABSTRACT SOLID FINITE ELEMENTS FOR SHEET METAL FORMING SIMULATION By Lorenzo M. Smith Two finite element formulations, suitable for sheet metal forming simulation, are proposed. The LNQS formulation features a linear and quadratic through thickness variation of the normal and shear strain, respectively. The CNQS formulation features a cubic and quadratic through thickness variation of the normal and shear stain, respectively. The LNQS model exhibits the ability to deform without the consequences of shear, Poisson’s or volume locking, while the CNQS model has been shown to exhibit some locking mechanisms due to an inconsistent shear strain field. The proposed elements exactly satisfy the shear strain (and stress) at the top and bottom of the element. These features, therefore, help qualify the elements to be used in sheet metal simulation procedures where shear strain is appreciable in the through-thickness direction. The LNQS in-plane strain due to both bending and membrane efi‘ects is accurately captured for plane strain metal forming cases with and without fi'iction. A system of evaluating and validating the elements, which involves the introduction of a FORTRAN subroutine within a commercial finite element software main program, has been established. In general, good correlation is found among the proposed model solutions and those found in literature. ACKNOWLEDGMENTS I begin by giving God thanks and praise for giving me everything that I needed to reach my goal. A special thanks goes to Dr. Ronald Averill, my advisor and friend. Thanks for your patience and guidance. I thank Dr. James Lucas, also an advisor and fiiend. Thank you for your words of encouragement and guidance. Dr. Srinivas Reddy, of MARC Analysis Corporation, very patiently helped me with technical issues. His dedication to customer support was remarkable. My friend, Dr. Ventateshwar Raoaitha, helped me in many ways. Thank you for your insights and helpful comments. Of course my best friend in the whole world, my beautiful and lovely wife of five years has been so patient and supportive. The only problem with finishing is that I can no longer use my research as an excuse not to do the dishes after dinner. Figure 2.1 Figure 2.2 Figure 2.3 Figure 2.4 Figure 2.5 Figure 2.6 Figure 3.1 Figure 6.1 Figure 6.2 Figure 6.3 Figure 6.4 Figure 6.5 Figure 6.6 Figure 6.7 Figure 6.8 Figure 6.9 Figure 6.10 Figure 6.11 Figure 6.12 LIST OF FIGURES AND TABLES Example of Typical Yield Surface in 2 Dimensional Principal Stress Space: von Mises Example of Isotropic and Kinematic Hardening Curves Elastic-Linear Hardening Model Graphical Description of Local Friction Geometry Idealizations of Two Major Friction Regimes at Die/Sheet Metal InterfacesIdealized Stribeck Curve Modified Friction Model with Cap Eight-noded Brick Element Topology Geometry for Axially Loaded Beam Load vs. Deflection for Axially Loaded Beam Geometry for Simply Supported Beam Locking Mechanisms Shear Strain on Single Element Cantilevered Section Under Transverse Tip Load Mesh Convergence Study: Refinement in X Mesh Convergence Study: Refinement in Z Load vs. Deflection for Transversely Loaded Beam (moderate bending) Load vs. Deflection for Transversely Loaded Beam (severe bending) Evolution of Material Constants in the Elasto-Plastic Material Matrix Material Property Sensitivity Study Page 35 36 37 38 39 4o 4s 87 88 89 91 92 93 94 95 96 97 98 Figure 6.13 Figure 6.14 Figure 6.15 Figure 6.16a Figure 6.16b Figure 6.17 Figure 6.18 Figure 6.19 Figure 6.20 Figure 6.21a Figure 6.21b Figure 6.22a Figure 6.22b Figure 6.22c Figure 6.23 Figure 6.24 Figure 6.25a Figure 6.25b Figure 6.26 Figure 6.27 Plasticity Solutions for Bending Gauss and Lobatto Points Integration Point Study for Plasticity Peak In-plane Strain for Elastic Case Peak In-plane Strain for Plastic Case Coulomb’s Friction Models: Stress or Force Based Nodal Friction Forces: Lee’s Benchmark Benchmark Geometry: Pinching Boundary Condition Transverse Coordinate vs. Transverse Normal Strain: Pinching Case Benchmark Geometry: Lee’s Stretch Forming Benchmark Final Mesh: Lee’s Stretch Forming In-plane Strain Solution: Lee’s Stretch Forming without Friction Benchmark Final Mesh: Modified Lee’s Stretch Forming In-plane Strain vs. In-planeCooediante: Modified Lee’s Adhesive and Cohesive Models Compared In-plane Strain Solution: Lee’s Stretch Forming with Friction Benchmark Geometry: Stoughton’s Stretch Forming Benchmark Final Mesh: Stoughton’s Stretch Forming In-plane Strain Solution: Stoughton’s Stretch Forming without Friction In-plane Strain Solution: Stoughton’s Stretch Forming with Friction 99 101 102 103 107 108 109 110 111 112 113 114 115 116 117 118 LIS'IWOFTABLES 2?: 5.1 Summary ofLNQS and CNQS Kinermtic Features 653'? (Iii/‘1'. . :3 , 1 , .. : . i ‘. I ' .r 2 CHAE'IF-i . > , ” ‘ ~ ' . . . , 31 '. ‘wr-mv "rtblmis W ' I} ‘ r M I ' 16 ” A as 1 ' I9 ’ W .2 2 1 £5 " " 26 CHAPTER ' i . ~ 38 3 - ‘ ‘i 39 . 1 ‘. ‘ . 4.; 'v CRAFTS! t “ A I f ‘3’ '3“ ' ,,' \fr : S; “ s. I \. .r. S; “if" murmur ” ‘- ' t . 57 ..,x , f 7 '1 1 I W ’3' Us Vi I; - CHAPTER 1 1.1 1.2 CHAPTER 2 2.1 2.2 2.3 2.4 2.5 2.6 CHAPTER3 3.1 3.2 3.3 CHAPTER 4 4.1 4.2 4.3 CHAPTER 5 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 TABLE OF CONTENTS Introduction Introduction Background Mathematical Preliminaries for Modeling Metal Forming Problems Implicit and Explicit Methods Stress/Strain Measures Newton-Raphson Solution Procedure Material Model 2.4.1 Yield Surface 2.4.1 Flow Rule 2.4.3 Hardening Law 2.4.4 Elasto-Plastic Material Matrix Contact Model Friction Models LNQS Element Kinematics Displacement Field Shear Strain Field Normal Strain Field CNQS Element Kinematics Displacement Field Shear Strain Field Normal Strain Field Finite Element Formulation Comparison of Features BL and BM Matrices Direct Stiffness Matrix Tangent Stiffness Matrix Finite Element General Solution Procedure Shear Traction Definitions Shear Locking Volume Locking Poisson’s Locking Page l3 14 16 17 18 19 21 25 26 38 39 43 49 54 55 57 57 58 60 6O 61 62 63 63 CHAPTER 6 Benchmark Results 6.1 Cantilevered Section Under In-Plane Tip Load 6.2 Simply Supported Section Under Uniform Transverse Load 6.3 Transverse Shear Strain on Single Element 6.4 Cantilevered Section Under Transverse Tip Load 6.5 Material Model Study 6.6 Preliminary Commercial Element Investigation 6.7 Lee’s Benchmark: Plane Strain Stretch with Zero Friction 6.8 Lee’s Benchmark: Plane Strain Stretch with Non-Zero Friction 6.9 Stoughton’s Benchmark: Plane Strain Stretch with Zero Friction 6.10 Stoughton’s Benchmark: Plane Strain Stretch with Non-Zero Frict 6.11 Thick Section Benchmark: Pinching Boundary Condition CHAPTER 8 Conclusions CHAPTER 9 Future Work BIBLIOGRAPHY APPENDIX A Commercial Code Utilization 120 125 126 132 Chapter 1 Introduction 1.1 Introduction Historically, the process of forming sheet metal into useful tools and objects has been carried out by highly experienced and accomplished craftsmen. The central aim of their task is quite straight forward; they must economically deform the metal into firnctional objects having an appearance pleasing to the eye. However, this process is not always straight forward. For centuries, craftsmen have successfully carried out this mission without a mathematical understanding of the mechanics associated with sheet metal forming. The spectacular work of these craftsmen was, therefore, attributable to their intuition based upon decades of trial and error. Because of this “sixth sense”, their work is more accurately and appropriately defined as an art rather than a science. Indeed, they were artists. Even today, elements of both the art and science of sheet metal forming coexist. In today’s automobile industry, for example, engineers draw upon the expertise of the tool and die makers to help direct the design of automobile components. For many years, the marriage of science and art in the sheet metal forming industry has been formally recognized by many as a necessary relationship. Where science fails, art must prevail; where art fails, science must prevail. Today, however, there is an enthusiastic campaign to rely exclusively upon the use of computer models, to determine if a particular sheet metal forming process can be successfirlly employed. Millions of dollars and decades of research have been devoted to this goal. Why? Although at first glance it may seem unreasonable to tamper with the success of the past, it is necessary. The middle 20th century industrial climate brought forth some harsh facts. First, there was less interest from young novices in the art of sheet metal forming. Thus, the pool of experienced craftsmen was not being replenished. Second, there was less lead time for production than there was in the past. Consequently, the artists had less time to find the optimum design parameters. Third, because of cost, weight and safety constraints, new materials were being introduced more ofien. The artist’s expertise in successfirlly forming a particular metal was acquired over years of trial and error. This was no longer practical in the increasingly fast paced industrial climate of that time. By the mid 1950’s, industry slowly began to usher aside its once seemingly immovable inertia of metal forming practice. A new era was born. Today, this new era is a witness to industry’s dependence upon computer simulation of sheet metal forming. Since the birth of this new era, there has been an intensifying effort by the research and industrial community to improve the capability of modeling the metal forming process. An understanding of the fundamental concepts of the micro and macromechanics associated with metal forming has been realized for some time. However, the computing tools, for many years, were not capable of efficiently accommodating the size and complexity of many models. It was only since the early 1970’s, with the advent of more sophisticated computers, that scientists were able to simulate relatively complex sheet metal forming processes in a more practical manner. Procedures such as the finite element method then became more common avenues to solutions. Researchers took full advantage of the available tools of the time to model the sheet metal forming process. Today, with the availability of workstations and supercomputers, new horizons have been opened to scientists for exploration. Nonetheless, the envelope of this new research domain has been pushed to its limit once again. As the pattern over the past several decades suggests, researchers find the complexity and size of their models being dictated by the availability of computing power. The real challenge to researchers is to develop a theory that can be cast into a computer model form that is both accurate and practical. 1.2 Background The finite element method is the most popular and arguably most powerfirl means of modeling the sheet metal forming process. Because of the kinematic assumptions used in the formulation of a finite element, it has been that finite elements should be classified as structural elements by definition. However, for the sake of discussion, the names of three categories of elements will be used in the vocabulary of this study; they are membrane, shell and three-dimensional (3D). Starting in the late 1970’s, membrane elements emerged as the first industrially usefirl elements for sheet metal forming simulation (Wang, 1978, Toh, 1985, Nakamachi, 1988, Huang, 1994). VVrth the exception of super-plasticity sheet metal forming (Argyris, 1984, Bellet, 1987, Bonet, 1990), the success of the membrane elements has been restricted to thin cross section applications because of the elements inability to model strain changes due to bending. In an effort to model bending effects while simultaneously preserving the simplicity of the membrane element, Yang (1995) proposed a bending- energy-augmented-membrane (BEAM) element. By introducing a rotational energy expression in terms of only translation degrees of freedom, Yang was able to introduce an element having the appearance of a membrane element, yet possessing the ability to simulate resistance to bending. The computational results given in Yang’s paper were used to compare difi‘erent solution algorithms; the BEAM element was essentially a “dummy" element. Accordingly, a rigorous assessment of the BEAM element’s performance is not available. In introducing a shell element formulation, Wang (1987) exposed the short comings of the membrane element with regards to sheet metal forming applications. Wang explained and demonstrated the membrane element’s inability, and the shell elements ability, to model variations in the in-plane strain with respect to the through-thickness coordinate. Many other shell element contributions have been made; those introduced by Belytschko (1981), Lee (1991) and Boubakar (1996) are just a few to be mentioned. These elements have proven to be quite useful in modeling a wide variety of sheet metal forming processes. In particular, Belyschko’s element has been very successful because of it’s combined efficiency, accuracy and robustness. In spite of the many advantages of the shell element, there are still many who are not satisfied with the shell for a variety of applications of sheet metal forming; the traditional shell cannot model double sided contact, for instance. Rebello (1990) and Gontier (1994) point out needs for explicit and accurate modeling of thinning in order to improve the workpiece/tool contact model. So, the relentless call for more improved element performance continued. In light of the double sided contact concern, one inclination may be to consider a 3D element. Although, the 3D element tends to feature more degrees of freedom than does the shell element, the 3D element advantages in many cases may justify the added expense of more degrees of fi'eedom. One challenge is to find a 3D element that comparably performs with the shell without using multiple 3D elements through the thickness of the moderately thick sheet metal domain. Using an enhanced assumed strain formulation, Sirno (1990,1992,l993) was able to overcome that challenge. Korelc (1995) introduced a 3D element for small strain and thin cross sections. Soon after, Wriggers (1996) unveiled a 3D element which shows great promise for thin cross sections and large strain conditions. A variety of other 3D elements proposed for sheet metal forming simulation has been suggested by others (Massoni, 1989, Onate, 1990, Oh, 1980, Butcher, 1994, Pian 1984,). It should be obvious that there are advantages and disadvantages to the use of shell or 3D elements. Shimizu (1991) may have been the first to attempt to combine many of the advantages of the membrane, shell and 3D element into one element suitable for sheet metal forming. The formulation features a four noded membrane element having the outward appearance of an eight-noded brick. The brick topology is introduced in order to enhance the contact model and to model thickness changes more accurately. The idea is very similar to that proposed in the work of Averill; that is, the structural element is disguised as a solid element. A trilinear interpolation is assigned to each transverse degree of fieedom. This most simple interpolation does not allow the element to model a quadratic shear stress variation through the thickness as accurately as some of the previously discussed elements. Furthermore, because the dilatational energy is integrated independently of the deviatoric energy, the incompressibility constraint is enforced by the penalty method. Therefore, to calculate the pressure required to maintain the incompressible behavior, an additional calculation involving the divergence of the velocity needs to be made. Finally, Shimizu’s solution is obtained via an explicit time integration scheme in conjunction with an updated geometry approach. Shimizu’s numerical results were in qualitative agreement with those from experiment. In the spirit of Shimizu’s intentions, a new thrust has been made in order to introduce an even more attractive element for sheet metal forming. A quest for improved shear stress and normal strain accuracy in a solid element suitable for metal forming using a total Lagrangian, implicit model is the primary goal of this thesis. A secondary focus is to obtain the primary goal by establishing a software infrastructure that can be used to study new finite elements for metal forming applications. To this end, an excursion into the area of finite elements for composite elements is taken. Not unlike finite elements for sheet metal forming, overcoming the challenges with finite elements for composite structural analysis is also a formidable task. When material delamination is to be modeled, a nonlinear material model must be developed for elements for composite analysis. For certain applications such as aircraft skin behavior under flight conditions kinematic nonlinearities must also be taken into account. Additionally, when transverse squashing occurs in composite structures, a nonlinear contact condition may possibly arise within the laminate structures of the body. For composite elements, accurately modeling the discontinuous in-plane normal stress and piece-wise continuous transverse shear stress, while maintaining some reasonable level of formulation simplicity, is a monumental challenge. These challenges have been addressed by many (DeSciuva 1987, Reddy 1984, Averill 1996, Cho 1997, Aitharaju 1997). One issue that is addressed in the work of Averill is that of Poisson’s locking. The remedy used to eliminate Poisson’s locking in composite elements is the ideal remedy to prevent volume locking for metal forming elements. Details of the locking behavior will be discussed in detail in Chapter 5. As indicated above, a critical capability for most elements in composite analysis is that which provides an accurate simulation of the shear stress through the thickness. By introducing shear traction degrees of fi'eedom at the nodes, Averill satisfied the shear traction exactly at the top and bottom of the element, while allowing a higher order variation through the thickness. This feature is remarkably well suited for metal forming applications because of the complex fiictional boundary conditions imposed upon the sheet metal by the forming tools. Furthermore, another consequence of the remedy used for Poisson’s locking is a form of the transverse normal strain which is useful in modeling the efi‘ects of squashing. Details of this feature are given in Chapter 5. It can be legitimately argued that the development of metal forming elements is a natural extension of the development of the type of composite elements developed by Averill. Averill’s work in composite elements has served as a catalyst for the introduction of a more accurate description of the transverse shear strain and transverse normal strain terms in the context of a 3D element for sheet metal forming applications. Two new finite element formulations will be proposed in this thesis. Both elements have only the standard translational degrees of fieedom at each node. Accordingly, the proposed elements lack some of the kinematic sophistication found in many shell elements, as discussed in the previous sections. In the early development of the proposed elements, a special form of shell kinematics was considered. However, the price for including shell kinematics was unacceptably high; additional rotation degrees of freedom had to be included, and some form of constraint (either penalty or Lagrange multiplier) used to relate the transverse displacements and rotation terms had to be imposed. In an effort to present the proposed elements in the most convenient, unintimidating and useful way, it was decided that the shell kinematics would not be included. The first proposed element features normal strain terms having a linear variation in the through-thickness direction and transverse shear strain terms having a quadratic variation in the through-thickness direction; this element will be referred to as the LNQS (linear normal, quadratic shear) element for the remainder of this thesis. The LNQS element is most suited for thin to moderately thick cross sections because of it’s linear through-thickness variation of the in-plane normal strain. It is worth noting that shell elements with the same through-thickness variation of the in-plane normal strain have been used extensively and successfirlly in automotive sheet metal forming simulation. The quadratic through-thickness variation of the transverse shear stress is necessary in order to exactly satisfy the shear stress at the top and bottom of the element. As will be discussed in Chapter 2, the shear stress is ofien included in the yield function. In thin to moderately thick sheet metal forming analysis, the shear stress is generally several orders of magnitude less than that of the normal in-plane stress; the exception is when the curvature is extreme or when double sided contact conditions impose severe shear tractions. In such cases, an accurate shear stress model is important. The second element features normal strain terms having a cubic variation in the through-thickness direction and transverse shear strain terms having a quadratic variation in the through-thickness direction. For the remainder of this thesis, the second element will be referred to as the CNQS (Cubic Normal Quadratic Shear) element. The derivation of the CNQS element formulation involves an initial assumption associated with shell kinematics. However, as Chapter 6 reveals, a critical omission of the rotation terms allows the CNQS formulation to take the form of a solid element with only translational degrees of freedom. With a through-thickness cubic variation of the in-plane strain, the CNQS element has a more realistic in-plane strain model than that of the LNQS element. Like the LNQS element, the CNQS element also can model the shear stress with more accuracy than the typical thin shell element. Some preliminary information needs to be presented prior to the introduction of the proposed element formulations. The three nonlinearities present in most sheet metal forming operations (ie. kinematic, material and contact) are discussed in a general sense in Chapter 2. Care is taken to review how the nonlinearities arise and how they can be handled in a finite element model. Details of the two major integration techniques (the implicit and explicit method) are discussed. Chapter 2 also includes a brief discussion of contact and fiiction models. The contact nonlinearity has two components. The first is the normal component. As the sheet metal slides over the tool, portions of the sheet metal come in and out of contact with the tool. This change in boundary condition obviously afl‘ects the normal stress, but also the shear stress if fiiction is non-zero. Therefore, the second component of the contact 10 nonlinearity arises; this second component is due to fiiction. Consideration of the appropriate models is made in Chapter 2. An element formulation is usefirl in the industrial environment when it is cast into a finite element model and used in finite element software package which offers a powerfirl pre and post processor. To this end, the LNQS and CNQS finite element models have been defined in a FORTRAN subroutine which is called by the commercial finite element software MARC. Miscellaneous documentation is included in Appendix A In Chapter 3 and 4, details of the LNQS and CNQS finite element derivations are provided. The elements share the same topology and degrees of freedom. Only the interpolation of the degrees of freedom is different among the two elements. In addition to defining the finite element models, Chapter 5 is used to highlight the issues of locking mechanisms. The root of the potential shear, Poisson’s and volume locking will be exposed and eliminated. A thorough review of locking is essential in most any new finite element formulation. Recall that the entire model is based upon assumed displacement (or stress/strain) fields within the element. These assumed fields are a manifestation of the necessary compromises that the scientist makes in order to mimic reality more precisely. In the end, the irnposture of reality will always be convicted of some weakness. Many locking weaknesses inherent to new element formulations may be avoided by either underintegrating the element or strategically manipulating certain strain fields. This is tantamount to mathematically “turning one’s heads to the problem”. Nonetheless, the literature has shown in many cases that this works quite well! ll The benchmark problems are introduced in Chapter 6. A deliberate and methodical assessment of the element is made by first considering simple academic cases, such as cantilevered sections under transverse and in-plane loads. Much information can be extracted fi'om such a study, as will be shown in the results. Some simple academic double sided contact cases are also studied. In these studies, the fiiction model is assessed. Finally, a plane strain sheet metal forming simulation with and without fiiction is carried out. The LNQS and CNQS solutions will be compared to MARC and experimental data. The thesis will be concluded with a summary of the results and future work in Chapter 7 and 8, respectively. Chapter 2 Mathematical Preliminaries In general, solving a set of nonlinear equations in the context of a finite element model requires three basic steps. The first is to, upon defining the equilibrium equations, specify the form of integration necessary to obtain the desirable solution (usually displacement). There are many questions that need to be answered when approaching this first step. Are the material properties deformation rate dependent? If so, can this dependence be ignored? Are inertial effects involved? Are there any heat transfer considerations that need to be made? How important is accuracy? How important is speed? Can the nonlinear terms be formulated in a practical manner? Answers to all of these questions can be used to help determine the appropriate technique of integration of the equilibrium equations. Direct integration procedures are common in many finite element programs. The basis of direct integration is two fold. First, a solution to the equilibrium equation is sought for a specific time interval. Therefore, within each time interval, the objective is simply to solve a static equilibrium set of equations which may include inertia and damping effects. Second, a variation of the displacement, velocity and acceleration within each time interval is assumed. The assumptions of such variations dictate accuracy, stability and efiiciency. Two types of direct integration methods will be discussed in subsequent sections. Once the integration of the equilibrium equations is established, a reference fi’arne must be defined. Again many questions need to be answered. Does the problem involve large strain or large rotations, or both? Does the deformation fit the mode of fluid or 13 solid? Once such questions are answered, a determination of the model formulation can be made. The equilibrium equations can be defined in the initial flame or current frame. These considerations are discussed in subsequent sections. Finally, once the integration and frame of the equilibrium equations are defined, then a solution procedure must be established. Since the equilibrium equations are generally nonlinear in nature for metal forming problems, an iterative sohrtion scheme usually is necessary. The firll Newton-Raphson solution scheme is a common method. This method is used in the currently proposed model. Chapter 2 will include more detailed explanations of direct integration, reference frames, Newton-Raphson iterative solution procedure and also the nonlinear aspects of the material and contact models. 2.1 Explicit and Implicit Methods For the explicit method, a solution to the displacement vector at time, t + At, is sought based upon the equilibrium conditions at time , t. Ifthe central difference method is assumed then equilibrium equations take on the following form: where: A7 , 5 , I? are the mass, damping and stiffiress matrices, respectively. (7 and E are the displacement and external load vectors, respectively. For the implicit method, the solution for the displacement vector at time, t + At, is sought based upon equilibrium conditions based upon time, t + At. If the Houbolt method is assumed then equilibrium equations take on the following form: (fifi+%6+fi)fi”m = Rt+Ar +[i21g+_3_@)0" _(—42—1t7+—3—-C‘)l7'—A‘ +[—i2—A~l +—l—C)l7t_2N 2.1-2 At A! A; 2A: Ar 3A! It is noted that in the implicit method, 2.1.2 can be utilized even when the damping and inertia matrices are neglected. Two schools of thought dominate the issue of which solution (from implicit or explicit means) is superior. One school of thought embraces the implicit form because, from an equilibrium standpoint, it is more reliable and rigorous at each step. However, convergence is not always guaranteed. The other school of thought supports the explicit form because, in spite of the fact that it is less rigorous from an equilibrium standpoint at each step, it has much more favorable convergence properties, provided that the appropriate time step is assumed. Ultimately, the nature of the application should drive the decision. Yang (1995) has shown that for various cases, both approaches provided comparable solutions. 2.2 Stress/Strain Measures As mentioned previously, there is an array of fi'ames from which the stress and strain (or equilibrium) can be defined. One side of the array is the spatial description, where the stress and strain measures are defined with respect to the current configuration. One advantage to the reference of the current frame is that because the geometry of the element is updated, some of the higher order terms in the strain measures nwd not be 15 included. One form of this approach is called the “updated Lagrangian method”. Additionally, the stress and strain terms correspond to the true stress and strain. The spatial description of the equilibrium equation is given as aw +f, = 0 2.2.1 where: 0,]. is the true or Cauchy stress acting on the deformed body under external traction. f, is the body forces vector The energy conjugate to the Cauchy stress is the Almonsi strain. On the other side of the array is the material description. In this description, the stress and strain measures are defined with respect to the original configuration. Solution accuracy can be achieved if the appropriate higher order terms in the strain tensors are defined. The stress and strain correspond to the engineering stress and strain. Typically this approach is called the “total Lagrangian” approach. Care must be taken when comparing the computed solutions from that of the total Lagrangian model to experimental data because of the lack of physical correspondence between the mathematics and practice. The equilibrium equation is given as 7w +f, =0 2.2.1 where: r, is the 2nd Piola Kirchhoff stress acting on the undeformed body The energy conjugate to the 2nd Piola-Kirchhoff stress is the Green-Lagrange strain. The total Lagrangian method is used for the current model. 2.3 Newton-Raphson Solution Procedure Defining the Newton-Raphson solution procedure for a nonlinear function with one independent variable is straight forward. However, the Newton-Raphson solution procedure, for a functional can be somewhat less straight forward. Therefore, it is worth the time to begin the introduction of the general Newton-Raphson method by considering the following functional. F(e)=3(ic'o+ er7)=0 2.3.1 where: 3 is an arbitrary fimctional such as an expression for the potential energy 2,, is an arbitrary position vector in space 17 is an arbitrary displacement vector in space e is an artificial parameter used as a vehicle for difi‘erentiation F is some function of e A Taylor series expansion of F about e=0 is given as dF 1 d’F F(E)= F(0)+fie=o E+§EL:O 62 +.... 2..32 Using 2.3.1 in 2.3.2 and preserving only the first order terms yields ~ ~ ~ d ~ ~ Sixo+ eu)-3(xo)~ed—J 2(xo+ eu) 2.3.4 5:0 The right hand side of 2.3.4 is called the directional derivative of 3 at x0 = 0 in the direction of i7 and is written as 0307 )[a]~ei{ :(r + at?) 2.3.5 o d e no ° Setting 6 equal to unity and using the notation of 2.3.5, 2.3.4 becomes 17 :(SEO +r7)z 21(20)+D3('fo)[z7] 2.3.6 Then fi'om 2.3.1, 2.3.6 becomes DSKYOXE] = {(350) 2.3.7 From the above equation, the general Newton-Raphson procedure can be expressed as DSKF,)[17,H] = —21(Y,,); in, = if, +17, 2.3.8 In the context of a finite element model, the terms on the lefi hand side of the above equation are identified as the tangent stiffness matrix and displacement vector. 1 2.4 Material Model For sheet metal forming analysis, the material model is nonlinear in nature. It may vary with stress/strain magnitude, strain path, temperature and many other factors. This section will introduce the three components of the material model used in the LNQS and CNQS element models; they are yield surface, hardening law and flow rule. 2.4.1 Yield Surface The yield surface is a fimction, commonly expressed in stress space, which defines the boundary between the elastic and plastic deformation modes in a material. Most yield surfaces are described by simple polynomials as shown in Figure 2.1. The appropriate order of the polynomial for sheet metal forming analysis is debatable. For an excellent account of how the shape of the yield surface afi‘ects the forrnability of sheet metal, the reader is referred to Barlat ( 1987). The most common yield surface used in industry for sheet metal forming simulation is Hill’s 1948 anisotropic surface (Hill, 1948). Although 18 there is little experimental evidence that the Hill 48 model is accurate (Hosford, 1993), it is used extensively in many finite element programs. In fact numerous studies (Barlat, 1991, Bramley, 1978) have led to the conclusion that the Hill 48 yield surface is inadequate for many sheet metal forming applications. Barlat, in fact, presents an argument which supports the use of a higher order yield surface, referencing solutions produced by polycrystal plasticity models. In contrast to such perspectives, Stoughton (1997) presents a formidable case which, indeed, supports the use of the Hill 48 model if (and only it) the yield surface is not assumed to be equivalent to the plastic potential. One of the more simple yield surfaces, the von Mises yield surface, will be used in this study. The von Mises yield surface, F, which is given below, is limited to isotropic material conditions. 1 F =|:%(0'1—02)2 +%(0‘2 —0'3)2 +%(0'3 —o,)2 +30: +30: +30%]: —a, 2.4.1.1 where subscripts 1,2,3 refer to the normal stress subscripts 4,5,6 refer to corresponding shear stress and a, is the yield strength in uniaxial tension 2.4.2 Flow Rule The “flow rule” is the expression that governs the relationship between the strain and the stress state during plastic deformation. It may be derived through plastic work considerations (Mendelson,1968) or by the introduction of a plastic potential function (Melan,1938). 69(3) f=—d1 2.4.2.1 g as 19 where: Q is a plastic potential fimction of stress, 5 d). is the magnitude of the plastic strain increment, d5 . Ifthe plastic potential and yield surface are equivalent, then the flow rule is characterized as “associated”. Bland (1957) showed that the plastic potential and yield surface must be the same function from a theoretical standpoint. However, Stoughton (1997) poignantly resurrects the old disclaimer put forth by Hill in 1948; there is little or no experimental evidence which shows that the yield surface and plastic potential are necessarily the same function! Yet, many inconsistencies between experimental data and analytical sohrtions based upon the Hill 48 theory have been attributed to the inadequacies of the Hill 48 model. Much of the blame has been directed towards the Hill 48 yield surface while little attention has been paid to the fact that an associated flow rule may not be appropriate. Stoughton lobbies that if a certain non-associated flow rule is used, then the alleged Hill 48 inconsistency would not be seen. Hence, credence to the Hill 48 model for sheet metal forming simulation can be established, according to Stoughton. What argument is one to accept? This issue is recognized. However, it is not the primary focus of this dissertation. For convenience the associated flow rule in conjunction with the Hill 48 (von Mises model is a special case of the Hill 48) model yield surface will be used as a basis in a plasticity model for the proposed element. 2.4.3 Hardening Law Of course, the mechanics of plastic flow can be precisely described by referring to the micro-mechanics or crystal mechanics of the metal. In this arena, such concepts as 20 crosspslip, latent hardening, twinning, and preferred slip systems come into play. A detailed discussion of polycrystal mechanics of sheet metal is beyond the scope of this study. Nonetheless, it is worth noting that the interaction among the crystals of the metal as plastic defamation occurs dictates not only the shape of the yield surface, but the evolution of the yield surface centroid as well. A description of the yield surface evolution is more commonly referred to as hardening of the yield surface. As illustrated in Figure 2.2a and 2.2b, there are two major types of hardening. The first type, kinematic hardening, involves a “rigid body type” of change in the yield surface in stress space. The second type, isotropic hardening, involves a “dilatational type” of change in the yield surface in stress space. Neither type is able to closely capture details of the actual hardening behavior of the metal. The true hardening description is likely some combination of the two as Dafalios (1982) points out. For many sheet metal forming problems, the elastic strain response is approximately one order of magnitude less than that of the plastic response. The following rough vahies for modulus and yield strength are given for example (Callister, 1997): Modulus (MPa) Yield Strength (MPa) Aluminum 69,000 400 Steel 207,000 1,500 From the above values, a simple calculation shows that the elastic strain for aluminum and steel is approximately 0.5%. This is a significant amount of strain when it is noted that the maximum strains generated in the current study are about 10%. One consequence of a rigid plastic model, therefore, is that it will tend to overestimate the compliance of the structure. The current proposed material model assumes an elastic-plastic response. 21 Strength increase due to hardening can be significant for many metals. Therefore, an overestimation of the structural compliance can also be attributed to the zero-hardening assumption if a rigid-perfectly plastic model were to be used. For this reason, a linear hardening model is assumed for the current material model. In theory, the slope of the tangent modulus can be varied from zero to the initial elastic modulus in numerical models. The transition from the elastic to plastic regime is, in general, smooth in stress strain space in metal forming practice. On the contrary, the proposed elastic-plastic numerical model imposes a rather sharp transition. This is particularly true when the tangent modulus is defined to be less than, say, one half that of the elastic modulus. In such cases, convergence or accuracy problems are likely. Numerous techniques have been proposed to circumvent such numerical dificulties. Cook (1989), for example, suggests a “corner rounding” technique where the current modulus is some function of the previous modulus. Such features have not been implemented into the current model. For this reason, some convergence difiiculties have been observed for cases of very sharp elastic- to-plastic transitions. Figure 2.3 describes a theoretical stress-strain model for a uniaxially loaded member. 2.4.4 Elasto-Plastic Relation Using the flow rule and yield surface, the stress may be related to the strain in a convenient form which resembles the linear stress-strain relation. The difference is that the components of the material matrix are fimctions of stress and hardening parameters 22 (Zienkiewicz, 1969). To derive the elasto-plastic matrix, the total strain increment is simply defined in terms of its elastic and plastic incremental components as follows: f, = d2, +dEp 2.4.4.1 Next using the elastic and plastic strains, the total strain is defined in terms of an inverted linear material matrix and the flow rule, respectively. d”, = are + gar 2.4.4.2 In 2.4.4.2, the plastic potential function, Q, is assumed to be the yield function, F. (This is a statement of the associated flow rule.) For plastic deformation, the yield function, F, is defined to be zero. Therefore, :5 rs equal to zero. Accordingly, we have in addition to 2.4.4.2, d? (T a7 a: ——d +——d +—d +...——d =0 50, 0'1 50,2 0'2 $3 03 6k, Kr where: x is a hardening parameter which represents a change in stress due to either kinematic or isotropic hardening. Ifno hardening is assumed, then x is set to zero. or a T {E} d{3} + AA = 0 2.4.4.3 67 l h :A=—d-— were 0hr). Equations 2.4.4.2 and 2.4.4.3 can be combined and written in matrix form as (for two dimensions) [ l .01- 50. d5, ~_I 01: d0, d8 D I a; do 2 = 2 z 2444 d8” I 57 dan 0 __ _ _ a"12 Z 2: or or A pa. an 50.2 J Inverting 2.4.4.4 yields d5=5¢dE 2.4.4.1 where: —1 .. ~ ~ a? a? ’ 5!? ’~ 6F C = -C— — A — C— 2.4.4.2 ... C {eiiaeifl +{65} {45” where: C is the elastic material matrix A is the plastic modulus of a uniaxial stress-strain curve For a two-dimensional case, 2.4.4.2 can be expanded as follows: C.. .2 o F... F... ' C.. C” o C C 0 F,I F,‘ C C 0 C" C" ° 5‘ n c F F 3' 32 c 60),: C21 C22 0 " 133 m m 33 2.4.4.3 o 0 C33 ’1: C11 C12 0 Fax 0 F F,z C2, C22 0 17,, +A F F 24 where: the subscripts x,z,and xz refer to differentiation with respect to the corresponding stress component. and for a two dimensional case without shear stress consideration in the yield function: 1 [an — 5(022 + 033)] sx_ , s;— afield [an—aw] 0' yield Expanding 2.4.4.3 yields: (p C121(F’x)2 +C122(Fu)2 +2C12C11F’r Fix C" =C"_ C F +C F +2C F F +A 2'4'4'4 11 ’x 22 ’2 12 ’x ’x ., C3.(F..)’ + C3.(F,.)’ + 26..C..F,. F... C” =C”_ C F +C F +2C F F +A 2'4'4'5 11 ’x 22 ’2 12 ’x ’z C11C12(F9x)2 +Cl21F’x Fax+C22C11sz Fix+C12C22(F!x)2 2 4 4 6 C" = C - ‘2 '2 C,,F,,+C,,F,,+2C,,F,, F,,+A 0;: = cg 2.4.4.7 cg = C,, 2.4.4.3 With the assumption of a von Mises yield criterion, Bathe (1982) conveniently expressed the same elastoplastic material matrix as follows: l—v l-V "flo'izi ‘flaiiazz 0 E 11-2v 1—12v ~ -v -V '=1__1_2v'”""°“ 3"”? ° ”4'9 l 0 O E‘flaiz 25 31 1 Where. fl— 5'0—32” Tm 1+ ~—— 3E—A E 2.5 Contact Model Perhaps the most important reason why the MARC subroutine procedure was implemented, is the availability of the contact model. MARC imposes a non-penetration constraint which is given in it’s general form as fl . a s D 2.5.1 where: (7 is the nodal displacement vector of the deformable body (sheet metal) ti is the normal vector of the of the rigid body (die or punch) D the distance between the node and the rigid surface. Within the context of a finite element model, this constraint can be imposed by several methods; three are the Lagrange multiplier method, the penalty method and the solver constraint method. MARC uses the solver constraint method. MARC will determine which deformable body nodes are close enough to the rigid body nodes to be considered candidates for penetration. Whenever the constraint is not satisfied for a given load step, MARC will impose a corrective displacement (not force) upon the deformable body node. Implementation of the proposed elements into MARC, in order that the contact model can be used is, one of the fruits of the current research. 26 2.6 Friction Models In the field of sheet metal forming, various compromises in modeling detail are made, ofien in the name of computational cost. Arguably, the most glaring example of such a compromise is the typical fiiction model used in commercial codes for industrial use. Two of the most prominent types of fiiction models noted in literature are the cohesive and adhesive models. The cohesive (Coulomb) model defines the tangential force to be a fi'action of the normal load. F, = We“, 2.6.1 where: u is the coefficient of friction (need not be constant) N is the normal force 2, is the direction of the force which is determined by the relative sliding velocities of the work piece and tool. The adhesive model defines the tangential force to be a fraction of a shear yield strength. f. =1"? 6 2.6.2 y t where: m is the coefficient of fiiction (need not be constant) I, is the shear yield strength of the material being formed In order to appreciate why an accurate and eflicient fiiction model has been so elusive, a study of the local contact mechanics is in order. From Figure 2.4, it is shown that there are two load carrying devices at the sheet metal/tool interface; they are the lubrication film and the metallic asperity peaks. ‘- 27 Accounting for the local contact conditions is not a trivial matter (Ronda, 1996). On one hand, the normal load can be totally carried by the lubrication film, if the film is relatively thick. On the other hand, the normal load can be totally carried by the asperity peaks if the peak heights are relatively large. Yet another scenario is when the normal load is carried by both the lubrication film and asperity peaks. If, indeed, the normal load is carried by the asperity peaks, then what is the slipping mechanism? In other words, will the peaks shear through or slide over each other? If shear is the mechanism, then which peaks will fail? To answer this question, both the material property and local geometry must be considered. If sliding is the mechanism, though, then the elastic response of the apserities and local normal forces may become more important. Another issue is that of relative velocity among the apserities of the sheet metal and tool. At low velocities, much of the load may be carried by the asperities. At high velocities, the load is more likely to be carried by the lubrication film. Carleer (1996) has suggested a dimensionless lubrication number that is useful in providing helpfiil insight. L = 1". 2.6.3 where: n, v, p and R are the dynamic viscosity of the lubricant, the relative sliding velocity, the mean contact pressure and the efi‘ective asperity height, respectively. From Figure 2.5, three main regimes which dictate the value of the fiiction coeficient are noticed. In region A and B, the normal load is carried by the asperity peaks and lubrication film, respectively. Region C is considered a mixed regime where both the asperity peaks and lubrication film carry the normal load. 28 In a finite element model, Liu (1994) introduced a variable coeflicient of fiction. Liu obtained good agreement with experimental data. Liu’s model fiirther suggests dramatic changes in the fiiction coefficient are produced, in response to changing boundary conditions. Schweizerhof (1991) has proposed a modified form of the classical Coulomb friction model (Figure 2.6). A limit coefficient has been established in order to more realistically model the local contact conditions of fiiction. Many types of models have been proposed (Ronda, 1996, Carleer, 1996, Wilson, 1988, Nagtegaal, 1988, Scheizerhof, 1991, Liu, 1994). Some models are based upon hydrodynamics, while others are based only upon shear yield strength. In spite of the wide array of approaches, nearly all proposed models were driven by the apparent inadequacies of a Coulomb’s model with a constant coeficient of friction. There seems to be no singular approach that works satisfactorily for all cases. Careful discernment of the mechanics of the particular forming problem should be made. In addition to the traditional Coulomb’s model, Nagtegaal has also endorsed a simple adhesive model. Additionally, Wilson has proposed that until a more universal and efficient model is found, a simple adhesive model is reasonable for general applications. Having obtained an approximation for the friction force at the element surface, the next step is to introduce the force into the nodal equivalent loads vector. Chandrasekaran (1987) proposes that if Coulomb’s model is used, then when the local and global coordinate systems are not aligned, a coordinate transformation must be performed on the force vector. Haber (1996) asserts that some difficulties may arise if the Total Lagrangian 29 (TL) method is used because the contact area is not known until punch contact is established. One key advantage to using the adhesive model in conjunction with the TL. method is that force transformations are not required because the fiiction force is not a function of the normal load; it is a fiinction of whether or not the punch is in contact. The fiiction forces will be treated as distributed nodal loads (Wertheimer, 1991) and introduced to the nodal equivalent loads vector at the element level. As Haruff (1995) explains, much care must be taken when determining the appropriate fiiction coeficient. Without great research into the local contact conditions of the metal forming problems used for benchmarks in this study, a realistic variable fiiction model coefficient is nearly impossible to attain. Therefore, a typical constant fiiction coefficient for the adhesive and cohesive model will be implemented into the proposed model. In the proposed model, the shear strain is partly a fimction of the shear stress due to the fiiction, which is defined via the fiiction model. A C'1 continuity constraint is allowed for the shear stress due to fiiction. This constraint is not ofi‘ensive because the shear stress due to friction, in theory, is not a degree of fieedom; it is a specified correction term for the shear strain field. The fiiction contribution to the shear strain field is not included in the shear strain energy. Rather it is simply utilized for post-processing purposes. However, as stated previously, the friction contribution to the externally applied load is included in the nodal equivalent loads vector. One obvious question that arises, then, is ‘Vvhat fiiction coefficient value is appropriate”? As illustrated in the previous paragraphs, the complexity of the fiiction force dynamics can be overwhelming. The cost of accounting for all of the significant 30 details is simply prohibitive with the current state of technology. It seems that the most efiicient method of determining a friction coefficient value is to assume a relatively simple model, then calibrate the model to match experimental data. Wang and Wenner (1978) took such an approach; they selected a coefficient of fiiction which best reproduced the experimentally observed data. Stoughton (1985) also selected a coefficient of fiiction by considering that coefficient which best represents experimental data. A sophisticated fiiction model is not the focus of this study. The use of fiiction values in the shear stress model is one of the foci of this study, though. Therefore, the fiiction coeficient for both the adhesive and cohesive fiiction models will be selected in a manner similar to that demonstrated in the published literature of Stoughton, Wang and Wenner. Unless otherwise stated, the LNQS results provided in the following sections for metal forming problems with fiiction are based upon an adhesive coeflicient of fiiction of approximately 0.8 and a constant cohesive (Coloumb) coefiicient of fiiction of 0.5. 31 83w $26 fiancee 3:07.356 N E .565 Bataw go; .33»... mm 2:9“. 33m .38.... new a .00... n! c ('Il scans Indlauud 381 32 8am 9.35m 3305;. 3:03:25 N 5 8.9.295... 032.3. tea 0.6505! “.0 oEEuxm ”New 2:9“. 08.5w 330:2...— can N F c 7 [Elia .n +0 ... n o n 0 4n 0 a o 4 D O D 0 4a 0 u o 4 n o n. 0 4n 0 n o 4 n o n o 4 a o o 4 n o n o 4 a o a F 4 a o u o announces. 4 4 4 4 4 4 4 u scans ludiouud 181 33 Stress ET Strain Figure 2.3: Elastic - Linear Hardening Constitutive Model. 34 AsperityPeak _\ Lubrication Film _\ Work Piece Tool Figure 2.4: Localized peaks and valleys at the workpiece/tool interface create load carrying mechanisms. The normal load can be carried by the lubricant, metal or both. Coefiicient of Friction 35 Figure 2.5 L (log) Coefiicient of fiiction is not constant in general for metal forming problems. The experimentally determined Stribeck curve indicates that the fiiction coefficient is a fimction of the parameter, L. where: L = fl pR and r], u , p and R are the dynamic coefficient of friction, relative sliding velocity, mean contact pressure and average asperity height, respectively. Region A and B identify regimes where load is carried predominantly by apserity peaks and lubrication film, respectively. Region C identifies a mixed regime. 36 Friction ‘ Force Fina L a Normal Force Figure 2.6: In order to better represent the fiiction behavior, the Coulomb fiiction model may be modified to feature a limit force. Chapter 3 LNQS Element Kinematics The outstanding features of the LNQS element are revealed not in it’s displacement interpolation. Rather, it is the improved transverse normal and shear strain models. It is the combination of the special kinematics and the application of the element in the commercial code via a user subroutine that makes the element attractive for many practical applications. In Chapter 3, the displacement field is first defined. A standard in-plane shear strain expression is given, followed by two assumed quadratic algebraic expressions for both of the transverse shear strains. Boundary conditions and Reissner’s principle are applied in order to explicitly define the transverse shear strain in terms of known variables. The transverse shear strain terms are carefully examined and liberated fi'om any field inconsistent terms which may cause shear locking under certain geometric and loading conditions. Along with the assumption that the transverse normal stress is constant in the through thickness direction, Reissner’s principle is again implemented to obtain an improved form of the transverse normal strain which accommodates isochoric deformation and prevents Poisson’s locking. The topology will be that of an eight-noded brick element as shown in Figure 3.1. The origin of the x-y-z triad axis is at the centroid of the element. 37 38 Displacement Field 3.1 A trilinear displacement field for the LNQS element is defined as follows: 3.1.1 3.1.2 313 i) + 2 {-11% 391% + (.1. H 1.4 We i) r—r 2X1 2 B 2 H + 2 59(1- 1T 1 2 3% iii 73-) l- 39 3.2 Shear Strain Field For each of the shear strain definitions, a, small shear strain assumption is imposed. Therefore, higher order terms are not included in the shear strain models. Notice that the transverse shear strains are defined to vary quadratically in the through thickness direction. 7:,» : Whyui' + Wi.xvr 3.2.1 in. = 1’, +20, +2219, 3.2.2 W... = ¢. +20. +223. 3.2.3 There are six unknowns in the algebraic expressions of the transverse shear strain. Four boundary conditions are now established; they are the shear stress at the top and bottom of the element. The shear stress terms can be defined in numerous ways. At this point in the derivation, they remain generally defined. The following equations were obtained by considering the value of z and the corresponding shear stress value at the top and bottom of the element. The stress and strain are assumed to be proportionally related by the shear modulus. Accordingly, we have t' H H2 €=¢y—-2—ay+7 y 3.2.4 r’ H H2 €=¢y+iay+7flr 3.2.5 I 2 %=¢,--§ia,+—’:—px 3.2.6 2 2 1’1:¢ +£a +fl—fl, 3.2.7 G‘2’4 Solve for B and or using equations 3.24-3.27. Substitute into assumed shear strain terms. . 422 —z 222 z 222 r... =¢’[I"IF)+T;‘(—GTI'+GH’)H:‘(_GF+GH’) 3'2'8 2 _ 2 2 7,, =¢,[l-fl,j+.;[_:+ 22,)..;[L..2_z_,) 3.2.9 H OH CH OH CH A form of the Hellinger-Reissner variational principle is now used to define the transverse shear strain. In a more traditional approach, the displacements and transverse shear strain are assumed to be related by a certain kinematic relationship. Here, however, the assumed transverse shear strain is defined as (initially) independent of the displacement field. The assumed form of the transverse shear strain is then equated to the traditionally defined kinematic form of the transverse shear strain. This equation is enforced by treating it as a constraint which is to be satisfied in the integral sense through the thickness. '--w|~ (n. -rf.)iz=0 3.2.10 h ”I Use 3.2.9 in 3.2.10. Carry out integration. Solve for 1b. Substitute into 3.2.9. The result is a shear strain model which, when used with the appropriate constitutive model, can satisfy the shear traction exactly at the top and bottom of the element. 7),, = Afw, +B,v, +Carya 3.2.11 7.. = Afw, +B.u. +C.,rm 3.2.12 4 \Il/ \l/ \Ifl/ \|_|3} \l/ \I../ \.|3/ \n) \l/ nl_.3} “J... \J 3 2 3 3 2 3 z 2 3 z_H 6_H 6_H 6_H &_H _H 6_H 6_H Arm—H 6_H 6_H &_H _ _ _ _ . 3H 3H 3_H 3_H 3H 3H 3_H 3H H 3_H 3_H 3H 2 2 , 2‘ I \ 2 2 2 2 3 2 . 2. 2 2 {H11 Not my mm Win at ml; (L) [1 my my; (I) 2 2 2 2 . . . . /|\ /l\ H_2 H_2 H_2 \l) \I) \l , ‘1 , \I/ \J \ / \Jl \I/ll \ / \ / )[ l_B 1_B __B __B 1_B l_B _b v,._b Y._b v,._b y_b v.70 /l\ {..l\ /|\ {|\ {l\ /|.\ _ _ _ _L .x_L, . _L, .x_L, . _L, . _L, 1.2. 1_2 1.2 , 111 . . 111 . 111 + _ + \l l \n) \.|/ \l I \l I \l) 1_2 1_2 1.2 1_2 1_2 1_2 _ L 1_L 1_L _ L __L 1_L {.I\ r\ {I\ /.|\ {I\ r\ /|\ (I.\ /.|\ /.|\ f.|\ {IK : _ __ : = = = : : I ”6 n1 2 36 A A A A A A A A A A A A 42 3__2_2_) 2H H3 2 %)(’1)( 2 55X 2 312 2412-6—59) 2.61:) 2H H3 M ..inZ 2 L 2 b l l 1 y) )(3 62’) ___ +1..__.____ 2 b( 2H H3 2 fl 1.. Mi???) 2’. 2 b 2 fll .1. J 322 2 4HH 43 3.3 Normal Strain Field The in-plane normal strain expressions are defined below. Here, the higher order von Karmon terms are included. 3.0: : Wigui +%(Wl,xui)2 +%(Wi,xvi )2 +%(Wi,xwi)2 3H3] 8)? : Whyvr' +é-(Vi,yvi)2 +£(Wi,yur‘)z +-;'(l//i.yW'-)2 3..32 Using the elastic material matrix, the transverse normal strain is assumed to be defined as s = —l—(a,, —C,,e',,r —— C235”) 3.3.3 33 _ Similar to the application of Reissner’s principle for the shear strain, the transverse normal strain is treated. Wit: (3,, -s:,)dz= 0 3.3.4 i "I where: 8; = mei +%(mei)2 +-;—(y/mu‘)2 +%(W’.xv,)2 Assuming that (321 is constant in z, substituting 3.3.3 into 3.3.4, and finally solving for on yields: 3.3.5 44 022 = C33{(W1,2W2)+‘;'(V2.2W2)2 +%(V 2,211, )2 +%(V1.2V2)2}+ u, +-é[(1),,fi%)2 +(1),..,,221,)2 .2-;—(D,.,u,)(1),.,u,)+-;—(D,,,u,)(1),,u,)]} + _C'} b N I i? r C13{% (022,321,!)2 + (07.1:er +% Dali’s Dr,xvr) + %(Dr,xvr DB’va {I} + r C13{% (DB'ow )2 +(DT'J‘WT)2 +%(D8.xw8 XDT.wa)+%(DT.wa XDBJWB ):l} + C23 {% Di'yv‘ + %[(DB-Yv3 )2 + (071va )2 + %(Da.yv8 XDT.va) + %(Dr.yvr XDB.va ):l} + F C,,{-‘6- (D,_,..,)’ +(D,,2,)’ +.;.(D,,2,)(D,,2,)2%(D,_,2,)(D,,2,)]} .2 C22 {.1 (wa, )2 +(D,,,w,)’ +%(D,.,w,XD,.,w,)+—;-(D,,,w,)(1),.,w,)]} . _ 1 _ 1 1 _ 1] Where' D’ ’ (2 L)(2 ' B 2112291111 2 L 2 B 2312:1112) 2 L 2 B .. 41-11112) 2 L 2 B \II/ \I) \l) \l} \J \I/ \I) \I/ 1 H ._._H 1_._H 1_._H 17 17 17 1_H _ 2 2 2 2 2 2 2 2 ( (1K {L [ {Il\ {\ {II\ ,III\ \I/ EJB £13 flaw \Mw \va BM \va y_B mum + . _ + + _ _ + + 1 {I\ (k [ {I\ {II\ [ {.|\ {.11\ ./||\ x _ _ 1 12 12 12 1_2 1_2 12 12 12 /|\ (K {I\ {L (I1\ [ {L (II\ {I\ __ _ 4 D E E E E E E E E .3.3. 3 Substitute 3.3.5 into 46 6:: : {(Wr :wr) '1' %(Wuwr) + %(Wi.zui) +§(Vi.xvi) }+ {,LE “,2, 2 (14(1) )2 2(13,,2, )2 2%(13,,2, )(13,,2,)2%(13,,2,)(15,,2, )1} 2 "g‘i{% r( “Bng )2 + (151',er )2 + %(Da.xva Xbmvr) + %(Dr.xvr XDB.va ):l} + 33 _ , %{% ( ‘19ow )2 + (137.3%): + £053.ow XDT.wa) + %(Dr.xwr XDABJWB ):l} + 33 L gig—12,2, 2%[(15,,2,)2 2 (13,,2, )2 ., %(13,,2,)(15,,2,)2 5(D,,,2,)(15,,2,)]} 2 %{% i( ,,2,)‘ +(zs,,2,)' 2;.(15,,2,)(15,,2,)2%(15,,,2,)(15,,2,)]} 2 33 .. E2213 *,,2,)’ 2(15,,2,)’ 2513,22, )(13,,w, )2§(15,,2,)(25,,2, )1} 33 - 3.3.6 where D 3: = Du V, x A Du: = D” ‘ Wu and B= 1,4 , T=S,8 It is worth noting that there are multiple ways of carrying out the algebra required in equation 3.3 .4. Consequently, the final form of the transverse normal strain is dependent upon the algebraic manipulations used in 3.3.4 prior to integration. For example, the gradient g— is defined as 9637’ + 26;} , where the subscripts T and B refer to top and bottom, respectively. As a result the square of % would feature only three terms. Namely 47 2 The above expression is much more manageable than that which would arise if (g) were to be expanded in terms of all eight nodal degrees of freedom; sixty-four terms could potentially arise! 48 u,v,w Figure 3.1: The topology of the proposed element makes it attractive to the applications oriented user. Chapter 4 CNQS Element Kinematics The LNQS element features normal strain models which vary linearly in the through thickness direction. Accuracy in simulating the normal strain values can be improved if the normal strain is allowed to vary cubically in the through thickness direction. This is apparent as bending becomes severe. The CNQS element features a normal strain that varies cubically in 2. To derive this model, a slightly difl‘erent approach is taken. Instead of beginning with a standard set of trilinearly interpolated displacement functions, a form which varies cubically in z is assumed a priori. Boundary conditions are set and a special assumption is made to define the in plane displacement in terms of known variables. Similar to the LNQS model, a constant transverse normal stress in the through- thickness direction is assumed. Aside from the higher order normal strains, the CNQS element is identical to the LNQS element. The topology will be that of an eight-noded brick element. The x-y axis is at the centroid of the element, while 2 = 0 at the bottom of the element. 4.1 Displacement Field First, cubic polynomials are assumed for the in plane displacements. u=uo+7t,z+,6,,zz+77,,z3 4.1.1 v = v0 +71'y2'5'21-flyz2 + 77,23 4.1.2 I W=Zyliw 4.1.3 - .— where: 152,132 and 112 are rotation terms and are functions of x. 49 1:28, and n, are rotation terms and are functions of y. u and v are the in-plane displacements u. and v. are the in-plane displacements at some reference line w is the transverse displacement w, are the trilinear shape filnction defined in Chapter 3. From equations 4.1. 1-4. 1.3, the strain field is given below. a! 0% 2 519 3 0"") = —°+ ——"+ —-"-+ —i‘- 4.1.4 8,, a: z a: z a: z 222 at e”=—l+z—1+z’@+z3@—’— 4.1.5 03’ 03' 4’ 03’ 3 £2: = ZVile 4..16 i=1 722 = 7r, +2213, + 32.21],‘ + V2.2”: (sum implied on i) 4_1,7 7,, = 7:, +225, +3227), + WWW, (sum implied on i) 4.1.8 It follows that the shear traction at the top and bottom surfaces are defined as 722(0) = 072(0) = If; 4.1.9 r20!) = 072.0!) = r; 4.1.10 92(0) = 072(0) = ti. 4.1.11 rn(h)=07y,(h)=r’y, 4.1.12 A pivotal assumption will now be made. If equations 4.1.7 and 4.1.8 are inserted into the above equations, the presence of a derivative of w, with respect to either x or y, will be introduced. Consequently, in order to develop a C0 continuous theory, rotation terms 51 would need to be introduced. One aim of this study is to explore the behavior of a CNQS element which features no rotational degrees of freedom. Accordingly, putting 4.1.7 and 4.1.8 into 4.1.9 - 4.1.12 and making the following assumption W222 ->0 4.1.13 w,,,,, —20 4.1.14 yields b .75.:n-x 4.1.15 1,! -(—’§-=7r,+2h,6,+3hzn, 4.1.16 b 1. —(’—;—=7ry 4.1.17 1,! y: _ 2 F-rt,+2hfly+3h 1], 4.1.18 Note further manipulations of 4.1.13 - 4.1.16 will result in the derivation of higher order in-plane strain corrections stemming fi'om tangential fiiction at the top and bottom of the element. Perhaps, some effectiveness in the formulation’s ability to model the influence of fiiction on the in-plane strain terms can be expected. However, the shear strain model is not compromised in any way by 4.1.13 and 4.1.14. Using 4.1.15 and 4.1.17 and solving for 17,. and r), in 4.1.16 and 4.1.18, and then introducing these expressions to 4.1.4 and 4.1.5 gives the following result. Tb 2 23 1" Tb = + .a + 2... _22__£__2h 4.1.19 u! ”a Z[G] 2 fix 3h2[G G fix] Tb 2 23 Tr Tb v =v +z —’“— +z + n - 22 —2h 4.1.20 t o [G] fly 3h2 G G fly Let u, and v., be defined as u, and vb, respectively. Then at z = h, 4.1.19 and 4.1.20 become 2'" h r' 1’" = + 4.2),! +_ 4-4-2}, 4.1.21 u! ”I: {0] fl: 3|:G G fix] Tb h T! Tb : + _’.'.'_+h2 +__:"1.__y_z__2h 4122 2, v” ’10] ’3’ 3[G G ’6’] Solving the above expressions for [3,, and Byand inserting into 4.1.19 and 4.1.20, yields a = A1": + A2142 + Clrf' + C22”; 4.1.23 v = A,v, + szz + Clrf' + C21? 4.1.24 w = 0,191 +D2w2 4.1.25 where: AI = ”(232—122 +(};2?1z3 3 2 A: ‘12—21" (171 a 212-12122 +1211 53 z D = — 2 h It is worth noting that the shape fimctions, C, , are based upon a linear shear strain assumption with no regards to rotation terms. Consequently, C, (and thus 2) may tend to lose some degree of validity under large shear strain deformation. In the in-plane directions, the degrees of fi'eedom ui, vi, 1:“ and r?” are assumed to vary linearly. Accordingly, u = aiu. + 1,1? 4.1.26 v = a,vn + 1,1,?" 4.1.27 1 x l W“ “1 457115—7114 2 D:— 2h 53 It is worth noting that the shape functions, C, , are based upon a linear shear strain assumption with no regards to rotation terms. Consequently, C, (and thus 17) may tend to lose some degree of validity under large shear strain deformation. In the in-plane directions, the degrees of freedom ui, vi, rf’ and If" are assumed to vary linearly. Accordingly, II = at". +11??? _ 7’ v— av +13, 4.1.26 4.1.27 54 l) 2 B + 2 i111 + 2.41 24 1) ZB 2 {-111 .1_ :)(1 £16, 2 L 2 B + 2.41 24 2.-.-( W £11 ac: 2 L 2 B + _1. ’1 2 B 2 £111 1. 4.1.28 V2”: Shear Strain Field 4.2 Using the displacement field, the shear strain field is given as 4.2.1 + Vixwi 1:: am”! + 11.: T! 7:: 4.2.2 y: “a": + 12.17.- + Whywr‘ 7,, 4.2.3 ai.yui + ar’ ’3: vi 71? 55 4.3 Normal Strain Field An identical procedure to that described in Section 3.3 for the assumed transverse normal strain derivation results in a similar expression. Note that, in spite of the assumptions of4.l. 10 and 4.1.11, the rotation terms are present in the shear strain terms of4. 1.4 - 4.1.5. Recall, that the presence of the shear terms in the in-plane displacement expressions was realized in an efi'ort to develop a suitable shear strain model. The target of this aim is not to improve the in-plane displacement model, but to enhance the shear strain model. To this end, the x shape fimctions will be assumed to be active only in the shear strain terms. Accordingly, we may re-wn'te the in-plane displacements as u = 62,11, 4.3.] v = am. 4.3.2 Then we have, a, = aura, +£(az,,,u,)2 + 522,392)2 + %(w,.,w,)’ 4.3.3 2,, = 2,2, 2.2102,», )2 259202,): +§(2,,2,)’ 4.3.4 56 F 2 2 2 1 2 2 2 (Duva) + (Dmvr) + %(DB.,vB)(D,Jvr) + %(Dmv,) 1r 2 1 2 1 1 2 2 1 2 2 3 (2,2,) 42,...) .,(2,,,.,) 2......) 2.2.(D.,..,)(2,,.,) . 2‘2- 5 D» II _D I Q i Chapter 5 Finite Element Models 5.1 Comparison of Features As pointed out in Chapter 4, the LNQS and CNQS elements are similar in form, yet different in function. In particular, the normal strain terms are of a difi‘erent order in 2. To highlight this relationship, Table 5.1 is given for study. A column is provided to describe the order of variation in particular directions for displacements and strains. From the Table 5.1, it can be seen that the displacement and strain fields for both elements are described in terms of the same degrees of freedom; only the interpolation of the degrees of freedom differs between the two elements in some cases. Therefore, for brevity, only the CNQS slmpe functions will be shown in the following section. 5.2 E, and EN, Matrices The B matrix, which when multiplied by the displacement vector, defines the element kinematics. The linear and nonlinear B matrices are provided below: P a,“ O 0 O a 0 B, = 2 C: 2 G: 5.2.4 Wit}, Wt! 0 0 an: Why _ at}: O Wi,x_j 5.2.5 57 58 a,“ O O ,y o o. PunIr 0 0 v” 0 w,x 0 0 - 02,2 0 0 o 24,, o o v,y o o w,, o 0 2.2 0 E _ o o u,, o o v,z o o w,, 0 a2, 0 "L " o o o o o o o o o 0 a... 0 o o o o o o o o o o o 1921‘)” o o o o o o o o o 2 C33 .. b d 0 19.313 i,x 2 C,, _ o E, 2 The above B matrices are further expanded in Appendix A. Note that the interpolation for the shear traction terms is not included in the B matrices; the contribution of the shear traction enters into the formulation through the force vector. As is shown in Appendix B, upon defining the shear traction, it is multiplied by the appropriate stifi'ness terms and brought into the force vector. 5.3 Direct Stiffness Matrix In Chapter 2, a general form of a linearized set of equations that can be manipulated in order to find a solution to a nonlinear problem via the Newton-Raphson method, was provided. This form is now made specific to finite element modeling. The total potential energy can be expressed as follows U = I{a,s,}dV—It,u,d4—Ib,u,dV 5.3.] V A V where: o and e are the stress and strain measures t is the externally applied traction b is the body force (assumed zero for this study) 59 The weak statement of the finite element model can be obtained by equating to zero the directiorml derivative of the potential energy in the direction of 5:4. a1=I{a,6£,}dV—It,&r,dA :0 5.3.2 V A Alternatively 5.3.2 may be expressed as, a]={ja,N,_,dV—jz,N,dA}ar=Rar=o 5.3.3 V A where: N is the vector of interpolation functions R is the residual Equation 5.3.3 is the finite element discretization of the pointwise equilibrium equation. Similar manipulation yields a more useful form of 5.3.2. u={j[[§, +%§N,)Z]T5[[§L +%§,,)Z]W+ [Zia/1} 5.3.4 V 6U={I[(§, +%§m)afi]rfi[(§, +-;—§N,)Z}/V+£763dA}=o 5.3.5 V The direct stifiiess matrix is defined below. 3:13 22 5.3.6 Note that the direct stifiress matrix is not necessarily symmetric. 5.4 Tangent Stiffness Matrix Upon taking the directional derivative of equation 5.3.5 with respect to an arbitrary A yields 151/ =[i&]a+az[ia] 5.4.1 dA dA dA =d‘i[(§ +BNL)&A]a+[(B +§m)&]D: —[(B,2 +—B NLJA] 5.4.2 = {13,342, +522)’5(f3'2 42m} 3.4.3 The tangent stifiress matrix is defined below. 1?, = 5;,57+§{5§L +3132)?” +§§L5§L +§[5§NL 5.4.4 5.5 Finite Element Model Using the notation used in previous sections, a standard description of the finite element model is given below. Efli =(7—ED)&'§ 5.5.1 5.6 Shear Traction Definitions As mentioned in Chapter 2, neither the adhesive or cohesive model is particularly accurate in modeling detailed aspects of fiiction. However, from a macro perspective, both models are adequate. One special feature of the proposed models is the ability to satisfy the tangential traction, exactly, at the top and bottom of the element for non-zero fiiction cases. This can be accomplished, in part, by accurately defining the shear traction. 61 The standard Coulombs fiiction model can be used to define the tangential forces in terms of the normal forces at each node. If Coulomb’s cohesive fiiction model is used, then knowledge of the global normal reaction load, and direction of the fiiction force need to be obtained. However, if an adhesive model is used, only the direction of the fiiction needs to be obtained. The shear yield strength and coeficient of fiiction can be defined to be constant. The contact forces (and direction of forces) can be determined fiom within the subroutine. An example of how the shear strain is expanded for the case where contact fiiction is present on the bottom surface of the element is given below. _. a: 7x2 - Wi,xwr‘ +ar’,zur + 11.31301 5.6.1 h . n _l 1:: + 3+ 2: + a: Not only do the shear strain terms need to be adjusted for fiiction, but the nodal equivalent loads vector must also be modified to reflect the nodal tangential nodal forces due to fiiction. From within the subroutine, the global normal reaction loads are used to approximate an effective normal force acting upon the element face. This efi‘ective normal force is, first, divided by the element face area, then second, multiplied by a coeficient of fiiction to yield an expression for the tangential force due to fi'iction. This tangential force is then appropriately distributed into the nodal equivalent force vector. (For some unknown reason, the MARC main program will not recognize non-zero user-defined nodal equivalent load vectors. Therefore, the program was “tricked”; the nodal equivalent loads where subtracted from the internal force vector which is recognized by the main program. Because the main program eventually adds the nodal equivalent load vector and the internal force vector, no modeling errors are realized.) 62 5.7 Shear Locking One of the most common complications associated with finite element research is that of locking. A locking mechanism may be defined as an increase in element stiffness due to strain field incompatibilities. Shear locking will now be considered. A two dimensional case is studied. All conclusions may be extended to three dimensions. For thin cross sections under bending, the shear strain tends towards zero. From Table 5.1, the following constraints imposed by a zero transverse shear strain condition are noted: From constant terms: a) 1+1: 0 dr 6% From x terms: b) 1:— = O Fromzzterms: c) év—+-a—‘=0 6% & Constraint a and c are completely compatible with the zero transverse shear strain condition. However, constraint b is not compatible with a zero transverse shear strain condition. Notice that for in-plane loads, constraint b is compatible because the in-plane displacement will not tend to vary in the through thickness direction. Bending, on the other hand, requires that the in-plane displacement changes with respect to the transverse coordinate. Consequently, constraint b will oppose bending action which results in looking behavior. The current model features field consistent transverse shear strain by eliminating constraint b. This is accomplished by removing all terms in the transverse shear strain associated with the x coordinates. 63 5.8 Volume Locking One of the challenges of developing three-dimensional elements for sheet metal forming is to model the isochoric nature of the deformation in the plastic mode. The restriction placed upon the element is that the sum of the normal strains must add to zero at each integration point used. Failure to do so likely will result in what is called numerical ‘Volume locking”. One method of avoiding volume locking is to underintegrate the element. In the case of the proposed element, this would require that the strain be evaluated only at the mid plane of the element. Such an integration plan would be very detrimental to the element because the bending efl‘ects would not be modeled. Fortunately, becausethetransversenorrnalstrainhasthesamevariationinzastheinplanestrain terms, a firll integration scheme through thickness can be implemented without the consequence of volume locking. This important feature is clearly seen when Table 5.1 is considered. From Table 5.1, the through thickness variation of the transverse normal strain is equal to that of the in-plane normal strain terms. 5.9 Poisson’s Locking To explain Poisson’s locking, the two dimensional case will be considered. Recall the form of the transverse strain. 2,, =g£_§ngn 5.9.1 C22 C22 If Poisson’s ratio is assumed to 0.3, then the following material constants arise for plane stress and plane strain. C,l = (1.13)E Pl Str : . . ane ess C12 = (0.38) E 5 9 2 , C = (1.50)E Pl Strain: " 5.9.3 m C,2 = (0.50)E For a very thin cross section, the transverse stress will approach zero in the absence of double-sided transverse loads. A zero transverse stress condition for a plane stress (or plane strain) cross section is tantamount to the following. = --1-8 5.9.4 3 The resulting constitutive relation for such a case is on = Es,“ for plane stress 5.9.5a on = 1.33E83 for plane strain 5.9.5b The finite element must satisfy 5.9.4. If the transverse strain field and in-plane strain field are constant and cubic (or linear) in 2, respectively, then 5.9.4 can only be satisfied when the transverse strain field becomes zero. The resulting constitutive relation is Plane Stress: on = (1.13)E£n 5.9.6 Plane Strain: (7,,r = (150115322. 5.9.7 Simple inspection 5.9.5 and 5.9.6-7 reveal approximately an 11.5% increase in stimress for both plane stress and plane strain. This efi‘ect is called Poisson’s stiffening efl‘ect (Prathap, 1994). The proposed formulation accommodates 5.9. 5 with no difficulty because the transverse and in-plane strain are both cubic (or linear) in the z direction. Therefore, Poisson’s stifl'ening effect is eliminated. 65 Table 5.1 Comparison of Features LNQS Field Order _ r r r u— Wt" x ,y ’z _ r r r v" Vivi x 3y 32 w - u! w x' y' z' - r r 3 3 as: = Wigui xos))lsz1 r o r 6», = Whyvi x 2y ,2 1 C .. l C . an -_— E,w 4.2—C13 Duu, +2_C23 Dwv, x',y',2:l 33 33 _ r 1 r 73y - Way”: + 1‘“,er x 3y ’7' _ I r 2 7y: _ Alflwr' +Blvt +Carya x sy ,2 7n = Ainwt +Biut +Catxa x‘,yl,22 CNQS u: aiuu x's.y1223 v: 0,11. xr’yr,zs w - u! w x' y' z' - t l s 2 ... 0 l 3 an: - auxui x 2” ,Z _ r o 3 l 13 2 1 C23 2 r r 3 5“ ‘2 EiW “FE-C—v—Dw‘u‘ +EE-Di'yvi x ,y ,2 33 33 _. t r r 729! ‘— Whyui + ngvi x 2y ,2 ._ l r 2 7y; " Waywi +auvr +£1.27.” x 3y 33 __ x: l l 2 7n “ Vaxwr +ar,zu1 +1127: x :y :2 Table 5. 1: A summary of the displacement and linear strain expressions for each of the two proposed elements is provided . Shape fimctions are defined in previous sections. fire superscripts in the comment column refer to the order of variation per coordinate. Chapter 6 Results and Discussion In order to verify the accuracy of the proposed elements, an array of academic- type studies are undertaken. Three major features of the elements are examined; they are the kinematic, material (plasticity) and fiiction models. When the kinematic model is verified, the linear material model is used instead of the plasticity material model. The use of the linear material model in the kinematics study allows the kinematic results to be studied without any additional kinematic influences introduced by the nonlinear material model. After the kinematic model is verified and validated, the nonlinear material model is then introduced and studied. This sequence of feature introductions assures that each feature can be studied most effectively. Following the kinematic and material model verifications, the fiiction model is studied. Like the material model, the LNQS fiiction model is based upon standard theory. (Recall that the CNQS element will not be evaluated for reasons described in Chapter 6.) It is appropriate to recall the main distinguishing features of the LNQS model. They are as follows: 1. Higher order approximation of transverse normal strain in the through-thickness direction. The benefit of this feature is the elimination of Poisson’s and volume locking. Additionally, threebdimensional efl‘ects are more closely modeled when a higher order transverse normal strain is included. 67 2. Ability to satisfy transverse shear traction exactly at the top and bottom of the element. This is perhaps most beneficial in cases of double (or single) sided contact for moderately thick sections. 3. A quadratic variation of the transverse shear strain in the through thickness direction. In cases where three dimensional effects become significant, as in very low R/t ratios, accurate modeling of the shear strain can become important. 4. A development of an efi'rcient programming environment (commercial program preprocessor, solver and post processor via FORTRAN subroutine). This network serves as a useful and convenient testbed tool for research and development. These four features will be evaluated in the context of several metal forming benchmark problems. The initial intention of this study was to use three dimensional elements for all of the benchmarks. Unfortunately, there exists some incompatibility with a user defined three-dimensional element and the contact algorithm. Afier many unsuccessful attempts to rectify the problem, it was decided that the two dimensional case would have to be implemented. Therefore, the metal forming verification problems are in two dimensions only. Unless otherwise stated, the units for length, mass and time are as follows: time: seconds length: millimeters (mm) mass: kilograms (kg) Therefore, force and stress will be expressed as follows: force: N 68 where N represents one Newton stress: _1Y_ = —N—2 = Pa x 106 mm2 (m x104) where Pa represents one Pascal 6.1 Cantilevered Section Under In-Plane Tip Load Figure 6.1 shows the boundary conditions of an axially loaded member subject to an axial load. The nodal displacements at the wall are set to zero, while a positive distributed force at the section tip is applied axially. The parameters of the problem are as follows: Length = l (Area)(Y. Modulus) = 1 Load = 1 Figure 6.2 shows the results for an axially loaded specimen. In this case, the nonlinear material model is not included. The exact linear and exact nonlinear solution are shown along with that of the LNQS element. Noting that the only difference between the LNQS and CNQS element is the through thickness interpolation of the in-plane strain terms, it is clear that LNQS and CNQS solution would be identical for this particular load condition. Accordingly, the CNQS solution is not shown. Excellent agreement is seen among the exact nonlinear solution and that of the LNQS solution. This trend demonstrates the ability of the element to accurately model in-plane response to in-plane loads. For most sheet metal fornring conditions, the maximum allowable in-plane strain is approximately 20% - 69 40%. Figure 6.2 shows that the in-plane kinematics easily accommodate such levels of strain without an introduction of significant error. 6.2 Simply Supported Section Under Uniform Transverse Load The parameters of the problems are as follows: Length = 20 Heights = 1 Y. Modulus = le+06 Load = 100/unit length Figure 6.3 illustrates the boundary conditions. Figure 6.4 quantifies the efi‘ect of shear and Poisson’s locking in the absence of the strain field correcting mechanism (higher order transverse nornral strain assumption). Displacements were taken from the center of the member for both the LNQS and analytical models. The results have been normalized to the analytical solution. The first column shows that without a consistent transverse shear strain field, shear locking can be very severe. The second column corresponds to a modeling condition with a consistent transverse shear strain field, but without a higher order transverse normal strain field. This solution is a result of Poisson’s locking. The locking is much less severe, although it still remains significant. The third column shows that with both the consistent transverse shear strain and transverse normal strain fields, the LNQS element is able to accurately model bending behavior; all locking mechanisms have been removed. In structural analysis, say for an aircraft wing section, the strain distribution is often driven by force distributions. However, in most sheet metal forming cases, the 70 strains are driven by displacement of the tools. The displacements of the nodes are almost entirely dependent upon the geometry of the tools. However, the strain distribution is heavily dependent upon both the tool geometry and strain models. Unrealistic strain models may result in unrealistic reaction loads because of locking effects. In the end, shear and Poisson’s locking can contribute to error-prone fiiction models (because of error- prone reaction loads). Consequently, the consistent strain fields of the LNQS element help promote more accurate metal forming results. 6.3 Transverse Shear Strain in Single Element Figure 6.5 illustrates the boundary conditions on a single element along with resulting transverse shear strain plotted out with respect to the z coordinate. Both solutions of the LNQS and standard bilinear continuum element of MARC Program are shown. The MARC.3 element is a four-noded quadrilateral element with no reduced integration. The difference in solutions is quite pronounced; where the MARC element fails to satisfy the transverse shear strain exactly at the top and bottom of the element, the LNQS element succeeds in doing so. Admittedly, the transverse shear strain (stress) is usually several order of magnitude less than that of the in-plane normal strain (stress) in many sheet metal forming problems. (Exceptions may be localized regions near draw beads or punch openings where the tool radius to sheet thickness is small.) However, the relatively small shear strain that is developed in may sheet metal parts during forming, does not suggest that the quadratic transverse shear model of the LNQS element is without noteworthy merit. First, any accuracy improvement that can be obtained without 71 compromising greatly the efficiency of the element should be promoted; the presence of an accurate shear strain model represents one less variable that is potentially a source of error. Second, as with most research elements, the LNQS element simultaneously rests upon the foundation of many elements of the past and upholds a portion of fiamework which many future elements may stand upon. So, although the quadratic shear strain model advantages may not appear to be very pronounced in the context of sheet metal forming, the LNQS quadratic shear strain model may help serve as a basis for research in metal forming simulation where the shear energy may dominate the mechanics. 6.4 Cantilevered Section Under Transverse Tip Load The parameters of the problem are as follows: Length = 20 Height = 1 Y. Modulus = 1e+06 Figure 6.6 shows the boundary conditions of a cantilevered section under transverse tip load. Before the kinematics can be studied for this case, a mesh refinement study was undertaken. Figure 6.7 shows the results corresponding to element discretizations which vary in the lengthwise direction. Complete convergence appears to be obtained at 50 elements in the lengthwise direction. However, convergence within reasonable limits is realized with about 20 elements in the lengthwise direction (about 5% difi‘erence from the 50 element solution). Figure 6.8 shows results for a mesh having 50 elements in the lengthwise direction and varying number of elements in the transverse direction. Complete 72 convergence is obtained with only two elements in the transverse direction. It is noted that the difference between the one and two element solutions is approximately 3%. Figure 6.9 provides solutions from the commercial finite element program, MARC, and those from the LNQS and CNQS elements. The purpose of Figure 6.9 is to identify the difl‘erence in solutions for a given mesh of 20 elements in X and 1 element in the Z direction. As stated previously, the MARC.3 elements are four-noded quadrilateral elements with no reduced integration. The MARC. 1 14 elements are four-noded elements with reduced integration. The MARC. 114 element, therefore, has an improved bending response (compared to the MARC.3 element). The most compliant elements were the MARC. 1 14 and LNQS elements. Very good agreement is seen among these two elements. As expected, the more stifl‘ elements were the MARC.3 and CNQS elements. It is worth noting that for smaller deflections, the MARC.3 element is in better agreement with the MARC. 1 l4 element than is the CNQS element. However, for larger deflections, the CNQS element is in better agreement with the MARC.114 element than is the MARC.3 element. The MARC. 1 l4 element is well verified and generally accepted fiom an accuracy standpoint when kinematics are considered. The results show that the LNQS element demonstrates a very favorable kinematic response not only for in-plane loads (as was shown in the previous section) but for out-of-plane loads as well. Although the CNQS element is quite accurate in modeling in-plane loads, it’s solution loses accuracy when out- of-plane loads are applied. The CNQS element is considerably more stifi‘ in bending than the LNQS element. The reason for this, of course, can be traced to the assumptions that were used in the derivation of the CNQS formulation. Recall the % terms of the shear 73 strain approximations were neglected during intermediate derivations of the formulation. This was deliberately done so that a cubic variation in z of the in-plane strain could be introduced without introducing rotation degrees of freedom. Because the CNQS element exhibits such dramatic locking behavior, it will no longer be considered through-out the remainder of this study. An extreme case, where the tip deflection is approximately equal to the length of the beam, is also studied for the LNQS and MARC.114 elements. From Figure 6.10, the kinematic accuracy of the LNQS element is further verified; even for extremely large deflections, the LNQS and MARC. 1 14 elements remain in very good agreement. 6.5 Material Model Study The LNQS material behavior is modeled via an elasto-plastic material model which features linear hardening. As described in the material model section of the thesis, the material constants are, in part, firnctions of stress. For fully three dimensional cases, all components of stress need to be included in the material model. However, for the typical plane strain sheet metal forming problem, the stress field is typically dominated by the in- plane normal stresses. In fact, according to Karafillis (1996), one may safely neglect the transverse normal stress when considering a material’s constitutive response for sheet metal forming problems, in general. It is not unrealistic to have the transverse normal strain about two orders of magnitude less than that of the in-plane normal strain. Figure 6.11 illustrates the evolution of the material constants as the in-plane and transverse normal stresses are increased (the transverse normal stress is two orders of magnitude less 74 than that of the in-plane normal stress). The expressions for the material constants are provided in the previous section and are plotted out with respect to the in-plane normal sum; the applied stress is normalized to the yield stress (assumed to be two orders of magnitude less than that of Young’s Modulus). The transverse normal stress is assumed to be two orders of magnitude less than the in-plane stress. From Figure 6.11, we see that C22 and C,2 are essentially constant for all applicable values of stress. However, Cll varies considerably. Five difi‘erent cases for Cu are considered. They correspond to tangent modulii of E, 0.2513, 0.55 and 0.75E and 1.013 (where E is Young’s modulus). The cantilevered section is next considered with nonlinear kinematics and a nonlinear material (plasticity) model. When studying the plasticity model for this case, carefial consideration of the number of integration points in the thickness direction is parasnount. The MARC and LNQS elements are compared for difi‘erent integration point and mesh conditions. Both plasticity models feature elasto-plastic behavior with linear hardening as described in Figure 2.3. The key material parameters are given below: Young’s Modulus, E = 1e+06 Tangent Modulus, E. = 0.5e+06 Yield Strength, Y = 1e+04 The MARC program does not allow the user to define the number of integration points for the elements in it’s library. Consequently, the only way to increase the number of integration points in the transverse direction is to increase the number of elements in the transverse direction. The LNQS user-defined element, of course allows the luxury of Specifying the integration point strategy. Because the MARC element discretization is not 75 the same as that of the LNQS, the comparison is not completely fair. For a given number of integration points, the MARC results are expected to be slightly more compliant; based upon the results shown in Figure 6.8, the difference should be approximately 4%. Indeed, Figure 6.13 shows that the MARC model, for 4 integration points in the transverse direction, is about 3.5%. A few interesting observations are worth noting. First, if the MARC curve is extrapolated back to two integration points, the LNQS and MARC solutions will match to within 1%. This excellent agreement is expected not only because the integration points are the same, but the number in elements in Z are the same as well. The second interesting observation is that the solution for the LNQS element becomes slightly more stifi‘ as the number of lobatto integration points increases fi'om 3 to 5. This behavior is rather uncommon, but nevertheless possible. Figure 6.14 is a description of the Gauss and Lobatto point distribution with respect to natural coordinate in Z. On the ordinate axis are labels identifying the location of the Gauss and Lobatto points. To the left and right of the ordinate axis are the Lobatto and Gauss weights, respectively. The thick solid line represents a stress distribution. Following the dotted lines extending from the elastic-to- plastic transition point of the stress distribution, one can observe the following: 3 point Lobatto Rule: 33% (1 of 3) of the points are in the elastic range 5 point Lobatto Rule: 60% (3 of 5) of the points are in the elastic range This qualitatively explains why it is possible for the LNQS solution (for the Lobatto rule) to become slightly more stifi‘ as the number of integration points is increased from 3 to 5; a 76 higher percentage of integration points remain elastic (in stifi‘er region). It is noted, though, that in the limit, the stimress will decrease as the number of Lobatto integration points is increased for plasticity analysis. Bathe (1982) provides a more rigorous explanation with a simple example. Figure 6.15 shows a two-noded bar element with a varying cross section. A state of stress is proposed such that the domain is in the plastic regime for the left side and elastic regime for the right side. If a one point integration rule is used, then the entire element is assumed to be in the plastic regime. This results in a very small stiffiress. If a two point rule is employed, then 50% of the element is assumed to be plastic regime. This results in a much more stifl‘ solution. If a three points rule is used, then 66% of the element is in the plastic regime. This results in a solution less stifi‘ than that of the two point rule. However, if a four point rule is used, then 50 % of the points are again in the plastic regime (location of the points are difi‘erent hour that ofthe two point rule). This results in a stifl‘er solution than that of the three point. The exact solution is provided at the bottom of Figure 6.15. Figure 6.14 and Figure 6.15 help verify that the results shown in Figure 6.13 are reasonable. Recall Figure 2.3 which describes the constitutive assumption for a uniaxially loaded member. Assume the parameters to be as follows: Length = 20 Area = 0.01 Y. Modulus = 100 Yield Stress = l 77 Tangent Modulus = 50 Load = 0.5 Then, the tip deflection based upon the analytical solution for the bilinear material model is 20. The LNQS finite element solution for the same case is 20.5. These results which are in good agreement, help to further demonstrate that the plasticity model used in conjunction with the LNQS element formulation is correct and accurate. 6.6 Preliminary Commercial Element Investigation The commercial finite element code, MARC, ofi‘ers a variety of elements from which to choose. The plane strain MARC 11 (MARC Eleven) element, in particular, may be used with a standard strain or an assumed strain fonnulation. The essential difi‘erence between the two is the variation of the transverse normal strain in the through thickness direction. In order to investigate the behavior of the elements in the context of a metal forming application, one element was isolated and studied. This element corresponds to the region of peak bending in Stougton’s square cylinder stretch problem which will be investigated in more detail in later sections. Figure 6.16a and Figure 6.16b show the transverse coordinate versus the in-plane normal strain solution for the following element formulations and meshes. . M. rest : MARC 11, refined mesh (5 elements in 2), standard strain This element features a strain field that is defined fi'om the bilinear displacement field. M.re.as: MARC 11, refined mesh (5 elements in z), assumed strain 78 This element features a transverse normal strain field that is not defined directly fiom the bilinear displacement field. Rather, it takes on an assumed form such that Poisson’s locking and volume locking are not possible. M.co.st : MARC 11, coarse mesh (1 element in 2), standard strain M.co.as: MARC 11, coarse mesh (1 element in z), assumed strain The results reveal that the elastic solution is, for most practical purposes, linear in z. The sohrtions using the assumed formulation tend to predict a slightly higher strain than that of the standard (kinematically correct) formulation. It is safe to conclude that the refined kinematically correct sohrtion most closely approximates the correct solution. Without specific information about the MARC element formulation, it is diflicult to explain the difference. However, comments on the LNQS results for this same test case may serve as useful insight into the behavior of the MARC assumed elements; these comments will be offered in the following sections. The plasticity results are very similar to the elastic results in that the assumed formulation solutions tend to show higher strains than those of the standard formulation. 6.7 Lee’s Benchmark: Plane Strain Stretch with Zero Friction In 1990, Lee introduced the benchmark geometry described in Figure 6.21a and 6.21b. The MARC.11 with assumed strain and LNQS solutions are compared. The same mesh discretization (one element through the thickness) is used for both the MARC and LNQS models. Both the LNQS and MARC results are based upon an elastic-plastic model 79 with linear hardening. The hardening is characterized by a tangent modulus of 0.5 times Young’s modulus. The material parameters and geometry for both finite element models are as follows: Young’s Modulus = l.0e+06 Yield Stress = l.0e+04 Poisson’s Ratio = 0.3 Sheet thickness, t = 1.0 Length, L = 59 Punch Radius, R, = 50.8 Unless otherwise stated, the material properties assumed for the MARC and LNQS models are the same as those specified above. Anisotropic material condition is assumed for all cases in this study. As mentioned previously, most sheet metal forming simulations are displacement driven. That is, the boundary conditions imposed upon the work piece are defined or constrained by the specified position of the tools. It is understood that some finite element contact algorithms apply a force imposed by a penalty or lagrangian constraint. But, these forces are defined ultimately such that the work piece satisfies some geometric contour. As a prelude to a closer study of the numerical results for this metal forming case, Figure 6.12 is considered. The results confirm that for the metal forming simulations considered in this study, key material properties (such as Young’s modulus and Tangent Modulus) have a minimal impact on the distribution of the strain. Of course, the reaction loads, as mentioned previously may be vastly different. It is important to point out that it is not 80 being suggested that the strain distributions are totally independent of Young’s modulus or the tangent modulus. For example, a material with zero hardening would show dramatically difi‘erent results. Because the material is unable to resist deformation at the “point” where yielding occurs, the other unyielded points will tend to remain unyielded. The result can be a highly nonhomogeneous strain distribution. Figure 6.22a shows a plot of the in-plane strain solution for the MARC 11 and LNQS model. Without fiiction, the problem is essentially reduced to geometric issues only. The membrane strain is constant in the x direction, as expected. (The top and bottom strains may be averaged to obtain the membrane strain.) At the end of the punch stroke, the sheet metal conforms closely to the contour of the punch fi'om x = 0 to x = 17. Since thepunchradiusisconstant,thepeakstrainisalsoconstantfiomx=0tox= 17. The bending action may be thought of as a mechanism that perturbs the in-plane membrane strain. On one hand, when the bending causes tension (top fiber), the in-plane strain will increase. On the other hand, when the bending causes compression (bottom fiber), the in- plane strain will decrease. Both elements demonstrate an ability to capture bending effects for this metal fornring case. Note that a very small difl‘erence in the top and bottom strain is noticed for most of the punch contact region among the MARC and LNQS solutions. This small difi'erence is attributed to the difi‘erence in contact conditions. When the sheet metal comes out of contact with the punch, the bending action is no longer imposed upon the sheet metal. Under such a boundary condition, the sheet metal is simply stretched between two points. Thus, the deformation mode which assumes the least amount of internal energy is purely membrane in nature. For this reason, the top and 81 bottom in-plane strain solutions tend to converge. The solutions depart somewhat as the metal comes out of contact with the punch (17 < x < 25). This departure of solutions suggests a difference in stifiress between the LNQS and MARC 11 models. Limited information is provided in the MARC manuals with regards to the exact element formulation for the assumed MARC 11 element. For this reason, it is difi'lcult to ascertain exactly why the models perform differently in this region. However, one obvious conclusion can be rrrade with regards to the LNQS solution in this region; the LNQS solution is incorrect. The top strain cannot exceed the bottom strain for this particular application. The author has extensively investigated this problem by varying the number of load steps, convergence criterion, mesh discretization, material ' property, element kinematics. No legitimate explanation was found. This problem remains to be an unexplained anomaly. This behavior is unique to this particular application. As will be shown in the Stoughton’s Benchmark problem, the anomaly disappears. Wlth the exception of the problem described above, the solutions are in good agreement. 6.8 Lee’s Benchmark: Plane Strain Stretch with N on-Zero Friction In the finite element simulation of sheet metal forming using Coulomb’s constant friction model, the normal reaction force may be recovered directly or indirectly. The direct method is to calculate the normal reaction loads at the beginning of each load step. An indirect method is to calculate the nOrmal reaction force fi'om the transverse normal stress at the boundary of the element. Figure 6.17 shows the MARC results using both the force based and stress based models along with the force based LNQS model. Contrary to 82 expectation, the LNQS forced based and MARC stress based model results are in very good agreement, while the MARC force based model results depart significantly from the other two results. The effect of fiiction seems to be much less apparent for the MARC force based model. The solutions can be justified via Figure 6.18 which reveals the actual nodal fi'iction forces that are applied to the sheet at the conclusion of the last step. Details of the exact method used in the commercial MARC program to calculate (which may include smoothing and redistributing) are not available in published literature. Without a complete understanding of the smoothing techniques employed by MARC in it’s commercial code, a fill] explanation of the difi'erences cannot be obtained. On one hand, the MARC stress based forces seem to be over estimated. The approximation shown below suggests that the nodal forces based upon the transverse normal stress should be about one order of magnitude less than that shown in Figure 6.18. Young’s modulus = 1e+06 Yield Strength z 1e+04 z 0'32. as: 0'“ z — 100 Let Friction Coeficient, v = 0.5 Friction Force or A, van , which is on the order of 10. On the other hand, the MARC force based forces appear to be somewhat low when compared to those of the LNQS. Literature published by MARC explains the normal reaction loads are to be used in the calculation of fiiction forces. Nonetheless, certain latitude within the fi'amework of the “standard” Coulomb’s model was likely taken (and maintained as a trade secret) by the developers. Because of the close correlation among 83 the MARC stress based model and the LN QS force based model, the MARC stress based model will be used as the benchmark cohesive model for the remainder of this study. W'rththeformofthecohesivemodels established,itisnecessaryto makea comparison with the adhesive model. As discussed in the previous sections, either a cohesive or adhesive model can be used for metal forming with fi'iction. Figure 6.23 shows the results for both models for Lee’s benchmark case. Very poor agreement is found among the two adhesive solutions. A study of the output file fiom the adhesive model runs revealed that the nodal fiiction forces applied in the MARC adhesive model were not constant. In light of the traditional definitions described in the previous sections, this finding is puzzling. Because of the correlation and lack of insight into the MARC adhesive fiiction model, the cohesive model is the preferred fiiction model to be employed in this study. Lee’s benchmark case is considered, now, with friction effects. The material parameters used by Lee are as follows: Hardening Coeficient, K = 589 MPa Elastic Strain, e. = 0.0001 N-value, n = 0.216 where: a = 1((2'o + 5p)" and e, is the plastic strain Coeficients of Friction, u = 0.3 for Lee’s solution rt = 0. 5 for LNQS cohesive solution u = 0.5 for MARC 11 cohesive solution 84 From Figure 6.24, it is shown that the MARC and LNQS solutions are in excellent agreement. As the punch displacement increases, so too does the normal reaction load which acts upon the sheet metal surface. Coulomb’s model then dictates an increase in fiiction force which acts tangent to the metal. In the absence of fiiction the in-plane strain will increase due to the stretching action imposed upon the sheet metal by the punch. With fiiction, the stretching is encumbered slightly as the fiictional forces counter act the stretching. All three models are successful in capturing this feature. As the metal comes out of contact with the punch, the fiiction forces no longer are able to restrain the stretching motion. For this reason, the maximum in-plane (membrane) strain is reached. Because Wang’s solution is based upon an implicit formulation with an updated Lagrangian solution scheme and Coulomb’s constant fiiction model, Wang’s fiiction model is not equivalent to that of the LNQS or MARC model. The solution difi‘erences are attributed to the difference in fi-iction models. 6.9 Stoughton’s Benchmark: Plane Strain Stretch with Zero Friction Figure 6.25a and 6.25b describe the geometry of the benchmark proposed by Stoughton. The main difl‘erence between Stoughton’s geometry and that of Lee’s is the flat shape of the punch on Stoughton’s model. The MARC and LNQS solution are provided in Figure 6.26. Refer to the previous section for material information. Unlike Lee’s case, though, the punch does not feature a constant radius. The Lee solutions showed that the peak strain was constant along the entire punch/work piece interface. However, in Stoughton’s benchmark solutions, the peak strain is confined only to the 85 punch surface with the highest curvature (or smallest radius). The sheet metal near the line of symmetry is essentially under a membrane load because there is no bending constraint. Likewise, in the region where the metal is not in contact with the metal, a membrane solution arises. Perhaps most noteworthy of all is that there is nearly perfect agreement (unlike the Lee solutions) among the two solutions as the metal comes out of contact with the die. Very good agreement is seen among the LNQS and MARC solutions. 6.10 Stoughton’s Benchmark: Plane Strain Stretch with Non-Zero Friction Figure 6.27 provides results for the MARC, LNQS and Wang model. Wang implemented an implicit, updated lagrangian model. Reasonable agreement is found among the three sohrtions. All three solutions are in nearly perfect agreement where the punch is not in contact with the metal. Yet, the strain is significantly lower for the Wang model than it is for the other two. This is possible because the contact region seems to be smaller for the Wang model. The result for Wang’s model is (with respect to the LNQS and MARC models) a lower strain from along the punch flat, and a higher strain along the circular region. These two effects apparently cancel each other out resulting in the strong agreement with the LNQS and MARC models from x > 11. The solutions are in reasonable agreement. 86 6.11 Thick Section Benchmark: Pinching Boundary Condition Up to this point, only relatively thin sections have been considered. The proposed element is primarily intended to be used for cases where three dimensional efi’ects do not prevail. However, for the sake of research, an extreme case, which stretches the element’s performance well beyond it’s intended limit is considered. Figure 6.19 describes the boundary conditions imposed. The length, thickness and modulus are 10, 1 and 1e+06, respectively. A single point load represents, perhaps, the most extreme case of single sided contact that a finite element domain can experience in sheet metal forming problems. Solutions from a MARC model with five elements in the thickness direction and a LNQS model with one element in the thickness direction are given in Figures 6.20. The region of study is that directly beneath the applied displacement of -0.04. The transverse efi‘ects are nicely captured by the refined MARC model. The transverse normal strain is relatively great near the top surface (~12%). However, it begins to decline prodigiously at the center of the element. The cubically varying solution finally settles to a small value of about - 0.75%. The LNQS solution, on the other hand, poorly represents the actual solution. The smallest strain is at the t0p while the greatest is at the bottom. The linearly varying solution bares witness to the fact that the LNQS model is not appropriate for use in cases where transverse effects, such as those shown in Figure 6.20, prevail. 87 1 b Figure 6.1: Cantilevered Section Under Axial Tip Load VV 88 .88 snow .82.. co. one: Eases oco .. u m< .. u ._ senses. 522m 23 883 2:32 5 co. coeoocoo 33. ...s 83 ”mo 2:9... A .3300 .2: 325.32 ”.22 I I 325 “.22. mOZJII DWI I?!" 89 e e A 65 Figure 6.3: Simply Supported Section Under Uniform Transverse Distributed Load .88 2.3m .592 92a $5820 on .3323 coal—6:23. ecu m.:omm.on_ :35 2:358 3.32 2 mgamscoo :ano Ea... ecu 23.5%.. 63. E6: .558 30:: 5.83 82% 0:2: 3:253 295m :8 5.80:3 3:3:8... ”to 8:9“. 20.". Eofimcoo 95.00.. 9.33.3. 55> 0:203 32a 5:5 - b . . o e .2 m6 .2 ed LI ”.0 .- ad _.._. uosnros mun leaMrIuv wanna 91 0.0% i I 8+0 rum .WIdao. 32m 80:: €2.20 20% :8 :o.Sn:Eu 58$ :35 835.5 50:25 ”me 230.". 38 53am 32.“ 032502». 0—0 3.0 «.20 ...0 00.0 00.0 3.0 No.0 41- .. ...0 .. No .. 0.0 .- 5.0 I I 1— z ‘rnaurprooo ossoasusu 92 L...___ h-J Figure 6.6: Cantilevered Section Under Transverse Tip Load 93 .~ :_ .:0E0_0 0:0 .Scpun. .0+0Pum ..n: .35 .5008 00.05550 3000 05.0 0000.0 :8 >030 0059350 :00: ”5.0 0.50.... . 30:20.2 I I .2.m..... 00 0V 0* x 0. ”co-:05 no 02.532 00 00 0N 0F or I'll f'll 'I'lf 0.0 0.. 0d 0.0 06 0.0 uorrosusa du, oasmuul. 94 x :. 0.5.5.0 00 .82"... 60.7050 .PuI .05.. .5003 0055.500 002.0 05.: 6000.0 .2 030 0052350 :00: H 0.0 059.... N 0. acoEofi 00 03:52 0 .8305 r ... ..2.m.u fl 0.0 0... 0.0 0.0 0.0 “0930“” sarcasm; jlllll! 008 2.020 8.82... 00050.5 5 52.0.0 P .522 2% 0.85% cu 8+2 n m .0 u < .855000 8.05..an 833 22260 .2 8.02.8 a: .9 83 a: 53 23E 5:00:00 023002... on 0.. 0.. 3 up 9. 0 0 v N o ”024' 3 v.01; I I p001 mum; 97 00000 0:0 00.3000... 000.95 2 000000.08 up 000 «N €209.00 00 02:02: 0.0::o> 0. m 80:30 20280000; .000 20 000.0 .0000 .003 s 0:385 E080. 0 2 .32 8 0:0 0m .00 .2. .5005 .3806 0.000.200_0 05 3 00.2858 ”:0 2:9“. Nwo ‘ I «No .I I 00.30% 0N._._.OIII' 00.....0alllu 05.30:! A0030 21> 05 8 3.05.55 one...» .32 60:02 _. 0.0 0.0 5.0 0.0 0.0 0.0 0.0 «.0 _ pl 4 1 l‘ HIIIIIIII-IIIIIIIIIIIIIII-IIIIIII IIIIIJICIIIIIII‘IIIIIII$OO Bum : on :9 100 .1 00—. 1. cup 00; ”I'A am I'm." 98 00000 ...... N300 P 0000 .Io.| 03".; .000~.0u=0.m .00'um ”0 0000 .00Fuo> .N00.ouc0.m 0.0an N 0000 .3+0Fu0> 60.700.05.20 .00+0Pum up 0000 .0000 000.5000. 0.00.. .8 >030 3.0020 .0020: H 3.0 059.... mm mm 0N «N 0m 0F 0? I. N— 0? 0 0 v 00 3 N0 (”DJ-Ins do» was OWN-Ill 96 400241.! 0.002..Ill Quasalol 0_. 90.000030 0.0; .mm.0n.0os_ 00... .000..ua .000nm .0_.u.. 3.0.00.0 5.3 000000 00050 000.0 00.26.0000 00.. >030 000090300 0.00 00.300020. . 9 050.". N 0. $0.00 000208... S .8052 0 0 h 0 0 0 r . .Ln ? . N 0. 0.09020 02“. N 0. 0009020 .50“. N 0. 0009020 020... N 0_ 00000.0 000 ./ N 0. 0009020 03... # SN 4 00.0 : 00.0w 00.9 uogaouoa mum; 100 0 2.000.. n 0 3.000.. x v 00:00 4 0 00:00 a N 00:00 0 2.00.0 0.0 PI 00020.2 ..0. 000390 00.5050... 002.0 0.00006 .N 0. 60000.. 5.3 0.00.03 000 0.0.00 00:208.... 2.000.. .000 00:00 .0 00.50505 ”I. 0.50.... 0 0.0 H...) 00.0 HS> 0. P. 00.005 0.0. 00.0 as _ I I I I I ' ' I I I . . 00.005 00000. a... .. m + 0.0- «0.002, 4 .... N.°I I0 . . . . ._ 00.005 £0.05 00 F .50 .- «.0 00.005 0 .- 0.0 00.05 0.0 00.002. 00.005 0.0 . 00.005 .. 0.20 H.2. 00.0 £5 101 Plastic. E/lOO Elastic. E CALCULATED STIFFNESS FOR BAR ELEMENT l —1 One point integration K = 0.0005El: 1 1 :l l —1 Two point integration K = 0.04164El: l 1 :l l —l Threepoint integration K = 0.027OOE[ l 1 J ' 1 —1 Four point integration K = 0.0402615 [ l 1 :I l —l Exact integration K = 0.03973E[ l l ] Figure 6.15 : For plasticity analysis of the two noded bar element, Bathe (1982) shows that for certain stress distributions, increasing the number of integration points does not always increase the accuracy of the plasticity model. 102 8.8.: ..x L 3.8.21' 8.2.5. .a n 3022* C .5228 032m .38 {2:550 9:02.032» be 581... goon 55 E9293 an 58% 0:29... so; 230E ...-ha 0:10.... i 0’ 0 0 h 0 r? Emu .gozaom , “aches: 5me BE: \ :26 89:8 \ . 82 \\ v.0 N0 0.0 #0 0.0 0.0 50 0.0 0.0 2 ‘ueutpiooo emoticon; 103 3.8.: .x I 360.2% 3.8.: __a I ~méé‘lolu 625.8 028$ .38 fiaezocon $520320 .8 £2.» goon 5.3 .362» .a 58% 23.9... ”no; 959“. 52..» 25...... «a 3 0.. a 0 h o m t N \ \d \ .5028 £90» .0 too \\ \ \ -.,ILEB: PE: \ 59% 38: a Lo 9 T R w . \ \ f \ \ \ \ \ \ \ llxll \ \ \ \ . s i \ \\ rd «.0 0.0 V0 0.0 0.0 . 5.0 0.0 0.0 2 “unugpiooo macaw; 12f 104 VN .mgame {258:3 M34 .326 38.6: 8835: no 33. or 3. NP or .35: 856 :8: 033 on has cozotn. £52.50 £03.: 502.“. 3 53858 N Z .0 2:3... . 83..moz._l gaudy—<2 . . O . . mmobmdmg l I I|.II.IT..I.I 1...... 0,0...- uo ui.oo.c inlll I‘ll- I 0' .C. Nw hummus) mus wild-u: 9s 105 .8283 {2.5.2.00 $3.. 3. 338E. 80.3.35... 3 8.8 5.8.... .802 “w. .0 2%.“. 00 . 092.01; .I I 09.0vaqu uEdemg I I 000 0310:] uonoyg 106 A or? 4'5 Figure 6.19: Section Under Double-sided Transverse Pinching Load 5.02883. .mOz; uca om; ..e N :. mEoEeo F 05 m. .3+opuo>.oo+om.0u:2m .oo+opum .Fu... .mu...uwo.. 055...... 8.05:8... .250 30.5 .3522 .8 58% .352 0232a; .m> 82.280 3.05:8... ”0N6 230.". 5g” _g—COZ agar-Eh. $ 0- 07 NT . _ n _ _ _ 108 Figure 6.21a: Geoméiry for Lee’s Benchmark Problem. 109 Inc : 1000 Tlmo:1.0000+00 lobt Figure 6 2 lb : Final Mesh for Lee’s Benchmark Plane Strain Stretch Problem 110 .3385... I .- _ao..:.om .00300ucafi .~ 0. E0820 P .x 0. 80060.0 cm 6.0an .00»: .0005 0.080.000 £0320 £3550 you. .2 8.05 5.; 05.50". 6.2.0 586 80E “3.0 200.". x 00.0w 00.0w 00.NN 00.0w 00.0.. 00.0.. 00.: 00...... 00.0“ 00.0 00.0 00... 00d 00.0 + i. . . . t. t. T i. i. n . 0.0 .. 0.0 r. 0.0 .. 0.0 .... °.N. 114 8..— ‘ I om(¢‘ 6 I 11 O.m 0024' : 0.0 .. 0.0.. i 0..... 0.N.. mans ouluqlmw Wild-"I % 115 "ii—at r' Figure 6.253: Geometry for Stoughton’s Benchmark Problem. 116 Figure 6.25b : Final Mesh for Stoughton’s Benchmark Plane Strain Stretch Problem 117 om<2 l I wOZJIl 3+3uo> .oo+om.onca.w .~ 5 52:20 . .x 5 8.3803 on €35 dd": .00+o..um .5039... {085:3 9:220:05 .5. 026...“. 505.3 5% 525 000?. ”0N0 230E 0p 0p 3 Ne 0.. 0 0 v "J7 “*1" 1 r L I ‘- O 1 I 0050:0va V 1 I 4 l l 1 U 1 I I ('0 N ‘- O '- 1- v- ‘- V F (”‘31 do” "I338 wild-III 96 118 95$. . a. . . om .00+mm.0uca.m .N 0. .coESu F .x E $000.20 0m .F.0nu._ .00"...— 60.70an .EoEEn. {menu—Em Mcofiusofi 8. 83:... 93.2230 59...; 222% £86 25E find 2:9". 9 0.. 3 n. N. ..F or 0 0 h 0 w v n N w ..4...c....¢...<...fl...<...<...< ugans euuqurew wind-m as Chapter 7 Conclusions As demonstrated in the previous chapters, the first step in validating a new finite element is the academic test, a series of very simple tests such as axially, or transversely loaded members under small loads is in order. The results from such simple tests can be compared with‘readily available analytical solutions. Eventually, large displacements are introduced which may require a previously accepted numerical solution from another model. However, for sheet metal forming, the element must be evaluated under sheet metal forming conditions for legitimate element validation. An established and proven contact model must be used. The solution algorithm has to be robust and reliable. Efi‘ective post-processing is essential. To this end, a very large portion of the research presented in this dissertation has been dedicated to simply establishing the software infrastructure necessary to evaluate the proposed element under realistic sheet metal forming boundary conditions. Reliable, useful and efiicient means of evaluating a new finite element formulation is often times more difiicult to obtain than is the actual new formulation itself. Indeed, such was the case for this study. The research described in this thesis has proven out a method to eficiently evaluate new finite element formulations. The LNQS model, with geometric and material nonlinearities, has been cast into a FORTRAN subroutine which is used interactively in a powerful commercial model. New formulations can be implemented quite eficiently by simply modifying the LNQS source code. 119 120 Rarely will all of the initial goals of an investigation be fulfilled as planned. The research carried out by the author was no exception. A fiill understanding of the interaction among the subroutine and the MARC main program is still yet to be realized. This problem has certainly made it more difiicult for the author to make meaningful contributions to the finite element technology base. Nonetheless, some contribution was made; the goal of establishing an infiastructure, limited as it may be, has been accomplished. Additional research in this area is now more approachable and inviting than ever before. The fi'uits of this research have made it possible for a new investigator to spend much less time developing infiastructure and more time evaluating their elements and applying creativity to make improvements. From this standpoint alone, this research has been a success. As outlined in the introduction, there are numerous finite elements that are suficiently accurate and practical for a limited range of sheet metal forming applications. The Belytschko-Tsai shell element is a good example. The LNQS and CNQS elements were not introduced to compete with such elements as the BelytschkooTsai element which has proven to be successful for very large models such as automotive door fiame stamping. The strength of the proposed elements is revealed in smaller models which can aflord more detailed descriptions of transverse normal strain and shear stress. For small deformations, the proposed elements are shown to be accurate. Shear and Poisson’s locking issues have been resolved. Because of the assumed transverse stress assumption, the transverse normal strain features a through-thickness interpolation which has the same order of variation as the in-plane normal strain. This relationship allows for 121 isochoric deformation. These features make the LNQS and CNQS elements field consistent. Thus, a reduced or selective integration procedure is not required. For large deformations, the LNQS and CNQS elements are shown to be accurate and reasonably accurate, respectively. The deformation imposed upon the element mesh for the large deflection studies were large compared to what the elements would typically experience in sheet metal forming applications. Only under severe deformations did the CNQS and LNQS solutions significantly depart. Both elements exhibit the ability to exactly satisfy the shear transverse shear strain (or stress) at the top and bottom of the element surface. At element surfaces, when fiction is zero or non-zero, the transverse shear strain is zero or non-zero, respectively. This is important when the order of shear strain approaches that of the normal strains. The shear strain model takes advantage of the MARC user subroutine to obtain appropriate values for the shear traction. This feature is unique to the LNQS and CNQS elements. The advantage of this feature is that not only is the nodal equivalent loads vector modified (in response to tangential fiiction forces), but the shear strain model is also modified in order to more accurately simulate the transverse shear strain approximation. Two plane strain metal forming problems were studied. The first was for and fiictionless case. Overall, good agreement was found among the LNQS and MARC models for frictionless cases. In particular, the proposed models were able to capture the membrane and bending effects quite well. Reasonable agreement was found among the LNQS, MARC and experimental results for cases with fiiction. When fi'iction was present, 122 the LNQS model was able to satisfy the shear traction exactly at the top and bottom of the element. Mth their shear strain models and unique application to the MARC program, the LNQS and CNQS elements can play an instrumental part in simulating various metal forming processes which involve boundary conditions with and without fiiction. An infi'astructure has been established in order that scientists may consider the MARC subroutine procedure a viable option for both evaluating academic-type elements and solving practical metal forming problems encountered in industry today. Chapter 8 Future Work Below is a short list of possible follow-up research activities that would be appropriate: 0 Continue development and refinement of MARC subroutine infrastructure. Determine methods of obtaining the global normal reaction loads directly from the main program. Determine how to implement a three dimensional version of the LNQS and CNQS element into the MARC subroutine. c From a material model standpoint, make a detailed study of the advantages of a more accurate strain field description within the element. Compare the LNQS, CNQS and various commercially available elements. 0 Make a detailed study of what the optimum definition of the shear traction at the top and bottom of the element should be (Coulombs model currently being used). 0 Compare the explicit version of the LNQS and CNQS element to the implicit version. 0 Study other types of metal forming such as forging or rolling. 0 Apply the LNQS and CNQS formulation to a time dependent case. 0 Investigate the issue of associated or non-associated flow rules that Stoughton addresses. 0 Apply a polycrytstal-based failure criterion to the LNQS element. 0 Apply a stress-based failure critertion to the LNQS element. 0 Apply an energy-based failure criterion to the LNQS element. 0 Investigate the contributions of bending to material failure in sheet metal forming. 0 Apply proposed post-process stress models to draw bead load boundary conditions. 123 124 Bibliography 125 Bibliography Aitharaju, V.R., “A C° Zig-Zag Finite Element for the Analysis of Laminated Composite Beams”, Accepted for publication in J oumal of Engineering Mechanics, ASCE Averill, R.C., “An Efficient Thick Beam Theory and Finite Element Model with Zig-Zag Sublaminate Approximations”, AIAA Journal, Vol. 34, pp. 1626-1632, 1996 Bathe, K, Finite Element Procedures in Engineering Analysis, Prentice-Hall, Inc, 1982 Barlat, F ., “A Six Component Yield Function for Anisotropic Materials”, Int. J. of Plasticity, pp. 693-712, 1991 Belytschko, T., “Explicit Algoritms for Nonlinear Dynamics of Shells”, ASME, Vol. 48, pp. 209-231, 1981 Belytschko, T., “Assumed Strain Stabilization Procedure for the 9-node Lagrange Shell Element”, International Journal of Numercial Methods in Engineering, Vol. 38, 1989 Bland, D.R., “The Associated Flow Rule of Plasticity”, Journal of the Mechanics and Physics of Solids, Vol. 6, pp. 71-79, 1957 Boubakar, M. “Finite Element Modeling of the Stamping of Anisotropic Sheet Metals”,Engineering Computations, Vol. 13, pp. 143-171, 1996 Buchter, N. “Three-Dimensional Extention of Nonlinear Shell Formulation Based on the Enhanced Assumed Strain Concept”, International Journal for Numerical Methods in Engineering. Vol. 37, pp. 2551-2568, 1994 Carleer, B. “Sheet Metyal Forming Simulations with a Frictional Model Based on Local Contact Conditions”, Numisheet 1996, pg. 40-46, 1996 Chandresekaran, N. “A Finite Element Solution Method for Contact Problems with Friction”, International Journal for Numerical Methods in Engineering, Vol. 24, pp. 477-495, 1987 Cho, Y.B., “An Improved Theory and Finite Element Model for Laminated and Sandwich Beams Using First Order Zig-Zag Sublaminate Approximations”, Composites and Structures, Vol. 37, pp. 281-298, 1997 Cook, R., Concepts and Applications of Finite Element Analysis, John Wiley & Sons, 1989 126 DeSciuva, M., “An Improved Shear Deformation Theory for Moderately Thick Multilayered Anisotropic Shells and Plates”, Journal of Applied Mechanics, ASME, Vol. 54, pp. 589-594, 1987 Gontier, C., “About the Numerical Simulation of Sheet Metal Stamping Process”, International”, International Journal for Numerical Methods in Engineering, Vol. 37, pp. 669-692, 1994 Haber, RB., “An Eulerian-Lagrangian Finite Element Approach to Large-Deformation Frictional Contact”, Computer and Structures, Vol. 20, pp. 193 - 201, 1985 Haruf, J. “The Sensitivities of Lab-Measured Friction Coeficients, JOM, July 1995 Hill, R, “Constitutive Modelling of Orthotropic Plasticity in Sheet Metals”, J. of Mech. Phys. Solids, Vol 38, pp. 405-417, 1990 Hill, R, “Constitutive Modelling of Orthotropic Plasticity in Sheet Metals”, J. of Mech. Phys. Solids, Vol 15, pp. 79-95, 1967 Huang, Y., Lu, Y., “An Analysis of the Axisymmetric and Nonaxisymmetric Sheet Stretching by a Hernispherical Punc ”, Computers and Structures, Vol. 51, 1994 Karafillis, AR, “Tooling and binder Design for Sheet Metal Forming Processes Compensating Springback Error”, Int. J. Mach. Tools Manu, Vol. 36, No. 4, pp. 503-526, 1996 Kobayashi, S. “Deformation Analysis of Axisymmetric Sheet Metal Forming Processes by the Rigid-Plasti Finite Element Meth ”, Symposium on Mechanics of Sheet Metal Forming”, General Motors Research Laboratories, Warren, MI, 1977 Lee, J .K., Cho, U.Y., Harnbrecht, “Recent Advances in Sheet Metal Forming Analysis”, ASME Advances in Finite Deformation Problems in Materials Processing and Structures, Vol. 125, 1991 Liu, W.K., “Finite Element Hydrodynamic Friction Model for Metal Forming”, International Journal for Numerical Methods in Engineering, Vol. 37, pp. 4015 - 403 7, 1994 Mellor, P.B. “Plasticity Analysis of Sheet Metal F orming”, Symposium on Mechanics of Sheet Metal Forming”, General Motors Research Laboratories, Warren, MI, 1977 Nagtegaal, J. “On the Development of A General Purpose Finite Element Program for Analysis of Forming Processes”, International Journal for Numerical Methods in Engineering, Vol. 25, pp. 113-131, 1988 127 Nakamachi, Eiji, “A Finite Element Simulation of the Sheet Metal Forming Process”, International Journal for Numerical Methods in Engineering, Vol. 25, 1988 Neale, K. ,W., Tugcu, P., “A Numerical Analysis of Wrinkle Formation Tendencies in Sheet Metals”, International Journal for Numerical Methods in Engineering, Vol. 30, 1990 Oh, S.I., Kobayashi, Shiro, “Finite Element Analysis of Plane Strain Sheet Bending”, International Journal of Mechanical Sciences, Vol. 22, 1980 Onate, E. “NUMISTAMP: A Research Project for Assessment of Finite Element Models for Stamping Processes”, Journal of Materials Proccessing Technology, Vol. 50, pp. 17-38, 1195 Pian, T., Sumihara, K, “Rational Approach for Assumed Strain Finite Elements”, International Journal of Numercial Methods in Engineering, Vol. 20, 1984 Pillinger, I, Hartley, P., et. al., “Finite Element Modeling of Metal Flow in Three- Dirnensional Forming”, International Method for Numerical Methods in Engineering, Vol. 25, 1988 Prathap, G. ,”The Displacement Type Finite Element Approach - From Art to Science”, National Aerospace Laboratories, Bangalore, India, Ref 03 76-0421(94)E0001-Z, 1994 Rama, S. “Numerical Modeling of 3-D Superplastic Sheet Foring Processes”, Symposium on Advances in Superplasticity and Superplastic Forming, The Minerals, Metals and Materials Society, 1993 Rebelo, N., Nagtegaal, J.C., “Finite Element Analysis of Sheet Forming Process”, Vol. 30, 1990 Reddy, J .N., “A Simple Higher Order Theory for Laminated Composite Plates”, Journal of Applied Mechanics, ASME, Vol. 51, pp. 745-752, 1984 Ronda, J. Colville, K. “Co-rotational Frictional Model for Metal F orming”, Numisheet ‘ 1996, pp. 47-54 Schweizerhof, K., “Improving Standard Shell Elements, Friction Models and Contact Algorithms for the Eflicient Solution of Sheet Metal Forming Problems with LS- DYNA 3D”, VDI (German Society of Engineering) FE Simulation of 3D Sheet Metal Forming Processes in Automotive Industry Conference, 1991 128 Shimizu, T., “A Modified 3-Dimensional Finite Element Model for Deep-Drawing Analysis”, ASME Advances in Finite Deformation Problems in Materials Processing and Structures, Vol. 125, 1991 Simo, J. “Improved Versions of Assumed Enhanced Strain Tri-linear Elements for 3D Deformation Problems”, Computer Methods in Applied Mechanics and Engineering, Vol 110, pp. 359-386, 1993 Simo, J. “A Class of Mixed Assumed Strain Methods and the Method of Incompatible Modes”, International Journal for Numerical Methods in Engineering, Vol. 29, pp. 1595-1638, 1990 Simo, J. “Geometrically Nonlinear Enhanced Strain Mixed Methods and the Method of Incompatible Modes”, International Journal for Numerical Methods in Engineering, Vol. 33, pp. 1413-1449, 1992 Stoughton, T.B., “Finite Element Modeling of 1008 AK Sheet Steel Stretched Over a Rectangular Punch with Bending Effects”, Computer Modeling of Sheet Metal Forming Processes, AIME, 1985 Torstenfelt, B., “Contact Problems with Friction in General Purpose Finite Element Computer Programs”, Computers and Structures, Vol. 16, pp. 487 - 493, 1983 Wang, N.M., Budiansky, 8., “Analysis of Sheet Metal Stamping by the Finite Element Meth ”, Journal of Applied Mechanics, Vol. 45, 1978 Wang, N.M., Tang, S., “Analysis of Bending Effects in Sheet Forming Operations”, International Journal for Numerical Methods in Engineering, Vol. 25, 1988 Wenner, M.L., “Elementary Solutions and Process Sensitivities for Plane Strain Sheet- Metal Forming”, Journal of Applied Mechanics, Vol. 114, 1992 Wertheimer, T., “Numerical Simulation of Metal Sheet Forming Processes”, VDI, Vol 894, pg. 517-548, 1991 Wilson, D. “Representation of Material Behaviour in Finite Element Methods of Modeling Sheet Forming Processes”, International Materials Reviews, Vol. 35, 1990 Wilson, W., Strategy for Friction Modeling in Computer Simulations of Metalforming”, 16th NAMRC (North American Manufacturing Research Conference) pg. 48-54, 1988 Wing, K. “Finite Element Hydrodynamic Friction Model for Metal Forming”, International Journal for Numerical Methods in Engineering, Vol. 37, pp. 4015-4037, 1994 129 Wriggers, P. “A Comparison of Three-Dimensional Continuum and Shell Elements for Finite Plasticity”, International Journal of Solids and Structures, Vol. 33, pp. 3309- 3326, 1996 Yang, D. “Comparitive Investigation into Implicit, Explicit and Iterative Implicit/Explicit Schemes for the Simulation of Sheet-Metal Foring Processes”, Journal of Materials Processing, Vol .50, pp. 39-53, 1995 Yip, Y.C., Averill, R.C., “A Simple Laminated Eight-Noded Brick Element with Zig-Zag Sublaminate Approximations for Composite Plate Problems”, .‘ ‘ You-Min, H, “Finite Element Analysis of Contact Problems for Sheet Metal Bending Process”, Computers and Structures, Vol. 57, pp. 15-27, 1995 You-Min, H. “An Analysis of the Axisymmetric and Nonaxisymmetric Sheet Stretching by a Hemispherical Punch”, Computers and Structures, Vol. 51, pp. 315-324, 1994 Zienkiewicz, O.C., “Elasto—Plastic Solutions of Engineering Problems ‘Initial Stress’, Finite Element Approach”, International Journal for Numerical Methods in Engineering, Vol. 1, pp. 75-100, 1969 130 Appendices 131 APPENDIX A All details regarding the subroutine implementation can be found in the MARC user manuals. The input decks are in the K62 version format. However, these input decks can be read into the newest version (K7 .2) without any errors. The following flow diagram is ofi‘ered to explain the general flow of information for the program. Wm Subroutine USELEM Global Displacement Vector ICalculatc Element Return K. R and F Tangent Stifiness (K), ‘ Int. Force Vector (R) and Nodal Equivelant Loads Vector (F) The source codes, “2d.f’ and “3d.f’, and also the input decks, “m1 11ee.dat” and “ml 1tangmod.dat” can be found in the following directory: lhome/d6/uc/smith/angela To execute compile and execute the program, the following command is required: launch -j (inputfile) -u (subroutine) where the inputfiledat and subroutinef are the input deck and subroutine names, respectively. and “launch” is a unix alias for “/path. . ./marck62” "llllllll'llllllll