Citation
Software improvement to simulate pure bending and stretching of sheet metal

Material Information

Title:
Software improvement to simulate pure bending and stretching of sheet metal
Creator:
Wang, Yang-Cheng
Publication Date:
Language:
English
Physical Description:
xi, 133 leaves : illustrations ; 29 cm

Subjects

Subjects / Keywords:
Sheet-metal -- Computer programs ( lcsh )
Genre:
theses ( marcgt )
non-fiction ( marcgt )

Notes

General Note:
Submitted in partial fulfillment of the requirements for the degree, Master of Science, Mechanical Engineering.
General Note:
Department of 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:
31250644 ( OCLC )
ocm31250644

Downloads

This item has the following downloads:


Full Text
I
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
\ f. 9 ij


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
Cat
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 i 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


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 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-


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).
(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
a = ~ "2)2 + (02 *~)Jg [(£l e2)2 + (e2 ~~ e3)2 + (e3 ~ £l)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
[oi + 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
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].
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.24,25


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 action along the direction s.
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 As1(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, Smj:
Asi(si,z) = Asmi(l + kiz), -^ < z < ^ (3-1)
Where
As ^ = ^(01 0i-i) (3-2)
-ft-i
The normal force resultant Fn; and the resultant moment Mj are given by:
rh/2 h h
F Ni=/-.h/2rdZ 2***2 (3~3)
/h/2
M; = / J-h/2
Where the resultant moment M; is also a function of the curvature K;:
Mj = M (Kj) (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 (Fivi = 0). The relation
of the forces can be expressed by equation 3-6.
^Fi = FNi = 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 i will be kept constant under
any arbitrary loadings. For the steady state case after one step increment:
Asmj_|-ih.;_|-i = ASjnjhj 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 (rjn) of the element i which
did not deform after one step increment:
As [si,z(77n)] = As [smi+^Si+i (rjn)] (3 8a)
Asmi [1 -(- K;Zj (f/n)] = ASny-j-i [1 k;+iZ;+i (^/n)] (3 86)
Where rjn 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 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.
* = Kr (3-9)
Where for the present plane strain case, the stress on the fiber will be given by:
1 + 7 W = Ca
v/r+2f
(3 10)
where:
7 : normal anisotropy parameter.
n 1+^
n/1+2?
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 q)
v(C,v) dC
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 (s;, 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:
0 < < <
(3 126)
Where the direction of the coordinate axis has been reversed, and a new coordi-
nates set (e*, a*) is now defined with the origin at the beginning of the unloading:
a* = cr\ e* = e + ea (3 12c)
The total effective strain ej (fj,Zj) of the material fiber j is not taken cumula-
tively as expressed in equation 3-11, but according to the active piecewise ( 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 deter-
mination of each ( The procedure consists of monitoring the concave and convex surface fiber of the


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
Astniiare now calculated using equations 3-3, and 3-4.
8. equilibrium equation (equation 3-6) must be satisfied; otherwise, another
trial Astinj 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.


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; Fn; is the normal
forces resultant of the element i and Fv, 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 Frj 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 boundary conditions. Figure 3-4 represents the angle i at the connection be-
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 \ at the Connection Between the Sheet Metal and the
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:
AFN; + KjAMi = 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:
Xi = XA + --- [Mi FR1 sin fa (Ya Y1)\ (3 15)
1? R1 cos Yi = Y;_1 + (cos 0i cos 0i_x) (3 16)
Where the parameters Fr\ and Xa represent the inlet conditions at the point A
and the abscissa of A, respectively.


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 Fri is constant in
both magnitude and direction along the whole region. Because the resultant force
Fri 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 Ki
2. a trial value for the length of a finite element i after one step deformation
AstIni 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 torn equation 3-1, and the location of the neutral fiber is
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.


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
Astmiare now calculated using equations 3-3, and 3-4.
8. equilibrium (equation 3-13) must be satisfied; otherwise, another trial Astmi
is assumed and steps (3-8) are repeated. Reiterate through the thickness
until equilibrium is fulfilled.
9. the angle for the element i for the corresponding trial value Fri (equa-
tion 3-12) is determined; (see Figure 3-4).
10. the length at the middle surface segment i (equation 3-2) is determined:
Asttmi = (^i $il)
Ai
11. compare the trial Astmi for which equilibrium was fulfilled (step 8) with the
calculated length Asttmi (step 10). The trial curvature K; is now readjusted
and the steps 2 through 11 are repeated until Astmj converges to the value
Asttmi-
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


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: Y; is the coordinate
of the sheet metal at i element in Y-direction; CBY is the coordinate at
the center of the pin in Y-direction; 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.
III-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.


36
Figure 3-5 The Flow Chart of Stretching Program with Wall-One-Pin Geometry


CHAPTER IV
NUMERICAL RESULTS AND DISCUSSIONS
IV-1 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-1 lists the comparison of
the CPU time for different numbers of the composed fiber.


38
Table 4-1: The Comparison of the CPU Time
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


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:
K\ : 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 (= quItc)
(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:
B =
0 n = 1
^[o.ISTSxCq ]
(4-3)
n1
n 1
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


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 users manuals.
Table 4-2: Material Properties and Initial Values
Material Initial
Properties Values
K (psi) 7833.4 Smo (in) 0.25
n 0.228 H0 (in) 0.2
7 0.177 Co 0.0
Ex (psi) 20E6 Deo 1.0E-6
Yield (psi) 20E3 N curve 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


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.


42
i
x>
s
o
3
-0.5 0 0.5 1 1.5 2 2.5 3 3.5 4
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,


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.
0 0.5 1 1.5 2 2.5
The Location Along the Horizontal Direction (in)
Figure 4-2 The Profile of the Sheet Metal.


44
The region between the two ends is unsupported area. The resultant force, Ffu,
is a constant in the unsupported region but the normal resultant force, Fn; 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 wall, 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.
The Normal Axial Force vs. The Horizontal Location of Sheet Metal
0 0.5 1 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


45
0.01 0.02 0.03 0.04
Curvature (1/in)
0.05
0.06
0.07
Figure 4-4 The Moment-Curvature Curve for Stretching Sheet Metal with
Wall-One-Pin Geometry
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 moment-
curvature 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


46
curve becomes nonlinear.
IV-4 Discussion
According to tbe 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.


CHAPTER V
CONCLUSIONS AND RECOMMENDATIONS
V-l 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.


48
3. No matrix manipulation is necessary.
In present method, the number of fibers along the thickness may have one hundred
but in finite 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 have 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 assumed 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.


APPENDIX A
NOTATIONS


50
NOTATIONS
Chapter I
e = strain
K = strength coefficient
n = strain hardening effect on the sheet due to the forming process
e; = 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


51
Chapter III
Ki Curvature of element i
Smi = Length of element i at mid-surface
FNi = Normal force resultant
Mi resultant moment
Fj the normal force acting on fiber j.
Vn The neutral axis corresponds to the location of the fiber of the element i
7 normal anisotropy parameter.
C 1+7 7i+2?
(.o') = the origin at the beginning of the unloading
*j(flj>zj) = The total effective strain of the material fiber j
xA the X coordinate of Fri at the point of contact with the wall
4>i the angle formed by Frj with the vertical


APPENDIX B
THE USERS MANUAL OF PURE BENDING PROGRAM


53
THE USERS MANUAL OF PURE BENDING PROGRAM
For this case, the input 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 parameters 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, ).
TITLE
TITLE : Title of this file (arbitrary characters up to 72 )
Example:
This is a example input data file of pure bending


CONTROL
CONTROL : There are six cases
PURE = pure bending moment (Only available now)
Example:
pure bending moment


55
MATE
COEF EN R EX YIELD
MATE: material properties (only the first four characters will be readied)
COEF : the coefficient, k of the equation a = k |e|en
EN : the power of the above equation
R : the normal anisotropy parameter
EX : Youngs Modulus
YIELD : yielding stress
Example:
material properties
506, 0.228, 0.177, 30E6, 36E3


INIT

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
D CURVO : 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


57
PARA
NETAO FNTOL SM2TOL NLOOP
PARA : parameters
NETAO : the total number of eta (77)
FNTOL : the tolerance of the summation of the stresses
SM2T0L : 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


58
END OF FILE
END OF FILE : end of file
Example:
end of file


APPENDIX C
PURE BENDING INPUT DATA FILE


60
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
initial values
0.40,0.1,0,0.0001, 30
parameters
200.0. 8, 0.02, 3
end of file


APPENDIX D
COMPUTER PROGRAM OF PURE BENDING


COMPUTER PROGRAM OF PURE BENDING
program main
C ===========================================================
C
C = THIS PROGRAM IS TO DETERMINE:
C (1). THE MATERIAL PROPERTIES RELATIONSHIP BETWEEN PURE
C BENDING MOMENT AND CURVATURE (PUREM)
C (2). 3333
C CREAT DATE : OCT. 16, 1993
C MODIFIED DATE
C AUTHOR : YANG-CHENG WANG
C
C ===========================================================
include 'control.1
include 'para.l'
C
C = OPEN AN OUTPUT DATA FILE WHOSE OUTPUT DATA CAN BE USED IN
C 3 MATLAB AND OPEN AN INPUT DATA FILE CALLED INPUT.DAT
C
open(unit=40.files' input.dat,status3'old')
C
open(unit=35,file='out.m' status3new')
open (unit=45,file=mk.dat,status3new')
open(unit=55,file3shistory.dat',status3'new)
C
C 3 CALL SUBROUTINE INPUT
C THE PURPOUSE:
C READ AND CALCULATE THE ORIGGINAL PARRAMETERS
C
call input
C
if((what .eq. 'pure') .or. (what .eq. 'PURE')) then


63
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 ((what .ne. pure) .or. (what .ne. 'none)
+ or.
+ (what .ne. 'PURE') .or. (what .ne. 'NONE'))
+ then
call inputerr
stop
end if
stop
end


64
C
C = DECK DF SUBROUTINE INPUT
C
subroutine input
include control.1*
c include peura.l
c = MPRO MATERIAL PROPERTIES
c = INVA INITIAL VALUES
c = PARA PARAMETERS
c = c c END END OF FILE
read(40,999)title
read(40,1000)what
write(45,999)title
if((what .eq. pure) .or. (what .eq. PURE)) then
write(45,5000)
end if
write(45,6000)
100 read(40,1000)which
if ((which .eq. mate)
read(40,*)coef, en, r,
write(45,1001)which
.or. (which
ex, yield
.eq. 'MATE1))then
write(45,*)' coef =, coef
write(45,*) en = en
write(45,*) r = '. r
write(45,*)' ex = \ ex
write(45,*) yield =>, yield
writ e(45,*)*======================================
go to 100
else if((which .eq. 1init)or.(which .eq. ,INIT,))then
read(40,*) smO, hO, curvO, dcurvO, ncurv
write(45,1001)which
write(45,*)' smO =*, smO
write(45,*)' hO =*, hO
write(45,*) curvO =, curvO
write(45,*) dcurvO =, dcurvO
write(45,*) ncurv =, ncurv


65
999
1000
1001
5000
2000
go to 100
else if((which .eq. para).or.(which .eq. PARA1))then
read(40,*) netaO, fntol, sm2tol, npathe
writ e(45,1001)which
write(45,*)' netaO =, netaO
write(45,*) fntol =, fntol
write(45,*)' sm2tol =', sm2tol
write(45,*)' npathe =', npathe
write(45,*)======================================'
go to 100
else if((which go to 2000 .eq. end ) .or.(which .eq. 'END ))then
else if((which .ne. 'mate') . and. (which .ne. init)
+ .and. (which .ne. 'MATE') .and. (which .ne. 'INIT')
+ . and. (which .ne. end ) .and. (which .ne. 'para')
+ + .and. (which ) then .ne. 'END ) .and. (which .ne. 'PARA')
. call inputerror
end if
format(a70)
format(a4)
format(lx,a4)
format(J
+
+
+
6000 format(
+
+
+
+
This is required to fine the relationship'
between the pure bending moment and'
curvature',//)
-----=====-----',//,
Verify the input data',//,

return
end


66
C
C = DECK OF SUBROUTINE INPUTERR
C
subroutine inputerr
print *,
print *,
write(45,1000)
1000 format(
+
return
end
Your option is NOT included.
Check that again !!!!!
Your option is NOT included.,/,
Check that againn!!!!!)


67
C
C = DECK OF SUBROUTINE PUREM
C
subroutine purem
C CREATED DATE : SEP. 20, 1993
C MODIFIED DATE : NOV. 7, 1993
C COMPUTER : PC VERSION
C AUTHOR : YANG-CHENG WANG
C =============================================================
C THIS PROGRAM IS DETERMINING THE RELATIONSHIP BETWEEN THE
C MOMMNT AND CURVATURE
C =============================================================
C =============================================================
C
C = ALL OF THE COMMON AND DIMENSION EXTERNAL FILE
C
C =============================================================
include peura.l'
C =============================================================
C
C CALCULATE THE STREE-STRAIN RELATIONSHIP AND THEIR
C COEFFICIENTS
C
C =============================================================
c ****** C = PARAMETER FOR THE RELATION OF STRESS AND
C STRAIN
C r = THE MATERIAL PROPERTIES
C =============================================================
C CALCULATE THE COEFFICIENT C
C =============================================================
c = (l.Od+OO + r) / (Cl.Od+OO + 2.0 d+00 r)**0.5)
C =============================================================


c
C CHANGE ALL OF THE INITIAL VALUES TO BE THE FIRST STEP DATA
C
c ****** JCOUNT = COUNT FOR THE NUMBER OF INTERATION
c ****** CURV1 = CURRENT CONFIGURATION OF CURVATURE
c ****** CURV2 = UPDATE CONFIGURATION OF CURVATURE
c ****** SMO = LENGTH OF THE NETURAL LINE OF ORIGINAL
c CONGIGURATION
c ****** SMI = LENGTH OF THE NETURAL LINE OF CURRENT
c CONGIGURATION
c ****** SM2 = LENGTH OF THE NETURAL LINE OF UPDATE
c CONGIGURATION
c ****** HO = THICKNESS OF THE OF ORIGINAL CONGIGURATION
c ****** HI = THICKNESS OF THE OF CURRENT CONGIGURATION
c ****** H2 = THICKNESS OF THE OF UPDATE CONGIGURATION
c ****** NETAO ' = NUMBER OF ETA OF ORIGINAL CONGIGURATION
c ****** NETA = NUMBER OF ETA OF CURRENT CONGIGURATION
C TRANNSFER THE INITIAL DADAS TO THE CUNNRENT CONFIGURATION
j count = 0
jcount1 = 1
curvl = curvO
sml = smO
nsm2 = 0
csm2 = sm2tol
hi = hO
neta = netaO
xneta = neta
deta = 1/xneta
C =============================================================
C CALCULATE ETA
do 100 i = 1, neta+1
etal(i) = (i-l)*deta
eta2(i) = (i-l)*deta


69
100 continue
C =============================================================
C INITIALIZE ALL OF THE STRESS AND MOMENT ARGUMENTS
C =============================================================
do 101 i = 1, neta
fdfn(i) = O.Od+OO
fdmn(i) = O.Od+OO
epsilon(i) = O.Od+OO
101 continue
fn = O.Od+OO
ffn = O.Od+OO
ftotalm = O.Od+OO
C =============================================================
C CALCULATE THE COEFFICIENT FOR THE CURVARURE
C =============================================================
C CURV2 = A * EXP( B*NCURV)
C 0 NCURV = 0, ===> CURV2 = DCURVO
C 0 NCURV = NCURV, ===> CURV2 = 1/0.1875
b = log(l/(0.1875*dcurv0))/(ncurv-l)
do 102 i = 1, neta+1
zO(i) = hi (i-1) deta hl/2.0
102 continue
C
C ****** CALL SUBROUTINE POSITIVE TO CALCULATE THE
C RELATIONSHIP
C BETWEEN MOMENT AND CUVATURE FOR POSITIVE CUVATURE
C
call positive
return
end


70
C ==
C
C =
C
c
c
c
c ==
DECK OF SUBROUTINE POSITIVE
THE PURPOSE:
CALCULATE THE RELATION BETWEEN THE CURVATURE AND MOMENT
FOR POSITIVE INCREMENT OF CURVATURE
subroutine positive
include para.l
C =============================================================
C
C = INCREASE CURVATURE
C
1111 jcount = jcount + 1
c curv2 = dcurvO + (1.0/1.875 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.
+ (j count .eq. (ncurv - 8)) .or.
+ (jcount .eq. (ncurv - 10)) .or.
+ (jcount .eq. (ncurv - 12)) .or.
+ (jcount .eq. (ncurv - 14)) .or.
+ (jcount .eq. (ncurv go to 1111 - 16)) ) then
else
curv2 = dcurvO*exp((jcount-i)*b)
end if
else if(ncurv .ge. 21) then
if ( (jcount .eq. (ncurv - 2)) .or.
+ (jcount .eq. (ncurv - 4)) .or.
+ (jcount .eq. (ncurv go to 1111 - 6)) ) then
else
curv2 = dcurvO*exp((j count-1)*b)
end if


else
curv2 = dcurvO*exp((jcount-l)*b)
end if
sm2 = sml
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
c ******
c
c ******
c ******
c ******
c ******
c ******
c ======
C = SET UP THE UPDAT THICKNESS h2 TO BE hi
C = AND CALCULATE Z1 AND Z2
C ===========================================================
2222 h2 = (hl*sml)/sm2
if (curvl .eq. O.O) then
do 200 i = 1, neta+1
zl(i) = hi* (i-1) deta hl/2
200 continue
else
do 300 i = 1, neta+1
zcom3a(i) = (i.0d+00-(0.50d+00*curvl*hl))**2
zcom4a(i) = sqrt(zcom3a(i) + 2.0d+00
+ * hi curvl etal(i))
zl(i) = (-1 + zcom4a(i)) / curvl
Z1 = THE THICKNESS AT EACH FIBER OF CURRENT
CONFIGURATION
Z2 = THE THICKNESS AT EACH FIBER OF UPDATE
CONFIGURATION
DEAT = INCREMENT OF ETA
ZC0M3A= PARAMETER FOR CALCULATION Z1
ZC0M4A= PARAMETER FOR CALCULATION Z1
ZC0M3 = PARAMETER FOR CALCULATION Z2
ZC0M4 = PARAMETER FOR CALCULATION Z2


300 continue
end if
C
C
C
if(curv2 .eq. 0.0) then
do 400 i = 1, neta +1
z2(i) = h2* (i1) 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 ******
c
c ******
c
c ******
c
c ******
c ******
c ******
c
c ******
c ******
c ******
c =======
C = CALCULATE THE LENGTH OF EACH FIBER
do 600 i = 1, neta+1
sl(i) = sml*(l+curvl*zl(i))
s2(i) = sm2*(l+curv2*z2(i))
check(i)= s2(i) sl(i)
51 = THE LENGTH OF EACH FIBER AT DIFFERENT LAYER
OF CURRENT CONFIGURATION
52 = THE LENGTH OF EACH FIBER AT DIFFERENT LAYER
OF UPDATE CONFIGURATION
EPSIL0N2 = STRAIN OF EACH LAYER OF UPDATE
CONFIGURATION
SG2 = STRESS OF UPDATE CONFIGURATION
YIELD = YIELD STRESS OF THE MATERIAL
ZNETURAL = THE LOCATION OF NETURAL AXIS
(AT THE FIBER OF STRESS EQUAL TO ZERO)
SM2NEV = PARAMETER FOR NEW SM2
DSM2 = DIFFERENCE BETWEEN NEW SM2 AND SM2
SM2T0L = TOLERANCE OF SM2


73
600 continue
do 700 i = 1, neta+1
C =================================
C = INITIALIZE THE STRAIN ARGUMENTS
C =================================
epsilon2(i) = 0.0
C =============================================================
C = CALCULATE THE STAINS EPSIL0N2
C =============================================================
epsilon2(i) = log(abs(s2(i)/sl(i)))
sig2(i) = c ex epsilon2(i)
C = CHECK THE ELASTICC RANGE
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 NEUTRAL AXIS
+
if(check(i-l) .It. 0.0 .and. check(i) .gt. 0.0)then
znetural = z2(i-l) + deta*abs(sig2(i-l)/
(sig2(i-l)-sig2(i) ) )
sm2old = sm2
sm2new = sm2*(l + curv2 znetural)
dsm2 = abs(sm2new sm2)
if (dsm2 .le. sm2tol) then
sm2 = sm2nev
else
sm2 = sm2nen


74
C = IF IT IS NOT CONVERGED THEN GOTO 2222 TO MODIFIED h2
C =============================================================

700
C ==
C
C
C
C
C
C
C
C
C
go to 2222
end if
else if (check(i) .eq. 0.0) then
sig2(i) '= 0.0
znetural = z2(i)
end if
continue
CALCULATE FM
THE PURPOSES:
(1) . DETERMINE THE FORCE FOR EACH LAYER
(2) . CALCULATE THE SUMMATION FORCE FOR THE STRUCTURE
(3) . DETERMINE THE MOMMENT ABOUT THE NETURAL FOR EACH LAYER
(4) . CALCULATE THE MOMENT ABOUT THE NETURAL AXIS FOR THE
STRUTURE
C
C
C
C
C
C
C
C
C
******
******
******
******
DFN
= THE STRESS OF EACH LAYER FOR UPDATE
CONFIGURATION
= THE TOTAL STRESS
= THE MOMENT OF EACH LAYER ABOUT ZNETURAL
FOR UPDATE CONFIGURATION
TOTLAM = THE TOTAL MOMENT ABOUT THE ZNETURAL
FN
DMN
THE EQUILIBRIUM EQUATION IS SUM(DFN) = FN = 0.0
C =============================================================
C = INITIALIZE THE fn AND totalm
C =============================================================
fnold = fn
fn = 0.0
totalm =0.0


C = CALCULATE STRESS(dfn) AND MOMENT(dinn) FOR EACH FIBER
C = AND TOTAL STRESS(fn) AND TOTAL MOMENT (totaim)
3333 do 800 i = 1, neta
dfn(i) = sig2(i) abs(z2(i) z2(i+l))
dmn(i) = dfn(i) * z2(i)
fn = fn + dfn(i)
totaim = totaim + dmn(i)
800 continue
C =============================================================
C
C = CKECK THE EQUILIBRIUM EQUATION SUM(FN)=?0
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. ((1.0/0.1875)-0.00001) .or.
+ jcount .ge. (ncurv 2)) then
work = fini*
end if
call writeout
C =============================================================
C = UPDATE ALL F THE PARAMETERS
C =============================================================
sm2 = smO
if (work eq. 'fini) then
C = CHECK THE LOADING AND UNLOADING CONDITIONS
C = NPATHE : NUMBER OF THE PATHES
C = NPATHE = 1 OR 0 (OR INDICATING NOTHING) THEN
C ONLY THE FIRST LOADIN PATH (POSITIVE


76
C MOMENT)
C = NPATHE = 2 POSITIVE MOMENT + 1 UNLOADING
C (NEGATIVE MOMENT)
C = NPATHE = 3 POSITIVE LOADING + NEGATIVE MOMENT +
C POSITIVE LOADING
if(npathe .eq. 2) then
call pathe2
else if(npathe .eq. 3) then
call pathe2
call pathe3
else if((npathe.eq.O).or.(npathe.eq.l))then
stop
else
stop
end if
else
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)
tol = fn/(fnold-fn)
else if(fn ,gt. fntol .and. fnold .gt. fntol) then
if (abs(tol) .le. 30.0) then
sm2 = sm2 tol (sm2old-sm2new)
else
sm2 = sm2 dsm2/2.0
end if
go to 2222
else if(-fn .gt. fntol .and. -fnold .gt. fntol) then
if (abs(tol) .le. 30.0) then
sm2 = sm2 tol (sm2old-sm2ne)


77
else
sm2 = sm2 + dsm2/2.0
end if
go to 2222
else if((fn .gt. fntol .and. -fnold .gt. fntol) .or.
+ (-fn.gt. fntol .and. fnold .gt. fntol))then
sm2 = (sm2new+sm2old)/2,0
go to 2222
else
go to 2222
end if
return
end


78
C =======================
C
C = SUBROUTINE INPUTERROR
C
C =======================:
subroutine inputerror
write(45,100)
write(*, 100)
100 format( Fatal Error ??????
+ //
+ * Your input data formation is not right
+ !!!!!,//,
+ Check that again, I am sorry!!!)
return
end


79
C
C = DECK OF SUBROUTINE WRITEOUT
C
subroutine writeout
include para.l
write(45,100)
write(45,200)totalm
write(45,300)curv2
write(45,400)znetural
write(45,500)sm2
100 format(1 ====================================== >)
101 format(v ======================================,
200 format(* The Moment is el4.7)
300 format( The curvature is ,2x, el4.7)
400 format( The netural axis is at ,2Xi el4.7)
500 format(' The length of netural axis is :,2x, el4.7)
C =============================================================
C HERE IS STRAIN HISTORY
if (jcount .eq. 1) then
write(55,600)
write(55,601)jcount
else
write(55,601)jcount
end if
write(55,101)
write(55,602)
do 1000 i 1, neta
write(55,603) i, dfn(i), dmn(i), s2(i)
1000 continue
write(55,604) fn, totaim
600 format ( This is the output data file for strain1
+ history.1,/)


601
format(l, This is the,lx,i4,lx,step iteration
+ of strain history.)
602 format( # Fiber ,1, Fn
+ M S2
+ /.
603 format(4x,i4, 2(lx,|, 2x, el7.9),lx,'|,2x,e20.12)
604 format( Total Fn = el4.7,
+ Total Moment = ,el4.7,///////)
C HERE IS CREATING THE DATA FILE FOR MATLAB TO PLOT H-K PLOT
if (jcount1 .eq. 1) then
write(35,703)
write(35,704)totalm, curv2
else if (work .ne. fini) then
write(35,705)totalm, curv2
else
write(35,706)totalm, curv2
end if
703 format(/, Total Moment ,
+ curvature ,/,
+ =====-------------------------------------=,
+ ==========================,//)
704 format(*mkl=[ 2(2x, e20.12))
705 format(7x, e20.12, 2x, e20.12)
706 format(7x, e20.12, 2x, e20.12,];)
C = SAVE THE YIELDING MOMENT(YMOMENT) AND YIELD CURVATURE
C (YCURVE)
if((abs( sig2(l) ) .ge. yield ) .or.
+ (abs(sig2(neta + 1) ) .ge. yield ) ) then
goto 212
else
ymoment = totalm
ycurve = curv2
nyield = jcount1
end if
212


81
c ===========================================
C = SAVE THE INFORMATION OF THE FIRST LOADING
C ===========================================
savem(jcountl) =
savec(jcountl) =
savez(jcount1) =
saves(jcountl) =
jmax = jcountl
jcountl = jcountl
return
end
totaim
curv2
znetural
sm2
+ 1


82

C
C
C
C
C
= DECK OF SUBROUTINE OF PATHE2
100
+
+
+
+
200
subroutine pathe2
include para.1
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)
continue
kkk = 1 + nyield
jjj = jmax + nyield
do 200 i = kkk, jjj
, p2curv2(i) = savec(jmax)
- savec(nyield)
p2totalm(i) = savem(jmax)
- savem(nyield)
p2znetu(i) = savez(jmax)
- savez(nyield)
p2sm2(i) = saves(jmax)
- saves(nyield)
continue
call patheout
return
end
savec(i-nyield)
savem(i-nyield)
savez(i-nyield)
saves(i-nyield)


83
C ========================================
C
C = DECK OF SUBROUTINE OF PATHE3
C
C ========================================
subroutine pathe3
include para.l
double precision p2curv2(maxstep),
double precision p2znetu(maxstep),
double precision p3curv2(maxstep),
double precision p3znetu(maxstep),
C
C
common/ path2/ p2curv2, p2totalm,
common/ path3/ p3curv2, p3totalm,
jjj = jmax + nyield
do 100 i = 1, nyield
p3curv2(i) = p2curv2(jjj) +
p3totalm(i) = p2totalm(jjj) +
p3znetu(i) = p2znetu(jjj) +
p3sm2(i) = p2sm2(jjj) +
100 continue
kkk = 1 + nyield
jjj = jmax + nyield
do 200 i = kkk, jjj
p3curv2(i) = p2curv2(jjj) +
+ + savec(nyield)
p3totalm(i) = p2totalm(jjj) +
+ + savem(nyield)
p3znetu(i) = p2znetu(jjj) +
+ + savez(nyield)
p3sm2(i) p2sm2(jjj) +
+ + saves(nyield)
continue
call patheoiit
return
end
p2totalm(maxstep)
p2sm2(maxstep)
p3t ot aim (max st ep)
p3sm2 (maxstep)
p2znetu, p2sm2
p3znetu, p3sm2
savec(i)
savem(i)
savez(i)
saves(i)
savec(i-nyield)
savem(i-nyield)
savez(i-nyield)
saves(i-nyield)


84
C
C = DECK QF SUBROUTINE PATHEOUT
C
subroutine patheout
include para.1
double precision p2curv2(maxstep), p2totalm(maxstep)
double precision p2znetu(maxstep), p2sm2(maxstep)
double precision p3curv2(maxstep), p3totalm(maxstep)
double precision p3znetu(maxstep), p3sm2(maxstep)
C
C
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 100 i = 1, jjj
if (i .eq. 1) then
write(35,703)
write(35,704)savem(jmax), savec(jmax)
write(35,705)p2totalm(i), p2curv2(i)
else if (i .ne. jmax+nyield) then
write(35,705)p2totalm(i), p2curv2(i)
else
write(35,706)p2totalm(i), p2curv2(i)
end if
100 continue
end if
if(npathe .eq. 3) then
write(35,*) *7, Should be continued!
do 200 i = 1, jjj
if (i .eq. 1) then
write(35,703)
write(35,7044)p2totalm(jjj),p2curv2(jjj)
write(35,705)p3totalm(i), p3curv2(i)


85
else if (i .ne. jjj) then
write(35,705)p3totalm(i), p3curv2(i)
else
write(35,706)p3totalm(i), p3curv2(i)
end if
200 continue
end if
if ((npathe .ne. 3) .and. (npathe .ne. 2)) then
write(35,*) Your path number is not right.
write(35,*) Check that again !
end if
c *************************************************************
703 format(/, Total Moment
+ curvature ,/,
+ -------------------------------==----===,
+ ,----------------------------======---==,//)
704 format(* 1^2=1 2(2x, e20.12))
7044 format(mk3=[ 2(2x, e20.12))
705 format(7x, e20.12, 2x, e20.12)
706 format(7x, e20.12, 2x, e20.12,];)
return
end


CONTROL.L
character*70 title
character*4 what, which
common/ control/title, what, which


87
C
C = PARAMETERS AND DIMENSION
C
parameter (maxstep = 1000, maxdim = 500)
character*4 work
double precision curvO, curvl, curv2, dcurvO, b, c
double precision coef, en, r, ex, yield
double precision smO, sml, sm2, dsm2, sm2new, sm2tol
double precision netaO, fntol, hO, hi, h.2
double precision fn, totaim, 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),
+ sigl (maxdim), sig2(maxdim),
+ zl(maxdim), z2(maxdim), zO(maxdim)
double precision etal(maxdim), eta2 (maxdim)
integer jcount, jcountl, jmax, nyield
integer neta, nerr, nsm2, ncurv, npathe
c c = f% C0MMM0NS
O common/ charac/ work
common/ mp/ coef, en, r, ex, yield
common/ ivl/ curvO, curvl, curv2, ncurv,
+ dcurvO,b,c
common/ iv2/ smO, sml, sm2, dsm2, sm2new,
+ sm2tol
common/ pal/ netaO, fntol, neta, nerr, jcount
common/ pa2/ fn, totalm, eta, deta, csm2,
+ sm2, ffn, ftotalm
common/ pa3/ znetural, hO, hi, h2, jcountl
common/ pa4/ zcom3, zcom4,zcom3a, zcom4a
common/ pa5/ si, s2
common/ exl/ epsilon2, check, epsilon


common/
common/
common/
common/
forcel/ dfn, dmn, sigl, sig2, zl, z2,
zO, fdfn, fdmn
twoeta/ etal, eta2
save/ savem, savec, savez, saves
paths/ jmax, npathe, ymoment, ycurve,
nyield


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

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 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.

PAGE 4

IV The computer programs are written based on a 486 IBM comparable computer, MAGTRON, with 55 MHz, 16 RAM memory and 1.4GB 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 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. 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 -1 Prologue . . . . . . . . . . . . . . 1 I-2 Historical Review . . . . . . . . . . . . 4 I-2-1 Influence of the Mechanical Properties of the Material on Sheet forming . . . . . . . . . . . 4 I-2-2 Influence of the Forming Geometry .. .. .. .. .. .. .. .. .. .. .. 7 I-2-3 Limits on the Formability of Sheets .......................... 8 I-2-4 Contact and Friction Effects ................................. 9 I-3 Problem Statement .............................................. 11 II. DESCRIPTION OF THE PROBLEM .. .. .. .. .. .. .. .. .. .. .. 13 II-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 Il-4-2 Theory of Beams and Shells ............................... 18 II-4-3 Numerical Methods ........................................ 19 Il-4-4 Finite Element Method .................................... 19 III. MODEL DEVELOPMENT . . . . . . . . . . 20 111-1 lnt:roduction ................................................... 20 111-2 Pure Bending .................................................. 21 111-2-1 Mathematical Formulation ................................ 22 111-2-2 Computational Algorithm .... : ............................ 26 111-2-3 Computational Implementation .................. :. . 27 III-3 Stretching for Wall-Pin Geometry . . . . . . . 29 111-3-1 Mathematical Formulation ................................ 31 III-3-2 Computational Algorithm ................................. 33 111-3-3 Computational 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 . . . . . . . . . . . . . 4 7 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

l.X 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 1 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

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 .............................................. 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

I-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. Eary1 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:re.4 pp.I-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) 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 2

PAGE 15

thickness in the fiat 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 a.ild 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.4GB 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). (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 Y'' 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 = u (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 11 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 and Framer .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 geometry has received more analytical attention in recent years. A number of articles10 1315 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. 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 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 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 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 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

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 Singh19 studied the springback effect in spherical caps formed from circular blanks without blank holder. They found that spring back 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 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 p, 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 .1. 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 strainrate 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."

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 DESCRIPTION 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 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 a.t the point of contact a.s .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 tlie 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 phis, the strip is

PAGE 26

14 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

B* .. .. .. Figure 2-3 Strip Deforming of Transient Case -----.;;;: i i+l Figure 2-4 Strip Deforming of Steady State Case 15 \ \ \ I I I I I

PAGE 28

16 11-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 long1tudinal :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 forining (plane strain deformation). 11-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 fiat 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 usually are based on the perfect elastic-plastic, as well as linear hardening 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. 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 a.n important a.dva.nta.ge in recent years. Some resea.rch22 26 proposed numerical procedures to do step-by-step incremental plasticity a.pproa.ch subjected to certain restric tions such a.s the boundary conditions a.nd 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 a.pproa.ch to sheet forming problems is based on the va.ria.tiona.l methods. It is the basic function of the finite element method. To analyze sheet forming, several types of nonlinea.rities in finite element method ha.ve to be involved such a.s geometric nonlinearity, ma.teria.l nonlinearity, a.s well a.s contact problems.

PAGE 32

III-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 111-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 hh and length Llsmi can be modeled by an arc of a circle of curvature ki 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 u(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 111-2-1 Mathematical Formulation Assuming cross section of the typical element to remain normal during the deformation, the deformed arc length .6.s1(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 .6.s = (B 8 1) rru Ki 1 1-(3-1) (3 -2) The normal force resultant FNi and the resultant moment Mi are given by: lh/2 FNi = crdz -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 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 (71n) of the element "i" which did not deform after one step increment: (3Sa) (3Sb) Where '1/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. u = Kf'l (3 -9) Where for the present plane strain case, the stress on the fiber will be given by: where: u= 1+1 u= Gu Jl+21 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: r 1 dv((,.,) f ( s' ., ) = J tXJ c v (''.,) d( d( (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 (sh 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 < < f* -J 0 (3-12a) (3-12b) Where the direction of the coordinate axis has been reversed, and a new coordinates set ( E"', o-*) is now defined with the origin at the beginning of the unloading: -0" = -O"j e* = -e+ fo (3-12c) The total effective strain fj (shzj) of the material :fiber "j" is not taken cumulatively as expressed in equation 3-11, but according to the active piecewise (0'*, "') 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 ( o-"', 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 mo ment 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 K1 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 As 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 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 now calculated using equations 3-3, and 3-4. 8. equilibrium equation (equation 3-6) must be satisfied; otherwise, another trial is assumed and steps 2 through 8 are repeated. Reiterate through the thickness until equilibrium is fulfilled. Ill-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.

PAGE 40

Change Tria 1 Sm Find Neutral Axis yes Printout Data yes urvatur Figure 3-2 The Flow Chart of Pure Bending Program 28

PAGE 41

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 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. FR1 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 FR1 at the point of contact with the wall, while 1 is the angle formed by FR1

PAGE 42

30 with the vertical. Note that, for equilibrium, the magnitude and direction of the resultant force FR1 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 III-3-1 Mathematical Formulation The unknown angle 1 is usually related to parameters defined by the boundary conditions. Figure 3-4 represents the angle 1 at the connection between the sheet metal and the retaining wall. Typical boundary conditions at the inlet are given by: @ A: FN = FNAi Fv = FvA where FNA1 FvA are the resultant normal and shear components of FRl 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 u( s, z) action along the direction s. The de:fini tions 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. I I Figure 3-4 The Angle 1 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: (} ,I. -1 ( FNi) i = 'f'l Sin FRl (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 Fn1 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 FR1 is constant in both magnitude and direction along the whole region. Because the resultant force FR1 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 K1 2. a trial value for the length of a finite element "i" after one step deformation .6-stmi is assumed 3. using constancy of volume (equation 3-7), the thickness of the element "i" for the assumed .6-stmi is determined. 4. the corresponding length .6-s 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 .6-sj, 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 .6.stm1are now calculated using equations 3-3, and 3-4. 8. equilibrium (equation 3-13) must be satisfied; otherwise, another trial.6.stmi is assumed and steps (3-8) are repeated. Reiterate through the thickness until equilibrium is fulfilled. 9. the angle 8i for the element "i" for the corresponding trial value FR1 ( 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 .6.Btmi for which equilibrium was fulfilled (step 8) with the calculated length .6.Sttmi (step 10). The trial curvature Ki is now readjusted and the steps 2 through 11 are repeated until .6.stmi 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 < CBYR where: Yi is the coordinate of the sheet metal at "i" element in Y-direction; CBY is the coordinate at the center of the pin in Y-direction; 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 Tria 1 Sm Change No < K1 > yes Find Angles yes Output ncr ease ria 1 Fa1 ecrease Trial Fa1 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-1 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-llists 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: B = l o.ult.c, n-1 n=1 (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 'f 0.177 Co 0.0 Ex (psi) 20E6 Dco l.OE-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 test file (n=lOO &. n=180) 140 120 100 80 j' 60 i _j_-:-I ---I / vi I / i 1/ .. 40 1':1 "' El 20 0 0 -20 -40 I v I /i i /VI I L-----l1----I I I -60 -0.5 0 0.5 1.5 2 2.5 3 3.5 4 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 1=1 0 :r Q cu .!:1 j:l ii "' :r ... cu > cu 1:1 0 < 1=1 0 :r "' Q .s cu ..c:: E6 5.9 5.8 5.7 5.6 5.5 5.4 5.3 5.2 0 The Profile of the Sheet Metal I I --I I I I i 0.5 1.5 2 The Location Along the Horizontal Direction (in) Figure 4-2 The Profile of the Sheet Metal. \ [\ I 2.5

PAGE 56

44 The between the two ends is unsupported area. The resultant force, FRb 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 wall, 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. .. Cl .. 0 r... ;; -;; < ... Cll s .. 0 z .. ..c: 1--40 20 0 -20 -40 -60 -80 -100 0 The Normal Axial Force vs. The Horizontal Location of Sheet Metal ; I \ I \ I \ I \ \ i i\ 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

d i b .... d Q) s 0 500 450 400 350 300 250 200 150 100 I 1 I I 50 1/ 0 0 Moment-Curvature Curve Along the Sheet Metal I J----r----1 v--!---" I I v I I I I I I l 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Curvature (1/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 -1 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 di:ft'erent 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 in finite 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 have 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 assumed 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 u = stress e =strain K = strength coefficient n = strain hardening effect on the sheet due to the forming process t:Ti = principal stress in direction i Ei = 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

51 Chapter III Ki Curvature of element "i" Smi -Length of element "i" at mid-surface FNi Normal force resultant Mi resultant moment Fj the normal force acting on fiber "j". 1/n The neutral axis corresponds to the location of the fiber of the element "i" normal anisotropy parameter. c ( e*, u*) the origin at the beginning of the u.nloading The total effective strain of the material fiber "j" the X coordinate of FR1 at the point of contact with the wall tP1 the angle formed by FR1 with the vertical

PAGE 64

APPENDIX B THE USER'S MANUAL OF PURE BENDING PROGRAM

PAGE 65

53 THE USER'S MANUAL OF PURE BENDING PROGRAM For this case, the input 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 u = k jjen 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 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 ( 11) 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 WHOSE OUTPUT DATA CAN BE USED IN C = MATLAB AND OPEN AN INPUT DATA FILE CALLED INPUT.DAT c c c open(unit=40,file='input.dat',status='old') open(uni t=35 ,file=' out .m'. status= '-new') open(unit::i45,file='mk.dat' ,status='new') open(unit=55,file='shistory.dat,status='new) C = CALL SUBROUTINE INPUT C THE PURPOUSE: C READ AND CALCULATE THE ORIGGINAL PARRAMETERS c call input c if((what .eq. 'pure') .or. (what .eq. 'PURE')) then 62

PAGE 75

c C = CALL SUBROtrl'INE 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 inputerr 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 vrite(45,999)title if ( ( vhat eq. ) or. ( vhat eq PURE )) then vrite(45,5000) end if vrite.(45,6000) 100 read(40,1000)vhich i't ( (vhich eq. 'mate') .or. (vhich eq. 'MATE' ))then read(40,)coef, en, r, ex, yield vrite(45,1001)vhich. vrite(45 ,*) coef =', coef vri t e( 45 ) en = en vrite(45,) r =', r vrite(45,) ex =, ex 'irrite(45,)' yield =' yield go to 100 else if((vhich .eq. 'init').or.(vhich .eq. 'INIT'))then read(40,*) smO, hO,curvO, dcurvO, ncurv vrite(45,1001)vhich vrite(45,)' smO =' smO vri t e ( 45 *) hO = hO vrite(46 ,*) curvO =', curvO vrite(45,)' dcurvO = dcurvO vri t e ( 46 *) 1 ncurv = ncurv vrite(45,*)'======================================' 64

PAGE 77

999 1000 1001 5000 ()000 2000 go to 100 else if((which .eq. 'para').or.(which .eq. 'PARA'))then read(40,*) netaO,_fntol, sm2tol, npathe write(46,1001)which write(45,*)'. write(45,*)' write(45,*)' netaO = 1 netaO fntol =' fntol sm2tol ='; sm2tol write ( 45, *) 1 npathe = 1 npathe write(45,*)'======================================1 go to 100 else if((which .eq. 'end 1 ) .or. (which .eq. 'END '))then go to 2000 else if( (which .ne; 'mate') (which .ne. 1 init 1 ) + .and. (which .ne. 'MATE') .and. (which .ne. 'INIT') + .and. (which .ne. 'end ') .and. (which ,ne. 'para') + (which .ne. 'END ') .and. (which .ne. 'PARA') + ). then + + + + + + + call inputerrQr end if format (a70). format(a4) format(1x,a4) format ( 1 This is required to fine the relationship' ,/, between pure bending moment and' curvature'.,//)' format ( 1 ===========================.;.==================== 1 ===============,;,.========='====== I I retUrn end Verify the input data',//, '===============================================' '====:===========================',//) 65

PAGE 78

c C = DECK OF SUBROUTINE INPUTERR c subroutine inputerr print print vrite(45,1000) 1000 format (, -+ return end Your option is NOT included., Check that again ! Your option is NOT included.,,/, Check that againn!!!! !1 ) 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+OO + r) I ((1.0d+OO.+ 2.0 d+OO r)**O.S) 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 ****** CURV1 = CURRENT CONFIGURATION OF CURVATURE c ****** CURV2 = UPDATE CONFIGURATION OF CURVATURE c ****** SMO = LENGTH OF THE NETURAL LINE OF ORIGINAL c CONGIGURATION c ****** SM1 = LENGTH OF THE NETURAL LINE OF CURRENT c CONGIGURATION c ****** SM2 = LENGTH OF THE NETURAL LINE OF UPDATE c CONGIGURATION c ****** HO = THICKNESS OF THE OF ORIGINAL CONGIGURATION c ****** H1 = THICKNESS OF THE OF CURRENT CONGIGURATION c ****** H2 = THICKNESS-OF THE OF UPDATE CONGIGURATION c ****** NETAO = NUMBER OF ETA OF ORIGINAL CONGIGURATION c ****** NETA = NUMBER OF ETA OF CURRENT CONGIGURATION c ======-=========-============================================== C TRANNSFER THE INITIAL DADAS TO THE CUNNRENT CONFIGURATION c ============================================================= jcount= 0 jcount1 = 1 curv1 = curvO sm1 = smO nsm2 = 0 csm2 = sm2tol h1 = hO neta = netaO xneta = neta deta = 1/xneta c ============================================================= C CALCULATE ETA c ============================================================= do 100 i = 1, neta+1 e.ta1 (:l.) = (i-1)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 = 0. 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.1875dcurv0))/(ncurv-1) do 102 i = 1, neta+1 zO(i) = h1 (i-1) detah1/2.0 102 continue c ============================================================= c C ****** CALL SUBROUTINE POSITIVE TO CALCULATE THE C RELATIONSHIP C BETWEEN MOMENT AND CUVATURE FOR POSITIVE CUVATURE c 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 if + + + + + + + jcount = jcount + 1 curv2 = dcurvO + (1.0/1.875-dcurvO)/xneta*(jcount-1) if(ncurv .ge. 31) then ( (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 = dcurvOexp((jcount-1)b) end if else if(ncurv .ge. 21) then if ( (jcount .eq. (ncurv -2)) .or. + (jcount .eq. (ncurv -4)) .or. + (j count eq. (ncurv -6)) ) then go to 1111 else curv2 = dcUrvOexp((jcount-1)*b) end if 70

PAGE 83

else curv2 = dcurvQexp((jcount-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 ****** Z1 = THE THICKNESS AT EACH FIBER OF CURRENT c CONFIGURATION c ****** Z2 = THE THICKNESS AT EACH FIBER OF UPDATE c CONFIGURATION c ****** DEAT = INCREMENT OF ETA c ****** ZCOM3A= PARAMETER FOR CALCULATION Z1 c ****** ZCOM4A= PARAMETER FOR CALCULATION Z1 c ****** ZCOM3 = PARAMETER FOR CALCULATION Z2 c ****** ZCOM4 = PARAMETER FOR CALCULATION Z2 c ============================================================= C = SET UP THE UPDAT THICKNESS h2 TO BE h1 C = AND CALCULATE Z1 AND Z2 c ============================================================= 2222 h2 = (h1sm1)/sm2 if (curv1 .eq. 0.0) then do.200 i = 1, neta+1 z1(i) = h1* (i-1) data h1/2 200 continue else do 300 i = 1, neta+1 zcom3a(i) = (1.0d+00-(0.50d+OOcurv1h1)) zcom4a(i) = sqrt(zcom3a(i) + 2.0d+OO + h1 curv1 eta1(i)) z1(i) = (-1 + zcom4a(i)) 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-1) deta h2/2 400 continue else do 500 i = 1, neta+1 zcom3(i) = (1.0d+OO (0.50d+OOcurv2*h2))**2 zcom4(i) = sqrt(zcom3(i) + 2.0d+OO + h2 curv2 eta2 (i)) z2(i) = (-1 + zcom4(i)) I curv2 500 continue end if c ============================================================= C CALCULATE THE LENGTH OF FIBER FOR ,EACH LAYER c ============================================================= c ****** c c ****** c c ****** c c ****** c ****** c ****** c c ****** S1 = THE LENGTH OF EACH FIBER AT DIFFERENT LAYER OF-CURRENT CONFIGURATION S2 : THE LENGTH OF EACH FIBER AT DIFFERENT LAYER OF UPDATE CONFIGURATION EPSILDN2 = STRAIN OF EACH LAYER OF UPDATE CONFIGURATION SG2 YIELD = STRESS OF UPDATE CONFIGURATION = YIELD STRESS OF THE MATERIAL ZNETURAL = THE LOCATION OF NETURAL AXIS (AT THE FIBER OF STRESS EQUAL TO ZERO) 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+curv1z1(i)) s2(i) = sm2(1+curv2z2(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) sig2(i) = log(abs(s2(i)/s1(i))) = c ex epsilon2(i) c ============================================================= C = CHECK THE ELASTICC RANGE c .lt. then sig2(i) = -c coe:f. 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 ============================================================= .lt. 0.0_ .and. check(i) .gt. O.O)then znetural = z2(i-1)+ detaabs(sig2(i-1)/ + (sig2(i-1)-sig2(i) ) ) sm2old = sm2 sm2new = sm2* (1 + curv2 znetural) dsm2 = abs (sm2new -sm2) if (dsm2 .le. sm2tol) then sm2 = sm2new else sm2 = sm2new 73

PAGE 86

c ============================================================= C =IF' IT IS NOT.CONVERGED THEN GOTO. 2222 TO MODIFIED h2 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 (1). DETERMINE THE FORCE FOR EACH LAYER C (2). CALCULATE THE SUMMATION FORCE FOR THE STRUCTURE C (3). DETERMINE THE MOMMENT ABOUT THE NETURAL FOR EACH LAYER C (4). CALCULATE THE MOMENT ABOUT THE NETURAL AXIS FOR THE C .STRUTURE c c ============================================================= 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 = fn fn 0.0 totalm = 0.0 c ============================================================= 74

PAGE 87

C = CALCULATE STRESS(dfn) AND MOMENT(dmn) FOR EACH FIBER 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 = 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; i875) c ============================================================= if( curv2 .ge. ((1.0/0.1875)-0.00001) .or. + j count ge. (ncurv 2)) then vork = 'fini' end if call vriteout c C = UPDATE ALL F THE PARAMETERS c ============================================================= sin2 = smO if (work .eq. 'fini') then C C =CHECK THE. LOADING AND UNLOADING CONDITIONS C = NPATHE : 'NUMBER OF THE PATHES C = NPATBE = 1 OR 0 (OR 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 if((npathe.eq.O).or.(npathe.eq.1))then 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 ============================================================= tol = fn/(fnold-fn) else if(fn .gt. fntol .and. fnold .gt. fntol) then if (abs(tol) .le. 30.0) then sm2 = sm2 -tol (sm2old-sm2nev) else sm2 = sm2 dsm2/2.0 end if go to 2222 else if(-fn .gt. fntol .and. -fnold .gt. fntol) then if (abs(tol) .le. 30.0) then sm2 = sm2 -tol (sm2old-sm2nev) 76

PAGE 89

else sm2 = sm2 + dsm2/2. 0 end if go to 2222 else if((fn .gt. fntol .and. -fnold .gt. fntol) .or. + (-fn.gt. fntol .and. fnold .gt. fntol))then else end if return end sm2 = (sm2nev+sm2old)/2.0 go to 2222 go to 2222 77

PAGE 90

c ============================================================= c C = SUBROUTINE INPUTERROR c c ============================================================= subroutine inputerror vrite(45,100) vrite(*, 100) 100 format(J Fatal Error ?????? + J,//, + + + return .. end Your input data formation is not rightJ !!!!!',//, Check that again, I am sorry! ) 78

PAGE 91

c ============================================================= c C = DECK OF SUBROUTINE WRITEOUT c c ============================================================= subroutine vriteout include 'para.l' vrite(45,10()) vrite(45,200)totalm vrite(45,300)curv2 vrite(45,400)znetural vrite(45,600)sm2 100 ormat(' 101 ormat(' + '===================================='.Ill) 200 ormat (' The Moment is : 2x 1 e14. 7) 300 ormat(' The curvature is :',2x, e14.7) 400 ormat (' The netural axis is at :' 2x t e14. 7) 500 ormat ( 1 The length o natural axis is :', 2x, e14. 7) c ============================================================= C HERE IS STRAIN HISTORY c ============================================================= i (jcount .eq. 1) then vrite(55,600) vrite(55,601)jcount else vrite(55,601)jcount end i 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) n, totalm 600 ormat(' This is the output data ile or strain' + history.',/) 79

PAGE 92

601 + 602 format('1',' This is the',1x,i4,1x,'step of strain history.') format('# Fiber 1 1 I',' M .I' Fn I iteration' I S2 + + + '==========================================='' + '===============================') 603 format(4x.i4, 2(1x,' I', 2x, e17.9),1x,'l' ,2x,e20.12) 604 format(' Total Fn = e14.7, + Total Moment= 1 ,e14.7,///////) c ============================================================= C HERE IS CREATING THE DATA FILE FOR MATLAB TO PLOT M-K PLOT c ============================================================= 703 + + 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 format(/, 'f.', 1 J Total Moment curvature I 'f.',' ========================' + '==========================',//) 704 format('mk1=[' 2(2x, e20:12)) 705 format(7x, e20.12, 2x, e20.12) 706 format(7x, e20.12, 2x, e20.1"2, 1]; ') c ============================================================= C = SAVE THE YIELDING.MOMENT(YMOMENT) AND YIELD CURVATURE C (YCURVE) c ============================================================= + 212 if ( (abs ( sig2 (1) (abs(sig2(neta + goto 212 else ymoment = totalm ycurve = curv2 nyield = jcount1 end if ) 1) ) .ge. yield ) .ge. yield ) .or. ) then 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) + -sa vee (nyield) p2totalm(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) = p.2curv2(jjj) + savec(i) p3totalm(i) = p2totalm(jjj) + savem(i) p3z:D.etu(i) = p2znetu(jjj) + savez(i) p3sm2(i) = p2sm2(jjj) + saves(i) 100 continue kkk = 1 + nyield jjj = jmax + nyield do 200 i = kkk, jjj p3curv2(i) = p2curv2(jjj) + savec(i-nyield) + + savec(nyield) p3totalm(i) = p2totalm(jjj) + savem(i-nyield) + + savem(nyield) = p2znetu(jjj) + savez(i-nyield) + + savez (nyield) p3sm2(i) = p2sm2(jjj) + saves(i-nyield) + + saves(nyield) 200 continue call patheolit return end 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 100 i = 1, JJJ 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(3S,*)'% Should be continued!' do 200 i = 1, jjj if (i .eq. 1) then vrite(36,703) vrite(36,7044)p2totalm(jjj),p2curv2(jjj) vrite(36,705)p3totalm(i), p3curv2(i) 84

PAGE 97

200 else if (i .ne. jjj) then write(36,705)p3totalm(i), p3curv2(i) else write(36,706)p3totalm(i), p3curv2(i) end if continue end if if ((npathe .ne. 3) write(35,*) write(36,*) end if .and. (npathe .ne. 2)) then Your path number is not right.' Check that again c ************************************************************* 703 format(/. 'Yo',' Total Moment curvature I + + Yo =======.================= + '==========================' ,//) 704 format ( 'mk2=P 2 (2x. e20 .12)) 7044 format ( ,-mk3= [' 2 (2x. e20 .12)) 706 e20.12, 2x, e20.12) 706 format(7x, e20.12, 2x, e20.12,'] ;') return end 85

PAGE 98

CONTROL.L character*70 title character*4 what, which common/ control/title, what, which 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, sm1, sm2, dsm2, sm2new, sm2tol double precision netaO, fntol, hO, h1, 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 s1(maxdim). s2(maxdim). + epsilon(maxdim). + epsilon2(maxdim). check(maxdim) double precision dfn(maxdim). dmn(maxdim). + fdfn(maxdim). fdmn(maxdim). + (maxdim) s ig2 (maxdim) + z1(maxdim). z2(maxdim). zO(maxdim) double precision eta1(maxdim). eta2(maxdim) integer jcount, jcount1, jmax, nyield integer neta, nerr, nsm2, ncurv, npathe C = COMMMONS c + + + common/ common/ common/ common/ common/ common/ common/ common/ common/ common/ charac/ work mp/ coef, en, r. ex, yield iv1/ curvO, curv1, curv2, ncurv, dcurvO,b,c iv2/ smO, sm1, sm2, dsm2, sm2new, sm2tol pa1/ netaO, fntol, nerr. jcount pa2/_ fn, totalm, eta; deta, csm2, sm2, ffn, ftotalm pa3/ znetural, hO, h1, h2. j count 1 pa4/ zcom3, zcom4,zcom3a, zcom4a paS/ s1, s2 ex1/ epsilon2, check, epsilon 87

PAGE 100

88 COJ11liJ.On/ orce1/ dn, dmn, sig1, sig2, z1, z2, + zO, dn, dmn common/ two eta/ eta1, eta2 common/ save/ savem, savec, savez, saves common/ paths/ jmax, npathe, ymoment, ycurve, + nyield

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 I MATE COEF EN EX YIELD MATE: material properties (only the first four characters will be readied) COEF : the coefficient, k of the equation u = k ieien 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 105

90 I INIT 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 originaJ. 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 I PARA NETAO PARA : parameters Example: parameters 100

PAGE 107

I GEOM Dl RPl GEOM : geometry Dl: the vertical distance between the wall and the top of the pin RPl: the radius of pin 1 Example: geometry 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

APPENDIXF STRETCHING FOR WALL-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 o.5, o.5, o.5 end of file

PAGE 112

APPENDIX G 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 BMETHODS WILL BE EXTENTED TO THREE DIMENSIONAL C CASE. c = C CREAT DATE C MODIFIED DAtE C FINAL MODIFIED DATE C AUTHOR C .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 11 OUTPUT. STR 11 TO STALL C THE STRAIN HISTORY ALONG THE SHEET METAL c open(unit=40, fil,e='input.s.tr' status='old') open(unit=36, file='output.aaa', status='nev', + carriagecontrol = 'list') open(unit=45, file='output.str', status='nev'.

PAGE 114

+ carri'agecontrol = 'list' ) open(unit=56, file='profile.m' status='new'. + carriagecontrol = 'list') open(unit=65, file='fm.m' status='new', + carriagecontrol = 'list') open(unit=75, file='output.m' status='new'. + carriagecontrol = 'list') c C = CALL SUBROUTINE INPUT. STR C THE PURPOSE OF THIS SUBROUTINE IS : C READ 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 wallpin else if ((geometry (geometry call.tvopili else if ((geometry (geometry call threepin else 'wall') .or. 'WALL')) then eq. 'twop') .or .eq. 'TWOP')) then eq. 'thre') .or .eq. 'THRE'))then 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 SUBROUTINE 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 vallpin include 'paramet.str' C =TRANSFER ALL OF THE DATA FROM THE INITIAL DATA TO THE C CURRENT CONFIGURATION DATA c call data vrite(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 vrite(45,2001)icount vrite(*,2001) icount ybold = yb xbold = xb call coord curvbold = curvb = 1.0d+OO I (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 (???) c if( ( yb .lt. rp1) .and. (curv2 .gt. 0.0) )then c C = MAKE SOME PARAMETERS IN ORDER TO BE USED IN C NEWTON-RAPHSON METaODS c fv0old2 = fvOold fvOold = fvO dtheold = dtheta the old = thetab thetab = acos(yb/rp1) dtheta = thetab -theta c C = IF THE ANGLES CONVERGENT INTO THE TOLERANCE, IT IS OK! C GO TO RETURN c + + vrite(45,2002) vrite(.2002) if (abs(dtheta) .lt. angletol) then vork = 'fini' vrite(46,2003) xb, yb, fnb2, fra1, fna1, total.m2. theta. curv2 vrite(,2003) xb, yb, fnb2, fra1, fna1, totalm2, theta, curv2 102

PAGE 118

write(65,5503) xb, yb go to 1099 c C = IF THE ANGLES ARE NOT CONVERGE, USE NEWTON-RAPHSON METHODS c + + else if(curvb .gt. curv2) then else write(45,2006) xb, yb, fnb2, fra1, fna1, totalm2, theta, curv2 write(*,2006) xb, yb, fnb2, fra1, fna1, totalm2, theta, curv2 fvO = fvO 1.15 call newdata write(45,2004) fvO write(55,5504) fvO write(*,2004) fvO go to 1002 if(((dtheta .gt. 0.0) .and. (dtheold .le. 0.0)) .or. 103 + + ((dtheta .le. 0.0) .and. (dtheold .gt. 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 = 0.5 (fv0old2 + fvOold) call newdata write(45,2006) fvO write(56 ,5605) fvO write(*,2005) fvO goto 1001

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 vrite(*,2005) fvO goto 1001 else fv01 = fvO fvO = fvO -(dtheta) (fvOold -fv0old2) I (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(45,2005)fv0 write(55,5505) fvO vrite(*,2005) fvO goto 1001 end if 104

PAGE 120

end if endif c C = IF THE yb NOT 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, + totalm2, theta, curv2 write(,2007) xb, yb, fnb2, fra1, fna1, + totalm2 theta, curv2 write(55,5502) xb, yb call trans goto 1001 end if C ******END THE'1ST LOOP 2001 format(' ========================================' 2002 2003 + + + The 1 I4,' segment has been finished! 1 ) format ( 1 Check the tangency! I) format( 1 The tangency has been reached. I 'I' Write out all of the data!',/, 105 + 6x, I xb =' E18.10, 5x, 'yb = E18.10,/, + 5x, 1 fnb2 =' E18.10, 5x, 'fra1 = E18 .10 ,/, + 5x, I fna1 =' E18.10, 5x, 'totalm2 = E18.10,/, + 5x, I theta =' E18.10, 5x, 1curv2 = I E18.10) 2004 format(' The trial vertical applied force is not', + reach the tangency.', I, + Change the new fvO = 'E18.10,1x, + 'and continue!',/, + + + + *****************************************' '*******************************',/, All of the results written previous 1 'are NOT good.',/,

PAGE 121

+ + *****************************************' '*******************************') 106 2005 format(' The trial vertical applied force is not', 2006 2007 + reach the tangency.', I, + Change the new FvO = 'E18.10,1x, + 'and continue!',/, + + + + + + + + + + + + + + + + *****************************************' '*******************************'' All of the results written previous 'are NOT good.',/, *****************************************' '*******************************') format (' The vertical 'top of pin and 5x, xb =' 5x, ,. fnb2 =' 5:i:, fna1 =' 5x, theta =' distance has reached the the data are:' .1. E18.10, 5x, 'yb E18.10, 5x, 'fra1 E18.10, 5x, 'totalm2 E18.10, 5x, 'curv2 = E18.10,/, = E18. 10 ./, =' E18.10,/, = E18.10) format(' The 'top of 5x, xb 5x, fnb2 vertical distance does not reach pin, yet. The data are:',/, the 5x,. 5x, fna1 theta =' E18.10, 5x, 'yb = =' =' =' !18.10, E18.10, E18.10, 5x, 5x, 5x, = 'fral 'totalm2 = 'curv2 = E18 .10 ,/, E18.10,/, E18.10 ./, E18.10) 5501 format(' Y. ----x -----==== y ==== ',/, 5502 5503 5504 + p=[' F15.5, ',', 3x, F15.5 ) + + + + + + + + format( 3x, F15.5, ',',.3:i:, F15.5 ) format( 3x, F15.5, '', 3x, F15.5 ];' ) format(' The trial vertical applied force is not', 'reach the tangency.',/, Change the new Sm1 = 'E18.10,1x, 'and continue!',/, *****************************************' '***************************'' All of the results written previous 'are NOT good.',/, *****************************************'

PAGE 122

+ 5505 + + + + + + + + + 1099 '***************************') format(' The trial vertical applied force is not', return end reach the tangency.',/, Change the nev FvO = 'E18.10,1x, 'and continue!',/, *****************************************' '***************************' ,/, All of the results vritten previous 'are NOT good.',/, *****************************************' '***************************') 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 Td BE THE FUNCTION OF h1 C ****** CALCULATE Z1 AND Z2 c ================================================================ c C ****** CALCULATE Z1

PAGE 124

c 1112 h2 = (h1 sm1) I sm2 if (curv1 .eq. 0.0) then do 200 i = 1, neta + 1 z1(i) = O.Od+OO z1(i) = h1 (i 1) deta-h1 l-2.0d+OO 200 continue else do 300 i = 1, neta + 1 zcom3a(i) = O.Od+OO zcom3a(i) = (1.0d+OO + (0. 5d+OO curv1 h1) ) **2. Od+OO zcom4a(i) = O.Od+OO zcom4a(i) = sqrt (zcom3a(i) + 2 .Od+OO h1 + curv1 eta1 (i)) z1(i) = O.Od+OO z1(i) = (-1 + zcom4a(i)) I curv1 300 continue endif c C ****** CALCULATE Z2 c if .eq. 0.0) then do 400 i = 1, neta + 1 z2(i) = O.Od+OO z2(i) = h2 (i 1) deta-h2 I 2.0d+OO 400 continue 500 else do 500 i = 1, neta + 1 zcom3(i) = O.Od+OO zcom3(i) = (1.0d+OO + (O.Sd+OO curv2 h2))**2.0d+OO + zcom4(i) = O.Od+OO zcom4(i) = sqrt(zcom3(i) + 2.0d+OO h2 z2(i) z2(i) continue curv2 eta2(i)) = 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 52 c ================================================================ do 600 i = 1, nata+ 1 s1(i) = O.Od+OO s1(i) = sm1 (1.0d+OO + curv1 z1(i) ) s2(i) = O.Od+OO s2(i) = sm2 (1.0d+OO + 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, nata + 1 c C ****** INITIALIZE THE STRAIN ARGUMENTS c epsilon2(i) = O.Od+OO c C ****** CALCULATE THE STRAINS EPSILON2 c c ???????????????????????? epsilon2(i) = log(abs(s2(i) I s1(i))) sig2(i) = c ex epsilon2(i) c C ****** CHECK THE ELASTIC RANGE c

PAGE 126

c if(sig2(i) .lt. -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) .lt. 0.0) .and. (check(i) .gt. 0.0)) then znetural = z2(i 1) + deta abs(sig2(i -1) I + (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 znetural) dsm2 = abs(sm2nev sm2old) 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

c C = THE FOURTH PURPOSE C (totalm) THE NETURAL AXIS FOR EACH LAYER (dmn) SUM THE MOMENT ALL OVER THE LAYERS 112 c ================================================================ c C ****** INITIALIZE fn AND totalm c c fn = O.Od+OO totalm = O.Od+OO C ****** CALCULATE STRESS (dfn) FOREACH 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) dmn(i) = dfn(i) fn = fn totalm = totalm 800 continue abs(z2(i) z2(i) + dfn(i) + dmn(i) -z2(i + 1)) c ================================================================ C =THE FIFTHPURPOSE : 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-1)] + Ki [M(i)-M(i-1)] = 0 c ================================================================ c C ****** CALCULATE THE PARAMETERS c fnb2 = fn totalm2 = totalm

PAGE 128

c ****** c ****** c diffn = fnb2 -fna1 difmn = totalm2 -totalm1 equiold2= equiold equiold = equi equi = diffn + curv2 difmn sm2old2 = sm2old sm2old = sm2 C ****** CHECK THE EQUILIBRIUM CONDITION c if (abs (equi) .le. equi tol) then c write(*,*) fnb2 =', fnb2 c write(*,*)' fra1 =', fra1 goto 1114 c ****** else if( ((equi .gt. 0.0) + .or. ((equi .le. 0.0) + then .and. .and. (equiold (equiold sm2 = 0.5d+OO (sm2old2 + sm2old) 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) (sm2old sm2old2) I (equi -equiold) go to 1112 end if C = CALCULATE THE ANGLE THETA c 1114 ratio1 = fnb2 I fra1 c if(ratio1 .gt. 1.0d+OO) then curv2 = 1.2 curv2 go to 1111 else if(-ratio1 .gt. 1.0d+OO) then CUrV2 = 0.8d+00 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) =? S(tmi) c sm2t = (theta -thetao) I curv2 sssold = sss sss = sm2t -sm2 curvold2 = curvold curvold = curv2 if(abs(sss) .le. sm2tol) then 114

PAGE 130

goto 1115 else if(((sss .gt. 0.0) .and. (sssold .le. 0.0)) .or. + ((sss .le. 0.0) .and. (sssold .gt. 0.0)) )then curv2 = 0.5 (curvold2 + curvold) sm1 = smO go to 1111 else if ( (abs(sss -sssold) .le. 0.00000001) .or. + (sssold .eq. 0.0 )) then + end if else curv2 = curv2 1.001 sm1 = smO go to 1111 curv2 = curv2 -(sss) (curvold -curvold2) I (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 Jparamet.strJ 116 c ================================================================ c C = ARGUMENTS c c c c c c c c c c c c c title vhat geom The title of this input data file The purpose of program (in this fileis 11STRETCH11 so far, in the futher may combine pure bending moment other capacities) case, or some Types of geometry: (1). Wall-pin (2). Wall-Tvo-pin (3). Wall-Three-pin (4). Wall-Tvo-pin (only function aviable nov) (in the further) (in the further) (in the further) C (5). Wall-N-pin (in the further) C (6). Three Dimensional (in the further) c ================================================================ c C = READ THE TITLE AND GEOMETRIC TYPE c read(40,1000)title read(40,2000)vhat read(40,2000)geometry if( vhat .ne. JstreJ) then call inputerror(1)

PAGE 132

else if( (geometry .ne. 'vall') .and. + (geometry .ne. 'WALL')) then call inputerror(2) endif c C =WRITE OUT THE CHARACTERS TO FILE OUTPUT.STR IN ORDERR TO C VERIFY THOSE DATA RIGHT. c c ****** c ****"'* vrite(45,*)title if((vhat .eq. 1stre1 ) .or. (vhat .eq. 'STRE') )then vrite(45,9001) else call inputerror(1) endif 117 if((geometry .eq. 'vall') .or. (geometry .eq. 1Wall'))then vrite(45,9002) else call inputerror(1) endif c ****** vri te ( 45, 9003) c C = READ THE REST OF THE DATA c 100 read(40,2000)vhich c ****** if( (which .eq. 'init') .or. (which .eq. 1INIT1)) then read(40, *) smO, hO, curveO, coefmu, dcurve, fvO vrite(45,4002) vrite(45,*)1 smO = smO vrite(45,)1 hO = hO vrite(46,*)' curveO = J c:urveO vrite(46,*)1 coefmu = I coefmu vrite(46,*)' dc:urve = dcurve vrite(45,*)1 fvO = fvO go to 100

PAGE 133

c ****** c ****** c ****** c ****** c ****** c ****** + + + + + 118 else if((vhich .eq. 'mate').or.(vhich .eq. 'MATE'))then read(40, ) coef, en, r, ex, yield .vri te (45 ,4001) vrite(45,*)' coef = coef vrite(45,*)' en = en vrite(45,*)' r = r vrite(45,*)' ex = ex vrite(45,*)' yield = yield go to 100 else if((vhich .eq. 'geom') .or. (vhich .eq. 'GEOM')) then read(40, *) d1, rp1 vrite(45,4003) vrite(45,*)' d1 =' d1 vrite(45,*)' rp1 =' rp1 go to 100 else if((vhich .eq. 'para').or.(vhich .eq. 'PARA')) then read(40, *) neta vrite(45,4004) vrite(46,*)' neta =' neta go to 100 else if((vhich .eq. 'tole').or.(vhich .eq. 'TOLE')) then read(40, ) sm2tol, angletol, equitol vrite(45 ,4005) vrite(45,*)' vrite(46,*)' vrite(45,*)' go to 100 sm2tol =' angletol =' equitol =', sm2tol angletol equitol else if((vhich .eq. 'end ').or. (vhich .eq. 'END ')) then go to 2200 else if((vhich ne. 'init').and.(vhich .ne 'INIT').and. (vhich ne. 'mate').and.(vhich .ne 'MATE').and. (vhich ne. 'geom') .and. (vhich .ne 'GEOM').and. (vhich ne. 'tole').and. (vhich .ne 'TOLE').and. (vhich ne. 'end ').and.(vhich .ne 'END ').and. (vhich .ne. 'para').and.(vhich .ne. 'PARA'))then

PAGE 134

call inputerror(3) c ****** endi c C = WRITE OUT THE INPUT DATA TO FILE OUTPUT. STR IN ORDERR TO C VERIFY THOSE DATA RIGHT. c 2200 vrite(46,9003) c C = FORMAT c 1000 ormat(a70) 2000 ormat(a4) 4001 format(' ****** The material properties are: 4002 4003 4004 4005 9001 9002 9003 + format(' format(' format(' format(' format(' format(' format(' ****** The initial values are: ****** The geometric properties ****** The parameters are: ****** 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 vrite(,1000) else if( n .eq. 2) then vrite(,2000) else if( n .eq. 3) then vrite(,3000) eildif format (' format (' format(' return end WARNING ! I 'I Your input data is WRONG on the fist line,, ,I, Check that again!') WARNING ! ,/ Your input data is WRONG on the second, I 1' I I 1ne. Check that again!1 ) WARNING ! I 'I' Your input data is WRONG after the third' I 1' I I 1ne. Check that again!1 )

PAGE 136

121 c ================================================================ c C = END OF 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 :xneta deta jcount icount kconut = (1.0d+OO + en)/( (1.0d+OO + 2.0d+OO en).5 ) = acos(O.Od+OO) = neta = 1.0d+OO I 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 cyp1 = d1 + rp1 h1 = hO sm1 = smO C -----COORDINATES xa = O.Od+OO ya = cyp1 0.5 hO yb = ya ybold = yb C -----FORCES AND MOMENTS fva1 = fvO fna1 = fvO coefmu fra1 = sqrt(fva1**2.0d+OO + fna1**2.0d+OO) totalm1 = O.Od+OO C -----ANGLES AND CURVATURES c phi1 thetao theta curp1 = atan(coefmu) = phi1 -asin(fna1lfra1) = thetao = 1.0d+OO I rp1 curv1 = curveO curvold.= O.Od+OO curvold2= O.Od+OO curvb = 1.0d+OO I (rp1 + 0.5d+OO hO) return end C = END OF FILE c 122

PAGE 138

c C = END OF SUBROUTINE DATA c c C = DECK OF SUBROUTINE NEWDATA c subroutine nevdata include 'paramet str.' C ------CHARACTERS control = 0 K c .. C = TRANSFER ALL OF 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 ETA2 do 100 i = 1, neta eta1(i) = (i.-1) deta eta2(i) = (i -1) deta 100 continue C -----LENGTH (DISTANCE) AND THICKNESS cyp1 = d1 + rp1 h1 = hO sm1 = smO C -----COORDINATES 123

PAGE 139

xa = ya = cyp1 0.5 hO yb = ya ybold = yb C FORCES AND MOMENTS fva1 = fvO fna1 = fvO coefmu fra1 = sqrt(fva1**2.0d+OO + fna1**2.0d+OO) totalm1 = O.Od+OO C -----ANGLES AND CURVATURES c phi1 thetao theta CUrV1 = atan(coefmu) = phi1 -asin(fna1/fra1) = thetao curveO curvold :::: O.Od+OO curvold2= O.Od+OO curvb = 1.0d+OO / (rp1 + .0.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 xb return end = yb = ybold + (l.Od+OO I curv2) ( cos(theta) -cos(thetao)) = xa + (totalm2 -(fra1 sin(phi1) (ya -yb)))l fra1 I 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+OOfnb2**2.0d+OO) fna1 = fnb2 C -----CURVATURES 126

PAGE 142

c curv1 = curv2 curv1 = curveO curvold = O.Od+OO C ------LENGTH AND THICKNESS sm1 = smO h1 hO C -----COORDINATES c c ybold = yb return end C = END OF SUBROUTINE STRANS c 127

PAGE 143

c C = PARAMETERS AND DIMENSIONS c (maxstep = 1000, maxdim = 500) c C = CHARACTERS c character 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 DCURVEFVO -. THE ORIGINAL LENGTH OF THE SEGMENNT THE ORIGINAL HEIGHTOF THE SEGMENT CURVATURE AT THE WALL THE INCREMENT OF CURVARURE THE APPLIED VERTICAL FORCE AT THE WALL (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, watch 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 = CYP1RP1) CURVATURE OF THE FIRST PIN (i.e., CURP1 = 1.0/RP1) c ------------------------------------------------------c C = MATERIAL PROPERTIES c double precision coef, en, r, ex, yield, coefmu c c = ARGUMENTS c c c c c c COEF EN R EX YIELD COEFMU ... .. k IN EQUATION (9a) OF REF. (1) n IN EQUATION (9a) .OF REF. (1) r IN. EQUATION (9c) OF REF. (1) YONG'S MODULUS ABOUT X-AXIS YIELD STRENGTH 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 = MEW PARAMETERS c double precision eta1(maxdim), eta2(maxdim) double precision c, spi, deta double precision sm1, sm2, h1, h2 double precision sm2nev, sm2old, sm2old2 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 va1, vb2, va1old, din double precision n, v0old, v0old2 double precision dn(maxdim) double precision totalm1, totalm2, dimn double precision dmn(maxdim) c ------------------------------------------------------c = ARGUMENTS C TOTALM1 C TOTALM2 C DIFMN THE RESULTANT MOMENT AT THE LEFT END THE RESULTANT MOMENT AT THE RIGHT END TOTALM2 TOTALM1 c ------------------------------------------------------c C = STRESSES AND STRAINS c double precision epsilon2(maxdim), sig2(maxdim) c -----------------------------------------------------c = ARGUMENTS c -----------------------------------------------------r c ============================================================= c C = COMHOMS c c ============================================================= c 131

PAGE 147

C = COMMOM BLOCK OF CHARACTERS c c common/ control1/ title, vhat, geometry, vhich, vork common/ control2/ control, checksm C = COMMON BLOCK OF INITIAL VALUES c c common/ common/ init1/sm0, hO, curveO, dcurve, fvO init2/xa, ya C = COMMON BLOCK OF GEOMETRY c common/ geom1/ cyp1, rp1, d1, xb, yb, ybold c C = COMMON BLOCK OF MATERIAL PROPERTIES c common/ materal1/ coef, en, r, ex, yield, coefmu c C = COMMON BLOCK OF PARAMETERS c c common/ common/ common/ param1/ neta, param2/ s1, param4/ equi jcount, icorint, kcount s2, check equiold C = COMMON BLOCK OF NEW PARAMETERS c common/ nevpara1/ eta1, eta2, data common/ nevpara2/ c. spi common/ nevpara3/ sm1 sm2, h1, h2 common/ nevpara4/ sm2nev, sm2old, sm2old2 common/ nevpara6/ z1, z2, zcom3a, zcom4a, zcom3, zcom4 132

PAGE 148

i_' c C = COMMON BLOCK OF TOLERANCE c common/ toler1/ sm2tol, angletol, equitol c C = COMMON BLOCK OF ANGLES AND CURVATURES c c common/ common/ common/ common/ angles1/ phi1, theta, thetao, thetab angles2/ dtheold, dtheta, theold curve1/ curv1, curv2, curvold, curvold2 curve2/ curp1 c = COMMON BLOCK OF FORCES AND MOMENTS c coumon/ force.1/ fva1, fvb2; fva1old, diffn coumon/ force2/ fn, dfn, fvOold, fv0old2 coumon/ force3/ fna1 fnb2, fra1 frb2 common/ moment1/ totalm1, totalm2, difmn coumon/ moment2/ dmn c c = COMMON BLOCK OF STRESSES AND STRAINS c common/ c c II END CF FILE c ss1/ epsilon2, sig2 133