Citation |

- Permanent Link:
- http://digital.auraria.edu/AA00005854/00001
## Material Information- Title:
- Software improvement to simulate pure bending and stretching of sheet metal
- Creator:
- Wang, Yang-Cheng
- Publication Date:
- 1994
- Language:
- English
- Physical Description:
- 1 online resource (xi, 133 leaves) : illustrations, tables ;
## Subjects- Subjects / Keywords:
- Sheet-metal -- Computer programs ( lcsh )
- Genre:
- theses ( marcgt )
non-fiction ( marcgt )
## Notes- Abstract:
- "An effective finite difference formulation and its computer implementation were derived by Professor Luis R. Sanchez to simulate sheet metal. His computer program was written for a main frame computer. In order to use the advantages of the improvement of the computer hardware in recent years, the computer program to simulate the pure bending and stretching of sheet metal is rewritten for a personal computer. The computer programs developed is this project are based on Dr. Sanchez's derivation of mathematical formulation and his computer algorithm..." - Abstract.
- General Note:
- Submitted in partial fulfillment of the requirements for the degree, Master of Science, Mechanical Engineering.
- Statement of Responsibility:
- by Yang-Cheng Wang.
## Record Information- Source Institution:
- University of Colorado Denver
- Holding Location:
- Auraria Library
- Rights Management:
- All applicable rights reserved by the source institution and holding location.
- Resource Identifier:
- 978292156 ( OCLC )
ocn978292156
## Auraria Membership |

Downloads |

## This item has the following downloads: |

Full Text |

/
SOFTWARE IMPROVEMENT TO SIMULATE PURE BENDING AND STRETCHING OF SHEET METAL by Yang-Cheng Wang B.S., Chinese Military Academy, 1980 M.S., University of Colorado, 1987 A thesis submitted to the Faculty of the Graduate School of the University of Colorado at Denver in partial fulfillment of the requirement for the degree of Master of Science Mechanical Engineering 1994 < o * M. ji This thesis for the Master of Science degree by Yang-Cheng Wang has been approved for the Department of Mechanical Engineering by James C. Gerdeen Da.t William H. Clohess Wang, Yang-Cheng (M.S., Mechanical Engineering) Software Improvement to Simulate Pure Bending and Stretching of Sheet Metal Thesis directed by Assistant Professor Luis R. Sanchez ABSTRACT An effective finite difference formulation and its computer implementation were derived by Professor Luis R. Sanchez to simulate sheet metal. His computer program was written for a main frame computer. In order to use the advantages of the improvement of the computer hardware in recent years, the computer program to simulate the pure bending and stretching of sheet metal is rewritten for a personal computer. The computer programs developed in this project are based on Dr. Sanchezs derivation of mathematical formulation and his computer algorithm. For pure bending of sheet metal, the bending and reversed bending have been considered. The strain history has been recorded and the moment-curvature has been presented in both tabular and graphic forms. For stretching of sheet metal, only the geometry of the wall with one pin has been studied in this project. The axial force, bending moment and the related coordinates along the sheet metal have been presented in tabular form. IV The computer programs are written based on a 486 IBM comparable com- puter, MAGTRON, with 55 MHz, 16 RAM memory and 1.4 GB hard disk space drivers. In this report, all of the listed CPU time is counted from this capable computer. This abstract accurately represents the content of the candidates thesis. I recommend its publication. Signed Luis R. Sanchez ACKNOWLEDGMENTS I would like to express my sincerest gratitude to my advisor, Professor Luis R. Sanchez, for his constant support and guidance throughout the course of this study. Gratitude is extended to Professor James C. Gerdeen, and Professor William H. Clohessy, for serving on my final examination committee. I would like to thank former superintendent of Chinese Military Academy, Taiwan, Republic of China, Lt. Gen. Jia-Chyi Hu, and former Vice Superinten- dent of Chinese Military Academy, Lt. Gen. Chung-Wu Ma for giving me the opportunity and financial support. Thanks is extended to the former Chancellor of Chinese Military Academy M. Gen. Hwa-Juh Gau for his encourage to study abroad. I would like to thank my wife, Shou-Chyn, my sons, George and William, for their love and patience. Finally, I would like to thank my parents, my sisters and my brother for their understanding and support. CONTENTS CHAPTER I. INTRODUCTION ................................................. 1 1-1 Prologue....................................................1 1-2 Historical Review ......................................... 4 1-2-1 Influence of the Mechanical Properties of the Material on Sheet forming ......................................... 4 1-2-2 Influence of the Forming Geometry .................... 7 1-2-3 Limits on the Formability of Sheets .................. 8 1-2-4 Contact and Friction Effects ...................... 9 I- 3 Problem Statement ....................................... 11 II. DESCRIPTION OF THE PROBLEM .................................. 13 II- l Introduction ........................................... 13 II-2 Idealization ............................................ 16 II-3 Definition of the Problem 16 VII II- 4 Methods of Approach A Review......................... 18 II-4-1 Analytical Methods ................................. 18 II-4-2 Theory of Beams and Shells ........................ 18 II-4-3 Numerical Methods................................... 19 II- 4-4 Finite Element Method ............................ 19 III. MODEL DEVELOPMENT ........................................... 20 III- 1 Introduction .......................................... 20 III-2 Pure Bending ........................................... 21 III- 2-1 Mathematical Formulation ......................... 22 III-2-2 Computational Algorithm .... ^.....................26 III-2-3 Computational Implementation .................... 27 III-3 Stretching for Wall-Pin Geometry ....................... 29 III-3-1 Mathematical Formulation ...........................31 III-3-2 Computational Algorithm.............................33 III-3-3 Computational Implementation ...................... 35 IV. NUMERICAL RESULTS AND DISCUSSIONS ........................... 37 vrn IV-1 Introduction .................................... 37 IV-2 Pure Bending ................................... 37 IV-3 Stretching for Wall-Pin Geometry ............... 42 IV- 4 Discussion ................................... 46 V. CONCLUSIONS AND RECOMMENDATIONS.......................47 V- l Conclusions ....................................47 V-2 Recommendations .............................. 48 APPENDIX: A. NOTATIONS ........................................ 49 B. USERS MENU OF THE PURE BENDING PROGRAM ............. 52 C. PURE BENDING INPUT DATA FILE..........................59 D. COMPUTER PROGRAM OF PURE BENDING .....................61 E. USERS MENU OF THE STRETCHING PROGRAM.................86 F. STRETCHING FOR WALL-PIN GEOMETRY INPUT DATA FILE 95 IX G. COMPUTER PROGRAM OF STRETCHING...............97 BIBLIOGRAPHY.....................................134 FIGURES Figure 2-1 Forming of the Sheet ........................................ 14 2-2 Mathematical Model .......................................... 14 2-3 Strip Deforming of Transient Case............................ 15 2- 4 Strip Deforming of Steady State Case ........................15 3- 1 The Free Body Diagram of An Element Subjected to Pure Bending Moment ....................................................... 21 3-2 The Flow Chart of Pure Bending Program ..................... 28 3-3 The Free Body Diagram of the Typical Sheet Element ........ 30 3-4 The Angle The Retaining Wall ........................................... 31 3- 5 The Flow Chart of Stretching Program with Wall-One-Pin Geometry 36 4- 1 The Moment-Curvature Curve for Pure Bending .................42 XI 4-2 The Profile of The Sheet Metal.................................. 43 4-3 The Normal Resultant Force Along the Sheet Metal................44 4-4 The Moment-Curvature Curve for Stretching Sheet Metal with Wall-One-Pin Geometry .........................................45 TABLES Table 4-1 The Comparison of the CPU Time ................................. 38 4-2 Material Properties and Initial Values ........................... 40 4-3 List Curvature and Moment for 100 Segments......................41 CHAPTER I INTRODUCTION 1-1 Prologue In recent years, formed sheet metal parts are widely used in the produc- tion of articles such as parts in aircraft, cars, refrigerators, cans, and machines, etc., which are produced in the order of millions of units per year in the United States. It is obvious that sheet metal forming is important in the worlds mod- ern economy. Any improvement of the metal forming process has an important economic impact. Eary1 and Keeler2,3 provided that a complex stamping is divided into a number of regions in which a specific type of forming operation takes place. The following basic forming operations are observed in the general change in shape:ref'4 ppA~4 1. Bending: The material on the outer surface of the bend will elongate during this type of operation and material on the inner side will shrink. No deformation occurs along the axis of the bending. 2. Bend-and-Straightening: After being bent, the material is straight- ened. The convex surface undergoes a successesive tensile (ten- sile phase) and compressive (straightening phase) sequence while the inner surface undergoes the reverse sequence. 3. Plane-Strain Stretch: A tensile strain component is added to the bending operation producing a tensile strain on the convex side and a smaller compressive strain (or a small tensile strain) on the concave side. The strain along the width remain zero. 4. Flanging: During this operation the radial zone is folded 90 to form the flange or wall, the line of bending is changed from a straight to a curved line. Two different types are identified: (a) Shrink flanging The flange length shrinks during forming, creating a compressive strain in the circumferential direction of the formed wall. This state has a tendency to produce buckles or wrinkles which are common for this type of oper- ation. (b) Stretch flanging The metal is stretched as it is flanged. A tensile strain state generates along the circumference of the wall and a minor strain ranging from zero at the bend line and a minimum negative value near the edge. This type of deformation applies to hole expansion operations as well. 5. Embossment: The metal is displaced without reducing the sheet 3 thickness in the flat section of the offset. It is performed in a localized areas such that the blank is large compared with the deformation zone. 6. Biaxial Stretch: In this forming operation, both the major and minor strain are tensile. 7. Drawing: As the material is pulled toward the die, the flange area of the blank are put into compression due to a reduction of its circumference. A radial elongation and thickening of the blank at the flanges results as the material is being drawn. Usually the basic operations interact with each other. Recently, height strength material has been more widely used instead of low-carbon mild steel. Therefore, the scientific understanding of the deform process is one of the key factors of the successful industrial forming processes. An effective finite difference formulation and its computer implementa- tion were derived by Professor Luis R. Sanchez4 to simulate sheet metal. The computer program of reference [4] was written in main frame computer. In order to use the advantages of the improvement of the computer hardware in recent years, the computer program to simulate the pure bending and stretching of sheet metal is rewritten for personal computer. The computer programs developed in this project are based on Dr. Sanchezs derivation of mathematical formulation and his computer algorithm. 4 For pure bending of sheet metal, the bending, unbending and reversed bending have been considered. The strain history has been recorded and the moment-curvature has been presented in both tabular and graphic forms. For stretching of sheet metal, only the geometry of wall with one pin has been studied in this project. The axial force, bending moment and the related coordinates along the sheet metal have been presented in graphic form. The computer programs are written based on a 486 IBM comparable com- puter, MAGTRON, with 55 MHz, 16 RAM memory and 1.4 GB hard disk space drivers. In this report, all of the listed CUP time are counted from this kind of capable computer. 1-2 Historical Review A large number of articles have discussed the factors influencing sheet metal behaviours. To survey recent developments in sheet metal forming, the following discussions can be divided into four primary categories. 1-2-1 Influence of the Mechanical Properties of the Material on Sheet Forming. Hollomon5 described the stress-strain curve obtained in simple tensile test 5 by using Equation (1-1). where: o : stress e : strain K : strength coefficient, n : strain hardening effect on the sheet due to the forming process. The constant n has been strongly considered related to formability (the ability of sheet metal to be processed into a desired shape at high speed and low cost)4 as it controls the necking strain and the strain distribution.6 Due to the strain rate which can cause increasing stress, the forming speed becomes another very important factor affecting the sheet metal behaviours. The Von Mises yield cri- terion is widely used6-9 to instigate the deformation other than the pure uniaxial tension case. The stress and strain can be expressed as the following equations W = r^2 + '3^2 + (<73 1^] U 2) 1 ~ ]jg [(Â£1 C2)2 + (Â£2 C3)2 + (^3 ~ el)2] (1_3) The effective stress-strain relationship a = 6 values obtained from uniaxial tension test (isotropic hardening criterion). Sev- eral yield criteria have been proposed10,11 but the criterion proposed by Hill7 is one of the most commonly used, possibly due to its relative simplicity. Hills approach assumes a uniformly distributed anisotropy along three mutually or- thogonal planes of symmetry. The stress-strain curve of simple tensile test can not well describe the stress-strain relationship of more complex geometry materials. Therefore, several new yield criteria have been proposed. Equation (1-4) is proposed by Hill and discussed by Mellor and Pramer.12 [ci + 0-2] + (T + 2r) 101 0-2 |m = 2 (1 + r) Fm (1 4) where: m : constant evaluated experimentally Y : the yield stress in uniaxial tension This equation has better agreement between the uniaxial and biaxial tests. In general no single criterion provided the best predictions for all combination of materials and tests. 7 1-2-2 Influence of the Forming Geometry The rotation of the principal strain axis in forming complex geometry has received more analytical attention in recent years. A number of articles1013-15 have discussed that the formability and characteristics of the deformation in some regions of the part are remarkably path dependent and highly influenced by the strain history. It also has been stated4 that the load-elongation curve, as obtained from the uniaxial tension test, is strongly dependent on the angle of the test specimens with respective to the rolling direction for prestrained sheet materials such as decarburize rimmed steel.15 The variations of the yield surface under conditions of reverse loading (Bauschinger effect) have received considerable attention in recent times. Con- sideration of this phenomenon is important since actual stamping, bending and reverse bending under tension cause a reverse state of stress and strain that has to be properly taken into account in the constitutive equations. E.Tanaka16 et al. performed a series of plastic strain controlled cyclic tests on 316 stainless steel tubes at room temperature and classified the obtained cyclic hardening curves according to the degree of hardening as: 1. simple or proportional paths, 2. cruciform and stellate paths, and 8 3. square and circular paths. They found that starting from virgin materials, the strain hardening in the proportional path is significant and the largest hardening is attained by the circular path. 1-2-3 Limits on the Formability of Sheets Limits on forming can be imposed by instabilities under compressive stress (wrinkling), defects due to lubrication (friction), and galling, inaccuracies of the formed shape due to buckling, surface wrap, surface deflection or springback, strength of the tools used in forming, etc. Instabilities are due to different factors such as localized necking fracture, fracture17, etc. Extensive bending operations can cause failure in sheet metal forming. Bending creates a strong strain gradient through the thickness ranging from a high tensile strain on the outer surface to a high compressive strain on the inner surface. These strains add to those caused by any eventual tensile force acting along the plane of the sheet. As deformation by bending is in plane strain along the width, failure will occur when the tensile strain reaches a certain limit. Springback is the result of elastic recovery that appears upon releasing the forming load. Therefore, it is important to have a method to predict the amount of springback. C. Monfort and A. Bragard18 presented a mathematical 9 model to predict springback. They claimed that friction, while significant in the calculation of the bending force, does not influence the total springback. Johnson and Singh19 studied the springback effect in spherical caps formed from circular blanks without blank holder. They found that springback decreases with the increasing of the polar angle defined by the shape. In general, as the complexity of the part shape increases, springback becomes more difficult to model. 1-2-4 Contact and Friction Effects The effect produced by the contact of the restraining surfaces (punch and die) with the sheet, and those produced by friction and lubrication at the interface are intimately coupled.15 The forces acting on the tooling and workpiece produce defections that could result in inappropriate geometric tolerances of the part, undesirable galling, scoring, etc., ultimately leading to tooling breakage or localized fracture of the sheet. The most commonly used form of accounting for friction in sheet metal forming is by the use of Amontons (Coulombs) law Fs = fiN, (1 5) where fi is the proportionality factor relating the shear and normal forces. In a study on friction and galling of sheet metals, it is obvious that when friction 10 decreased, sliding velocity increased; when friction increased with surface rough- ness in the same longitudinal direction using good lubricants were required; an increase in surface pressure appeared to decrease the friction coefficient, which seemed unaffected by the sliding length. A. Curnier20 proposed a variational formulation for the contact problem of two deformable bodies. He proposed a generalization of Coulombs law where the force of friction is dependent on: 1. the normal loading resulting in a nonlinear friction effect, 2. the isotropic or anisotropic initial rugosity of the surfaces in contact, and 3. the subsequent wear of these particles. His theory is limited to small slips or to flat surfaces of contact. J.L. Duncan, R. Sowerby, and E. Chu21 stated the influence of friction on the strain distribution between the tool and sheet as well as the effects of variations in the hold-down pressure may be more significant than factors such as plastic anisotropy or strain- rate sensitivity. Professor Sanchez4 also stated: it is obvious that an appropriate approxi- mation to contact friction problems does not depend as much on the mathematical or computational complexity of the proposed model as it does on the experimental substantiation of parameters such as force and displacement fields. 11 1-3 Problem Statement In this project, we have improved, modified and verified the existing finite difference program available from Dr. Sanchezs doctoral dissertation4 in order to perform the following tasks: Task 1 : Pure Bending In order to investigate the stress-strain relationship of the material sub- jected to the loading, unloading and reverse loading, the pure bending model has been considered in this project. The computer program simulates a hypothetical case where loading occurs along the curve f(e), and the unloading and reverse loading along the curve f*(e). The moment-curvature curve has been presented in both tabular and graphical forms and parts of the strain history has been listed in tabular form. From the moment-curvature curve, the stress-strain relationship can be easily read. Task 2 : Restraining Wall and One Pin For the sheet forming procedure, the geometry of tools can be arbitrary. The device of the sheet forming procedure can be modeled as a numbers of pins. In this project, the only geometry considered is the restraining wall and one pin. 12 The program simulates the connection between the sheet metal and the wall as well as the pin. The location of the sheet metal and the axial forces, moment along the sheet metal have been also presented in both tabular and graphical forms. CHAPTER II DESCRIPTION OF THE PROBLEM II-1 Introduction The model can be best described with referring to Figure 2-1 and Figure 2-2. Figure 2-1 shows the forming of the sheet and the Figure 2-2 shows the mathematical model. Any arbitrary die geometry is simplified by a wall or by cylindrical pins of equivalent die curvature at the point of contact. The sheet strip geometry is represented by a reference mid-surface composed of circular segments whose curvature depends on the process variables. Since the calculations are referred to the middle fiber of the sheet, the radii of the pins are increased by half the thickness of the strip at the point of contact as shown in Figure 2-3 and Figure 2-4. Figure 2-3 presents the strip deforming of transient case and Figure 2-4 presents the strip deforming of steady state case. Depending on the relative position of the wall and pins, the strip may theoretically establish line contact or completely conform to the curvature of the pin. The points of the contact are initially unknown, but the curvature of any assumed pin boundary is the same as the curvature of the tool geometry accounted for half the thickness of the strip. The pins can be remodeled at different locations if required, but usually their convenient location is indicated by the geometry of the tool. Since a general die profile is conveniently modeled by pins, the strip is 14 considered to deform within pins of known curvature and the influence of each specific pin on the deformation of the strip is determined. In order to achieve the tasks of this project, the cases described in the section 1-3 have been considered. These cases are considered to provide the basis for a further generalization of an n pin system. As will be shown, the number of iterations may increases significantly as a new pin is added depending on its effects on the preceeding pins. The present model will be thoroughly developed for the steady state case. PUNCH A Figure 2-2 Mathematical Model Figure 2-3 Strip Deforming of Transient Case Figure 2-4 Strip Deforming of Steady State Case 16 II-2 Idealization An exact formulation was made within the limitations of the following assumptions: Members are initially straight. In elastic range, the moduli of elasticity E in tension and compression are equal. The transverse plane section of the members remained plane and normal to the longitudinal fibers before and after bending. The effect of residual stresses is negligibly small. The contact surfaces are assumed rigid and friction. The width of the sheet reminds unchanged during forming (plane strain deformation). II-3 Definition of the Problem The present work is based on the assumptions described in the previous section. The geometry of the devices is arbitrary as shown in Figure 2-1. In the sheet forming procedure, the deformation starts from the flat sheet. The sheet 17 metal keeps forming until it is in contact with the tools. As the sheet is formed, the sheet is subjected to bending moment, and stretching force. The magnitude changed is related to the position of the geometry corresponding to the sheet metal and the tools. Due to the various magnitude of the forces coupled with bending moment along the sheet, the thickness of the sheet metal may be varied. According to the volume constitution, the strain history of each fiber should be traced. For pure bending moment analysis, a hypothetical case has been consid- ered in the present work. A piece of the sheet metal is considered subjected to pure bending moment. The sheet starts bending at the flat. As the sheet is bent, the corresponding curvature can be found and the corresponding bending moment can be calculated as well. Since this is a pure bending, the sheet metal should not be subjected to any pulling or stretching forces. Therefore, the equi- librium condition is that the summation of the axial forces must be zero (or in the tolerance range). As the equilibrium condition is reached, the bending moment increases to perform the new position and the corresponding curvature can be found. Unloading and reversed loading can follow the same procedure. Finally, the strain history and the moment-curvature curve can be determined. 18 II-4 Methods of Approach A Review A number of methods have been published to study the behaviors of sheet forming under bending as well as coupled with stretching. The following discus- sions are based on theoretical and numerical methods and these discussions refer to reference [4]. II-4-1 Analytical Methods Analytical methods usually are based on the perfect elastic-plastic, as well as linear hardening materials such as Swift22, Chakrabarty,23 etc. The type of analysis does not take into account any effects either due to cyclic loading on the yield strength of the materials or due to large strains. The strain history is also disregarded by this type of analysis. II-4-2 Theory of Beams and Shells The theory of beams and shells is based on the constitutive law. The ratio of the length and the thickness is another key factor of this kind of method. Several methods have been published.2425 19 II-4-3 Numerical Methods Using a digital computer to do the endless iterative procedure is an impor- tant advantage in recent years. Some research2226 proposed numerical procedures to do step-by-step incremental plasticity approach subjected to certain restric- tions such as the boundary conditions and equilibrium conditions. II-4-4 Finite Element Method The finite element method is one of the preferred methods for researches in sheet forming procedure. The most common approach to sheet forming problems is based on the variational methods. It is the basic function of the finite element method. To analyze sheet forming, several types of nonlinearities in finite element method have to be involved such as geometric nonlinearity, material nonlinearity, as well as contact problems. CHAPTER III MODEL DEVELOPMENT III-l Introduction In this chapter, the model of both pure bending and stretching coupled with bending is discussed. Basically the model of the sheet metal and its math- ematical derivation as well as its computer algorithm are based on Professor Sanchezs previous work.4 In calculations of the stresses and strains of an ma- terial element, the mathematical formulation is only expressed in pure bending case since the formulation is the same for both cases. The computer programs have been improved and written for IBM compa- rable personal computer. The computer programs have been separated by each case, pure bending and stretching, respectively. Users can use the input forma- tion described in users manuals to create the input data files. Any data in the input data file can easily be changed. It is much easier to use and it has a great advantage to users doing parametric study. Both of the programs, pure bending and stretching, have their own users manual, and they are listed in appendix B, and appendix E, respectively. In order to illustrate the user manuals, example files of these both cases are also listed in appendix C and appendix F, respectively. 21 III-2 Pure Bending In order to investigate the relationship between the moment and curvature of materials, a simple case is considered. The deformation of a wide thin strip being bent over an arbitrary radius has been considered and shown in Figure 3-1. Figure 3-1 represents the free body diagram of an element subjected to pure bending moment. Defining the mid-surface of the strip as the reference surface, a typical element i of thickness h;, and length As^ can be modeled by an arc of a circle of curvature k,- as shown in the free body diagram in Figure 3-1. (s,z) are the convective coordinates of the deformed element defined along the mid-surface arc length s, and along the thickness, z, respectively. The material element i is assumed to be composed of a number of fibers through the thickness and they are able to be stretched or compressed under the application of a stress Figure 3-1 The Free Body Diagram of An Element Subjected to Pure Bending Moment 22 III-2-1 Mathematical Formulation Assuming cross section of the typical element to remain normal during the deformation, the deformed arc length As^s, z) of any fiber can be linearly related to the curvature K; and to the arc length of the reference surface for the element i, S^: Asi(s;,z) ASnu(l -f- k;z), 2 ^ 2 (3 1) Where AsBi = i-(fli-fli_1) (3-2) -ft-i The normal force resultant Fn; and the resultant moment Mj are given by: /*/2 h h Fm=LlSiz 2***2 (3'3) fh/2 Mi = / Where the resultant moment M; is also a function of the curvature K;: Mi = M(Ki) (3-5) Since the element i is only subjected to pure bending moment, the equilibrium 23 condition is that the normal force resultant equals to zero (Fuji = 0). The relation of the forces can be expressed by equation 3-6. 5]Fj = FN, = 0 (3-6) j where Fj is the normal force acting on fiber j. Since any fiber within a typical element i stretches of compresses along the s direction, the same fiber will compress or stretch along the thickness in the some proportion due to the plane strain condition, and the total column of element iw will be kept constant under any arbitrary loadings. For the steady state case after one step increment: As>nu_j-ih;_|-i = ASmjhj As0h0 (3 7) Where, for this case, the size of the arc length interval As0 corresponds to the size of the step increment as defined at the inlet of the deformation process. The neutral axis corresponds to the location of the fiber (7]n) of the element i which did not deform after one step increment: As [si,z(77n)] = As [Smi+ijSi+i (7/n)] (3 8a) Asmi [1 -f- KjZ; (?7n)] = ^^nu+l [1 H" ^i+l^i+l (^/n)] (3 86) Where T)n is the material coordinate along the z direction corresponding to the 24 neutral fiber. 7}n can occupy a different position along the thickness depending on the loading conditions. To relate the stress is involved due to the building-unbending and reverse bending nature of the process. The effective true stress-strain curve has been provided Wang26 as shown in equation 3-9. a = Ke (3-9) Where for the present plane strain case, the stress on the fiber will be given by: _ 1 + 7 V^+2f a = Ccr (3 10) where: 7 : normal anisotropy parameter. n i+t7 ^ x/1+27 and where anisotropy has been taken into account using Hill4,s anisotropic theory. For the steady state case, Wang26 also proposed an expression for a cumulative effective strain that accounts for the strain history in a fiber: e 1 d C c (3 -11) 25 Where v (Â£, prefers to the velocity of the fiber tj located at ( along the arc length. Considering a fiber j located at (si, zj) on an element i, the cyclic loading conditions acting on the fiber j can be represented in a piecewise fashion as follows: Loading: Unloading and reverse loading: 0 < 4 < < (3 126) Where the direction of the coordinate axis has been reversed, and a new coordi- nates set (e*, o-*) is now defined with the origin at the beginning of the unloading: tively as expressed in equation 3-11, but according to the active piecewise ( history of the fiber has to be accounted for. This procedure requires the deter- mination of each ( 26 strip while it is subjected to a pure bending moment applied to its ends. III-2-2 Computational Algorithm In this procedure, the hypothetical case is subjected to pure bending mo- ment. The bending moment and the curvature of the element are unknowns. The curvature depends on the boundary conditions; that means the intial curvature K; is a trial value. Once the curvature is assumed, the corresponding resultant force Fn and the moment M; can be determined. Since this is pure bending mo- ment case, the equilibrium condition is that the normal resultant force is equal to zero. The iterative procedure starts when the curvature is assumed and it stops when the equilibrium condition is reached. The iterative procedure can be described in the following steps: 1. a trial curvature Ki 2. a trial value for the length of a finite element i after one step deformation Astmi is assumed 3. using constancy of volume (equation 3-7), the thickness of the element i for the assumed Astmi is determined. 4. the corresponding length As for each of the fibers forming the element i is determined from equation 3-1, and the location of the neutral fiber is 27 found comparing the current value Asj, with respective to the length of the same fiber one step before. 5. the natural strains can now be calculated using equation 3-12. 6. the stresses acting on each fiber can be determined from equation 3-10 from a power hardening law, or from equation 3-12 for any other type of constitutive law. 7. the resultant forces and moment acting on element i for the assumed trial Astmiare now calculated using equations 3-3, and 3-4. 8. equilibrium equation (equation 3-6) must be satisfied; otherwise, another trial Astini is assumed and steps 2 through 8 are repeated. Reiterate through the thickness until equilibrium is fulfilled. III-2-3 Computational Implementation Based on the mathematical formulation and the computer algorithm, the computer programs are written in FORTRAN to simulate the procedure of sheet metal subjected to pure bending moment. The flow chart of the pure bending program is presented in Figure 3-2. 28 Figure 3-2 The Flow Chart of Pure Bending Program 29 III-3 Stretching for Wall-Pin Geometry The model of the sheet metal under stretching for wall-pin geometry can be described in Figures from 2-1 through 2-4. The geometry of the forming tools can be arbitrary as shown in Figure 2-1. The tool can be modeled as a number of cylindrical pins as shown in Figure 2-2. The pins can be remodeled at different locations if it is necessary to indicate the location of the geometry of the tools. In this section, the only geometry considered is the retaining wall and one pin. The geometry can be seen in Figures 2-3 and 2-4. The deformation of a wide thin strip being pulled and bent over an arbitrary radius has been considered and been modeled in this section. Figure 3-3 represents the free body diagram of the typical sheet element located in the unsupported region. During the steady state deformation the element i occupies the element i+l along the curve of the unsupported region AB after a one step increment. The definitions of the geometric properties in the typical element i are the same as those of pure bending case shown in Figure 3-1. Fri is the resultant force acting on the unsupported area AB shown in Figure 3-3; Fuji is the normal forces resultant of the element i and Fvi is the shear force resultant. The (X,Y) system corresponds to a fixed Cartesian system, and (s,z) are also the same definition of the coordinate system before. represents the X coordinate of Fri at the point of contact with the wall, while 30 with the vertical. Note that, for equilibrium, the magnitude and direction of the resultant force Fpti must be the same at any point of the unsupported area AB. Figure 3-3 The Free Body Diagram of the Typical Sheet Element Located in the Unsupported Region 31 III-3-1 Mathematical Formulation The unknown angle tween the sheet metal and the retaining wall. Typical boundary conditions at the inlet are given by: @ A: Fn = Fna;Fv = Fva where F^AiFva are the resultant normal and shear components of Fri. Fna may be the result of previous loading and/or the result of friction due to the clamping on the binding. The material element i is also assumed to be composed by fibers through the thickness and they are able to be stretched or compressed under the applica- tion of a stress cr(s, z) action along the direction s. The definitions of deformed arc length As;(s, z), the normal force resultant Fn; and the results moment M; are the same as those in pure bending case. Figure 3-4 The Angle Retaining Wall 32 The equation of equilibrium of the moments with respective to the instantaneous center of curvature for the element i located at s; is: AFNi + KiAMi = 0 (3 13) The mathematical formulation of the arc length, the relation between stress and strain has been described in section Â§111-2-1. It must be remarked that given the iterative character of the present method, any type of known constitutive law can be applied to the model. From Figure 3-4: 0i = fa- sin-1 (|T^) (3 14) From Figure 3-4, the Cartesian coordinates of a segment i can be related to the resultant force acting on the unsupported area AB and the bending moment acting on i as following: Xs = XA + --- [M; FR1 sin fa (Ya Yi)] (3 15) 1? R1 cos
Yi = Y;_i + (cos 6i cos 0i_!) (3 16) |

Full Text |

PAGE 1 SOFTWARE IMPROVEMENT TO SIMULATE PURE BENDING AND STRETCHING OF SHEET METAL by' Yang-Cheng Wang B.S., Chinese Military Academy, 1980 M.S., University of Colorado, 1987 A thesis submitted to the Faculty of the Graduate School of the University of Colorado at Denver in partial fulfillment of the requirement for the degree of Master of Science Mechanical Engineering 1994 PAGE 2 This thesis for the Master of Science degree by Yang-Cheng Wang has been approved for the Department of Mechanical Engineering by Luis R. Sanchez James C. Gerdeen Date AA.#14 / 1994 PAGE 3 Wang, Yang-Cheng (M.S., Mechanical Engineering) Software Improvement to Simulate Pure Bending and Stretching of Sheet Metal Thesis directed by Assistant Professor Luis R. Sanchez ABSTRACT An effective finite difference formulation and its computer implementation were derived by Professor Luis R. Sanchez to simulate sheet metal. His computer program was written for a main frame computer. In order to use the advantages of the improvement of the computer hardware in recent years, the computer program to simulate the pure bending and stretching of sheet metal is rewritten for a personal computer. The computer programs developed in this project are based on Dr. Sanchez's derivation of mathematical formulation and his computer algorithm. For pure bending of sheet met. al, the bending and reversed bending have been considered. The strain history has been recorded and the moment-curvature has been presented in both tabular and graphic forms. For stretching of sheet metal, only the geometry of the wall with one pin has beeIi. studied in this project. The axial force, bending moment and the related coordinates along the sheet metal have been presented in tabular form. PAGE 4 IV The computer programs are written based on a 486 IBM comparable computer, MAGTRON, with 55 MHz, 16 RAM memory and 1.4 GB hard disk space drivers. In this report, all of the listed CPU time is counted from this capable computer. This abstract accurately represents the content of the candidate's thesis. I recommend its publication. Signed Luis R. Sanchez PAGE 5 ACKNOWLEDGMENTS I would like to express my sincerest gratitude to my advisor, Professor Luis R. Sanchez, for his constant support and guidance throughout the course of this 'Gratitude is extended to Professor James C. Gerdeen, and Professor William H. Clohessy, for serving on my final examination committee. I would like to thank former superintendent of Chinese Military Academy, Taiwan, Republic of China, Lt. Gen. Jia-Chyi Hu, and former Vice Superinten dent of Chinese Military Academy, Lt. Gen. Chung-Wu Ma for giving me the opportunity and financial support. Thanks is extended to the former Chancellor of Chinese Military Academy M. Gen. Gau for his encourage to study abroad. I would like to thank my wife, Shou-Chyn, my sons, George and William, for their love and patience Finally, I would like to thank my parents, my sisters and my brother for their understanding and support. PAGE 6 CONTENTS CHAPTER I. INTRODUCTION ................................................... 1 I-I Prologue .......................................................... 1 1-2 Historical Review ................................................. 4 1-2-1 Influence of the Mechanical Properties of the Material on Sheet forming ....................... 0........................ 4 1-2-2 Influence of the Forming Geometry .. ....................... 7 1-2-3 Limits on the Formability of Sheets .......................... 8 1-2-4 Contact and Friction Effects ................................. 9 1-30 Problem Statement .............................................. 11 II. DESCRIPTION OF THE PROBLEM ............................... 13 11-1 Introduction .................................................... 13 11-2 Idealization ..................................................... 16 II-3 Definition of the Problem ....................................... 16 PAGE 7 Vll II-4 Methods of Approach -A Review ........................... ... 18 II-4-1 Analytical Methods ........................................ 18 II-4-2 Theory of Beams and Shells ............................... 18 II-4-3 Numerical Methods ........................................ 19 II-4-4 Finite Element Method ... ................................ 19 III. MODEL DEVELOPMENT ......................................... 20 111-1 Int:roduction .................................................... 20 111-2 Pure Bending .................................................. 21 III-2-1 Mathematical Formulation ................................ 22 1II-2-2 Computational Algorithm .... : ............................ 26 III-2-3 Computational Implementation .................. :........ 27 III-3 Stretching for Wall-Pin Geometry .............................. 29 III-3-1Mathematical Formulation ................................ 31 III-3-2 Computational Algorithm ................................. 33 III-3-3 Computa.tional Implementation ........................... 35 IV. NUMERICAL RESULTS AND DISCUSSIONS ...................... 37 PAGE 8 viii IV-1 Introduction ..................................................... 37 IV-2 Pure Bending .................................................. 37 IV-3 Stretching for Wall-Pin Geometry .............................. 42 IV -4 Discussion ;..................................................... 46 V. CONCLUSIONS AND RECOMMENDATIONS ...................... 47 V -1 Conclusions ..................................................... 47 V-2 Recommendations .............................. ; ............ ; .. 48 APPENDIX: A. NOTATIONS ...................... ................................. 49 B. USER'S MENU OF THE PURE BENDING PROGRAM ............ 52 C. PURE BENDING INPUT DATA FILE .............................. 59 D. COMPUTER PROGRAM OF PURE BENDING .................... 61 E. USER'S MENU OF THE STRETCHING PROGRAM ............... 86 F. STRETCHING FOR WALL-PIN GEOMETRY INPUT DATA FILE. 95 PAGE 9 lX G. COMPUTER PROGRAM OF STRETCHING ....................... 97 BIBLIOGRAPHY ........................................................ 134 PAGE 10 FIGURES Figure 2-1 Forming of the Sheet ................................................ 14 2-2 Mathematical Model ................................................ 14 2-3 Strip Deforming of Transient Case ................................... 15 2-4 Strip Deforming of Steady State Case ................................ 15 3-1 The Free Body Diagram of An Element Subjected to Pure Bending Moment ............................................................ 21 3-2 The Flow Chart of Pure Bending Program ........ .................. 28 3-3 The Free Body Diagram of the Typical Sheet Element ............... 30 3-4 The Angle cPl at the Connection Between the Sheet Metal and The Retaining Wall ................................................. 31 3-5 The Flow Chart of Stretching Program with Wall-One-Pin Geometry. 36 4-1 The Moment-Curvature Curve for Pure Bending ..................... 42 PAGE 11 Xl 4-2 The Profile of The Sheet Metal ...................................... 43 4-3 The Normal Resultant Force Along the Sheet Metal .................. 44 4-4 The Moment-Curvature Curve for Stretching Sheet Metal with Wall-One-Pin .............................................. 45 PAGE 12 TABLES Table 4-1 The Comparison of the CPU Time .................................. 38 4-2 Material Properties and Initial Values ............................... 40 4-3 List Curvature and Moment for 100 Segments ........................ 41 PAGE 13 1-1 Prologue CHAPTER I INTRODUCTION In recent years, formed sheet metal parts are widely used in the produc tion of articles such as parts in aircraft, cars, refrigerators, cans, and machines, etc., which are produced in the order of millions of units per year in the United States. It is obvious that sheet metal forming is important in the world's mod ern economy. Any improvement of the metal forming process has an important economic impact. Earyl and Keeler2 ,3 provided that a complex stamping is divided into a number of in which a specific type of forming operation takes place. The following basic forming operations are observed in the general change in shape:ref.4 pp.1-4 1. Bending: The material on the outer surface of the bend will elongate during this type of operation and material on the inner side will shrink. No deformation occurs along the axis of the bending. 2. Bend-and-Straightening: After being bent, the material is straight ened. The convex surface undergoes a successesive tensile (ten- PAGE 14 sile phase) and compressive (straightening phase) sequence while the inner surface undergoes the reverse sequence. 3. Plane-Strain Stretch: A tensile strain component is added to the bending operation producing a tensile strain on the convex side and a smaller compressive strain (or a small tensile strain) on the concave side. The strain along the width remain zero. 4. Flanging: During this operation the radial zone is folded 90 to form the flange or wall, the line of bending is changed from a straight to a curved line. Two different types are identified: (a) ShriIik flanging -The flange length shrinks during forming, creating a compressive strain in the circumferential direction of the formed wall. This state has a tendency to produce buckles or wrinkles which are common for this type of oper ation. (b) Stretch flanging -The metal is stretched as it is flanged. A tensile strain state generates along the circumference of the wall and a minor strain ranging from zero at the bend line and a minimum negative value near the edge. This type of deformation applies to hole expansion operations as well. 5. Embossment: The metal is displaced without reducing the sheet 2 PAGE 15 thickness in the flat section of the offset. It is performed in a localized areas such that the blank is large compared with the deformation zone. 6. Biaxial Stretch: In this forming operation, both the major and minor strain are tensile. 7. Drawing: As the material is pulled toward the die, the flange area of the blank are put into compression due to a reduction of its circumference. A radial elongation and thickening of the blank at the flanges results as the material is being drawn. 3 Usually the basic operations interact with each other. Recently, height strength material has been more widely used instead of low-carbon mild steel. Therefore, the scientific understanding of the. deform process is one of the key factors of the successful industrial forming processes. An effective finite difference formulation and its computer implementa tion were derived by Professor Luis R. Sanchez4 to simulate sheet metal. The computer program of reference [4] was written in main frame computer. In order to use the advantages of the improvement of the computer hardware in recent years, the computer program to the pure bending and stretching of sheet metal is rewritten for personal computer. The computer programs developed in this project are based on Dr. Sanchez's derivation of mathematical formulation and his computer algorithm. PAGE 16 4 For pure bending of sheet metal, the bending, unbending and reversed bending have been considered. The strain history has been recorded and the moment-curvature has been presented in both tabular and graphic forms. For stretching of sheet metal, only the geometry of wall with one pin has been studied in this project. The axial force, bending moment and the related coordinates along the sheet metal have been presented in graphic form. The computer programs are written based on a 486 IBM comparable computer, MAGTRON, with 55 MHz, 16 RAM memory and 1.4 GB hard disk space drivers. In this report, all of the listed CUP time are counted from this kind of capable computer. 1-2 Historical Review A large number of articles have discussed the factors influencing sheet metal behaviours. To survey recent developments in sheet metal forming, the following discussions can be divided into four primary categories. 1-2-1 Influence of the Mechanical Properties of the Material on Sheet Forming. Hollomon5 described the stress-strain curve obtained in, simple tensile test PAGE 17 5 by using Equation (1-1). u=Kf? (1 1) where: u : stress f : strain K : strength coefficient, n : strain hardening effect on the sheet due to the forming process. The constant n has .been strongly considered related to formability (the ability of sheet metal to be processed into a desired shape at high speed and low cost)4 as it controls the necking strain and the strain distribution.6 Due to the strain rate which can cause increasing stress, the forming speed becomes another very important factor affecting the sheet metal behaviours. The Von Mises yield criterion is widely used6 9 to instigate the deformation other than the pure uniaxial tension case. The stress and strain can be expressed as the following equations (1 -2) (1-3) The effective stress-strain relationship u = (j (e) is very often described by the PAGE 18 6 values obtained from uniaxial tension test (isotropic hardening criterion). Sev eral yield criteria have been proposed10,1l but the criterion proposed by Hill7 is one of the most commonly used, possibly due to its relative simplicity. Hill's approach assumes a uniformly distributed anisotropy along three mutually or thogonal planes of symmetry. The stress-strain curve of simple tensile test can not well describe the stress-strain relationship of more complex geometry materials. Therefore, several new yield criteria have been proposed. Equation (1-4} is proposed by Hill and discussed by Mellor Pramer. 12 (1 -4) where: m : constant evaluated experimentally Y : the yield stress in' uniaxial tension This equation has better agreement between the uniaxial and biaxial tests. In general no single criterion provided the best predictions for all combination of materials and tests. PAGE 19 7 1-2-2 Influence of the Forming Geometry The rotation of the principal strain axis in forming complex geom .etry has received more analytical attention in recent years. A number of articles10,13-15 have discussed that the formability and characteristics of the deformation in some regions of the part are remarkably path dependent and highly influenced by the strain history. It' also has been stated4 that the load-elongation curve, as obtained from the uniaxial tension test, is strongly dependent on the angle of the test specimens with respective to the rolling direction for prest rained sheet materials such as decarburize rimmed steel.15 The variations of the yield surface under conditions of reverse loading (Bauschinger effect) have received considerable attention in recent times. Con sideration of this phenomenon is important' since actual stamping, bending and reverse bending under tension cause a reverse state of stress and strain that has to be properly taken into account in the constitutive equations. E. Tanaka 16 et al. performed a series of strain controlled cyclic tests on 316 stainless steel tubes at room temperature and classified the obtained cyclic hardening curves according to the degree of hardening as: 1. simple or proportional paths, 2. cruciform and stellate paths, and PAGE 20 8 3. square and circular paths. They found that starting from virgin materials, the strain hardening in the proportional path i& significant and the largest hardening is attained by the circular path. 1-2-3 Limits on the Formability of Sheets Limits on forming can be imposed by instabilities under compressive stress (wrinkling), defects due to lubrication (friction), and galling, inaccuracies of the formed shape due to buckling, surface wrap, surface deflection or spring back, strength of the tools used in forming, etc. Instabilities are due to different factors such as localized necklng fracture, fracture17, etc. Extensive bending operations can cause failure in sheet metal forming. Bending creates a strong strain gradient through the thickness ranging from a high tensile strain on the outer surface tQ a high compressive strain on the inner surface. These strains add to those caused by any eventual tensile force acting along the plane of the sheet. As deformation by bending is in plane strain along the width, failure will occur when the tensile strain reaches a certain limit. Springback is the result of elastic recovery that appears upon releasing the forming load. Therefore, it is importii.nt to have a method to predict the amount of springback. C. Monfort and A. Bragard18 presented a mathematical PAGE 21 9 model to predict springback. They claimed that friction, while significant in the calculation of the bending force, does not influence the total springback. Johnson and Singh 19 studied the springback effect in spherical caps formed from circular blanks without blank holder. They found that springback decreases with the increasing of the polar angle defined by the shape. In general, as the complexity of the part shape increases, springback becomes more difficult to model. 1-2-4 Contact and Friction Effects The effect produced by the contact of the restraining surfaces (punch and die) with the sheet, and those produced by friction and lubrication at the interface are intimately coupled. IS The forces acting on the tooling and workpiece produce defections that could result in inappropriate tolerances of the part, undesirable galling, scoring, etc., ultimately leading to tooling breakage or localized fracture of the sheet. The most commonly used form of accounting for friction in sheet metal forming is by the use of Amonton's (Coulomb's) law Fs = p,N, (1 -5) where 1 is the proportionality factor relating the shear and normal forces. In a study on friction and galling of sheet metals, it is obvious that when friction PAGE 22 10 decreased, sliding velocity increased; when friction increased with surface rough ness in the same longitudinal direction using good lubricants were required; an increase in surface pressure appeared to decrease the friction coefficient, which seemed unaffected by the sliding length. A. Curnier20 proposed a variational formulation for the contact problem of two deformable bodies. He proposed a generalization of Coulomb's law where the force of friction is dependent on: 1. the normal loading resulting in a nonlinear friction effect, 2. the isotropic or anisotropic initial rugosity of the surfaces in contact, and 3. the subsequent wear of these particles. His theory is limited to small slips or to flat surfaces of contact. J .L. Duncan, R. Sowerby, and E. ChU21 stated the influence of friction on the strain distribution between the tool and sheet as well as the effects of variations in the hold-down pressure may he more significant than factors such as plastic anisotropy or strainrate sensitivity. Professor Sanchez4 also stated: "it is obvious that an appropriate approximation to contact friction problems does not depend a.s much on the mathematical or computational complexity of the proposed model as it does on the experimental substantiation of parameters such as force and displacement fields." PAGE 23 11 1-3 P roblem Statement In this project, we have improved, modified and verified the existing finite difference program available from Dr. Sanchez's doctoral dissertation4 in order to perform the following tasks: Task 1 : Pure Bending In order to investigate the stress-strain relationship of the material sub jected to the loading, unloading and reverse loading, the pure bending model has been considered in this project. The computer program simulates a hypothetical case where loading occurs along the curve f("), and the unloading and reverse loading along the curve f*(e). The moment-curvature curve has been presented in both tabular and graphical forms and parts of the strain history has been listed in tabular form. From the moment-curvature curve, the stress-strain relationship can be easily read. Task 2 : Restraining Wall and One Pin For the sheet forming procedure, the geometry of tools can be arbitrary. The device of the sheet forming procedure can be modeled as a numbers of pins. In this project, the only geometry considered is the restraining wall and one pin. PAGE 24 12 The program simulates the connection between the sheet metal and the wall as well as the pin. The location of the sheet metal and the axial forces, moment along the sheet metal have been also presented in both tabular and graphical forms. PAGE 25 CHAPTER II DESCRlPTION OF THE PROBLEM 11-1 Introduction The model can be best described with referring to Figure 2-1 and Figure 2-2. Figure 2-1 shows the forming of sheet and the Figure 2 2 shows the mathematical model. Any arbitrary die geometry is simplified by a wall or by cylindrical pins of equivalent die curvature at the point of contact. The sheet striP geometry is represented by a reference mid-surface composed of circular segments whose curvature depends on the process variables. Since the calculations are referred to the middle fiber of the sheet, the radii of the pins are increased by half the thickness of the strip at the point of contact as shown in Figure 2-3 and Figure 2-4. Figure 2-3 presents the strip deforming of transient case and Figure 2-4 presents the strip deforming of steady state case. Depending on the relative position of the wall and pins, the strip may theoretically establish line contact or completely conform to the curvature of the pin. The points of the contact are initially unknown, but the curvature of any assumed pin boundary is the same as the curvature of the tool geometry accounted for half the thickness of the strip. The pins can be remodeled at different locations if required, but usually their convenient location is indicated by the geometry of the tool. Since a general die profile is conveniently modeled by piIis, the strip is PAGE 26 14 considered to deform within pins of known curvature and the influence of each specific pin on the deformation of the strip is determined. In order to achieve the tasks of this project, the cases described in the section 1-3 have been considered. These cases are considered to provide the basis for a further generalization of an "n" pin system. As will be shown, the number of iterations may increases significantly as a new pin is added depending on its effects on the preceeding pms. The present model will be thoroughly developed for the steady state case. PUNCH FORCE I SHEET FLOW Figure 2-1 of the Sheet Figure 2-2 Mathematical Model PAGE 27 A* .. ----..1------._ R Figure 2-3 Strip Deforming of Transient Case I ---i i+l Figure 2-4 Strip Deforming of Steady State Case B* .. ... ... ... 15 \ \ \ \ PAGE 28 16 II-2 Idealization An exact formulation was made within the limitations of the following assumptions: Members are initially straight. In elastic range, the moduli of elasticity E in tension and compression are equal. The transverse plane section of the members remained plane and normal to the longltudinal fibers before and after bending. The effect of residual stresses is negligibly small. The contact surfaces are assumed rigid and friction. The width of the sheet reminds unchanged during forzmng (plane strain deformation ). II-3 Definition of the Problem The present work is based on the assumptions described in the previous section. The geometry of the devices is arbitrary as shown in Figure 2-1. In the sheet forming procedure, the deformation starts from the flat sheet. The sheet PAGE 29 17 metal keeps forming until it is in contact with the tools. As the sheet is formed, the sheet is subjected to bending moment, and stretching force. The magnitude changed is related to the position of the geometry corresponding to the sheet metal and the tools. Due to the various magnitude of the forces coupled with bending moment along the sheet, the thickness of the sheet metal may be varied. According to the volume constitution, the strain history of each fiber should be traced. For pure bending moment analysis, a hypothetical case has been consid ered in the present work. A piece of the sheet metal is considered subjected to pure bending moment. The sheet starts bending at the fiat. As the sheet is bent, the corresponding curvature can be found and the corresponding bending moment can be calculated as well. Since this is a pure bending, the sheet metal should not be subjected to any pulling or stretching forces. Therefore, the equi librium condition is that the summation of the axial forces must be zero (or in the tolerance range). As the equilibrium condition is reached, the bending moment increases to perform the new position and the corresponding curvature can be found. Unloading and reversed loading can follow the same procedure. Finally, the strain history and the moment-curvature curve can be determined. PAGE 30 18 11-4 Methods of Approach -A Review A number of methods have been published to study the behaviors of sheet forming under bending as well as coupled with stretching. The following discus sions are based on theoretical and numerical methods and these discussions refer to reference [4]. 11-4-1 Analytical Methods Analytical methods usua.lly are based on the perfect elastic-plastic, as well as linear hardening such as Swift 22 Chakrabarty,23 etc. The type of analysis does not take into account any effects either due to cyclic loading on the yield strength of the materials or due to large strains. The strain history is also disregarded by this type of analysis. 11-4-2 Theory of Beams and Shells The theory of beams and shells is based on the constitutive law. The ratio of the length and the thickness is another key factor of this kind of method. Several methods have been published.24,25 PAGE 31 19 11-4-3 Numerical Methods Using a digital computer to do the endless iterative procedure is an important advantage in recent years. Some research22 ,26 proposed numerical procedures to do step-by-step incremental plasticity approach subjected to certain restric tions such as the boundary conditions and equilibrium conditions. 11-4-4 Finite Element Method The finite element method is one of the preferred methods for researches in sheet forming procedure. The most common approach to sheet forming problems is based on the variational methods. It is the basic function of the finite element method. To analyze sheet forming, several types of nonlinearities in finite element method have to be involved such as geometric nonlinearity, material nonlinearity, as well as contact problems. PAGE 32 111-1 Introduction CHAPTER III MODEL DEVELOPMENT In this chapter, the model of both pure bending and stretching coupled with bending is discussed. Basically the model of the sheet metal and its math ematical derivation as well as its computer algorithm are based on Professor Sanchez's previous work.4 In calculations of the stresses and strains of an ma terial element, the mathematical formulation is only expressed in pure bending case since the formulation is the same for both cases. The computer programs have been improved and written for IBM compa rable personal computer. The computer programs have been separated by each case, pure bending and stretching, respectively. Users can use the input forma tion described in user's manuals to create the input data files. Any data in the input data file can easily be changed. It is much easier to use and it has a great advantage to users doing parametric study. Both of the programs, pure bending and stretching, have their own user's manual, and they are listed in appendix B, and appendix E, respectively. In order to illustrate the user manuals, example files of these both cases are also listed in appendix C and appendix F, respectively. PAGE 33 21 III-2 Pure Bending In order to investigate the relationship between the moment and curvature of materials, a simple case is considered. The deformation of a wide thin strip being bent over an arbitrary radius has been considered and shown in Figure 3-1. Figure 3-1 represents the free body diagram of an element subjected to pure bending moment Defining the mid-surface of the strip as the reference surface, a typical element "i" of thickness hi, and length Llsmi can be modeled by an arc of a circle of curvature kj as shown in the free body diagram in Figure 3-1. (s,z) are the convective coordinates of the deformed element defined along the mid-surface arc length s, and along the thickness, z, respectively. The material element "i" is assumed to be composed of a number of fibers through the thickness and they are able to be stretched or compressed under the application of a stress O"(s, z) action along the direction s. Figure 3-1 The Free Body Diagram of An Element Subjected to Pure Bending Moment PAGE 34 22 III-2-1 Mathematical Formulation Assuming cross section of the typical element to remain normal during the deformation, the deformed arc length .lsI (s, z) of any fiber can be linearly related to the curvature Ki and to the arc length of the reference surface for the element "i", Smi: Where h h --< z <-2 --2 1 .ls = -((j. -(j. 1) IIU Ki 1 1-(3 -1) (3 -2) The normal force resultant FNi and the resultant moment Mi are given by: jh/2 FNi = udz -h/2 h h -- PAGE 35 23 condition is that the normal force resultant equals to zero (FNi = 0). The relation of the forces can be expressed by equation 3-6. (3 -6) where Fj is the normal force acting on fiber "j". Since any fiber within a typical element "i" stretches of compresses along the s direction, the same fiber will compress or stretch along the thickness in the some proportion due to the plane strain condition, and the total column of element "i" will be kept constant under any arbitrary loadings. For the steady state case after one step increment: (3 -7) Where, for this case, the size of the arc length interval b.so corresponds to the size of the step increment as defined at the inlet of the deformation process. The neutral axis corresponds to the location of the fiber (7]n) of the element "i" which did not deform after one step increment: (3 -Sa) (3 -8b) Where 7]n is the material coordinate along the z direction corresponding to the PAGE 36 24 neutral fiber. 1/n can occupy a different position along the thickness depending on the loading conditions. To relate the stress u to the strains resulting from the displacement fields shown in equation 3-1, a constitutive equation is required. This relationship is involved due to the building-unbending and reverse bending nature of the process. The effective true stress-strain curve has been provided Wang26 as shown in equation 3-9. (3 -9) Where for the present plane strain case, the stress on the fiber will be given by: where: u = 1 +1 u = Cu Jl + 'Vi 1 : normal anisotropy parameter. (3 10) and where anisotropy has been taken into account using Hill4 's anisotropic theory. For the steady state case, Wang26 also proposed an expression for a cumulative effective strain that accounts for the strain history in a fiber: "(s,1])= 'oos c 1 dV((,1/) dC j" v (C, 1/) dC (3 11) PAGE 37 25 Where v ("", )refers to the velocity of the fiber", located at along the arc length. Considering a fiber "j" located at (Si, Zj) on an element "i", the cyclic loading conditions acting on the fiber "j" can be represented in a piecewise fashion as follows: Loading: (3 -12a) Unloading and reverse loading: (3 -12b) Where the direction of the coordinate axis has been reversed, and a new coordinates set (f*, 0"*) is now defined with the origin at the beginning of the unloading: --0" = -0"; e* = -e+eo (3 -12c) The total effective strain ej (sj,Zj) of the material fiber "j" is not taken cumula-tivelyas expressed in equation 3-11, but according to the active piecewise (0'*, f*) curve for that specific fiber. Therefore, on the general case the entire strain history of the fiber has to be accounted for. This procedure requires the determination of each (0"*, e*) curve corresponding to each reversible of the loading. The procedure consists of monitoring the concave and convex surface fiber of the ----------------_._. -- PAGE 38 26 strip while it is subjected to a pure bending moment applied to its ends 111-2-2 Computational Algorithm In this procedure, the hypothetical case is subjected to pure bending mo ment. The bending moment and the curvature of the element are unknowns. The curvature depends on the boundary conditions; that means the intial curvature Ki is a trial value. Once the curvature is assumed, the corresponding resultant force FN and the moment Mi can be determined. Since this is pure bending moment case, the equilibrium condition is that the normal resultant force is equal to zero. The iterati"ve procedure starts when the curvature is assumed and it stops when the equilibrium condition is reached. The iterative procedure can be described in the following steps: 1. a trial curvature Kl 2. a trial value for the length of a finite element "i" after one step deformation is assumed 3. using constancy of volume (equation 3-7), the thickness of the element "i" for the assumed is determined. 4. the corresponding length for each of the fibers forming the element "i" is determined from equation 3-1, and the location of the neutral fiber is PAGE 39 27 found comparing the current value ASj, with respective to the length of the same fiber one step before. 5. the natural strains can now be calculated using equation 3-12. 6. the stresses acting on each fiber can be determined from equation 3-10 from a power hardening law, or from equation 3-12 for any other type of constitutive law. 7. the resultant forces and moment acting on element "i" for the assumed trial AStmiare now calculated using equations 3-3, and 3-4. 8. equilibrium equation (equation 3-6) must be satisfied; otherwise, another trial AStmi is assumed and steps 2 through 8 are repeated. Reiterate through the thickness until equilibrium is fulfilled. III-2-3 Computational Implementation Based on the mathematical formulation and the computer algorithm, the computer programs are written in FORTRAN to simulate the procedure of sheet metal subjected to pure bending moment. The How chart of the pure bending program is presented in Figure 3-2. PAGE 40 Change Tr i a 1 Sm Find Neutral Axis yes Printout Data yes urvatur Figure 3-2 The Flow Chart of Pure Bending Program 28 PAGE 41 29 111-3 Stretching for Wall-Pin Geometry The model of the sheet metal under stretching for wall-pin geometry can be described in Figures from 2-1 through 2-4. The geometry of the forming tools can be arbitrary as shown in Figure 2-1. The tool can be modeled as a number of cylindrical pins as shown in Figure 2-2. The pins can be remodeled at different locations if it is necessary to indicate the location of the geometry of the tools. In this section, the only geometry considered is the retaining wall and one pin. The geometry can be seen in Figures 2-3 and 2-4. The deformation of a wide thin strip being pulled and bent over an arbitrary radius has been considered and been modeled in this section. Figure represents the free body diagram of the typical sheet element located in the unsupported region. During the steady state deformation the element "i" occupies the element "i+l" along the curve of the unsupported region AB after a one step increment. The definitions of the geometric properties in the typical element "i" are the same as those of pure bending case shown in Figure 3-1. FRl is the resultant force acting on the unsupported area AB shown in Figure 3-3; FNi is the normal forces resultant of the element "i" and FVi is the shear force resultant. The (X,Y) system corresponds to a fixed Cartesian system, and (s,z) are also the same definition of the coordinate system before. XA represents the X coordinate of FRl at the point of contact with the wall, while l is the angle formed by FR! PAGE 42 30 with the vertical. Note that, for equilibrium, the magnitude and direction of the resultant force FRI must be the same at any point of the unsupported area AB. Figure 3-3 The Free Body Diagram of the Typical Sheet Element Located in the Unsupported Region PAGE 43 31 1II-3-1 Mathematical Formulation The unknown angle cPI is usually related to parameters defined by the boundary conditions. Figure 3-4 represents the angle cPI at the connection between the sheet metal and the retaining wall. Typical boundary conditions at the inlet are given by: @ A: FN = FNA; Fv = FVA where FNA, FVA are the resultant normal and shear components of FRI. FNA may be the result of previous loading and/or the result of friction due to the clamping on the binding. The material element "i" is also assumed to be composed by fibers through the thickness and they are able to be stretched or compressed under the application of a stress 0'( S, z) action along the direction s. The definitions of deformed arc length ASi(S,Z), the normal force resultant FNi and the results moment Mi are the same as those in pure bending case. ------1L------_-j{ I I Figure 3-4 The Angle cPI at the Connection Between the Sheet Metal and the Retaining Wall PAGE 44 32 The equation of equilibrium of the moments with respective to the instantaneous center of curvature for the element "i" located at Si is: (3 13) The mathematical formulation of the arc length, the relation between stress and strain has been described in section III-2-1. It must be remarked that given the iterative character of the present method, any type of known constitutive law can be applied to the model. From Figure 3-4: (3 14) From Figure 3-4, the Cartesian coordinates of a segment "i" can be related to the resultant force acting on the unsupported area AB and the bending moment acting on "i" as following: (3 15) (3 16) Where the parameters FRl and XA represent the inlet conditions at the point A and the abscissa of A, respectively. PAGE 45 33 III-3-2 Computational Algorithm In this case, the deformation of the unsupported region has been con cerned. Since the region is unsupported, the resultant force FRl is constant in both magnitude and direction along the whole region. Because the resultant force FRl and the coordinates of A in horizontal direction are unknowns, therefore, the intial values are assumed and iterative procedures are taken. The curvature is depended on the boundary conditions; that means the intial curvature Ki is also a trial values. The iterative procedure can be described in the following steps: 1. a trial curvature Kl 2. a trial value for the length of a finite element "i" after one step deformation Llstmi is assumed 3. using constancy of volume (equation 3-7), the thickness of the element "i" for the assumed Llstmi is determined. 4. the corresponding length Lls for each of the fibers forming the element "i" is determined from equation 3-1, and the location of the neutral fiber is found comparing the current value Llsj, with respective to the length of the same fiber one step before. 5. the natural strains can now be calculated using equation 3-12. PAGE 46 34 6. the stresses acting on each fiber can be determined from equation 3-10 from a power hardening law, or from equation 3-12 for any other type of constitutive law. 7. the resultant forces and moment acting on element "i" for the assumed trial now calculated using equations 3-3, and 3-4. 8. equilibrium (equation 3-13) must be satisfied; otherwise, another trial is assumed and steps (3-8) are repeated. Reiterate through the thickness until equilibrium is fulfilled. 9. the angle (Jj for the element "i" for the corresponding trial value FRI (equation 3-12) is determined; (see Figure 3-4). 10. the length at the middle surface segment "i" (equation 3-2) is determined: 11. compare the trial for which equilibrium was fulfilled (step 8) with the calculated length (step 10). The trial curvature Kj is now readjusted and the steps 2 through 11 are repeated until converges to the value 12. The curvature of the sheet metal is determined by the trial value of the resultant force FRI' If the sheet metal contacts with the pin, the curvature --'--'----------------- PAGE 47 35 of the sheet should be the same .as that of the pin. As the sheet metal contacts with the pin, the coordinate of the sheet metal in Y-direction is below the top of the pin, i.e., Yi < CBY -R where: Yi is the coordinate of the sheet metal at "i" element in Y-directionj CBY is the coordinate at the center of the pin in Y-directionj R is the radius of the pin. In present procedure, once the sheet metal contacts with the pin, the iteration starts. The iterative procedure will not stop until the difference between these two curvatures, the curvature of the sheet metal and the curvature of the pin, is less than tolerance. 111-3-3 Computational Implementation Based on the mathematical formulation and the computer algorithm, the computer program is written in FORTRAN to simulate the procedure of stretch ing the sheet metal with the geometry of retaining wall and one pin. The flow chart of the stretching program is illustrated in Figure 3-5. PAGE 48 Change Tr i a 1 Sm Change No Curvaturet--..::...;....;"-< (Ki ) yes Find Angles yes Output ncr ease ria 1 FRI ecrease Tr i a 1 FRI es 36 Figure 3-5 The Flow Chart of Stretching Program with Wall-One-Pin Geometry PAGE 49 CHAPTER IV NUMERICAL RESULTS AND DISCUSSIONS IV-I Introduction Based on the mathematical formulation and its computer algorithm as well as its computer implementation, both cases, the pure bending case and stretching case, have been studied. The typical material element has been considered as composed of "n" fibers. In this iterative procedure, several parameters need to be assumed as initial values. The intial values always affect the speed of the convergence; they may even affect the numerical stabilities. According to parametric study, the better intial guess has been recommend by the author for using these programs. The numerical results have also been presented in both graphical form and tabular forms. IV -2 Pure Bending In oder to examine the speed of the convergence under different numbers of composed fibers for the normal cross section, Table 4-11ists the comparison of the CPU time for different numbers of the composed fiber. PAGE 50 38 II Table 4-1: The Comparison of the CPU Time II n m Time(sec) 200 30 over 300 180 30 40 160 30 30 140 30 35 120 30 45 100 30 30 80 30 60 60 30 75 40 30 over 300 n : Number of Fibers m : Number of Curvature Increments Time : CPU Time in sec The first column of Table 4-1 represents the number of the fibers along the normal cross section shown in Figure 3-1. The second column of Table 4-1 is the number of the curvature increment in which the increment is the same The iterative procedure starts at the curvature equal to zero and it stops at curvature equal to the maximum which is 1/0.1875. The curvature increment being used PAGE 51 39 here is not an equal space. It is an exponential function. The function is expressed as: Ki = A X exp (B x N) (4 -1) where: Ki curvature at ith step A : constant determined by initial and final conditions B : constant determined by initial and final conditions N : the total number of increment The initial and final conditions: { @N = 0 Ki = initial curvature (Co = 0) @N = n Ki = maximum curvature (= (4 -2) The coefficient of A and B can be determined by solving equations 4-1 and 4-2 in which A is 1.0 and B is: n=l (4-3) Table 4-1 shows that it is not necessary to get less CPU time if the user increases the number of fibers along the cross section. According Table 4-1, the number of PAGE 52 40 fibers ranging from 80 to 180 is reasonable. One hundred fibers is recommend by the author for using pure bending program. The material properties being used in this case are listed in Table 4-2. All of the symbolic in Table 4-2 are referred to the mathematical formulations and Figures in the previous definitions as well as user's manuals. Table 4-2: Material Properties and Initial Values Material Initial Properties Values K (psi) 7833.4 SmO (in) 0.25 n 0.228 Ho (in) 0.2 '1 0 177 Co 0.0 Ex (psi) 20E6 Dco 1.0E-6 Yield (psi) 20E3 Ncurve 30 Another very important result of this case is the relationship between the moment and its corresponding curvature. By taking the results of Tables 4-1 and 4-2, the material properties in Table 4-2 have been used in this case study and the number of fibers has been taken to be 100. Table 4-3 lists the numerical results. The first column of Table 4-3 is the number of fibers which is 100. The second column of Table 4 3 is CPU time (in seconds) which is counted by using a 486 IBM compatible personal computer with 50 MHz, 16 RAM memory and 1.4 GB hard PAGE 53 41 drivers. The third column represent the curvatures of the typical element. The last column of Table 4-3 is the corresponding moment to the curvature listed in the third column. Table 4-3: List Curvature and Moment for 100 fibers n Fibers CPU time(sec) Curvature Moment 0.0001 0.274 0.0006 1.791 0.0014 3.795 0.0062 17.030 0.0130 36.078 100 30 0.0590 56.975 0.1250 66.289 0.2640 77.062 0.8170 96.492 1. 7290 112.202 3.6640 129.838 In graphic form, Figure 4-1 represents the moment-curvature curve. In Figure 4-1, the horizontal axis is curvature and the vertical axis is moment. Figure 4-1 shows when the element is unloading and reverse loading, the magnitude of the linear range is twice as much as the first loading. PAGE 54 42 This is a lest file (n=100 &: n=180) 140 120 100 80 'i:i' i 60 I I J--I -I I ---I I L i I I ... 40 d Q) a 20 0 ::s 0 -20 -40 I I I V I / I I Y ---I --I I i -60 -0.5 o 0.5 1.5 2 2.5 3 3.5 Curvature Figure 4-1 The Moment-Curvature Curve for Pure Bending IV -3 Stretching for Wall-Pin Geometry The same material listed in Table 4-2 is being used here again but the initial guesses are not completely the same as those listed in Table 4-2. More initial guesses need to be used in this case. The trial values are listed in appendix F. According to the geometric and material properties as well as the initial values, --'---'-------------------- PAGE 55 43 the iterative procedure has been completed. Figure 4-2 represents the profile of the sheet metal after forming. In Figure 4 -2, the horizontal axis is the x-axis and the vertical axis is the y-axis of the Cartesian coordinate system. At the left end, the sheet metal contacts with the retaining wall. At the right end, the sheet metal contacts with the pin. g I=l 0 :J 0 QI 1: Ai I>-:::l .. QI > QI :l I>CI I=l 0 :;;: I=l 0 :::l III CI 0 ..:l QI .d Eo< 6 5.9 5.8 5.7 5.6 5.5 5.4 5.3 5.2 o The Profile of the Sheel Metal I I ----"" I I I I I I I 0 5 1 1.5 2 The Location Along the Horizontal Direction (in) Figure 4-2 The Profile of the Sheet Metal. \. 1\ I I I 2.5 PAGE 56 44 The between the two ends is unsupported area. The resultant force, FR1, is a constant in the unsupported region but the normal resultant force, FNi along the sheet metal is various. Figure 4-3 represents the normal resultant force along the sheet metal in x-direction. Figure 4-3 shows that when the sheet connects with the retaining wail, the normal resultant force is the friction force; when the sheet is close to the pin, the normal resultant force decreases; the normal resultant force even becomes negative, and the relationship is nonlinear. ., CJ .. o r.. 4 0 20 o OJ -20 >4 < OJ -40 S .. o Z ., -60 t: -80 -100 o The Normal Axial Force VB. The Horizontal Location of Sheet Metal I i I I \ I \ I \ I 1\ I I 1 \ I I i I I 0.5 1.5 2 2 5 3 The Horizontal Location Along the Sheet Metal (in) Figure 4-3 The Normal Resultant Force Along the Sheet Metal PAGE 57 '1 i .. 1:1 4> e 0 ::I't 500 450 400 350 300 250 200 150 100 50 o o Moment-Curvature Curve Along the Sheet Metal I I I I J---t-I -----I / / ( I / I I I I I / I I I I 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Curvature (l/in) Figure 4-4 The Moment-Curvature Curve for Stretching Sheet Metal with Wall-One-Pin Geometry 45 As discussed in the case of pure bending moment, the moment-curvature cure is an important factor for sheet forming behavior. Figure 4-4 represents the momentcurvature curve for the sheet metal along the unsupported area. Figure 4-4 shows that when the sheet metal is closer to the retaining wall, the moment-curvature keeps linear; when the sheet metal is closet to the pin, the moment-curvature PAGE 58 46 curve becomes nonlinear. IV -4 Discussion According to the numerical results, the most critical point of this proce dure is convergence. The initial guesses also cause the speed of convergence and numerical scheme stability. Another important factor designed by users is the tolerances. Some tolerances does not very seriously affect the accuracy but they do seriously delay the CPU time, such as the tolerances of the equilibrium of both pure bending and stretching cases. Therefore, the user needs to be familiar with the programs. Then, he can much more effectively use the advantages of the procedure to simulate this type of work. Since this procedure has several parameters, the effective initial guesses and stable numerical schemes are also important to deal with this type of work. PAGE 59 CHAPTER V CONCLUSIONS AND RECOMMENDATIONS V-I Conclusions Based on the work of this project, the following conclusions can be made: The numerical iteration is the most critical of the present procedure. An arbitrary geometry can be modeled by the numbers of pins and be solved by the present procedure. Compared to finite element methods, this present procedure is much more effective to solve this type of problem. In nonlinear finite element analysis, the nonlinearities can be divided into four primary categories: geometric nonlinearity, material nonlinearity, force nonlinearity and boundary non linearity (contact problem). In this problem, at least three of the four categories of nonlinearities are involved. It is a challenge to solve three kinds of nonlinearities simultaneously. Although both finite element meth ods and present procedure have to make numerical iteration, the present procedure still has the following advantages: 1. No meshing definition is needed explicitly. 2. No solving of explicitly different types of nonlinearities simultaneously is required. PAGE 60 48 3. No matrix manipulation is necessary. In present method, the number of fibers along the thickness may have one hundred but infinite element methods, the number of elements along the thickness may have only a few elements (maybe less than five elements). V -2 Recommendations Because the geometry of sheet metal device can be arbitrary, this type of analysis also can be extended to an arbitrary number of pins to simulate the behavior of the stretching of the sheet metal. The mathematical formulation and computational algorithm ha.ve been derived by Professor Luis R. Sanchez.4 In further research, in order to run the program much more effectively by us ing personal computer it seems necessary to improve the software4 for arbitrary geometry (numbers of pins and wall) .. In this project, the sheet metal was a.ssumed as plane strain. In further re search work, the three dimensional analysis seems an effectively way to reduce the idealization (assumptions). In three dimensional analysis, all of the pins can be modeled as either cylinders or spheres. The expanse analysis of the mathematical formulation is also necessary. PAGE 61 APPENDIX A NOTATIONS PAGE 62 50 NOTATIONS Chapter I (j = stress E = strain K = strength coefficient n = strain hardening effect on the sheet due to the forming process C7"j = principal stress in direction i Ej = principal strain in direction i m = contact evaluated experimentally Y = the yield stress in uniaxial tension v = proportionality factor relating the shear and normal forces PAGE 63 Chapter III Kj smi 1/n 'f C (e*, u*) fj (Sj, Zj) XA -l 51 Curvature of element "i" Length of element "i" at mid-surface Normal force resultant resultant moment the normal force acting on fiber "j". The neutral axis corresponds to the location of the fiber of the element "i" normal anisotropy parameter. the origin at the beginning of the unloading The total effective strain of the material fiber "j" the X coordinate of FRI at the point of contact with the wall the angle formed by FRI with the vertical PAGE 64 APPENDIX B THE USER'S MANUAL OF PURE BENDING PROGRAM PAGE 65 53 THE USER'S MANUAL OF PURE B .ENDING PROGRAM For this case, the inp1,lt data file name is called "INPUT .DAT". All of the input format of this program is "free format". The user can use comma(,) or space to separate two sequential parameters or data. At the first line of each control card, there is the parameter's description and only the first four characters will be read. All of the input data shoUld be associated with their descriptions. For every control card, it is NOT necessary to follow the menu sequence (i.e., on can specify like: material properties, initial values, "', or he can specify like: initial values, material properties, ... ). I TITLE TITLE : Title of this file (arbitrary characters up to 72 ) Example: This is a example input data file of pure bending PAGE 66 54 I CONTROL CONTROL: There are six cases PURE = pure bending moment (Only available now) Example: pure bending moment PAGE 67 55 I MATE COEF EN R EX YIELD MATE: material properties (only the first four characters will be readied) COEF : the coefficient, k of the equation (j = k Iflen EN : the power of the above equation R : the normal anisotropy parameter EX : Young's Modulus YIELD : yielding stress Example: material properties 506, 0.228, 0.177, 30E6, 36E3 PAGE 68 56 ----_______________ 1 SMO HO CURVO DCURVO NCURV INIT : .initial values of the parameters SMO : the original length of central fiber of the beam HO : the original thickness of the beam CURVO : the original curvature of the beam DCURVO : the original increment of curvature NCURV : the total number of the increment of curvatUre Example: initial values 0.25, 0.2, 0, 0.000001, 50 PAGE 69 57 I PARA NETAO FNTOL SM2TOL NLOOP PARA : parameters NETAO : -the total number of eta (1/) FNTOL : the tolerance of the summation of the stresses SM2TOL : the tolerance of the length of the central fiber NLOOP : the number of the loops of the moment-curvature curve Example: parameters 100, 0.001, 0.001, 3 PAGE 70 58 I END OF FILE END OF FILE: end of file Example: end of file PAGE 71 APPENDIX C PURE BENDING INPUT DATA FILE PAGE 72 PURE BENDING INPUT DATA FILE This is only a test for the program. pure bending moment material propeties 73386.51,0.2,0.7,30000000;25000 inj,.tial values 0.40,0.1,0,0.0001, 30 parameters 200,0.8, 0.02, 3 end of file 60 PAGE 73 APPENDIX D COMPUTER PROGRAM OF PURE BENDING PAGE 74 COMPUTER PROGRAM OF PURE BENDING program main C ============================================================= C C = THIS PROGRAM IS TO DETERMINE: C (1). MATERIAL PROPERTIES RELATIONSHIP BETWEEN PURE C BENDING MOMENT AND CURVATURE (PUREM) C (2) C CREAT DATE OCT. 16, 1993 C MODIFIED DATE : C AUTHOR YANG-CHENG WANG C C ============================================================= C include control. l' include 'para.l' C = OPEN AN OUTPUT DATA FILE OUTPUT DATA CAN BE USED IN C = MATLAB AND OPEN AN INPUT DATA FILE INPUT.DAT C C C open(unit=40,:file='input.dat',status='old') open(uni t=35 ,:file=' out .m' ,status= 'nell') open (unit.=45 ,:file='mk. dat status='nell') open(unit=55,:file='shistory.dat',status='nell') C = CALL SUBROUTINE INPUT C THE PURPOUSE: C READ AND CALCULATE THE ORIGGINAL PARRAMETERS C call input C i:f((llhat .eq. 'pure') .or. (llhat .eq. 'PURE' then 62 PAGE 75 C C = CALL SUBROUTINE PUREM C THE PURPOSE: C DETERMINE THE MATERIAL PROPERTIES RELATIONSHIP BETWEEN C THE MOMENT AND CURVATURE C + + + call purem else if vhat .ne. 'pure') .or. (what .ne. 'none') .or. (vhat ne PURE') or (vhat ne NONE' ) ) then call input err stop end if stop end 63 PAGE 76 C C = C C C = C = C = C = C C DECK OF SUBROUTINE INPUT subroutine input include 'control.l' include 'para.l' MPRO MATERIAL PROPERTIES INVA INITIAL VALUES PARA PARAMETERS END END OF FILE read(40,999)title read(40,1000)vhat write(45,999)title if ( (what eq ) or (vhat eq. PURE' then vrite(45,5000) end if vrite.(45,6000) 100 read(40,1000)vhich vhich .e9,. 'mate') .or. (vhich .eq. 'MATE'then read(40,*)coef, en, r, ex, yield vrite(45,1001)vhich vri tee 45, *)J coef = J, coef vri t e ( 45 *) en =', en vrite(45.*)' r ='. r vrite(45,*)' .. ex =' ex vrite(45.*)' yield =' yield go to 100 .. else ifvhich .eq. 'iriit').or.(which .eq. 'INIT'then read(40,*) smO, hO." curvO. dcurvO, ncurv vrite(45.1001)vhich vrite(45.*)' smO smO vrite(45.*)' hO =' hO vrite(45.*)' curvO ='. curvO vrite(45.*)' dcurvO =', dcurvO write (45,*) ncurv =', ncurv vrite(45,*)J======================================J 64 PAGE 77 999 1000 1001 5000 6000 2000 go to 100 else ifwhich .eq. 'para') .or.(which .eq. 'PARA'then read (40, *) netaO" fntol, sm2tol, npathe write(45,1001)which writ, e (45, *) netaO =', netaO write(45,*)' fntol =' fntol write(45,*)' sm2tol =" sm2tol wri te (45, *) npathe =' npathe write(45,*)'======================================' go to 100 else if which eq. end ,). or. (which eq. "END ,) )then go to 2000 else if( (which .ne; 'mate')' (which .ne. init') + and. (which .ne. 'MATE') '.and. - PAGE 78 C C = DECK OF SUBROUTINE INPUTERR C subroutine input err print *. print *. vrite(45.1000) 1000 format(' + return end Your option is NOT included.' Check that again !!!!!' Your option is NOT included.' .1. Check that againn!!!! !') 66 PAGE 79 C C = DECK OF SUBROUTINE PUREM C subroutine purem C =====================================================.======== C CREATED DATE C MODIFIED DATE C COMPUTER C AUTHOR SEP. 20, 1993 NOV. 7, 1993 PC VERSION YANG-CHENG WANG C ============================================================= C ============================================================= C THIS PROGRAM IS DETERMINING THE RELATHIONSHIP BETWEEN THE C MOMMNT AND CURVATURE C ============================================================= C C C = ALL OF THE COMMON AND DIMENSION EXTERNAL FILE C C ============================================================= include 'para.l' C ============================================================= C C CALCULATE THE STREE-STRAIN RELATIONSHIP AND THEIR C COEFFICIENTS C C ============================================================= C ****** C C C r = PARAMETER.FOR THE RELATION OF STRESS.AND STRAIN = THE MATERIAL PROPERTIES C ============================================================= C CALCULATE THE COEFFICIENT C C ============================================================= c = (1.0d+00 + r) / .0d+00+ 2.0 d+OO r)**0.5) C ============================================================= 67 PAGE 80 C C CHANGE ALL OF THE INITIAL VALUES TO BE THE FIRST STEP DATA C C ============================================================= C ****** JCOUNT = COUNT FOR THE NUMBER OF INTERATION C ****** CURVi = CURRENT CONFIGURATION OF CURVATURE C ****** CURV2 = UPDATE CONFIGURATION OF CURVATURE C ****** SMO = LENGTH OF THE NETURAL LINE OF ORIGINAL C CONGIGURATION C ****** SMi C = LENGTH OF THE NETURAL LINE OF CURRENT CONGIGURATION C ****** SM2 C = LENGTH OF THE NETURAL LINE OF UPDATE CONGIGURATION C ****** C ****** C ****** C ****** C ****** HO Hi H2 NETAO NETA = = = = = THICKNESS OF THE OF ORIGINAL CONGIGURATION THICKNESS OF THE OF CURRENT CONGIGURATION THICKNESS-OF THE OF UPDATE CONGIGURATION NUMBER OF ETA OF ORIGINAL CONGIGURATION NUMBER OF ETA OF CURRENT CONGIGURATION C ======_=========-============================================== C TRANNSFER THE INITIAL DAD AS TO THE CUNNRENT CONFIGURATION C ============================================================= jcount= 0 jcounti = 1 curvi = curvO sm1 = smO nsm2 = 0 csm2 =sm2tol hi = hO neta = netaO xneta = neta deta = i/xneta C =============:================================================= C CALCULATE ETA_. C ===============a============================================= do 100 i = 1, neta+1 e .ta1 (1) = (i-i)*deta eta2(i) = (i-1)*deta 68 PAGE 81 100 continue C ============================================================= C INITIALIZE ALL OF THE STRESS AND MOMENT ARGUMENTS C ============================================================= do 101 i = 1, neta dn(i) = O.Od+OO dmn(i) = O.Od+OO epsilon(i) = O.Od+OO 101 continue n = O.Od+OO n = O.Od+OO totalm = O.Od+OO C ============================================================= C CALCULATE THE COEFFICIENT FOR THE CURVARURE C ============================================================= C CURV2 = A EXP( B*NCURV) C NCURV = 0, ===> CURV2 = DCURVO C NCURV = NCURV, ===> CURV2 1/0.1875 C ============================================================= b = log(1/(0.1875*dcurvO/(ncurv-1) do 102 i = 1, neta+1 zO(i) = hi (i-i) deta -h1/2.0 102 continue C ============================================================= C C ****** C C C CALL SUBROUTINE POSITIVE TO CALCULATE THE RELATIONSHIP BETWEEN MOMENT AND CUVATURE FOR POSITIVE CUVATURE C ============================================================= call positive return end 69 PAGE 82 C ============================================================= C C = DECK OF SUBROUTINE POSITIVE C THE PURPOSE: C CALCULATE THE RELATION BETWEEN THE CURVATURE AND MOMENT C FOR POSITIVE INCREMENT OF CURVATURE C C ============================================================= subroutine positive include 'para.l' C ============================================================= C C = INCREASE CURVATURE C C ============================================================= 1111 c jcount = jcount + 1 curv2 = dcurvO + (1.0/1.876 -dcurvO)/xneta*(jcount-1) if(ncurv ge 31) then if ( (jcount eq. (ncurv -2 or + (jcount .eq. (ncurv .4 .or. + (jcount eq. (ncurv -6 .or + (jcount .eq. (ncurv -8 .or. + (jcount eq. (ncurv -10 .or + (jcount eq (ncurv -12 or + (jcount eq. (ncurv -14 .or + (jcount .eq. (ncurv 16 ) then go to 1111 else curv2 = dcurvO*expjcount-1)*b) end if else if(ncurv .ge. 21) then if ( (jcount .eq. (ncurv -2 .or. + (jcount .eq. (ncurv -4 .or. + (jcount .eq. (ncurv -6 ) then go to 1111 else curv2 = dCUrvO*expjcount-1)*b) end if-70 PAGE 83 else curv2 = dcurvQ*expjcount-1)*b) end if sm2 = sm1 C ============================================================= C C = CALCULATE SM C THE PURPOSES: C (1). FIND THE Z(ETA) C (2). DETERMINE THE LENGTH OF EACH LAYER FIBER C (3). DETERMINE THE STRAINS FOR EACH LAYER C (4). CALCULATE THE STRESS FOR EACH LAYER C C ============================================================= C ****** Zi = THE THICKNESS AT. EACH FIBER OF CURRENT C CONFIGURATION C ****** Z2 C = THE THICKNESS AT EACH FIBER OF UPDATE CONFIGURATION C ****** DEAT = INCREMENT OF ETA C ****** ZCOM3A= PARAMETER FOR CALCULATION Z1 C ****** ZCOM4A= PARAMETER FOR CALCULATION Zi C ****** ZCOM3 = PARAMETER FOR CALCULATION Z2 C ****** ZCOM4= PARAMETER FOR CALCULATION Z2 C ============================================================= C = SET UP THE UPDAT THICKNESS h2 TO BE hi C = AND CALCULATE Zi AND Z2 C ============================================================= 2222 h2 = (h1*sm1)/sm2 if (curvi .eq. 0.0) then do 200 i = 1, neta+i z1(i) = h1* (i-i) deta h1/2 200 continue else do 300 i = 1, neta+1 zcom3a(i) = (1.0d+00-(0.50d+00*curv1*h1**2 zcom4a(i) = sqrt(zcom3a(i) + 2.0d+00 + hi curv1 eta1(i z1(i) = (-1 + zcom4a(i / curv1 71 PAGE 84 300 continue C C C end if if(curv2 .eq. 0 .. 0) then do 490 i = 1, neta +1 z2(i) = h2* (i-i) deta h2/2 400 continue else do 500 i = 1, neta+1 zcom3(i) = .(1.0d+00 (0.50d+00*curv2*h2**2 zcom4(i) = sqrt(zcom3(i) + 2.0d+00 +* h2 curv2 eta2 (i z2(i) = (-1 + zcom4(i / curv2 500 continue end if C ============================================================= C CALCULATE THE LENGTH OF FIBER FOR, EACH LAYER C ==============="============================================== C ****** S1 = THE LENGTH OF EACH FIBER AT DIFFERENT LAYER C OF CURRENT CONFIGURATION C ****** S2 = THE LENGTH OF EACH FIBER AT DIFFERENT LAYER C OF UPDATE CONFIGURATION C ****** EPSILON2 = STRAIN OF EACH LAYER OF ,UPDATE C C ****** SG2 C ****** YIELD CONFIGURATION = STRESS OF UPDATE CONFIGURATION = YIELD STRESS OF THE MATERIAL C ****** ZNETURAL = THE LOCAtION OF NETURAL AXIS C (AT THE FIBER OF STRESS EQUAL TO ZERO) C ****** SM2NEW = PARAMETER FOR 'NEW SM2 C ****** DSM2 = DIFFERENCE BETWEEN NEW SM2 AND SM2 C ****** SM2TOL = TOLERANCE OF SM2 C ============================================================= C = CALCULATE THE LENGTH OF EACH FIBER C ============================================================= do 600 i = 1 neta+1 s1(i) = sm1*(1+curv1*z1(i s2(i) = sm2*(1+curV2*z2(i check(i)= s2(i) -s1(i) 72 PAGE 85 600 continue do 700 i = 1, neta+1 C ============================================================= C = INITIALIZE THE STRAIN ARGUMENTS C ============================================================= epsilon2(i) = 0.0 C ======================================='====================== C = CALCULATE THE STAINS EPSILON2 C ============================================================= epsilon2(i) = log(abs(s2(i)/s1(i) sig2(i) = c ex* epsilon2(i) C C = CHECK THE ELASTICC RANGE C .It. sig2 (i) = -c coet abs (epsilon2 (i **en else if(sig2(i) .gt. yield) then sig2(i) = c coef abs (epsilon2 (i **en end if C ============================================================= -. C = CALCULATE THE NEUTRAL AXIS C it(check(i-:1) .It. 0 0 and. check(i) .gt. O.O)then znetural == z2(i-1). + deta*abs(sig2(i-1)/ + (sig2 (i-1)-sig2 (i) ) ) sm20ld = 8m2 sm2nev = sm2*(1 + curv2 znetural) dsm2 = abs(sm2nev-8m2) if (dsm2 .le sm2tol) then sm2 = sm2nev else 8m2 = 8m2nev 73 PAGE 86 C ======================.======================================= C = IF IT IS nOT. CONVERGED THEN GOTO. 2222 TO MODIFIEDh2 C ============================================================= go to 2222 end if else if (check(i) .eq. 0.0) then sig2 (i) =0.0 znetural = z2(i) end if 700 continue C ============================================================= C C = CALCULATE FM C THE PURPOSES: C(l). DETERMINE THE FORCE FOR EACH LAYER C (2). CALCULATE THE SUMMATION FORCE FOR THE STRUCTURE C (3) DETEltMINE THE MOMMENT ABOUT THE NETURAL FOR EACH LAYER C (4). CALCULATE THE MOMENT ABOUT THE NETURAL AXIS FOR THE C .. STRUTURE C C ================:d::============================================ C ****** DFN. = THE STRESS OF EACH LAYER FOR.UPDATE C CONFIGURATION C FN = THE TOTAL STRESS C ****** DMN = THE MOMENT OF EACH LAYER ABOUT ZNETURAL C FOR UPDATE CONFIGURATION C ****** TOTLAM = THE TOTAL MOMENT ABOUT THE ZNETURAL C ====,;==============,;.===================,,;,===================== C THE EQUILIBRIUM EQUATION IS SUM(DFN) = FN = 0.0 C .'. C ============================================================= C = INITIALIZE.THE fn AND totalm C ============================================================= fnold =n fn 0.0 totalm = 0.0 c ============================================================= 74. PAGE 87 C = CALCULATE STRESS(dfn) AND MOMENT(dmn) FOR.EACHFIBER C = AND TOTAL STRESS(fn) AND TOTAL MOMENT (totalm) C ============================================================= 3333 do 800 i =1. neta dfn(i) = sig2(i) abs(z2(i) z2(i+1 dmn(i) = dfn(i) z2(i) fn = fn + dfn(i) = totalm + dmn(i) 800 continue C ============================================================= C C = CKECK THE EQUILIBRIUM EQUATION SUM(FN)=?O C IF YES, WRITE IT OUT AND INCREASE (OR DECREASE THE C CURVATURE) C IF NOT, SUM(FN) < 0.0, INCREASE SM C IF NOT, SUM(FN) > 0.0, DECREASE SM C C ============================================================= C ****** FNTOL = THE TOLERRANCE OF THE FN if (abs(fn) .le. fntol) then C ============================================================= C = THE MAXIMUM CURVATURE IS (1/0:1875) C ============================================================= if ( curv2 ge. % .1875) -0.00001) or. + j count ge. (ncurv -2 then vork = 'fini' end if call vriteout C C = UPDATE ALL F THE PARAMETERS C ============================================================= 8m2 = smO if (vork.eq. 'fini') then C C = C .HECK THE LOADING AND UNLOADING CONDITIONS C = NPATHE : 'NUMBER OF THE PATHES C = NPATHE = 1 DR 0 (DR INDICATING NOTHING) THEN C ONLY THE FIRST ,LOADIN PATH (POSITIVE. 75 PAGE 88 C MOMENT) C = NPATHE = 2 POSITIVE MOMENT + 1 UNLOADING C (NEGATIVE MOMENT) C = NPATHE = 3 POSITIVE LOADING + NEGATIVE MOMENT + C POSITIVE LOADING C ===================================='========================= else if (npathe eq. ,2) then call pathe2 else if(npathe .eq. 3) then call pathe2 call pathe3 else ifnpathe.eq.0).or.(npathe.eq.1then stop else stop end if C ============================================================= C = GOTO 1111 TO INCREASE CURVATURE C ============================================================= go to 1111 end if C ============================================================= C = IF FN NOT SATISFY fn TOLERANCE THE MODIFY sm2 C = fn SUPPOSE TO BE ZERO (EQUILIBRIUM EQUATION) C to1, = fn/(fno1d-fn) else if(fn .gt. fntol .and. fno1d .gt. fnto1) then if (abs(to1) .le. 30.0) then sm2 = sm2 -to1 (sm201d-sm2new) else, sm2= sm2 dsm2/2.0 end if go to 2222 else if(-fn .gt. fnto1 .and. -fno1d .gt. fnto1) then if (abs(to1) .1e. 30.0) then sm2 = sm2 -to1 (sm201d-sm2new) 76 PAGE 89 else sm2 = sm2 + dsm2/2. 0 end if go to 2222 else iffn .gt. fntol .and. -fnold.gt. fntol) .or. + (-fn.gt. fntol and. fnold .gt. fntolthen else end if return end sm2 = (sm2new+sm2old)/2.0 go to 2222 go to 2222 77 PAGE 90 c ================.============================================= c C = SUBROUTINE INPUTERROR C C ============================================================= 100 + + + + subroutine input error write(45.100) write(*. 100) format (, I I. Fatal Error 111111 Your input data formation is not right' 1 !!!!!'.II. return end Check that again. I am sorry!! 1 ) 78 PAGE 91 C ============================================================= C C = DECK OF SUBROUTINE WRITEOU! C C ============================================================= 100 101 200 300 400 500 subroutine .vriteout include 'para.l' vrite(45.100) vrite(45.200)totalm vrite(45.300)curv2 vrite(45.400)znetural vrite(45.500)sm2 format(' format(' + '====================================' .111) format(' The Moment is :'.2x. e14.7) format ( The curvature is : '. 2x. e14. 7) format ( The netural axis is at : 2x; e14. 7) format (' The length of netural axis is :'. 2x. e14. 7) C ============================================================= C HERE IS STRAIN HISTORY C ============================================================= if (jcount .eq. 1) then vrite(55.600) vrite(55.601)jcount else vrite(55.601)jcount end if vrite(55.101) vrite(55.602) do 1000 i = 1. net. a vrite(55.603) i. dfn.(i). dmn(i). s2(i) 1000 continue vrite(55.604) fn, totalm 600 format(' This is the output data file for strain' + history.'. /) 79 PAGE 92 601 + 602 format('1',' This is the',1x,i4,1x,'step of strain history.') formate' # Fiber ','1',' M I Fn I' iteration' I S2 + + + '===========================================' + '===============================') 603 format (4x,i4, 2(1x,' I', 2x, e17.9),1x,'I' ,2x,e20.12) 604 format ( Total Fn = e14. 7, + Total Moment = ,e14.7,1111111) C ============================================================= C HERE IS CREATING THE DATA FILE.FOR MATLAB TO PLOT M-K PLOT C ============================================================= if (jcount1 .eq. 1) then vrite(35,703) vrite(35,J04)totalm, curv2 else if (vork .ne. 'fini') then vrite(35,705)totalm, curv2 else vrite(35,706)totalm, curv2 end if 703 format(/, 'Ye',' Total Moment + curvature ,I, + 'Yo',' ========================' + '==========================' ,II) 704 format ('mk1=[' 2(2x, e20:12 705 format (7x, e20.12, 2x, e20.12) 706 format (7x, e20.12, 2x, e20.1"2,'];')C ============================================================= C = SAVE THE YIELDING-MOMENT(YHOMENT) AND YIELD CURVATURE C (YCURVE) C if( (abs( sig2(1) ) .ge. yield ) .or. + (abs(sig2(neta + 1) ) .ge. yield) ) then goto 212 else ymoment = totalm ycurve = curv2 nyield = jcount1 212 end if 80 PAGE 93 c ============================================================= C = SAVE THE INFORMATION OF THE FIRST LOADING C ============================================================= savem(jcount1) = totalm savec(jcount1) = curv2 savez(jcount1) = znetural saves (jcount1) = sm2 jmax = jcount1 jcount1 = jcount1 + 1 return end 81 PAGE 94 C ============================================================= C C = DECK OF. SUBROUTINE OF PATHE2 C C ============================================================= subroutine pathe2 include 'para.l' double precision p2curv2(maxstep), p2totalm(maxstep) double precision p2znetu(maxstep), p2sm2(maxstep) common/ path2/ p2curv2, p2totalm, p2znetu, p2sm2 do 100 i = 1, nyield p2curv2(i) = savec(jmax) -savec(i) p2totalm(i) =savem(jmax) -savem(i) p2znetu(i) = savez(jmax) -savez(i) p2sm2 (i) = saves (jmax) -saves (i) 100 continue kkk = 1 + nyield jjj = jmax + nyield do 200 i = kkk, jjj p2curv2(i) = savec (jmax) -savec(i-nyield) + -savec(nyield) p2to1;alm(i) = savem(jmax) -savem(i-nyield) + -savem(nyield) p2znetu(i) = savez(jmax) -savez(i-nyield) + -savez(nyield) p2sm2(i) = saves (jmax) -saves(i-nyield) + -saves(nyield) 200 continue call patheout return end 82 PAGE 95 C ============================================================= C C = DECK OF SUBROUTINE OF PATHE3 C C ============================================================= C C subroutine pathe3 include 'para.l' double precision p2curv2(maxstep), p2totalm(maxstep) double precision p2znetu(maxstep), p2sm2(maxstep) double precision p3curv2(maxstep), p3totalm(maxstep) double precision p3znetu(maxstep), p3sm2(maxstep) common/ path2/p2curv2, p2totalm, p2znetu, p2sm2 common/' path3/ p3curv2, p3totalm, p3znetu, p3sm2 jjj = jmax + nyield do 100 i = 1, nyield. p3curv2(i) = p2curv2(jjj) + savec(i) p3totalm(i) = p2totalm(jjj) + savem(i) p3znetu(i) =p2znetu(jjj) + savez(i) p3sm2(i) = p2sm2(jjj) + saves(i) 100 continue .+ + + + kkk= 1 + nyield jjj = jmax + nyield do 200 i = kkk, p3curv2(i) p3totalm(i) p3znetu(i) p3sm2(i) jjj = p2curv2(jjj) + savec(nyield) = p2totalm(jjj) + savem(nyield) = p2znetu(jjj) +savez(nyield) = p2sm2(jjj) + saves(nyield) 200 continue call patheotit return end + savec(i-nyield) + savem(i-nyield) + savez(i-nyield) + saves(i-nyield) 83 PAGE 96 C ============================================================= C C = DECK OF SUBROUTINE PATHEOUT C C ============================================================= C C subroutine patheout include 'para.l' double precision p2curv2(maxstep). p2totalm(maxstep) double precision p2znetu(maxstep). p2sm2(maxstep) double precision p3curv2(maxstep). p3totalm(maxstep) double precision p3znetu(maxstep). p3sm2(maxstep) common/ path2/ p2curv2. p2totalm. p2znetu. p2sm2 common/ path3/ p3curv2. p3totalm. p3znetu. p3sm2 C ============================================================= C HERE IS CREATING THE DATA FILE FOR MATLAB TO PLOT M-K PLOT C============================================================= jjj = jmax + nyield if npathe eq 2) or. (npathe eq. 3 then do 1 00 i = 1. j j j if (i .eq. 1) then vrite(35,703) vrite(35,704)savem(jmax). savec(jmax) vrite(35.705)p2totalm(i), p2curv2(i) else if(i .ne. jmax+nyield) then vrite(35,705)p2totalm(i), p2curv2(i) else vrite(35.706)p2totalm(i), p2curv2(i) end if 100 continue end if if(npathe .eq. 3) then vrite(35,*)'% Should be continued!' do 200 i = 1, jjj if (i .eq. 1) then vrite(35,703) vrite(35,7044)p2totalm(jjj),p2curv2(jjj) vrite(35,705)p3totalm(i), p3curv2(i) 84 PAGE 97 200 else if (i ,ne, jjj) then write(36.705)p3totalm(i). p3curv2(i) else write(35.706)p3totalm(i). p3curv2(i) end if continue end if if ((npathe ,ne, 3) write(36.*) write(36.*) end if ,and, (npathe ,ne, 2 then Your path number is not right,' Check that again !' c ************************************************************* 703 704 7044 706 706 + + + format(/. 'Y.' Total Moment curvature / Y.' .' =======.=================' '==========================' .11) format ('mk2=I' 2 (2x .e20 ,12 format (fInk3= [' 2 (2x, e20 ,12 e20,12, 2x, e20,12) format (7x, e20,12, 2x, e20,12,'] ;') return end 85 PAGE 98 CONTROL.L character*70 title character*4 vhat, vhich common/ control/title, vhat, vhich 86 PAGE 99 C C = PARAMETERS AND DIMENSION C C parameter (maxstep = 1000. maxdim = 500) character*4 work double precision curvO. curv1. curv2. dcurvO. b. c double precision coef. en. r. ex. yield double precision smO. smi. sm2. dsm2. sm2new. sm2tol double precision netaO. fntol. hOt hi. h2 double precision fn. totalm. ffn. ftotalm. eta. deta double precision znetural. csm2, ymoment, ycurve double precision savem(maxstep), savec(maxstep) double precision savez(maxstep), saves(maxstep) double precision zcom3(maxdim), zcom4(maxdim) double precision zcom3a(maxdim). zcom4a(maxdim) double precision si(maxdim) s2(maxdim), + epsilon(maxdim) + epsilon2(maxdim), check(maxdim) double precision dfn(maxdim), dmn(maxdim) + fdfn(maxdim), fdmn(maxdim), sig2 (maxdim) + zi(maxdim), z2(maxdim), zO(maxdim) double precision etai(maxdim) eta2(maxdim) integer jcount, jcounti. jmax, nyield integer neta, nerr, nsma, ncurv, npathe C = COMHMONS C + + + common/ common/ common/ common/ common/ common/ common/ common/ common/ common/ charac/worlt mp/ coef, en, r, ex, yield ivi/ curvO, curvi. curv2. ncurv. dcurvO.b.c iv2/ smO. 8mi. sm2. dsm2, sm2new. sm2tol pai/ netaO, fntol, nerr. jcount pa2/ fn. totalm. eta; deta. csm2. 8m2. ffn. ftotalm pa3/ znetural. hO, hi. h2. jcounti pa4/zcom3. zcom4.zcom3a, zcom4a paS/ 81, s2 ex1/epsilon2. check. epsilon 87 PAGE 100 + + common/ common/ common/ common/ orce1/ dn, dmn, sig1, sig2, z1, z2, dn, dmn twoeta/ eta1, eta2 savel savem, savec, savez, saves paths/ jmax, npathe, ymoment, ycurve, nyield 88 PAGE 101 APPENDIX E THE USER'S MANUAL OF STRETCHING PAGE 102 87 THE USER'S MANUAL OF STRETCHING For this case, the input data file name is called "INPUT .STR". All of the input format of this program is "free format". The user can use comma(,) or space to separate sequential parameters or data. At the first line of each control card, there is the parameter's description and only the first four characters will be read All of the input data should be associated with their descriptions. For every control card, it is NOT necessary to follow the menu sequence (i.e., on can specify like: material properties, initial values, ... or he can specify like: initial values, material properties, ... ). I TITLE TITLE: Title of this file (arbitrary characters up to 72 ) Example: This is a example input data file of stretching PAGE 103 I CONTROL I GEOMETRY CONTROL: There are six cases STRE = stretching GEOMETRY : There are six cases WALL = wall-one-pin (Only available now) TWOP = wall-two-pin (Not available now) THRE == wall-three-pin (Not available now) Example: stretching wall 88 PAGE 104 89 1 MATE COEF EN EX YIELD MATE: material properties (only the first four characters will be readied) COEF : the coefficient, k of the equation (T = k IElen EN : the power of the above equation R : the normal anisotropy parameter EX : Young's Modulus YIELD : yielding stress Example: material properties 506, 0.228, 0.177 j 30E6, 36E3 PAGE 105 90 IINIT SMO HO CURVO COEFMU DCURVE FVO INIT : initial values of the parameters SMO : the original length of central fiber of the beam HO : the original thickness of the beam CURVO: the original curvature of the beam COEFMU: the coefficient of friction DCURVE : the origin8.l increment of curvature FVO : the trial value of the vertical load at the wall Example: initial values 0.25, 0.2, 0, 0.15, 0.000001, 200.0 PAGE 106 91 [PARA NETAO PARA : parameters Example: parameters 100 PAGE 107 I GEOM Dl RPI GEOM : geometry Dl: the vertical distance between the wall and the top of the pin RPl: the radius of pin 1 Example: geomet ry 7.0, 5.6 92 PAGE 108 93 I TOLE SM2TOL ANGLETOL EQUITOL TOLE : tolerance SM2TOL : tolerance of segment length ANGLETOL : tolerance of angles EUITOL : tolerance of equilibrium condition Example: tolerance 0.001, 0.001, 0.1 PAGE 109 94 I END OF FILE END OF FILE: end of file Example: end of file PAGE 110 APPENDIX F STRETCHING FOR WA.LL-PIN GEOMETRY INPUT DATA FILE PAGE 111 96 STRETCHING FOR WALL-PIN GEOMETRY INPUT DATA FILE This is just a test file stretching wall-one-pin initial values 0.25. 0.0 0.150. 0.0000001. 25.0 material properties 73386.51.0.228. 1.77. 29000000.20000 geometry 0.41, 5.6 parameter 100 tolerance 0.5. 0.5. '0.5 end of file PAGE 112 APPENDIXG COMPUTER PROGRAM OF STRETCHING PAGE 113 98 COMPUTER PROGRAM OF STRETCHING program stretch C ================================================================ C C = THE PURPOSE OF THIS PROGRAM C THIS PROGRAM IS TO SIMULATE THE STRETCHING OF SHEET METAL. C THE WALL-PIN GEOMETRY IS ONLY CONSIDERED IN THIS PROJECT. C IN THE FUTHER,-THE FOLLOWING GEOMETRIES WILL BE COVERED: C (1). WALL-PIN C (2): TWO-PIN C (3). C (4). N-PIN C -FINALY, THIS NMETHODS WILL BE EXTENTED TO THREE DIMENSIONAL CCASE. C == C C C C C CREAT DATE MODIFIED DAtE FINAL MODIFIED DATE AUTHOR _ADVISOR March 17, 1994 YANG-CHENG WANG LUIS R. SANCHEZ C ================================================================ include 'paramet.str' c -. C = OPEN AN INPUT DATA FILE CALLED "INPUT.STR" C = OPEN AN-OUTPUT DATA FILE CALLED "OUTPUT.M" AND ITS DATA CAN C BE USED IN 'MATLAB' SOFTWARE C = OPEN ANOTHER OUTPUT DATA FILE CALLED "OUTPUT.STR" TO STALL C THE STRAIN HISTORY ALONG THE SHEET METAL C open(unit=40, fil,.e='l.nput.str' status='old') open (unit=35 ,file='output .aaa', status='nev', + carriage control = li st ) open(unit=45, file='output.str', status='nev', PAGE 114 + carri'agecontrol = 'list') file='profile.m' status='new', + carriagecontrol = 'list') open (unit=65, file='fm.m' status='new', + carriagecoritrol = 'list') open (unit=75, file='output.m' ,status='new', + carriagecontrol = 'list', ) C C = CALL SUBROUTINE INPUT. STR C THE PURPOSE OF THIS SUBROutINE IS: CREAD AND CALCULATE THE INITIAL DATA FILE. C CALL INPUT C C = CALL SUBROUTINE STRETCH C 'THE PURPOSE OF THE SUBROUTINE IS: C CALCULATE AND DETERMINE THE SHEET METAL BEHAVIOR AND C PROPETERIES UNDER STRETCHING. C C C C C C C + if geometry .eq. (geometry eq call vallpin else if geometry (geometry calltvopin else if geometry (geometry call threepin else 'vall') .or. 'WALL' then eq. 'twop' ) .or .eq. 'TWOP' then ,.eq. 'thre') .or. .eq. 'THRE') )thEtn call inputerror(2) 'endif C C = END OF THE PROGRAM C 99 PAGE 115 stop end 100 PAGE 116 101 C C = DECK OF 'SUBROUTINE WALL-PIN C THE PURPOSE OF THIS StmROUTINE IS TO CALCULATE THE STRAIN C HISTORY THE AXIAL FORCE, MOMENT AND POSITION ALONG THE SHEET C METAL WITH WALL-PIN GEOMETRY. C C subroutine wallpin include 'paramet.str' C = TRANSFER ALL OF THE DATA FROM THE INITIAL'DATA TO THE C CURRENT CONFIGURATION DATA C call data write(55,5501)xa,ya C ================================================================ C = CALL SUBROUTINE SEGMENT C PURPOSES: C (1). CALCULATE EQUILIBRIUM CONDITION FOR EACH SEGMENT C (2). WRITE DOWN THE INFORMATION FROM SEGMENT C ================================================================ 1001 icount = icount + 1 1002 call segment write(45,2001)icount write(*,2001) icount ybold = yb xbold = xb call coord curvbold = curvb = 1.0d+00 / (rp1 + 0.5 h2) PAGE 117 C C = IF THE SHEET (yb) REACH THE TOP OF THE FIRST PIN, C AND THETA_b > THETA_i > 0, CHECK THE CONVERGENCE C THETA_b THETA_i = 0 (111) C if ( (yb .It. rpi) .and. ,(curv2 .gt. 0.0) )then C C = MAKE SOME PARAMETERS IN ORDER TO BE, USED IN C NEWTON-RAPHSON METaODS C fvOold2 = fvOold fvOold = fvO dtheold = dtheta theold = thetab thetab = acos (yb/rpi) dtheta = thetab -theta C C = IF THE ANGLES CONVERGENT INTO THE TOLERANCE, IT IS OK! C GO TO RETURN C + + write(45 ,2002) write(*,2002) if (abs(dtheta) .It. angletol) then work = 'fini' write(45,2003) xb,yb, fnb2, frai, fnai, totalm2, theta, curv2 write(*,2003) xb, yb, fnb2, frai, fnai. totalm2, theta, curv2 102 PAGE 118 xb. yb go to 1099 C C = IF THE ANGLES ARE NOT CONVERGE. USE NEWTON-RAPHSON METHODS C else if(curvb .gt. curv2) then write(46.2006) xb. yb. fnb2. fra1. fna1. + totalm2. theta. curv2 + + + + + else write(*.2006) xb. yb. fnb2. fra1. fna1. totalm2. theta. curv2 fvO = fvO 1.15 call newdata write(46.2004) fvO write(56.5504) fvO wri te (* 2004) fvO go to 1002 if (dtheta gt. 0.0) and (dtheold .le 0.0 .or. dtheta .. le. 0.0) .and. (dtheold .gt. 0.0 ) write(46.2006) xb. yb. fnb2. fra1. fna1. totalm2. theta. curv2 write(*.2006) xb. yb .fnb2. fra1. fna1. totalm2. theta. curv2 fvO = 0.5 (fvOold2 + fvOold) call newdata write(45.2005) fvO fvO write(*.2005) fvO goto 1001 103 then PAGE 119 + + + + + + + else if dtheta eq. dtheold) or. (dtheold .eq. 0.0 then write(45,2006) xb, yb, fnb2, fra1, fna1, totalm2. theta. curv2 write(*.2006) xb, yb. fnb2. fra1. fna1. totalm2. theta. curv2 fvO = fvO 1. 1 call newdata write(45.2005) fvO write(55,5505) fvO write(*,2005) fvO goto 1001 else fv01 = fvO fvO = fvO -(dtheta) (fvOold fvOold2) / (dtheta -dtheold) if(fvO.lt. 0.0) then fvO = fv01 1.1 end if call newdata write(45,2006) xb yb, fnb2. fra1, fna1, totalm2. theta, curv2 write(*,2006) xb, yb. fnb2. fra1. fna1. totalm2, theta. curv2 write(46,2005)fvO write(66,6505) fvO write(*.2006) fvO goto 1001 end if 104 PAGE 120 end if endif c C = IF THE yb NO. T RREACH THE TOP OF THE FIRST PIN YET, C = CONTINUE TO DO THE NEXT SEGMENT C else write(45,2007) xb, yb. fnb2. fra1. fna1. + tota1m2. theta, curv2 write(*.2007) xb. yb. fnb2. fra1. fna1. + totalm2 theta. curv2 wri t 8 (55,5502) xb. yb call trans goto 1001 end if C ****** END THE1STLOOP 2001 formate' ========================================' 2002 2003 + + + The '. I4,' segment has been finished ,) format ( Check the tangency! ) format (, The tangency has been reached. ,I, Write out all of the data!' .1. 105 + 5x, xb -, -, Ela.l0. 5x, 'yb = E18.10,!, + 5x. fnb2 =' E18.l0. 5x. 'fral = E18 .10,!. + 5x, fnal =' E18.10, 5x, 'totalm2 = E18 .10,!, + 5x, theta = E18.10, 5x. 'curv2 = E18.10) 2004 formate' The trial vertical applied force is not'. + reach the tangency.', I. + Change the new fvO = 'E18.10.1x, + 'and continue !',I. + + + + *****************************************' '*******************************',1. All of the results written previous 'are NOT good.',I, PAGE 121 + + 2005 + + + + + + + + + 2006 + + + + + 2007 + + + + + 5501 + 5502 5503 5504 + + + + + + + + *****************************************' '*******************************') formate' The trial vertical applied force is not'. reach the tangency. '.1. Change the new FvO = 'E18.10.1x. 'and continue !'.I. *****************************************' '*******************************'.1. All of the results written previous 'are NOT good.'.I. *****************************************' ,*******************************,) format ('. The vertical distance has reached the 'top of pin and the data are:' .1. 106 5x. xb =' E18.10. 5x. 'yb = E18.10.l. 5x. fnb2 = E18.10. 5x. 'fra1 = 5x. fna1 =' E18.10. 5x. 'tota1m2 = 5x. theta =' E18.10. 5x. 'curv2 = format ( The vertical distance does not reach 'top of pin. yet. The data are: I. 5x. xb =' E18.10. 5x. 'yb 5x. fnb2 =' E18.10. 5x. 'frat 5x fna1 =' E18.10. 5x. 'totalm2 5x. theta =' E18.10. 5x. curv2 format ( 1. ----x ---------, p=[' F15.5. '.'. 3x. F15.5 format ( 3x. F15.5. '.', 3x. F15.5 y ) ) = = = = E18 .10.1. E18.10.1. E18.10) the E18 .10.1. E18 .10.1. E18.tO .I. E18.10) I format ( 3x. F15.5. '.'. 3x. F15.5 formate' The trial vertical applied force ];' ) is not'. reach the tangency. '.1. Change the new Sm1 = 'E18.10.1x. 'and continue !'.I. .*** **.*.* '**'.1. All of the results written previous 'are NOT good.'.I. *.******* *.**.' PAGE 122 + 5505 + + + + + + + + + 1099 '***************************,) formate' The trial vertical applied force is not', reach the tangency.' ,I, return end Change the new FvO = 'E18.10,1x, 'and continue !' ,I, *****************************************' '***************************' ,I, All of the results written previous 'are NOT good.' ,I, *****************************************' '***************************,) 107 PAGE 123 C C = END OF SUBROUTINE WALL-PIN C 108 C ================================================================ C C = DECK OF SUBROUTINE SEGMENT C THE PURPOSE OF THIS SUBROUTINE IS TO CALCULATE THE EQUILIBRIUM C CONDITION FOR SEGMENTS C C C INITIALIZE ALL OF THE STRESS AND MOMENT ARGUMENTS C ================================================================ subroutine segment include 'paramet.str' curv2 = curv1 + dcurve C ================================================================ C = INCREASE CURVATURE C ======================"========================================== 1111 sm2 = sm1 C ================================================================ C = CALCULATE SM C PURPOSES:' C (1). FIND Z(ETA) C (2). DETERMINE THE LENGTH OF EACH OF FIBER C (3). DETERMINE THE STRAIN OF EACH LAYER C (4). CALCULATE THE STRESS OF EACH LAYER C ================================================================ C = THE FIRST PURPOSE:" FINE Z(ETA) C ****** SET UP THE UPDATE THICKNESS h2 TO BE THE FUNCTION OF h1 C ******CALCULATE Z1 AND Z2 C ================================================================ C C ****** CALCULATE Z1 PAGE 124 C 1112 h2 = (hi sm1) / sm2 if (curv1 .eq. 0.0) then do 200 i = 1. neta + 1 z1(i) = O.Od+OO z1(i) = hi (i 1) deta -hi / 2.0d+00 200 continue else do 300 i = 1, neta + 1 zcom3a(i) = O.Od+OO zcom3a(i) = (1.0d+00 + (0 5d+00 curv1 hi) )**2.0d+00 zcom4a(i) = O.Od+OO zcom4a(i) = sqrt (zcom3a(i) + 2.'Od+OO hi + curv1 eta1(i z1(i) = O.Od+OO z1(i) = (-1 + zcom4a(i / curv1 300 continue endif C C ****** CALCULATE Z2 C if (cury2 .eq. 0.0) then do 400 i = 1. neta + 1 z2(i) = O.Od+OO z2(i) = h2 (i 1) .deta h2 / 2.0d+00 400 continue 500 else do 500 i=1. neta +1 zcom3(i) = O.Od+OO zcom3(i) = (1.0d+OO + (0.5d+OO curv2 h2**2.0d+00 + zcom4(i) = zcom4(i) = O.Od+OO sqrt(zcom3(i) + 2.0d+00 h2 curv2 z2(i) z2(i) continue = O.Od+OO = (-1 + zcom4(i I curv2 109 PAGE 125 110 endif C ================================================================ C = THE SECOND PURPOSE: CALCULATE THE LENGTH OF FIBER FOR EACH C FIBER C ****** FIND S1 AND S2 C do 600 i = 1, neta + 1 s1(i) = O.Od+OO s1(i) = sm1 (1.0d+00 + curv1 z1(i) ) s2(i) = O.Od+OO s2(i) = sm2 (1.0d+00+ curv2 z2(i) ) check(i) = O.Od+OO check(i) = s2(i) -s1(i) 600 continue C ================================================================ C = THE THIRD PURPOSE: DETERMINE STRAIN OF EACH FIBER C ****** FIND EPSILON2 C ================================================================ do 700 i = 1, neta + 1 C C ****** INITIALIZE THE STRAIN ARGUMENTS C epsilon2(i) = O.Od+OO C C ****** CALCULATE THE STRAINS EPSILON2 C C 111111111111111111111111 epsilon2(i) = log(abs(s2(i) / s1(i) sig2(i) = c ex epsilon2(i) C C ****** CHECK THE ELASTIC RANGE C PAGE 126 C if(sig2(i) .It. -yield) then sig2(i) =-c coef abs(epsilon2(i**en else if(sig2(i) .gt. yield) then sig2(i) = c coef abs(epsilon2(i**en end if C ****** CALCULATE THE NETUAL AXIS C 111 if check(i-1) .It. 0.0) .and. (check(i) .gt. 0.0 then znetural = z2(i 1) + deta abs(sig2(i -1) / + (sig2(i 1) -sig2(i) ) ) else if (check(i) .eq. 0.0) then sig2(i) = O.Od+OO znetural= z2(i) end if 700 continue sm2new -sm2 *(1 +curv2 znettiral) dsm2 = abs(sm2new -sm201d) C ====================================================================== C = CALCULATE FN C PURPOSES: C (1). CALCULATE FORCE FOR EACH LAYER C (2). SUM THE FORCES ALL OVER THE LAYERS C (3). CALCULATE THE MOMENT ABOUT THE NETURAL AXIS FOR EACH C LAYER C (4); SUM THE MOMENT ALL OVER THE LAYERS C (5). CHECK THE EQUILIBRIUM CONDITION C ================================================================ C = THE FIRST PURPOSE CALCULATE FORCE FOR EACH LAYER (dfn) C = THE SECOND PURPOSE SUM THE FORCES ALL OVER THE LAYERS (fn) C = THE THIRD PURPOSE CALCULATE THE MOMENT ABOUT PAGE 127 112 C THE NETURAL AXIS FOR EACH LAYER (dmn) C = THE FOURTH PURPOSE SUM THE MOMENT ALL OVER THE LAYERS C (totalm) C C C ****** INITIALIZE fn AND totalm C C fn = O.Od+OO totalm = O.Od+OO C ****** CALCULATE STRESS (dfn) FOR EACH LAYER C ****** CALCULATE MOMENT (dmn) FOR EACH LAYER C ****** SUM THE STRESSES ALL OVER TO BE TOTAL STRESSES (fn) C ****** SUM THE MOMENT ALL OVER TO BE TOTAL MMOMENT (totalm) C 1113 do 800 i = 1. neta dfn(i) = sig2(i) abs(z2(i) dmn(i) = dfn(i) z2(i) fn = fn + dfn(i) totalm = totalm + dmn(i) 800 continue -z2(i + 1 C ================================================================ C = THE FIFTH PURPOSE : CHECK EQUIVILIENT CONDITION C ****** THE EQUIVILIENT CONDITION IS: C (EQUATION (6) OF REFERENCE 1) C ****** DF(Ni) + Ki DMi = 0 C ****** i.e [FN(i)-FN(i-l)] + Ki [M(i) -M(i-l)] = 0 C ================================================================ C C ****** CALCULATE THE PARAMETERS C fnb2 = fn totalm2 = totalm PAGE 128 C ****** diffn = fnb2 -fna1 difmn = totalm2 -totalm1 C ****** equiold2= .equiold equiold = equi equi = diffn + curv2 difmn sm201d2 = sm201d sm201d = sm2 C C ****** CHECK THE EQUILIBRIUM CONDITION C if(abs(equi) .le. equitol) then c write(*,*)' fnb2 =' fnb2 cwrite(*.*)' fra1 ='. fra1 goto 1114 C ****** else if( equi .gt. 0.0) + .or. equi .le. 0.0) + then .and. .and. (equiold (equiold sm2 = 0.5d+00 (sm201d2 + sm201d) go to 1112 else .le. 0.0 .gt. 0.0) if ( (abs(equi -equiold) .le. 0.000000001) .or. + (equiold .eq. 0.0 then sm2 =sm2 1.01 go to 1112 else 113 PAGE 129 + endif C sm2 = sm2 -(equi) (sm201d sm201d2) / (equi -equiold) go to 1112 end if C = CALCULATE THE ANGLE THETA C 1114 C ratio1 = fnb2 / fra1 if (ratio1 .gt. 1.0d+OO) then curv2 = 1;2 curv2 go to 1111 else if(-ratio1 .gt. 1.0d+OO) then cUrv2 = O.8d+OO curv2 go to 1111 end if theta = phi1 -asin(fnb2 I fra1) C = IF THERE IS A REVERSE CURVATURE, C *** THETA = PHI1 + ASIN(FNB2 /FRA1) C c C = CHECK THE CONVERGENCE OF S(ttmi) =1 S(tmi) C sm2t = (theta -thetao) / curv2 sssold = sss sss = sm2t -sm2 curvold2 = curvold curvold = curv2 if(abs(sss) .le. sm2tol) then 114 PAGE 130 goto 1115 elseif(((sss .gt. 0.0) .and.(sssold .1e.0.0 .or. + ((sss .le. 0.0) .and. (sssold .gt. 0.0 )then curv2 = 0.5 (curvold2 + curvold) sml = srnO go to 1111 else if ( (abs(sss -sssold) .le. 0.00000001) .or. + (sssold .eq. 0.0) then + endif else curv2 = curv2 1.001 sml = smO go to 1111 curv2 = curv2 -(sss) (curvold -curvold2) / (sss -sssold) sml -= smO go to 1111 end if 1115 return end 115 PAGE 131 C C = END OF SUBROUTINE SEGMENT C C C = DECK OF THE SUBROUTINE INPUT C PURPOSE: C READ ALL OF THE INFORMATION FROM USERS AND WRITE THE C INPUT DATA OUT IN ORDER TO VERIFY THEY AR CORRECT! C subroutine input include 'paramet.str' 116 C ================================================================ C C = ARGUMENTS C C C C C C C C C C C C C C C title what geom The title of this input data file The purpose of this program (in this fileis "STRETCH" so far, in the futher may combine pure bending moment other capacities) case, or some Types of geometry: (1). Wall-pin (2). Wall-Two-pin (3). Wall-Three-pin (4). Wall-Two-pin (5). Wall-N-pin (6). Three Dimensional (only function aviable now) (in the further) (in the further) (in the further) (in. the further) (in the further) C ================================================================ C C = READ THE TITLE AND GEOMETRIC TYPE C read(40,1000)title read(40.2000)what read(40,2000)geometry if( what .ne. 'stre') then call input error (1) PAGE 132 else if( (geometry .ne. 'wall') .and. + (geometry .ne. 'WALL' then call input error (2) 'endif C C = WRITE OUT THE CHARACTERS TO FILE OUTPUT.STR IN ORDERR TO C VERIFY THOSE DATA RIGHT. C C ****** C ****** write(45.*)title ifwhat .eq. 'stre') .or. (what .eq. 'STRE') )then write(45.9001) else call inputerror(1) endif 117 if ( (geometry eq wall') or. (geomet ry eq Wall' then write(45.9002) C ****** else call inputerror(1) endif write (45.9003) C C = READ THE REST OF THE DATA C 100 read(40.2000)which C ****** if( (which :eq. 'init') read(40. *) smO. hOI write(45.4002) write(45.*)' write(45.*)' vrite(45.*)' write(45.*)' write(45 )' write(45.*)' go to 100 .or. (which .eq. IN!T' then curveO. coefmu. dcurve. fvO smO = smO hO = hO curveO = curveO coefmu = coefmu dcurve = dcurve fvO = fvO PAGE 133 c ****** C ****** C ****** C ****** C ****** C ****** + + + + + 118 else ifwhich .eq. 'mate').or.(which .eq. 'MATE'then read(40, *) coef, en, r, ex, yield write (45,4001) write(45,*)' coef = coef write(45,*)' en = en write(45,*)' r = r write(45,*)' ex = ex write(45,*)' yield = yield go to 100 else ifwhich .eq. 'geom').or.(which .eq. 'GEOM' then read(40, *) d1, rp1 write(45,4003) write(45,*)' d1 =' dt write(46,*)' rp1 =' rp1 go to 100 else ifwhich .eq. 'para').or.(which .eq. 'PARA' then read(40, *) neta write(45,4004) write(45,*)' neta =' neta go to 100 else .eq. 'tole').or.(which .eq. 'TOLE' then read(40, *) sm2tol, angletol, equitol lirite(46,4006) write(45,*)' write(46,*)' vrite(45,*)' go to 100 sm2tol ='. =', equitol =', sm2tol angletol equitol else ifwhich .eq. 'end ').or. (which .eq. 'END then go to 2200 else ifwhich ne. 'init').and.(which .ne 'INIT').and. (which ne. 'mate').and.(which .ne (which ne. 'geom') .and. (which .ne 'GEOM').and. (which ne. 'tole').and. (which .ne 'TOLE').and. (which ne. 'end ').and.(which .ne 'END ') .and. (which ne. 'para').and.(which .ne 'PARA'then PAGE 134 call inputerror(3) C ****** endi:f C C = WRITE OUT THE INPUT DATA TO FILE OUTPUT .STR IN ORDERR TO C VERIFY THOSE DATA RIGHT. C 2200 write(45,9003) C C = FORMAT C 1000 :format (a70) 2000 :format (a4) 4001 format (, ****** The material properties are: 4002 4003 4004 4005 9001 9002 9003 + foriDat (, :format (, :format (, :format (, :format (, :format (, :format'" ****** The initial values are: ****** The geometric properties ****** The parameters ar.e: ****** The tolerance are: This :file is doning stretching') and the geometry is Wall-Pin.') are: '===================================') return end ****** ****** ****** ****** ****** 119 ) ) ) ) ) PAGE 135 120 c C C = END OF SUBROUTINE INPUT C C ================================================================ c C = DECK OF SUBROUTINE INPUTERROR C THE PURPOSE OF THIS SUBROUTINE C PRINT OUT WARNING TO USERS. C 1000 + + + 2000 + + + 3000 + + + subroutine inputerror(n) if( n .eq. 1) then write(*.1000) else if( n.eq. 2) then write(*.2000) else if( n .eq. 3) then write(*.3000) endif format (, format ( format (, return end WARNING !'. I Your input data is WRONG on the fist line.' .1. Check that again!') WARNING !'.l. Your input data is WRONG on thesecond' l' I 1ne. Check that again!') WARNING !'. I Your input data is WRONG after the third' line.'. I Check that again!') PAGE 136 121 C ================================================================ C C = END OF SUBROUTINE INPUTERROR C C ================================================================ C C = DECK OF SUBROUTINE DATA C subroutine data include 'paramet.str' C ------CHARACTERS C control = 'O.K.' checksm = 'O.K.' C = TRANSFER ALL OF THE DATA FROM THE INITIAL DATA TO THE C CURRENT CONFIGURATION DATA C C ------PARAMETERS c spi meta deta jcount icount kconut = (1.0d+00 +en)/( (1.0d+00 + 2.0d+00 en)**0.5 ) = acos(O.Od+OO) = neta = 1.0d+00 / xneta = 0 = 0 = 0 C ****** CALCULATE ETA1 AND ETA2 do 100 i = 1, neta eta1(i) = (i 1) deta eta2(i) = (i 1) deta 100 continue PAGE 137 C ------LENGTH (DISTANCE) AND THICKNESS cypi = di + rpi hi = hO smi = smO C ------COORDINATES xa = O.Od+OO ya = cypi -0.5 hO yb = ya ybold = yb C ------FORCES AND MOMENTS fvai = fvO fnal = fvO coefmu fral = sqrt(fval**2.0d+00 + fnal**2.0d+00) totalml = O.Od+OO C ------ANGLES AND CURVATURES C phil = atan(coefmu) thetao = phil -asin(fnal/fral) theta curpl = thetao = 1.Od+00 / rpl curvl = curveO curvold. = O.Od+OO curvold2= O.Od+OO curvb= 1.0d+00 / (rpl + hO) return end C = END OF FILE C 122 PAGE 138 C C = END OF SUBROUTINE DATA C C C = DECK OF SUBROUTINE NEiDATA C subroutine newdata include 'paramet.str' C ------CHARACTERS control = '0. K. C C = TRANSFER ALL ot THE DATA FROM THE INITIAL.DATA TO THE C CURRENT CONFIGURATION DATA C C ------PARAMETERS jcount = 0 icount == 0 kcount = 0 C ******CALCULATE ETA1 AND do 100 i = 1, neta eta1 (i) = (i -1) deta eta2(i) = (i -1) deta 100 continue C ------LENGTH (DISTANCE) AND THICKNESS cyp1 = di + rp1 hi = hO sm1 = smO C ------COORDINATES 123 PAGE 139 xa = O.Od+OO ya = cyp1 -:' 0.5 hO yb = ya ybold = yb C FORCES AND MOMENTS fva1 = fvO fna1 = fvO co efmu fra1 = sqrt(fva1**2.0d+OO + fna1**2.0d+OO) totalm1 = O.Od+OO C ------ANGLES AND CURVATURES C phi 1 = atan(coefmu) thetao = phi1 -asin(fna1/fra1) theta = thetao cUrv1 curveO curvold == O.Od+OO curvold2= O.Od+OO curvb = 1.0d+OO / (rp1 + .O.5d+OO hO) return end C = END OF FILE C 124 PAGE 140 C C = END OF SUBROUTINE DATA C 125 C ==========================================?===================== C C = DECK OF SUBROUTINE COORD C PURPOSES: C (1). C (2). C ================================================================ + + subroutine coord include 'paramet.str' ybold = yb o yb -ybold + (1.0d+OO / curv2) ( cos (theta) xb return end = cos (thetao xa + (tota1m2 -(fra1 II< sin (phil) (ya -yb)/ofra1 0 / cos(phi1) PAGE 141 C C == DECK OF SUBROUTINE TRANS C C subroutine trans include 'paramet.str' C = TRANSFER ALL OF THE.DATA FROM THE INITIAL DATA TO THE C CURRENT CONFIGURATION DATA C C CHARACTERS control = 'N.G.' C PARAMETERS jcount = 0 sss -O.Od+OO sssold = O.Od+OO C ------ANGLES thetao = thetao + theta C ------MOMENTS totalm1 = -totalm2 C ------FORCES' C C ****** FRA1 IS NEVER CHANGED C ****** FNB2 IS KNOWN FROM SEGMENT.FOR C fva1 = sqrt(fra1**2.0d+OO fnb2**2.0d+OO) fna1 = fnb2 C CURVATURES 126 PAGE 142 c curvi = curv2 curvi. = curve O curvold =O.Od+OO C ------LENGTH AND THICKNESS smi = smO hi -hO C ------COORDINATES c C ybold, = yb return end C = END OF SUBROUTINE STRANS C 127 PAGE 143 C C = PARAMETERS AND DIMENSIONS C parameter (maxstep = 1000, max dim = 500) C C = CHARACTERS C character*70 title character*4 vhat,geometry, vhich, vork character*4 control, checksm C ------------------------------------------------------C = ARGUMENTS C ------------------------------------------------------C C = INITIAL VALUES C double precision smO. hO, curveO. dcurve. fvO doubie precision xa. ya C ------------------------------------------------------C = ARGUMENTS C C C C C C SMO HO CURVEO DCURVE -FVO : THE ORIGINAt LENGTH OF THE SEGMENNT THE ORIGINAL HEIGHT OF THE SEGMENT CURVATURE AT THE VALL THE INCREMENT OF CURVARURE THE APPLIED VERTICAL FORCE AT THE VALL (INITIAL GUESS FORCE) C ------------------------------------------------------C C = GEOMETRY C double precision cyp1. rp1, d1. xb. yb. ybold C ------------------------------------------------------128 PAGE 144 C = ARGUMENTS C CYP1 C C C THE Y-COORDINATE OF THE CENTER OF THE FIRST PIN (The original point is at the vall. Therefore, vatch out the positive or negative signs.) THE RADIUS OF THE FIRST PIN C RP1 C D1 C .THE VERTICAL DISTANCE BETWEEN THE WALL AND THE TOP OF THE FIRST PIN C C CURP1 C (i.e., D1 = CYP1 -RP1) CURVATURE OF THE FIRST PIN (i.e., CURP1 = 1.0/RP1) C ------------------------------------------------------C C = MATERIAL PROPERTIES C double precision coef, r, ex, yield, coefmu C ------------------------------------------------------C = ARGUMENTS C COEF k IN EQUATION (9a) OF REF. (1) C EN n IN, EQUATION (9a) ,OF REF. (1) C R : r IN' EQUATION (9c) OF REF. (1) C EX YONGJS MODULUS ABOUT X-AXIS C YIELD YIELD STRENGTH C COEFMU FRICTION COEFFICIENT C ----------------------------------------.--------:...-----C C C = PARAMETERS C integer neta" jcount. icount, kcount double precision s1(maxdim). s2(maxdim), check(maxdim) double precision equi, equiold C ------------------------------------------------------C = ARGUMENTS C ------------------------------------------------------129 PAGE 145 C C = NEW PARAMETERS C double precision eta1(maxdim) eta2(maxdim) double precision c, spi, deta double precision sm1, sm2, hi, h2 double precision sm2new, sm201d, sm201d2 double precision z1(maxdim) ,z2 (maxdim), zcom3a(maxdim), + zcom4a(maxdim), zcom3(maxdim), zcom4(maxdim) C ------------------------------------------------------C = ARGUMENTS C C C SPI C IN EQUATION (9c) OF REF. (1) PI c ______ C C = TOLERANCE C double precision sm2tol, angletol, equitol C---------------------------------------------________ C = ARGUMENTS C ----------------------------------------------_______ C C = ANGLES AND CURVATURES C double precision phi1, theta, thetao, thetab double precision dtheold, dtheta, theold double precision curv1, curv2, curvold. curvold2 double precision curp1 double precision fna1, fnb2, fra1. frb2 130 PAGE 146 c ------------------------------------------------------C = ARGUMENTS C ------------------------------------------------------C C = FORCES AND MOMENTS C double precision val, vb2, valold, din double precision n, vOold, vOold2 double precisiondn(maxdim) double precisiontotalml, totalm2, dimn double precision dmn(maxdim) c ------------------------------------------------------C = ARGUMENTS C TOTALMl C TOTALM2 C DIFMN THE RESULTANT MOMENT AT THE LEFT END THE RESULTANT MOMENT AT THE RIGHT END TOTALM2 -TOTALMl C ------------------------------------------------------C C = STRESSES AND STRAINS C double precision epsilon2(maxdim), sig2(maxdim) C ------------------------------------------------------C = ARGUMENTS C ___________________ 0 __________________________________ C ============================================================= C C = COMMOMS C C c 131 PAGE 147 C = COMMOM BLOCK OF CHARACTERS C common/ control1/ title, vhat, geometry, vhich, vork common/ control2/ control, checksm C C = COMMON BLOCK OF INITIAL VALUES C common/ common/ hO, curveO, dcurve, fvO init2/xa, ya C C = COMMON BLOCK OF GEOMETRY C common/ geom1/ cyp1, rp1, d1, xb, yh, ybold C C = COMMON BLOCK OF MATERIAL PROPERTIES' C common/ materal1/ coef" en, r, ex, yield, coefmu c C = COMMON BLOCK OF PARAMETERS C commonl common/ common/ param1/ neta, param2/ s1, param4/ equi, jcount, icoUnt, kcount s2, check equiold C C = COMMON BLOCK OF NEW PARAMETERS C -common/ nevpara1/ eta1, eta2 t deta common/ nevpara2/ c, spi common/ nevpara3/ sm1, sm2, hi, h2 common/ nevpara4/ sm2nev, sm201d, sm2old2 common/ nevpara6/ z1, z2, zcom3a, zcom4a, zcom3, zcom4 132 PAGE 148 C C = COMMON BLOCK OF TOLERANCE C common/ toler1/ sm2tol. angletol. equitol C C = COMMON BLOCK OF ANGLES AND CURVATURES C C commonl common/ commonl common/ angles1/ phi1. theta. thetao. thetab angles2/ dtheold. dtheta. theold curve1l curv1. curv2. curvold. curvold2 curve2/ curp1 C = COMMON BLOCK OF FORCES AND MOMENTS C commonl force"i/ fva1." fyb2; fva1old. diffn common/ force2/ fn. dfn. fvOold. common/ force3/ fnai. fnb2.fra1. frb2 common/ momenti/ totalm1." tota1m2. difmn common/. moment2/ dum C C = COMMON BLOCK OF STRESSES AND STRAINS C common/ C C = END OF FILE C ss1/ epsilon2. sig2 133 |