‘

[wag

2;;-

‘ fig‘gc‘gfiag
%* 4

._ ‘ . ’_ qf'r‘zfii“ -
v ‘1 : V @%§g§;5$*=‘#=
‘ ‘ ‘ i‘t

 

Fit
‘3 v ‘2 '2
s" e“ 1
fififikf

KP”
LA
I . . . “ha 5 ~: > ‘  '; .J.‘ . '2, “f. - 4 2:.
' .Ai. " 7. ' . .t‘égfiflt"9§ . I ' -‘ ‘61 I , I rfiZp» V f ' , “fi-‘fit‘égggi E
t 7 ‘ ~. " W . . ., " ‘ h. 3f. I , “3 ' ‘ -~ '1 . 4' “' '5 " fl . .u '1 . ‘ '43 ‘r-‘
-. H ”7. ,. : .1 , . - , é-rk _ ‘ 1% . hi; . «wig, N
2 . *‘ ‘ -‘ - ., .. '1‘ . . .. . , . f, a
W _ gfiflg .' ’ ‘ ~ > .V 13? ', ' , 5w; ,‘ f“ 2.. j. Wfifi‘é‘ \
..
I <.
.8
:3
7 ‘ E ' h!
l '1 .' I I” {3.5; “ ~ ‘
‘v gl‘:%;k§§? A ‘ . . ~. ‘
1' ‘ ."iau‘u."‘f: 'i" ' ‘I‘I' ‘ A ‘ , ”1 .
“5. ~ .
«g ‘

‘A
. «at; A" I
, w 33",“
l'iaygfififi
mi ‘zizftw
% r
:ss’t‘ ;

 

n‘. W; ‘

',-""'1'5.‘T‘f m

:g‘ _-;-$§:;;Yg,, 5',
{22.15% .

J -.A-:w ..,.:v‘\..‘..- .

J3?

flak

1%. f 2 ”dirk:
{inuc ' a;
r‘

v I ‘ 31‘-
318:5

" 'Wf-‘u
. 5%

M3

 

k 2. ‘ : ' e
. L :, '~ ‘. j {'1 ,.' .* aha;
. ‘ ‘ 1 ~ ~ r’ , ~. 4w
3 , , ~‘ ‘ ‘73 1.". i1. Eb 33“.:
".5 ~. , J: : , - '» L: ‘ T ‘-
~ .. fig; «. I "22%
_ g Jig, ' g1. .v V‘JIIS'
‘. w . , * ~ -« «gm
A, 6 ”‘3‘ "I: ' a?!”
. ‘ A ‘ 1‘17 .
. . . ,’ “535.) ‘ ”1‘
4 . . _ . L ‘ 5:,» . u »‘
. > V 2' ' .. - ' ‘1

 

1
I

3'
. 4“ C~
":éh‘fibv ‘2} ,7 I
téfifirfiu"

 

 

MtCHIGAN STATE UNTEF‘SITY LIBRARIES

NW \“HWWHH um NM W

i
H Hi, i
3 1293 01399 4417

\
LIBRARY
Mlchlgan $tate
University

 

 

 

THESS

 

 

 

 

This is to certify that the
thesis entitled
FINITE ELEMENT SIMULATION OF SWAGING IN
THIN-WALLED CYLINDERS

presented by

Eric James Leaman

has been accepted towards fulfillment
of the requirements for

Master's degree in Mechanics

 

 

Mam

Major professor

Date March 28, 1995

 

0-7639 MS U is an Affirmative Action/Equal Opportunity Institution

 

 

CE II RETURN BOX to mow-rhi- chockoutm your mood.
Motown duo.

FLA
To AVOID FtNES mum on or

   

 

 

 

 

 

 

 

W‘W

MSU |IM W WHO
p.301

FINITE ELEMENT SIMULATION OF SWAGIN G IN
THIN-WALLED CYLINDERS

BY

Eric James Leaman

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Materials Science and Mechanics

1995

ABSTRACT

FINITE ELEMENT SIMULATION OF SWAGIN G IN
THIN-WALLED CYLINDERS

By

Eric James Leaman

The sheet metal forming operation of swaging can produce wrinkles in a thin-
walled cylinder. In order to alleviate wrinkling, an understanding of the deformations and
stress state was necessary to isolate those variables controlling the formation of wrinkles.
In this investigation, a finite element model for simulating the swaging process was devel-
oped using the commercial finite element code MARC. Revealed in the analysis was the
presence of tensile circumferential stresses on the inner surface of the cylinder in the wrin-
kling region. These stresses, in conjunction with the transcending of critical circumferen-
tial strain levels were hypothesized to be responsible for wrinkling. A proposal was made
to increase the cylinder wall’s effective moment of inertia by the introduction of a circum-

ferential stiffening rib.

ACKNOWLEDGEMENTS

I would like to acknowledge the Faculty and Students of the Department of Mate-
rials Science and Mechanics for their guidance in completing this project. In particular, I
would like to thank Mr. Mike Wood and Dr. James Lucas for their contribution of the
mechanical properties, and also Mr. Yue (Jeffery) Qiu for drawing Figure 1.1 of the text. I
would like to acknowledge the financial support of the United States Can Company under
grant 61-8502. 1 would like to thank the personnel at MARC Analysis Research Corpora-
tion for their assistance in using their product. I must, also extend my gratitude to my
research advisor, Dr. Ronald C. Averill, for his guidance and support throughout this
project. Although simple words are inadequate, I must say thank you to my wife Laurie
for all of her encouragement, support, and sacrifices that she has made during the comple-

tion of this work.

Table of Contents

List of Tables vii
List of Figures x
Chapter 1 Introduction 1
1.0 Introduction 1

1.1 Problem Statement 5

1.2 Literature Review 7

1.3 The Present Study 25
Chapter 2 Modelling of the Swaging Operation 28
2.0 Introduction 28

2.1 Geometric Model of Cylinder and Forming Tools 29

2.1.1 Preliminary Element Selection 33

2.1.2 Evaluation of Element Performance for a Linear Problem 35

2.1.3 Discussion of Results for a Linear Problem 37

2.1.4 Nonlinear Evaluation of Element Performance 43

2.2 Boundary Conditions and Contact Constraints 45

2.2.1 Boundary Conditions 45

2.2.2 Initial Stresses 47

2.2.3 Contact Conditions 50

iv

2.3 Material Model
2.3.1 Mechanical Properties
2.3.2 Plasticity Model
2.4 Solution Parameters and Analysis Options
2.4.1 Solution Method
2.4.2 Analysis Options
Chapter 3 Numerical Simulation of Swaging
3.0 Introduction
3.1 Simulation of the Swaging Process
3.2 Final Element Selection
3.2.1 Forming Analysis Results Using Element 1
3.2.2 Forming Analysis Results Using Element 116
3.3 Convergence of Forming Problem
3.3.1 Mesh Refinement Study
3.3.2 Discussion of Mesh Convergence Study Data

3.3.3 Comparison of Stress Trends Resulting from Mesh
Refinement

3.3.4 Convergence Study on Number of Load Increments

3.4 Results from Forming Simulation

3.4.1 Strain and Stress States in the Critical Region

3.4.2 Assessment of Strain Hardening and Initial Stress
Effects on Strain and Stress States

Chapter 4 Discussion of Wrinkling

4.0 Introduction

52

52

54

57

57

59

60

60

61

66

67

67

69

71

73

81

91

91

91

93

102

102

4.1 Existence of Critical Strain
4.2 Bending Effects
4.3 Introduction of Stiffening Rib
Chapter 5 Summary and Future Work
5.1 Summary
5.2 Future Work
APPENDIX
A] Mesh Refinement Study Data
A.2 Analysis Results without Strain Hardening

REFERENCES

vi

103

106

108

122

122

123

128

137

142

Table 2.1

Table 2.2

Table 2.3

Table 2.4

Table 2.5

Table 2.6

Table 2.7

Table 3.1

Table 3.2

List of Tables

MARC Elements

Predicted maximum deflection due to a circumferential line

load in a simply supported cylinder of radius=l.0, thickness=0.l,
and length=100.0. All values are normalized to the closed form
solution.

Predicted maximum deflection due to a circumferential line
load in a simply supported cylinder of radius=5.0, thick-
ness=0.05, and length=300.0. All values are normalized to the

closed form solution.

Predicted maximum deflection due to a circumferential line load in a
simply supported cylinder of radius=100.0, thickness=l.0, and
length=200.0. All values are normalized to the closed form solu-
tion.

Nonlinear test case results corresponding to the predicted
maximum deflection due to a circumferential line load in a simply
supported cylinder of radius: 100.0, thickness=1.0, and
length=100.0. All values are normalized to the closed form solu-
tion.

Mechanical Properties of Double-Reduced Tinplate Steel, Type L.

Stress versus Plastic Strain for Double-Reduced Tinplate Sheet
Steel, Type L.

Percentage Change in Axial Strain Values in the Critical Region.
(%)

Percentage Change in Circumferential Strain Values in the Critical
Region. (%)

vii

34

38

39

42

53

56

76

76

Table 3.3

Table 3.4

Table 3.5

Table 3.6

Table A.1

Table A2

Table A3

Table A4

Table A5

Table A6

Table A7

Table A8

Table A9

Table A. 10

Table A. 11

Stress Values for 250 vs. 500 incs for the First Stage in the
Critical Region.

Stress Values for 250 vs. 500 incs for the Second Stage in the
Critical Region.

Effects of Strain Hardening on the Global Minimum and Maximum
Strain and Stress Values.

Effects of Initial Stresses on the Global Minimum and Maximum
Strain and Stress Values.

Axial Strain on the Outer Surface in the First Bending Region.
(x10’3)

Axial Strain on the Mid-Plane Surface in the First Bending Region.
(x10'3)

Axial Strain on the Inner Surface in the First Bending Region.
(x10’3)

Circumferential Strain on the Outer Surface in the First Bending
Region. (x10‘3)

Circumferential Strain on the Mid-Plane Surface in the First Bending
Region. (x10’3)

Circumferential Strain on the Inner Surface in the First Bending
Region. (x10’3)

Axial Stress on the Outer Surface in the First Bending Region.
(szi)

Axial Stress on the Mid-Plane Surface in the First Bending
Region. (szi)

Axial Stress on the Inner Surface in the First Bending Region.
(szi)

Circumferential Stress on the Outer Surface in the First Bending
Region. (szi)

Circumferential Stress on the Mid-Plane Surface in the First
Bending Region. (szi)

viii

92

92

94

95

128

128

129

129

129

130

130

130

131

131

131

Table A. 12

Table A. 13

Table A. 14

Table A. 15

Table A. 16

Table A. 17

Table A. 18

Table A. 19

Table A.20

Table A.21

Table A.22

Table A.23

Table A.24

Circumferential Stress on the Inner Surface in the First Bending
Region. (szi)

Axial Strains on the Outer Surface in the Critical Region. (x10'3)

Axial Strains on the Mid-Plane Surface in the Critical Region.
(x10'3)

Axial Strains on the Inner Surface in the Critical Region. (x10’3)

Circumferential Strains on the Outer Surface in the Critical
Region. (x10‘3)

Hoop Strains on the Mid-Plane Surface in the Critical Region. (x10’3)
Hoop Strains on the Inner Surface in the Critical Region. (x10—3)
Axial Stresses on the Outer Surface in the Critical Region. (szi)

Axial Stresses on the Mid-Plane Surface in the Critical Region.
(szi)

Axial Stresses on the Inner Surface in the Critical Region. (szi)

Circumferential Stresses on the Outer Surface in the Critical
Region. (szi)

Circumferential Stresses on the Mid-Plane Surface in the Critical
Region. (szi)

Circumferential Stresses on the Inner Surface in the Critical
Region. (szi)

132

132

132

133

133

133

134

134

134

135

135

135

136

Figure 1.1
Figure 2.1

Figure 2.2

Figure 2.3
Figure 2.4

Figure 2.5

Figure 2.6

Figure 2.7

List of Figures

Wrinkling Behavior Exhibited in Swaged Cylinder.

Relationship Between Modelled Surface and Full Cylinder.

Geometry of Forming Tools with Reference to the Cylinder Axis.

(Not to scale)
MENTAT Model of Cylinder and Forming Tools.
Geometric Model with Applied Boundary Conditions.

Cross Section of Cylinder with Stress Resulting from
Pure Bending.

Stress versus Plastic Strain for Double Reduced Tinplate
Sheet Steel.

Newton-Raphson Solution Method.

Figure 3.1(a) First Stage of Swaging at the Zero‘h Increment.

Figure 3.1(b) First Stage of Swaging at the 50th Increment.

Figure 3.1(c) First Stage of Swaging at the 100th Increment.

Figure 3.1(d) First Stage of Swaging at the 250th Increment.

Figure 3.1(e) Release of the First Stage Tool.

Figure 3.1(0 Tool Moved into Position in First Increment of Second Stage.

Figure 3.1(g) Second Stage of Swaging at the 352nd Increment.

Figure 3.1(h) Final Increment of the Second Stage.

30

31

32

46

48

56

58

62

62

63

63

65

65

Figure 3.1(i) Release of Second Stage Tool.

Figure 3.2
Figure 3.3
Figure 3.4

Figure 3.5

Figure 3.6

Figure 3.7

Figure 3.8

Figure 3.9

Figure 3.10

Figure 3.11

Figure 3.12

Figure 3.13

Figure 3.14

Figure 3.15

Figure 3.16

Figure 3.17

Figure 3.18

Swaging Analysis Using Element 116 with a 150x2 Mesh.

Sampling Locations for Convergence Study of Forming Analysis.

Axial Strain on the Outer Surface in the Critical Region vs. Mesh Size

Circumferential Stress on the Outer Surface in the Critical Region vs.
Mesh Size.

Axial Strain on the Inner Surface in the Critical Region vs. Mesh Size.

Circumferential Strain on the Inner Surface in the Critical Region vs.
Mesh Size.

Axial Stress on the Outer Surface in the Critical Region vs. Mesh Size.

Circumferential Stress on the Outer Surface in the Critical Region
vs. Mesh Size.

Axial Stress on the Inner Surface in the Critical Region vs. Mesh Size.

Circumferential Stress on Inner Surface in the Critical Region vs.
Mesh Size.

Computational Time vs. Mesh Size.

Axial Stress in the Critical Region for the Second Stage of the
150x8 Mesh.

Circumferential Stress in the Critical Region for the Second Stage
of the 150x8 Mesh.

Axial Stress in the Critical Region for the Second Stage of the
300x8 Mesh.

Circumferential Stress in the Critical Region for the Second Stage of
the 300x8 Mesh.

Axial Stress in the Critical Region for the Second Stage of the
150x2 Mesh.

Circumferential Stress in the Critical Region for the Second Stage of
the 150x2 Mesh

xi

66

68

74

74

75

75

78

78

79

79

80

83

84

85

86

87

Figure 3.19
Figure 3.20
Figure 3.21

Figure 3.22
Figure 3.23
Figure 3.24
Figure 3.25

Figure 3.26

Figure 3.27

Figure 3.28

Figure 4.1
Figure 4.2
Figure 4.3
Figure 4.4
Figure 4.5
Figure 4.6
Figure 4.7

Figure 4.8

Figure 4.9

Figure 4.10

Deformations of the 800x2 Mesh at Increment 303.
Deformations of the 150x2 Mesh at Increment 303.
Deformations for the 800x2 Mesh at Increment 373.

Deformations for the 150x2 Mesh at Increment 373.
Axial Strains in the Critical Region at the 503rd Increment.
Circumferential Strains in the Critical Region at the 503rd Increment.

Axial Stresses in the Critical Region at the 503rd Increment.

Circumferential Stresses in the Critical Region at the 503rd Increment.

Axial Stresses in Critical Region at 503rd Increment, Calculated
without Initial Stresses.

Circumferential Stresses in the Critical Region at the 503rd
Increment, Calculated without Initial Stresses.

Euler Buckling Column.

Shell Element Subject to Bending and Resultant Moments.
Circumferential Stress versus Strains for the First Stage of Forming.
Simple Beam Subjected to a Bending Moment.

Plate Subjected to a Bending Moment.

Initial Geometry of Cylinder with Added Stiffening Rib.

Formed Shape of the Cylinder with Added Stiffening Rib.

Final Shape of the Cylinder with Added Stiffening Rib after
Elastic Snapback.

Axial Strains in the Critical Region of the Cylinder with Added
Stiffening Rib.

Circumferential Strains in the Critical Region of the Cylinder with
Added Stiffening Rib.

xii

89

89

90

9O

96

97

98

99

100

101

104

104

109

114

114

115

116

117

118

119

Figure 4.11

Figure 4.12

Figure A.1
Figure A.2
Figure A.3

Figure A.4

Axial Stresses in the Critical Region for the Cylinder with
Added Stiffening Rib.

Circumferential Stresses in the Critical Region for the Cylinder
with Added Stiffening Rib.

Axial Strains in the Critical Region Without Work Hardening.
Circumferential Strains in the Critical Region Without Work
Hardening.

Axial Stresses in the Critical Region Without Work Hardening.

Circumferential Stresses in the Critical Region Without Work
Hardening.

xiii

120

121

138

139

140

141

Chapter 1

Introduction

1.0 Introduction

The current thrust in the manufacturing community is toward the conservation of
resources and energy, while at the same time maintaining product quality and responding
to the demands of the consumer. This philosophy has its roots in the latter days of the
stone age and has evolved over the next six thousand years to its present state. The pri-
mary concern of early man was simply to obtain food to eat. The tools necessary to kill
were fashioned by chipping flint into weapons. Eventually flint tools began to take on spe-
cific roles, such as axes, knives, or borers. Toward the end of the Stone Age man discov-
ered the ease with which some materials could be worked into shape and metal forming

was born. A couple of millennia later casting was found to make the working of metals

2

more efficient by eliminating some of the hammering necessary, but, because the castings
were of copper, a soft material, it was still advantageous to use stone. Bronze improved
the utility of cast tools, but shortcomings were still present. For instance, bronze Chisels
could not cut stone. Thus a new material was needed and iron was found to satisfy those
needs, giving us the iron age (Parr, 1967). In each of these cases, societal needs for
improved tools and weapons required new materials and methods of production to be
developed. Many discoveries were made by pure luck and others were the result of invest-
ments in manpower, raw material, and time. The experiences of those working in metal
forming built up over time and eventually a group of artisans was established.

The needs of the United States government in the late 18th century made it neces-
sary to obtain a large stand of arms. The time and the gunsmiths necessary to fill the need
was scarce. Eli Whitney was driven by these two circumstances to develop a new system
of manufacturing whereby he would “form the tools so the tools themselves shall fashion
the work and give to every part its just proportion - which when accomplished will give
expedition, uniformity, and exactness to the whole” (Green, 1956). The process Whitney
was proposing was currently being used to produce goods, but not on the scale of com-
plexity that the musket required. The lock on a musket required an exact fit of the parts in
order to work. This being the case, the standard of the day called for crude parts to be
made and then assembled by a craftsman who filed and fitted the parts together to produce
the working lock, a very time consuming task. Whitney’s plan would eliminate the need
for filing, and introduce the system of interchangeable parts. This evolution in manufac-
turing was consistent with the goals of speeding up the manufacturing process, reducing

the labor necessary, and improving the quality in musket production. But the system still

3

required craftsmen to assemble the parts and tools needed to do the work. These two diffi-
culties prevented the delivery of the arms on time, but the system was successful nonethe-
less.

Finally, Henry Ford brought mass production to its present form. He recognized
the need for a low cost, high quality vehicle that almost anyone could afford. To do so
required the marriage of five components: precision, standardization, interchangeability,
synchronization, and continuity. While all of these existed individually, Ford brought them
all together in the moving assembly line where the work moved to the worker who per-
formed a single operation. Production costs were once again driven downward by this new
manufacturing concept and the joys of driving were delivered to a wider audience. Thus
the advent of mass production had arrived on the scene, with its roots stemming from
Whitney’s system of interchangeability. However, elaborate planning, expensive tooling,
and exact synchronization of the parts moving through the plant were still required.

While each of these evolutionary steps has brought about greater efficiency and
reduced cost in manufacturing, societal requirements are for even greater cost reductions,
quicker product development times, and environmental impact awareness. To accomplish
these mandates manufacturers have turned to a concept known as simultaneous or concur-
rent engineering. This represents a fundamental change in the manner in which products
are designed and manufactured. Previously a sequential process was followed. A need was
defined, a part designed to meet the need, a means of manufacturing the part devised, and
the process was tested. Rarely, though, is the cycle correct the first time. One or more of
the components in the design cycle may have to be altered, thus forcing a retrial of the

process. Trial and error is expensive and not consistent with the current thrust of conserva-

tion.

Concurrent engineering is “a more co-operative, simultaneous, computer-aided”
approach to product development (Butman, 1991). Design and manufacturing engineers
are voicing their concerns to each other at the initial stages. This parallel approach gives
rise to a more efficient design cycle, but methods of analyzing the integrity of the part and
the manufacturing process are also essential. Although trial and error methods served this
function in the past, the decline of skilled tradesmen and the need for rapid results render
this scheme insufficient. What is needed is' a means to do the analysis without the physical
product or process. Numerical methods are readily adept at performing this function.

Numerical methods consist of building mathematical models that are representa-
tive of the product and tools of the process. In this way testing is done on the computer
rather than on the plant floor.The finite element method is one such numerical scheme
available for the analysis of continuous structures that are currently being manufactured.
The approach of the method is to model the structure by discretizing the geometry with a
group of elements. The equations that govern the behavior of each of these individual ele-

ments are then assembled to produce a system of algebraic equations of the form
Ki = F (1.1)

where E is the stiffness matrix of the entire structure and is a function of the elemental
deformation modes and the material properties, 35 is a vector of the displacements and rota-
tions (or the degrees of freedom, DOF) that describe the deformations of the elements, and
[—7 is a vector of the forces at discrete points within each element. By substituting the

applied loading along with the DOE which are fixed on the boundary into Equation 1.1,

5

the unknown displacements, rotations, and forces may be solved for in each element and
thus the global response of the structure may be predicted. The results obtained, however,
are rarely exact because of the approximations necessary to manipulate the governing
equations of the elements into the form of Equation 1.1. Even so, solutions accurate for
engineering purposes can usually be obtained, though the large number of equations nec-
essary to obtain theses results may demand the use of large scale computing machinery.
In the above paragraph it was assumed that the stiffness matrix, 7?, and the force
vector, 7‘, were independent of the displacement vector, 3. In many practical situations this
is not true and the analysis becomes nonlinear. Such is the case when analyzing a manu-
facturing process. One such operation is sheet metal forming where a set of rigid tools
introduces plastic deformations into a relatively simple geometry to produce the more
complex configuration of the final product. The analysis of such an operation is nonlinear
because of the large deformations required, the relationship between stress and strain, and
the contact of the rigid tools and the blank. Not only does this add to the complexity of the
finite element analysis, but it also renders nearly impossible any analytical solution.
Therefore, the means to analyze sheet metal forming are either the trial and error method,
which has been shown to be costly and contrary to the requirements of simultaneous engi-
neering, or the development of computer models where part and process evaluation may

be done for relatively little cost on the computer.

1.1 Problem Statement

A typical example of a metal forming operation is swaging or necking wherein the

diameter of a cylinder is reduced by the action of a set of rigid dies. When using this pro-

 

 

 

 

 

 

 

 

 

 

Figure 1.1 Wrinkling Behavior Exhibited in Swaged Cylinder.

cess on thin-walled cylinders the possibility then exists for the formation of wrinkles in
the necked-in region, as shown in Figure 1.1. The wrinkles are the result of a build up of
compressive circumferential stresses in the cylinder. When these stresses exceed some
critical value, any disturbance, whether due to material imperfections or the manufactur-
ing process, will cause the cylinder to buckle (or bifurcate) to an adjacent equilibrium
position (the wrinkled configuration). One possible solution to the problem is to increase
the thickness of the cylinder wall. However, arbitrary use of thicker material increases the
cost associated with the purchasing of raw materials which in the end is passed on to the
consumer or absorbed by the manufacturer, neither of which is desirable. In addition,
other solutions may exist that will eliminate the misuse of resources. In other words, the
increased wall thickness places greater demands on the environment both in terms of raw

materials for steel production and waste management at the end of product life. However

7

an investigation of the stress state and kinematics of the cylinder during the swaging pro-
cess may lead to an alternate solution that is consistent with the goal of conserving
resources and energy. Trial and error techniques coupled with the die makers experience
will not produce the desired information, and this also inherently conflicts with the current
thrust in manufacturing as already explained. No analytical means of evaluating the pro-
cess is available. The nonlinearities due to large deflections, plasticity, and contact
between the cylinder and tools associated with the. problem render the solution analyti-
cally intractable. Thus, a finite element simulation of the necking process is sought in
order to achieve the stated objectives. With the stress states and deformations of the cylin-
der in hand, changes, whether in the material, the forming path, and/or the final geometry,

with an engineering basis may be proposed to eliminate the wrinkling phenomena.

1.2 Literature Review

Early attempts to predict the success of a forming operation centered on the use of
forming limit diagrams (FLD). The formation of the FLD is through experimental proce-
dures wherein the sheet blank is etched with a grid pattern of circles. The blank is then
stretched over an unlubricated punch and the deformation of the circles in the region of
failure are observed. The shape of the etchings will have changed from circular to ellipti-
cal. The major strain is calculated by measuring the major axis of the ellipse and the minor
strain from the minor axis. After a series of tests the delineation between a successful for-
mation and one that fails is established (Kalpakjian, 1991). While experimental construc-
tion of FLD’s has been successful, numerical prediction of the curves requires a failure

criterion. Toh (1989) proposed a method of predicting the forming limit curves of sheet

8

materials using an incremental rigid-plastic finite element method and a new limit crite-
rion. In the analysis of a hemispherical punch stretching a circular blank with various cir-
cular cutoffs and various friction conditions at the tool-sheet interface, the material was
assumed to obey Hill’s anisotropic yield criterion and its associated flow rule. The result-
ing set of load versus displacement curves were then used to derive a critical slope condi-
tion which could subsequently be used as a criterion to determine the corresponding
critical major and minor strains. This approach has been shown to give reasonable predic-
tion of the FLD’s, but this is only a method of predicting the failure of the material. To
understand the underlying behavior of a structure, simulation of the forming process is
necessary.

Simulation requires a sophisticated finite element formulation, especially concem-
ing those components giving rise to non-linearities. These include the large deformations
that are usually introduced to the workpiece. The changes in the configuration of the
workpiece should be permanent and as such the material experiences plastic deformations,
another source of nonlinearity. Finally, the deformations are produced by a tool or set of
tools that come into contact with the workpiece, thus introducing nonlinear boundary con-
ditions into the analysis. While many considerations are necessary in the analysis of metal
forming operations, these components draw the most attention.

Keck, et a1. (1990) summarized results obtained for various sheet metal forming
operations that were simulated using an implicit finite element code. In general, a variety
of element types are used, but all are formulated based on the principle of virtual displace-
ments and an updated Lagrangian description for the incremental calculation of the defor-

mation process.

9

Use was made of the rate form of the elastic-plastic material model. However with
the assumption of small strains, a hypoelastic constitutive law may be implied, because
the strains may be decomposed into their elastic and plastic components. This required a
yield function for the material which, in general, was assumed to conform to von Mises
yield criterion. But in the case of a membrane element formulation, Hill’s quadratic yield
criterion was used, thereby taking into account the normal anisotropy. The rate form
requires the integration of the constitutive relationship which was performed using either a
radial return method or a return mapping algorithm.

The treatment of the changing contact conditions of the problem were handled in
the following manner: The dies were treated as rigid bodies and in order to fulfill the con-
tact conditions a penalty formulation was adopted. Penetration of the tool was checked
only at nodes that were specified as contact nodes. Then a penalty factor was calculated
for each possible contact node from the stiffness of the workpiece perpendicular to the
tool surface and scaled by a user defined factor. The resulting contact friction was mod-
elled with a modified Coulomb friction model. This overcame the problem of singularities
arising in the first derivative of the Coulomb friction model at vanishing relative motions
and enables the linearization of the friction law in connection with the Newton-Raphson
iteration scheme.

The equilibrium of each incremental loading of the structure was obtained by a
dynamic solution scheme. In general, the N ewton-Raphson iteration method was used, but
in cases where the rate of convergence reached some predetermined value, the so called
BFGS method (named for Broyden Fletcher, Goldfarb, and Shanno, see for example,
Zienkiewicz, 1989) was employed. The Newton-Raphson method was reintroduced upon

a deceleration of the convergence rate.

10
The finite element formulation proposed was applied to the analysis of hydrostatic

bulging of a circular sheet of annealed aluminum, hydrostatic bulging of rectangular
sheets of anisotropic mild sheet steel, an axisymmetric stretch forming process, and the
deep drawing of a rectangular cup. In each case the results showed good agreement with
previous experimental work. The authors point out, however, that the formulation pre-
sented is limited to materials which have accurately identified constitutive models and
material parameters. Emphasis is thus placed on the accurate description of the material
behavior through experimental techniques before proceeding with more complex analyses
than those performed.

A more complicated analysis was conducted by Lee, et a1. (1991) on a stretch-
drawing process particularly used in the automotive industry. The aim was to develop a
finite element method which resulted in acceptable computational times, while retaining
the ability to model the key features of a metal forming analysis. Chief among these being
the contact algorithm which is generally regarded as a major issue for numerical stability.
As such the authors implemented a membrane line element consistent with a shell formu-
lation for the contact formulation. Additionally, to aid in the reduction of computational
time, the number of degrees of freedom within the shell elements was reduced.

The model was formulated using an updated Lagrangian description of motion. In
addition, the shells used to describe the sheet metal were assumed to be thin and therefore,
the through the thickness shear deformations could be neglected. “With a full bending the-
ory being considered, the strain expressions contained the second derivative of the incre-
mental in-plane displacements. The degrees of freedom therefore consisted of the change

in axial displacement, the change in transverse deflection, and their first derivatives. But,

11

by assuming a linear variation of the incremental axial displacement, its derivative was no
longer required, resulting in a net reduction in the number of DOF in the model. There-
fore, the interpolation functions used in the finite element formulation were linear for
inplane displacements and Hermite cubic for the transverse deflection. These interpolation
functions were then used in the principle of virtual work, which when linearized along
with a geometric constraint for tool contact resulted in the stiffness matrices and force
vectors due to inertia, contact forces, and internal reactions. When evaluating the resulting
algebraic equations the through the thickness integration was reduced to single point
quadrature rule in order to reduce computational time, while retaining the key features of
the bending model.

In regard to the material model of the formulation, the response was assumed to be
rigid visco-plastic so that incompressibility could be applied using small natural strains. In
addition, normal anisotropy was taken into account through the use of Hill’s yield crite-

rion. Finally, a power law hardening relation with strain rate sensitivity of the form

a = oo+K(a0+eo+A£)"(é/y)m (1.2)

was implemented, where 00 is the yield stress, K is the material strength parameter, a0 is
the pre-strain, n is the strain hardening exponent, m is the strain rate sensitivity index, and
Y is the base strain rate. Using this expression a wider variety of sheet metals could be
considered.

Finally, the detection of contact between the tool and the workpiece was indicated

by satisfaction of the following,

12
gnEn-( —x3) = 0 , (13)

where n is the normal between the tool and the workpiece and x" and .15s are position vec—

tors of a point on the sheet and it nearest point on the tool. However, because 1:" and x‘ are
unknowns the equation must be linearized and solved iteratively. This was a general gap
condition. In the event that the tool moved only in the vertical direction, a vertical gap
parameter would be used. In either case, when the tool and sheet metal were in contact
external forces would be transmitted through the tool to the sheet metal and therefore a
consistent contact formulation was required.

The methods used above show good agreement with other numerical results as
well as with experimental data. However, it was found that when the strain reaches the n
value of the strain hardening equation (1.2), the model failed to converge. This difficulty
was due to a loss of definiteness in the stiffness matrix and is termed softening.

Rebolo, et a1. (1990) are members of a team developing a general robust finite ele-
ment code for the analysis of sheet metal forming problems. As such the techniques which
they employ have been proven to be effective in the solution of these types of analyses.
The focus of the discussion was again toward the element types used, the modelling of the
contact condition, and the frictional model used in the metal forming analyses. Little con-
sideration was given to the constitutive model as it appeared to the authors that the con-
ventional elasto—plastic models together with small strain linear elasticity, a smooth flow

surface, and isotropic hardening provided adequate representation of the material behavior

13

during forming.

In discussing the types of elements which are available, the authors point out that
the use of continuum elements is prohibitively expensive. These computational difficulties
can be overcome by the use of structural elements, such as membrane or shell elements,
but these elements introduce the inability to handle two-sided contact due to their lack of
ability to accommodate thickness stresses. In contrast, the introduction of assumed strain
and stress elements allow for the advantages of structural elements, while at the same time
providing for two-sided contact.

The formulation presented was based on a membrane element which employed
finite strain plasticity. Generally this type of formulation would require the introduction of
an additional rotation tensor to account for the stresses and strains. However, a more direct
method was offered by calculating the strain measures from the polar decomposition of
the incremental deformation gradient, the incremental rotation vector, and the incremental
plastic strain. The stress was then a function of the elastic strain. The advantage of this for-
mulation was the ease in which it may be implemented in an updated Lagrangian frame-
work.

A simplified approach to the detection of contact and the use Lagrange multiplier
techniques to impose the no penetration constraint were employed. A measure of overclo-
sure, quite similar to that proposed by Lee, et al., was formulated by calculating the dis-
tance along the normal from the closest point on the rigid surface to the contact node on
the deformable body. When the overclosure is greater than or equal to zero, contact is
detected and a Lagrange multiplier was introduced which provided the contact pressure.

Finally, when contact was identified a Coulomb friction model was employed to

14

account for the sticking or slipping stresses generated within the model. Again this was
done by introducing Lagrange multipliers. In this situation, however, the parameters could
take on different roles. In the event that sheet metal and the die stuck together the
Lagrange multiplier was used to enforce the constraint that the shear strain was zero. Con-
versely for a sliding condition, the Lagrangian took on the value of the frictional stress.

Within each of these formulations the technique used for the detection of contact
was carried out in similar fashions. A normal vectorfrom the tool(s) to the workpiece was
constructed and when the magnitude of the vector goes to zero contact was assumed. Once
contact has been realized, a constraint must be added to the governing equations to pre-
vent penetration of the tool into the material. In general this constraint was enforced in one
of two manners, either a Lagrange multiplier was used or a penalty function parameter
was used. The Lagrange multiplier enforces the contact condition exactly, but at the
expense of additional DOF in the solution. Also, a singular stiffness matrix can result from
this formulation. Both of these difficulties can be avoided by using the penalty function,
but satisfaction of the contact condition is not exact and the solution is highly sensitive to
the choice of the penalty factor. In each case the contact constraint is enforced in a “point-
wise” sense.

A method of enforcing the contact constraint in an average sense over the element
boundary was proposed by Simo and Taylor (1985). This work was conducted within the
context of bilinear isoparametric interpolation of the displacement field, thereby allowing
the assumption that the contact pressure was constant over each contact segment, and dis-
continuous across segments. The contact constraint was therefore enforced in an average

sense over the contact segment, rather than in a pointwise manner. This resulted in the

15

average gap being a kinematic variable.

The introduction of the average contact constraint was accomplished by using a
perturbed Lagrangian formulation of the total potential energy functional. By defining the
gap function in the usual manner, a modified total potential energy functional for the con-

strained system may be written as
1 2 2 A A 1 2
”8(9 ,LI ,A) = A2111 (9 )‘i’Ag-z—EA (1.4)

where the first term on the right hand side is the potential energy, the second term is the
contribution of the Lagrangian constraint, and the last term has the form of a penalty term.
As this penalty term goes to infinity, 8 -* co, the standard Lagrangian formulation is recov-
ered. The real value in adding the additional ‘penalty term’ to the functional was that the
possibility of a non-positive definite stiffness matrix was eliminated. The perturbed
Lagrangian functional could then, after taking the first variation, be used to formulate a
finite element solution to the contact problem.

Yet the perturbed Lagrangian formulation, as with the other classical approaches,
is not without faults. In each method only those constraints that were active during the
iteration contributed to the incremental equations. Undesirable consequences may arise
from this approach. For example, the quadratic convergence of the Newton-Raphson
method may be seriously affected if there are frequent changes in the active set.

By introducing a method in which all the inequality constraints are enforced
directly by means of additional equations together with the equilibrium equations to form

a system of nonlinear equations in the displacements and the contact tractions, the New-

16

ton-Raphson method may be used to overcome these difficulties. Just such an algorithm
was introduced by Eterovic and Bathe (1991). The goal was to derive a set of finite ele-
ment equations that would explicitly account for the presence of the contact conditions.
This was done by formulating the virtual work for a continuous body in a total Lagrangian
fashion. To this statement was added the virtual work done by the contact tractions over
the virtual relative displacements between the contacting bodies. A second governing
equation was derived from the contact conditions. Thus, using the two equations in con-
junction resulted in a complete variational formulation of the motion. These equations
were then cast into a finite element model which contained n+4m unknowns, where n is
the number of nodal dof and m is the number of contact nodes. But the elegance of the
model was that all the inequality constraints associated with the contact conditions were
contained within the finite element equations. But this gives rise to a highly nonlinear set
of equations, and therefore a robust solution technique was paramount.

Riccobono (1992) presented a comparison of two methods which are used to ban-
dle the nonlinearities presented by the plasticity model in metal forming problems. The
first was a matrix method proposed by Prager and Hodge in which the actual strain rate

field of the analysis was the one that satisfied the following:

2 1. .
J = 7—560! /§eijeinV— ITividS (1.5)
v s,

where S is the surface and T,- are components of the surface traction. In this formulation,

the incompressibility condition is enforced using Lagrange multipliers. Thus the problem

17

is defined as finding the minimum to a functional which combines (1.5) and the
Lagrangian constraint of the incompressibility condition.

The second method used the infinitesimal strain increment form of (1.5). Lineariz-
ing the yield surface by a set of tangent planes and using Gaussian numerical integration,

the functional was written as,

4 n, n,
J = Jic 2 2115.5]. 2 (1).,- 2 Fidui (1.6)

j=I i=1 i=1

where for the axisymmetric state, 21cerJ. is the volume associated with the integration

point j, and u,- the number of the nodes on which the forces F ,- are applied. Thus, using this

new functional along with two sets of restraints (one which states that the infinitesimal
strain increments derived by the velocity field equals the one derived by imposing the nor-
mality law at each integration point and the other stating that the die never penetrates the
blank) the linearization of the yield function is complete.

A comparison of the two methods revealed that the first formulation is based on an
iterative process and uses a non-linear system of equations, while, on the other hand, the
second formulation is given in a closed form and thus eliminates any chance of conver-
gence problems. This also meant that the matrix method may allow the tool to penetrate
the blank during a solution step, but the linearization technique enforces the non-penetra-
tion restraint explicitly. Therefore the linearization method may provide superior perfor-
mance to the matrix method, but the increased number of unknowns associated with the
formulation may indeed cancel any advantages gained.

Finite element models are usually formulated for fully three dimensional or possi-

18

bly two dimensional spaces when assumptions of plane stress or plane strain are made.
This may lead to a prohibitively expensive analysis in the case of a nonlinear metal form-
ing problem. In cases when the geometry, properties, and loading are constant in one
direction, for example an axisymmetric problem, the model may be simplified by model-
ling the mechanical fields in that coordinate with a traditional, trigonometric Fourier
series. In general this method fails to satisfy the anisotropic form of the governing equa-
tions for laminated composite structures. Padovan (1974) developed a quasi-analytic finite
element procedure for axisymmetric anisotropic shells and solids by using the exact func-
tional representation for the circumferential variable dependency.

In shell theory for anisotropic materials the appearance of the shearing terms in the
stiffness matrix cause the traditional, trigonometric Fourier series approach to break down
when trying to obtain a solution for the equilibrium equations. One possible resolution is
to transform the equilibrium equations into an orthotropic form, but the analytical solution
is intractable due to the distortion of the boundary conditions. Therefore, it was proposed
to cast the displacement vector into a complex Fourier series expansion. Applying the new
formulation for the displacement and force vectors in the equilibrium equations yielded a
complex differential equation with the real and imaginary parts of the displacement vector
being coupled for the fully anisotrOpic case. The difficulty in obtaining an analytical solu-
tion lead to the development of a finite element solution. The formulation of the finite ele-
ment equations was the same as a traditional approach with the exception that the
displacement and force vectors were in the complex, Fourier expansion form. This lead to
a complex stiffness matrix for solution.

The procedure outlined was tested against problems with known solutions and it

was found that the complex form of the equations gave very good approximations in all

19

cases. Given the accuracy of the results and the difficulties in using the traditional Fourier
series, the complex Fourier series representation appeared to have some advantages, espe-
cially in that the significant effects of material anisotropy were revealed by its use.

Two constraints on the use of the semi-analytic, or finite strip method are that the
geometry of solids of revolution be axisymmetric and that the material properties in the
circumferential direction be constant. To broaden the scope of the method Sedaghat and
Hermann (1983) proposed relaxing the first constraint and removing the second altogether.
The formulation started with the removal of the circumferential or theta dependence of the
total potential energy of the system. This was done by expanding the displacements and
body forces with a trigonometric Fourier expansion. However the elastic coefficient
matrix may also be theta dependent. This difficulty was handled by decomposing this
matrix into two components called the base properties and the deviation properties respec-
tively. The base properties are theta independent and are in most cases found by averaging
the material properties in theta coordinate. The deviation properties are simply the devia-
tion from the base values.

From the decomposed elastic matrix the stress-strain relationship was in the form
of

(IN = De”+be”” (1.7)

where the first term was independent of theta and yielded the usual uncoupled equations
and the second term was theta dependent and contributed to the load matrix in a manner

similar to initial stresses. The finite element model could then be constructed from the

20
total potential energy of the system.

The variation of the material properties aided in the formulation of the semi-ana-
lytic method for non-axisymmetric geometries. By idealizing the geometry and then vary-
ing the material properties a semi-analytic formulation could be brought about. For
example a semi-infinite plate with a hole loaded by internal pressure was idealized by a
circular plate of radius p, where p was chosen such that the value of stress in the x direc-
tion was small. If the hole is near the top edge of the plate, the edge effects would be
accounted for by varying the properties within the idealized plate. Results from a finite
element model formulated using the idealized geometry and the varying material proper-
ties were in good agreement with the closed form solution for the problem. Thus it appears

that the scope of the Fourier series formulation had been broadened.

Attention will now be turned to analyses involving the formation of wrinkles in
sheet metals. Currently the most popular method of predicting wrinkle failures is the
Yoshida Buckling test. However, this method, as well as other empirical techniques, have
proved to be inadequate for observed trends. In particular the local curvature of sheet
metal during forming has been shown to have a significant effect on the conditions for
wrinkling, but the empirical methods fail to account for this parameter.

Neale and Tugcu ( 1990) concentrated on some of the basic features of wrinkle for-
mation, in particular those occurring within the context of the plastic buckling theory for
thin plates and shells. Wrinkling can be viewed as a plastic buckling process in which a
wavelength of the buckles in one direction is very short. This buckling is local and depen-
dent on the local curvature, the thickness of the sheet, on the material properties, and the

local stress state.

21

The basic element considered had constant radii of curvature R1 and R2 and con-
stant thickness t, the state of stress was assumed to be a uniform membrane state (no bend-
ing), and the investigation was limited to those regions of the sheet that were not in
contact with the die. Simplifications arose by exploiting the fact that the short-wavelength
wrinkling modes were shallow and could be analyzed by the Donnell—Mushtari-Vlasov
(DMV) shallow shell theory, which restricted the analysis to modes where the characteris-
tic wavelength was large in comparison to the thickness, but small compared to the local
radii of curvature.

To determine the critical stress state for buckling, a so called “bifurcation func-
tional” was developed. For all admissable displacement fields the condition that the sec-
ond variation of the bifurcation functional is greater than zero ensures buckling does not
occur. However, buckling becomes possible when the bifurcation functional is zero for
some non-zero field. The analysis was done by considering three fields in the bifurcation
functional and then integrating over some local region, S, which was separated from the
rest of the sheet. Because the integration is carried out over a local region, the boundary or
continuity conditions are relatively unimportant.

The bifurcation model may also be written as a function of the buckling displace-
ments and a coefficient matrix, M. When the determinant of the matrix M goes to zero,
buckling was again predicted. Also, by minimizing the determinant with respect to some
waveform parameters and setting it equal to zero, the critical stresses associated with
buckling could be found. Because of the complexity of the equations, closed-form solu-
tions to the problem are difficult to obtain, therefore numerical solutions are generated by
means of the Newton—Raphson technique. This gave rise to an understanding of how the

material flowed in the plastic regime. The most widely used constitutive relations in buck—

22
ling are based on the J2 flow theory or the J2 deformation theory of plasticity. In addition,

the material behavior has been assumed to be isotropic and the power-law hardening rule
used.

A parametric study was done to determine the influence of material properties and
geometry of the critical stress state for local buckling. Material constant values of mild
steel and low-hardening steel were used and the results were plotted in terms of principal
compressive stresses and angle between the principal stress axes and the principal axes of
curvature. When comparing the flow theory results to the deformation theory results, the
deformation theory gave a more conservative prediction for the onset of wrinkling than
that of the flow theory and as a result was the preferred theory. Also, the results showed
that the critical stresses for the onset of wrinkling decreased as the strain hardening expo-
nent decreased and the critical stresses were also shown to decrease with the thickness of
the sheet.

Chan (1993) presented a method for locating the bifurcation points on the general-
ized force-displacement curve. This method was to be carried out in conjunction with any

technique which could traverse the limit point using the following equation,

[AFf + Alf/AF] = [KT] f[Auf + Afon] (1.3)

. . k .
where A)»: rs the load parameter vector, AF: rs the vector of unbalanced forces, Au 1. is

the residual displacement vector, AF is the reference load vector parallel to the applied

load vector, and A12 is the vector of conjugate displacement. The superscript is the load

increment and the subscript is the iteration within the load increment. The load parameter,

23

Alf, was chosen such that the equilibrium error was minimized and the load increment

was chosen using the constant arc-length method. This procedure has shown to have good
convergence and reliability against divergence.

A change in equilibrium state was detected by checking,

De:(K’;)-Der(l<’;")<o . (1.9)

When the above equation was satisfied, a critical point was found to lie between the solu-
tions at k and k+1. The location of the critical point was within the displacement degree of

freedom with the largest increase during the load increment. The secant method was then

used to predict the displacement increment, Aum, which caused the determinant of the

tangent stiffness matrix to go to zero. The iterative procedure outlined satisfied simulta-
neously the equilibrium and semi-definite conditions for the structural system and was
more efficient than satisfying them separately. Finally, equilibrium was assumed when
both the Euclidean norms of the residual displacements and of unbalanced forces were

less than 0.1% of their respective total or accumulated normal. The semi-definite condi-
tion was fulfilled when the change of correction load factor, Alf, was less than 1x10"5 of

the first load.

Yu and Zhang (1988) investigated the plastic wrinkling of an circular plate being
transformed into a cup by means of a deep drawing operation. To do this they attempted to
develop a criterion for combinations of geometry and material properties that gave rise to

wrinkling on the upper flange of the cup. The geometry of the investigation was simplified

24

by assuming an annular ring with a uniform tension applied on the inner edge. By non-
dimensionalizing the wrinkling equation, the Kantorovich method was used to produce an
ordinary differential equation (ODE). However before the Kantorovich method was
applied it was assumed that the transverse deflection of the ring was the product of two
functions, f, a function of the non-dimensionalized radius, and g, a function of the angle
theta and the number of waves. The resulting ODE was then a function of the non-dimen-
sionalized radius and was not readily solvable. Therefore the Galerkin method was applied
to obtain an approximate solution. The result of the Galerkin analysis yielded two equa-
tions which could be used to establish a criterion to show which combinations of geometry
and material properties would result in wrinkling. The authors then used these equations
to develop curves showing the faulty combinations and thus gave the limits for safe pro-
duction.

Finally, Adams (1993) investigated the wrinkling of a pre-tensioned, circular plate
subject to a concentrated load at the center. Using the von Karman plate equations, analy-
sis of the axisymmetric response was performed followed by an algorithm to determine if
wrinkling resulted from the calculated displacements. In finding the deflections and stress
state due to the load, the static, axisymmetric governing equations were cast into finite dif-
ference form. The resulting set of N nonlinear, banded algebraic equations were solved by
an iterative solution routine. This reference configuration must then be checked to see if

wrinkling was induced. Accomplishing this goal required returning to the full von Karman
equations and replacing the displacements by perturbation expansions in 8. Collecting

terms would reveal that the zeroth order terms reflect the axisymmetric results while the

first order terms gave a new set of governing equations which were coupled to the refer-

25

ence configuration by the transverse deflection and the circumferential stresses. The first
order perturbation displacements were then taken to be the product of a function of radius
and circumference and a harmonic function of time. Dependence on the circumferential
direction was removed by an exponential Fourier series expansion and using the orthogo-
nal property. This resulted in a set of 111 ordinary differential equations, m = 0, 1,..., 00. For
each Fourier number, m, a set of linear algebraic equations was generated by using a finite
difference approximation. Combining the boundary conditions with the finite difference
equations resulted in a linear, matrix eigenvalue problem. The eigenvalues found corre-
sponded to the square of the natural frequency, and when a zero eigenvalue resulted, wrin-
kling was assumed to have occurred. An equivalent method of prediction was to compute
a zero determinant for the finite difference coefficients. Adams’ results showed that the
axisymmetric transverse displacements could give rise to compressive circumferential
stresses which become large enough to induce wrinkling. The use of nonlinear membrane

theory would give wrinkling when the circumferential stress goes to zero.

1.3 The Present Study

The information in the literature has shown that the finite element method has been
successfully applied to the analysis of metal forming operations. In addition, the method
has also been used in the prediction of wrinkling in circular plates. Therefore, rather than
attempting to use trial and error methods or develop a closed form solution, the finite ele-
ment method will be used in the present work as a means for determining the stress state
and the deformation modes of the cylinder during the swaging process.

Before building the mathematical models necessary to simulate the necking pro-

26

cess with the finite element method, it would be useful to make some initial simplifying
assumptions about the operation. Thus an examination of the physical process of swaging
a cylinder was initiated to determine which parameters have the greatest impact on the for-
mation of wrinkles. These variables will be classified as primary parameters, as their
omission from computational analysis would render any results meaningless. A second
classification would contain those parameters which contribute to the deformation
response of the cylinder as it is formed, but do not directly contribute to the wrinkling phe-
nomena. That is the wrinkling will occur for any value of the secondary parameters within
the range of interest, and their omission from any analysis will not alter the numerical
results or the conclusions based thereon.
In the current study, the primary parameters are identified as:

1) geometry of the can - diameter, thickness, and final shape;

2) material properties - elastic and yield moduli, anisotropy, and plastic flow;
3) tool geometry - contact conditions;

4) forming steps - path dependency of the response.

The secondary parameters are:

1) rate of the forming process - material strain rate dependency;
2) friction - tangential interaction between the tool and the cylinder;

3) lack of axisymmetry - presence of weld seam and microstructural imper-
fections;
4) elastic snap-back - release of stored energy after the forming process.
In addition to the above assumptions, an analysis of the swaging process requires a
computational model that accurately accounts for large deformations, large strains and
rotations, nonlinear material behavior, rate-independent elastic-plastic material model,

and contact conditions, an algorithm for tracing the contact conditions between the tool

and the cylinder at all stages of the forming process. The complexity involved in formulat-

27

ing a fully operational program for each of these situations leads to the use of a commer-
cial finite element code for the analysis of the forming problem. Two programs from
MARC Analysis Research Corporation were well suited for this task, Mentat II and
MARC version K6.l (Mentat II, 1994 and MARC, 1994)

Using Mentat II, the geometry and parameters for the mathematical models of all
the occurring phenomena can be generated. These models are then translated into input for
the analysis program, MARC. MARC reads the input created by Mentat II, processes the
information, and calculates results for the unknown variables. These solutions are then
read back into Mentat H wherein they may be plotted in several manners. Most often the
results will be plotted as bands of color on the geometric model of the cylinder, with a
range of values being assigned a different color. In this way it will be possible to visualize
the levels of stress and strain within the cylinder as it a passes through the forming pro-

cess.

Chapter 2

Modelling of the Swaging Operation

2.0 Introduction

Simulation of the swaging process involves the creation of a computational model
to approximate the behavior of the cylinder during deformation. Therefore the domain of
the problem is defined by constructing a geometric model of both the cylinder and the
forming tools. The continuous region constituting the cylinder will then be discretized into
a set of elements and the deformation of each element approximated using the finite ele-
ment method. The selection of the element to be used for this purpose is the subject of
some scrutiny. While several elements are available to model the cylinder, not all will pro-
vide accurate results for a reasonable level of computational cost. Once an element has

been chosen the mechanical properties must by applied to the model. This includes not

28

29
only the elastic properties, but also the definition of material behavior after yielding. In

addition to geometric and material properties, the conditions which exist in the cylinder
prior to swaging (the initial stresses due to rolling the flat sheet into a cylinder) and the
known conditions at the boundaries must also be applied to the model in order to calculate
adequate results. Finally, the contact conditions of the problem, such as which bodies may
come into contact and the motion of the forming tools, are required as well. The accuracy
of the solution obtained is dependent on the quality of the definition of each of these com-
ponents because each contributes to the final result.

Once the model has been completed, the nonlinear nature of the problem forces the
need for a special solution algorithm. Usually the load-deflection relationship is linearized
and solved using the Newton-Raphson algorithm. The resolution of the linearized curve is
determined by the number of loading steps requested. Additionally, the accuracy of each
step in the solution will be tested with a convergence criteria and specification of the
degree of accuracy required. Finally, the large deformations involved in the necking pro-
cess necessitate the definition of a coordinate system to locate the material points as the

model of the cylinder is deformed.

2.1 Geometric Model of Cylinder and Forming Tools
The initial step in modelling the problem was to establish the geometry of the cyl-
inder. For the purposes of this study the physical dimensions used for the cylinder were
given as:
[D = 2.581 in. h = 0.007 in. , (2.1)

where ID is the inside diameter of the cylinder and h is the thickness of the material. How-

‘\
\
\\\\\‘

Figure 2.1 Relationship Between Modelled Surface and Full Cylinder.

ever, it was not necessary to model the fully three dimensional cylinder. Because the
geometry, the applied load, and the deformation response were all axisymmetric, the cyl-
inder and the subsequent models for the forming tools were represented by simple geo-
metric entities. In the case of the cylinder, the geometry consisted of a surface created by
slicing the cylinder along its length and using the cross-sectional plane which made up the
wall thickness. The cross hatched region in Figure 2.1 shows the modelled surface and its
relationship to the fully three dimensional cylinder. Only a small portion of the length was
modelled because the region in which wrinkling occurs is rather small in comparison to
the entire length of the cylinder. Thus, the deformations due to the forming occur only at
some finite distance from the necked region, which is approximately 0.7 inches in this
case.

The desired diameter reduction was to be achieved in two successive forming

stages by similar operations. In each stage the tool moves toward the cylinder and upon

31

contact between the tool and the cylinder, the leading edge of the cylinder was pushed
down to a new diameter, contained in a region called the “alley”. While each tool operates
in an analogous fashion, different levels of diameter reduction result from them.

The ramp of the first stage tool was established at 30 degrees from horizontal,
measured clockwise, and the diameter of the cylinder was reduced 5.73%, or until ID:
2.4330 in. The final diameter of 2.367 in. was realized through the use of a 35 degree ramp
from horizontal. Each of these rigid dies is represented as a set of bold lines in Figure 2.2,
the first stage being on the left and the second on the right. Also shown in Figure 2.2 for
reference are the cylinder, the cross hatched region, and the centerline of the geometry, the
dashed lined. The actual model as created by the MENT AT program is displayed in Figure
2.3. Here the surface representing the cylinder is in black and the dies are presented in red
with the first stage tool above the second stage tool. Again the center line of the geome-

tries is given for reference.

 

 

 

 

 

 

 

The “alley”
600
K :%5
1.2975 in. . /
1.2165 1n.
1.1835 in.
.11.- .1---__----__.---.__----_-_.__---_..--__--_--_-_-.--._

Figure 2.2 Geometry of Forming Tools with Reference to the Cylinder Axis. (Not to scale)

32

 

 

 

\f
‘\ I‘M..- v
-_ .. -- - ‘
\._
\ _ ‘ .. -..‘
,/'

 

Figure 2.3 MENTAT Model of Cylinder and Forming Tools.

 

33
2.1.1 Preliminary Element Selection

An important step in the modelling process was the selection of an appropriate
finite element which would allow for the most accurate results at an efficient
computational cost. Inspection of the geometry and the deformation modes of the process
helped in the identification of a limited group of appropriate elements. The reliability and
efficiency of each was then estimated by using each of the appropriate elements to model
problems with known solutions. The problem of reducing the diameter or necking of a
thin-walled cylinder requires an axisymmetric element which can be subjected to bending.
Given this criterion two groups of MARC elements were tested in a series of validation
problems. Table 2.1 shows a simple schematic and gives relevant information for both
groups of elements under consideration.

The first group contained axisymmetric, quadrilateral isoparametric elements. The
kinematics of the elements are based on a continuum theory. Only two DOF are required
to describe the deformations, axial displacement along the generator axis and radial
displacement perpendicular to the generator axis. The behavior of each element in the
group is governed by the interpolation of the displacement field and the number of
integration or Gauss points used for numerical integration. Element 10 consists of four
corner nodes and four integration points, which gives two point integration through the
thickness of the element. In contrast, element 28 has nine integration points, three through
the thickness, and assumes a quadratic displacement field. In addition, because element 28
consists of four corner nodes and four midside nodes, it has twice the number of degrees
of freedom per element, 16, when compared to element 10, which contains 8 dof. Finally,

element 116, like element 10, uses a linear displacement assumption, but the number of

34

Table 2.1 MARC Elements

 

Element

Topology / DOF

Attributes

 

 

10

- Axisymmetric

- Bilinear displacement
interpolation

- 2 x 2 Gaussian integration

 

28

U,V

- Axisymmetric

- Biquadratic displacement
interpolation

- 3 x 3 Gaussian integration

 

116

1:1 {:1 1:1

U,V

- Axisymmetric

- Bilinear displacement
interpolation

- 1 point Gaussian integration with
hourglass control

 

I

U, V,<I>

- Axisymmetric

- Linear displacement
interpolation

- 1 point Gaussian integration

- Simpson’s rule integration
through the thickness

 

 

89

 

l

U, V,<I>

 

- Axisymmetric

- Quadratic displacement
interpolation

- 1 point Gaussian integration

- Simpson’s rule integration
through the thickness

 

 

35

integration points has been reduced to one and hour glass control has been added. (More
information on the meaning of these two aspects of element 116 will follow.)

The second group of elements included element 1 and element 89. Both elements
are also axisymmetric and isoparametric, but, in contrast to the previous group, the kine-
matics are based on first order shear deformation shell theory (FSDT). The FSDT assumes
that planes originally straight and normal to the midplane remain straight but not necessar-
ily normal during deformation. This assumption allowed for inclusion of the shear defor-
mations which may arise. The displacement field was then be approximated using three
DOF at each node, axial and radial displacements and a right hand rotation in the plane.
Again the behavior of the elements is governed by the displacement field interpolation and
number of integration points along the length. Element 1 uses a linear shape function with
one integration point and element 89 uses a quadratic interpolation and a two point inte-
gration rule. For both of these elements the integration through the thickness is carried out
using Simpson’s rule based on a user specified number of points or shell layers determined

by the type of problem to be analyzed. At this time, the number of layers was set to three.

2.1.2 Evaluation of Element Performance for a Linear Problem

Given this set of eligible elements, the validity of each may be verified through the
comparison of results obtained for simple problems using these elements to those from
closed form solutions to the same problems. Therefore, each element was used to model a
test case of an infinitely long cylinder subjected to a circumferentially distributed line load
of unit magnitude. The data required for modelling the test case includes Young’s modulus

and Poisson’s ratio of the material and the radius and the thickness. Because some finite

36

length must be declared in order to use the finite element method, a finite cylinder was
modelled with sufficient length that the displacements near the ends were neglible. Addi-
tionally, it was appropriate to apply simply supported boundary conditions to the end
points of the finite cylinder. Two sets of geometric data were used to generate the model,

one being a moderately thick-walled cylinder:

R = 1.0 in. h = 0.1 in. L = 100.0 in. (2.2)

where R was the radius, h was the thickness, and L was the length of the cylinder. The

second set of data was representative of a thin-walled cylinder with the following

geometry:

R = 5.0 in. h = 0.05 in. L = 300.0 in. (2.3)

The material properties were the same for each of the cases:

E = 6.895x1010 Psi v = 0.3 (2.4)

where E was Young’s modulus and v is Poisson’s ratio. After obtaining approximate linear
solutions from MARC, the results were normalized to the closed form solution developed

for a thin-walled cylinder based on FSDT (Reddy, 1984) given by

37

QR B (2.5)

u3 (max) = - 25h

 

where u 3(max) was the maximum radial deflection, Q was the magnitude of the circumfer-
entially distributed line load, R was the radius of the cylinder, E was Young’s Modulus, h

was the thickness, and B was a parameter defined by

2 1/4
(3 = [ill—SJ] . (2.6)
R h

Using a load of 1.0/21: , the maximum deflection for the thick—walled cylinder was
calculated to be -4.8072x10‘11 in., while the thin-walled cylinder had a maximum

deflection of -2.9745x10'10 in. The normalized MARC results for the thick cylinder are
shown in Table 2.2 and the results for the thin cylinder are in Tables 2.3. In all cases the

advantage of symmetry was taken to reduce computational times.

2.1.3 Discussion Results for a Linear Problem

A discussion of the results obtained from the MARC analysis of the test problem
requires an examination of the assumptions that were involved in the formulation of the
elements and the closed form solution. The closed form solution was based on FSDT.
Thus any result obtained using Equation 2.5 was approximate. Since the normalizing
solution was based on shell theory, one would expect elements 1 and 89 to have better
agreement with the closed form solution because of the consistency in formulation.

Elements 10, 28, and 116, being based on a continuum theory may not converge exactly to

Table 2.2 Predicted maximum deflection due to a circumferential line load in a
simply supported cylinder of radius = 1.01m, thickness = 0.1 in., and length = 100.0

in. All values are normalized to the closed form solution.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Mesh. Size Element
(wrth Element 10 Element 28 116 Element 1 Element 89
symmetry)
50 x 1 0.4458 0.6689 0.6555 1.0321 1.0285
50 x 2 0.4462 0.6683 ------------------------------
100 x 1 0.6020 0.8499 0.8782 1.0321 1.0285
200 x 1 0.7666 0.9810 0.9776 1.0321 1.0285
200 x 2 0.6027 0.9817 --------------------------
400 x 1 0.8861 1.0063 1.0021 1.0321 1.0285
600 x 1 0.9208 1.0103 1.0009 1.0321 1.0285
600 x 2 0.8559 1.0215 ------------------------------
750 x 1 0.9312 1.0142 1.0193 1.0321 1.0285
1000 x 1 0.9387 1.0180 1.0301 1.0321 1.0285

 

the results given by Equation 2.5 for thick shells, but should approach the same results for

thin shells.

Examination of the results for elements 1 and 89 in Table 2.2 showed good

correlation with the closed form solution, with the error for each being approximately 3

percent. In comparison, substantial mesh refinement was required for the continuum

elements to converge. Using the 50 x l mesh size as a baseline, elements 28 and 116

required an increase of about four times as many elements, while twenty times the number

of elements were necessary for convergence of element 10. Despite the need for an

 

39

Table 2.3 Predicted maximum deflection due to a circumferential line load in a
simply supported cylinder of radius = 5.01m, thickness = 0.05 in., and length = 300.0

in. All values are normalized to the closed form solution.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Mesh Size Element
(wrth Element 10 Element 28 116 Element 1 Element 89
Symmetry)
50 x 1 0.1833 0.2786 0.4101 1.0032 1.0032
50 x 2 0.1833 0.2786 ------------------------------
100x1 0.2499 0.4365 0.6695 1.0032 1.0032
200 x 1 0.3438 0.6986 0.8760 1.0032 1.0032
400 x 1 0.4762 0.9116 0.9659 1.0032 1.0032
750 x 1 0.6235 0.9839 0.9921 1.0032 1.0032
750 x 2 0.6267 0.9850 ---------- ---------- ---------
1250 x 1 0.7494 0.9981 0.9989 1.0032 1.0032
1500 x 1 0.7896 0.9998 1.0001 1.0032 1.0032
2000 x 1 0.8396 0.9996 1.0012 1.0032 1.0032
2000 x 2 0.8562 ---------- --------- -------- ----------

 

 

increased number of elements, 10, 28, and 116 did provide acceptable solutions to the test
problem. Examination of the results for the thin-walled case should provide a better
understanding of the behavior of the elements for the diameter reduction problem as it too
is thin-walled problem.

In Table 2.3 the shell elements again quickly converged to the result of Equation
2.5. Also as before, the mesh had to be refined in order for elements 10, 28, and 89 to
approach convergence, with elements 28 and 116 continuing to yield results nearly

corresponding to the analytical result. On the other hand, element 10 required a far greater

40

number of elements to converge and the final results obtained (200x2 mesh) were still in
error by approximately 15 percent. The large error may be a result of the element
experiencing locking and/or the size aspect ratio of the element.

Locking may be understood by considering the case of a thin beam subject to pure
bending. For a situation such as this the transverse shear and membrane energies should be
negligible throughout the element. However, the element approximation may not be able
to satisfy the kinematic constraints of vanishing shear and membrane strains and also
adequately represent the bending strains. As a result, the element may be too stiff and not
converge. The most prevalent method for eliminating locking is to use fewer Gauss
integration points than the minimum required for exact integration of the energy
quantities, thereby relaxing the constraints in the element and allowing the bending energy
to control the deformation (Averill and Reddy, 1990). One drawback to this solution
technique is that a singular or nearly singular stiffness matrix may result which will cause
large deformations to appear without any change in strain energy, a so-called zero energy
or hourglass mode. Control of these modes may be attained by adding a small amount of
stiffness, akin to adding springs, to their contributions in the coefficient matrix of the
model. However, insufficient stiffness will not prevent the spurious energy modes, while
too much will prevent viable modes from being found. Element 116 has been derived
using both reduced integration and hourglass control, thereby presenting a viable
alternative to element 10 which exhibits locking. While element 116 mandated more
elements for convergence than the standard mesh, by a factor of 8, the increase in mesh
size was manageable and the error for a mesh of this size was less than four percent.

Another solution to the problem of element locking is to increase the number of

41

DOF within the model. Two ways of doing this were available: a) increase the elements in
the mesh or b) use an element with a higher order interpolation scheme. As can be seen
from the results in Tables 2.2 and 2.3, the mesh refinement technique seemed to be
working, though a converged solution was not yet attained using element 10 for these two
cases. Use of element 28 increased the DOF in the analysis by increasing the interpolation
order. Once again the table shows that a greater number of elements were needed to model
the cylinder (about fifteen times more) than the benchmark, but the converged result was
less than two percent in variance with the normalizing result.

The size aspect ratio, the length of the element divided by the height, also
contributed to the failure of element 10 to converge. In the thick-walled case the 50 x l
mesh resulted in an aspect ratio of 10, meaning the element length was 10 times longer
than its height. Thus, to achieve the desired deflection the element was required to bend
while experiencing very little shear. This was in contrast to the case of the 100x] mesh,
where the minimum aspect ratio was one half and shearing was the dominant mode of
deformation. So, the mesh of 1000 x 1 converged to the closed form solution because the
shear was not forced to vanish. The maximum and minimum aspect ratios in the thin case
were 60.0 and 1.5 respectively. In both cases the bending energy was governing the
deformation of the elements and, thus, locking was a consequence. Keeping the aspect
ratio close to one or less, the phenomena of locking may be relieved. In order to confirm
this hypothesis a third test case was run.

In this final case the same material properties as before were used and the

geometry was as follows:

42
R = 100.0 in. h = 1.0 in. L = 200.0 in.. (2.7)

This geometry allowed greater flexibility in creating meshes with the aspect ratio of one.

The analytical solution was again calculated using Equation 2.5 and found to be

1.4835x10'll in. Because elements 1 and 89 did not exhibit any difficulties in converging,
only elements 10, 28, and 116 were used in modelling this second thin walled case. Of
course both elements 28 and 116 did converge to the normalizing solution, but the hope
was that by maintaining the aspect ratio as close to one as possible, these elements would
converge faster. The results from the MARC analysis are given in Table 2.4. As expected
element 10 did converge to the normalizing solution. Also, as was expected, elements 28
and 116 both converged to the analytical solution at a faster rate. Actually, they converged
with the initial mesh. The results obtained indicated that the effect of the aspect ratio was

significant and should be kept to a value of 1.0 or less.

Table 2.4 Predicted maximum deflection due to a circumferential line load in a
simply supported cylinder of radius = 100.01n., thickness = 1.0 in., and length = 200.0
in. All values are normalized to the closed form solution.

 

 

 

 

 

 

 

 

 

Mesh Element 10 Element 28 E16111“: nt
== _ = =
100 x 1 0.8989 1.0041 1.0046
200 x 1 0.9717 1.0049 1.0043
400 x 4 0.9961 1.0050 1.0048
800 x 4 1.0000 1.0050 1.0048

 

 

 

 

 

 

43
2.1.4 Nonlinear Evaluation of Element Performance

The ability of element 10 to represent the kinematics of a thin structure was still in
question and the decision was made to test this element in a nonlinear situation. So a
nonlinear test case was devised using the same geometry as given in Equation 2.7, but

with I;100.0 and the following material properties:

1511 = 3.230x107 Psi 1522 = 2.72x107 Psi 1533 = 2.72x1o7 Psi

v12 = 0.3 v13 = 0.3 v23 = 0.3
7 7 7 (2.8)
G12 = 1.242x10 Psi Gl3 = 1.242x10 Psi 023 = 1.242x10 Psi

0), = 7.785x10“ Psi

where Em are the Young’s moduli, Gm are the shear modulus and 0'), is the yield stress.

The subscripts given refer to the direction in which the property is valid. The 1 direction
was defined to be parallel to the axis of the cylinder, 2 was in the circumferential
coordinate, and 3 was parallel to the radial direction. In order to keep the computational

times to a minimum plasticity was not included in this test case. A compressive

circumferential line load of 4.5x106 pounds was applied to one end of the cylinder in the
model. This load was chosen as it would provide for approximately the same amount of
diameter reduction as in the swaging process under study. Also, because this problem was
to be solved nonlinearly, Equation 2.5 was no longer valid and no other closed form
solution was readily available. Hence the converged result of element 28 was chosen to
serve as the normalizing solution for this case due to its previously demonstrated ability to

converge. Table 2.5 contains the values obtained for this validation problem.

Table 2.5 Nonlinear test case results corresponding to the predicted maximum
deflection due to a circumferential line load in a simply supported cylinder of radius
= 100.01n., thickness = 1.0 in., and length =100.0 in. All values are normalized to the

converged solution of element 28.

 

 

 

 

 

 

Mesh Element 10 Element 28
100 x 1 0.3078 0.9451
200 x 2 0.6322 0.9811
400 x 4 0.9498 1.0000
800 x 8 0.9537 -----

 

 

 

 

 

The data attests to the fact that element 10 will indeed converge in a nonlinear
analysis. This, coupled with earlier information from the other test cases, indicates that
any of the elements tested would be suitable for use in the forming process, but other
information helped to eliminate some elements from further consideration.

One of the major features of the MARC finite element code is the ease in which the
user can define contact bodies and conditions. Potential contact bodies are identified by
the user, and the code internally defines the contact constraints. However, this algorithm
does not detect contact at midside nodes, thus eliminating higher order elements from
contact analysis. Elements 28 and 89 were, therefore, not eligible to be used in the
analysis of the forming problem as the topology of these elements includes midside nodes.
This left elements 10, 116, and 1 as the remaining possibilities. Each of these elements had
capabilities as well as insufficiencies for use in forming problems. One drawback to

element 116 was the reduced integration technique which provides only one Gauss point

45

within the element. Because the strains and stresses (and, hence, the elasto-plastic material
properties) will vary from the outer surfaces inward, a single integration point through the
thickness would not adequately predict the variation of these values through the thickness
of the cylinder. It then becomes necessary to increase the number of elements in the
transverse direction to gain more Gauss points. In contrast, Element 1 allows the user to
input the number of points to use for Simpson’s integration rule through the thickness,
eliminating the need for multiple elements in the radial direction. The final decision of

which element to use was left until each could be used in the necking analysis.

2.2 Boundary Conditions and Contact Constraints

2.2.1 Boundary Conditions

Some values of the primary degrees of freedom must be known on the boundaries
of the cylinder in order to constrain rigid body motions. Therefore, the next step in the
modelling process was to establish boundary conditions for the model. For the solution of
the diameter reduction problem, it was appropriate to fix the axial displacements on the
end opposite the necked region to zero. In doing so the continuity requirements between
the modelled and unmodelled portions of the cylinder were satisfied. Since transverse dis-
placements are constrained by the axisymmetric response of the cylinder and none of the
elements under consideration have the ability to deform in the circumferential direction,
this was the only boundary condition needed for the solution of the forming problem. The
applied boundary conditions are illustrated in Figure 2.4 by the arrows located on the left

hand end of the cylinder.

46

 

 

 

aunt:

 

Figure 2.4 Geometric Model with Applied Boundary Conditions.

 

47
2.2.2 Initial Stresses

Prior to the actual forming process, the cylinder had a pre-existing state of stress
due to the rolling of a flat sheet into a cylinder. These effects were also included in the
model. It was assumed that the bending was the result of equal yet opposite moments
being applied to the edges of the plate resulting in a state of pure bending and, hence, an
exact elastic-plastic solution could be obtained. The effects of anticlastic curvature were,
however, neglected. For elastic or plastic deformations, the circumferential strain in the

cylinder after bending was obtained from geometric considerations as:
599 = g (2.9)

where R=1 .294 in. was the mean radius of curvature of the cylinder and 2 was the distance

from the midplane of the cylinder. Of course, the maximum strain occurred at 2 = i;

and was equal to 0.0027, which exceeded the yield at:

O’
.39 = 0.00227 = am (2.10)
E00

where (type was the yield stress in the circumferential direction, E99 was the elastic mod-
uli in the same direction, and Eype is the yield strain. While some plasticity would be evi-

dentin the cylinder, the entire cross-section would not yield due to bending. Rather, there

would exist an elastic core near the midplane of the flat sheet and yielding would occur

48

near the edges of the plate, with the plastic zone increasing toward the midplane as curva-

ture was increased. In the final configuration, the elastic core extended a distance 2e, away

from the midsurface found by:

O'
YPB .
ze, = 3.3-R = 0.00294 in. (2.11)

and illustrated in Figure 2.5 in which the cross section of the cylinder is shown with the
distribution of stresses resulting from bending the flat sheet.

Determination of the limits of the elastic core was important in that MARC does

Plastic Region

N I

 

 

 

 

Plastic Region

Figure 2.5 Cross Section of Cylinder with Stress Resulting from Cylindrical Bending.

49
not allow the initial stresses to exceed the elastic yield limit of the material. Being that the

initial stresses must be input at the Gauss points of each element, it was important to
ensure that the stress at these points did not exceed the yield stress. Thus, the coordinates
of the Gauss points in reference to the inner diameter of the cylinder needed to be deter-

mined. Using this information the } coordinate of each was also found with
2 = y(GP) -- ° (2.12)

where y(GP) was the coordinate of the Gauss point measured from the inner diameter of
the cylinder and h was the cylinder thickness.

Initial stresses were determined in the following manner. From the bending strain

found in Equation 2.9, the stresses at 2 = i; were found using the one dimensional

Hooke’s Law,

91090 = 500300 (2°13)

where fine was the calculated yield stress in the circumferential direction. Finally,

assuming a linear distribution of stress through the thickness, and using the results from

Equations 2.12 and 2.13, the stress at the Gauss points was found using:

6
a (GP) = #292 (2.14)

50
where 0(GP) is the stress at the Gauss point. In the event that the Gauss point fell within

the plastic zone, the initial stress value was assumed to take on the value of stress at the
elastic/plastic interface.

An alternate method of calculating the initial stresses would have been to assume
a linear distribution through the elastic region of the cylinder and having the stress in the
plastic zone maintain the stress value at the elastic/plastic interface. Employing this
scheme would have resulted in a more accurate model of the initial stresses and a 16 per-
cent increase in their magnitude. The results to be presented later will show that either
approach does not significantly impact the final results of the model. In either case, the
input of this data into the model results in circumferential stresses being present at the out-

set of analysis.

2.2.3 Contact Conditions

Using MARC, the user is able to define a deformable body through the definition of a set
of finite elements and a rigid body is defined as a set of geometrical entities, such as lines,
circles, splines, or surfaces.

In many commercial finite element codes, contact is defined by placing gap-
friction elements between the two bodies that may come into contact during the analysis.
The difficulty in this method is that the user is required to know a priori which bodies
come into contact. While this approach is valid in MARC, another option is also available
to automate the contact definition procedure. The basic thrust behind the option is the
definition of bodies rather than gap elements. All the required information to enforce non-
penetration is contained on the surfaces of these boundaries. In the case of a deformable

body, which is defined by a set of finite elements, all the nodes at the boundary become a

51

set of candidate nodes for contact. For rigid bodies only the surfaces of potential contact
need to be defined by a set of geometric entities, such as lines. As the analysis is
performed, the contact condition is checked by calculating the vector product of the nodal
displacements and the normal of the rigid body. If this result is less than or equal to the
distance between the bodies contact is detected and constraints are applied to the contact
nodes to prevent penetration. Because the distance between the deformable and rigid
bodies is a calculated quantity, some round off error will be expected and lead to some
penetration. To account for this phenomenon, a contact tolerance is provided to allow
nodes to go an incremental distance below a surface and yet still be considered in contact.
This may be a user defined quantity or assigned by the algorithm. For a continuum
element, the default is 1/20 the smallest element edge and for a shell element it is 1/2 the
element thickness. Finally the rigid body motion of the dies is given in the contact option
by defining an instantaneous velocity. By explicitly integrating the velocity over time the
motion path of the tools is provided (MARC, 1994).

For the swaging analysis the elements constituting the cylinder were classified as
the deformable body and the forming tools declared as separate rigid bodies. The action of
the forming process is such that the tools move toward the cylinder, so each rigid body
defined was given a series of velocities to simulate their displacements. At the present
time, displacements rather than velocities will be discussed with assumption that the time
integration limits have been chosen properly to achieve the correct velocities. The first
stage forming process was defined by applying an axial displacement of -0.3729 inches to
the body defining the first stage tool. This die was then pulled off the cylinder by using a

release option provided in MARC to decouple the nodes of the deformable body from the

52

rigid die in conjunction with a second rigid body translation of 0.4271 inches in the axial
direction. Simulation of the swaging process was completed by the motion of the second
tool. Initially the position of the tool relative to the cylinder was quite significant and
would prevent the analysis from finding a converged solution. Therefore the first
movement of the second stage was to bring the tool to a point just before contact would be
detected. The required displacement was -0.1060 inches. Forming was done by then
moving the die -0.2558 inches axially. This amount of movement was sufficient to bring
the ramp of the tool into a position parallel with the previously formed angle in the
cylinder. The process was completed by again using the release option with the 0.4271

inch translation of the tool.

2.3 Material Model

The response of the cylinder to the action of the forming is also dependent on the
mechanical properties of the material used to construct the cylinder. Among these
quantities are properties to describe the elastic response of the cylinder, but more
important are the characteristics which dictate the behavior after the elastic limit has been

exceeded.

2.3.1 Mechanical Properties

An elasto-plastic model was developed to describe the behavior of double-reduced
tinplate sheet steel conforming to A.S.T.M. A-623 Type L, the material of choice in this
analysis. The model began with the specification of the mechanical properties which are

given in Table 2.6 (Wood, 1994). The in-plane values for Young’s moduli and yield stress

53

were determined by averaging the results obtained from five samples used in standard

tensile testing procedures. Given the thickness of the material, current experimental

procedures were inadequate to measure the transverse material properties.

Table 2.6 Mechanical Properties of Double-Reduced Tinplate Steel, Type L.

 

 

 

 

 

 

 

 

 

 

 

 

Property Axial, x Circumferential, 0 Transverse, z
Young’s Modulus. 3.2311107 in. 2.721(107 in. 2.72x107 in.
E (psi)
Poisson’s ratio, v 0.3 0.3 0.3
Shear Modulus. 1.242x107 in. 1.24211107 in. 1.24211107 in.
G(psi)
Yield Stress. 7.785x104 in" 6.544x104 in. 6.544x104 in.
oyp(psr)

Therefore, the shear moduli were calculated using,
E

with E being given by the axial value of Young’s modulus. This resulted in isotropic

values for the shear moduli, in addition to the given isotropic values of Poisson’s ratio.

This assumption of isotropy in the shearing values will not lead to significant error. In fact,

transverse shear effects are usually neglected in a structure with a radius to thickness ratio

of greater than 20 (Dym, 1974). In this case the radius to thickness ratio is in excess of

180.

54

2.3.2 Plasticity Model

The definition of the plasticity model requires the yield stresses, a method for
determining the onset of yielding (a yield surface), a definition of the hardening curve, and
finally the effects of hardening on the yield criteria. The yield stresses are given in Table
2.6 and were found using the same tensile tests as used in finding the elastic properties.
Since the cylinder was in a state of multi-axial loading during the simulation, a multi-axial
yield criterion was required. Several models exist for this purpose, most notably among
them is the von Mises yield criterion wherein the material is assumed to have yielded
when the distortion energy is equal to the measured uniaxial yield stress (Mendelson,
1968). For an anisotropic material, the relationship of distortion energy to yield stress is

given as

2 2 2 2

2 2 2
+ 304123 + 3a5‘t31 + 3a6t12

where am corresponds to the axial yield stress, oij are the calculated stresses referred to

the coordinate axes of the cylinder (11 - axial, 22 - circumferential, and 33 - transverse), ti]-
are the calculated shear stresses in the cylinder during analysis. The coefficients aa, 01:1 to

6, account for directional variations in the yield stress from the axial and may be found by
defining the Hill’s yield ratios YRDIRl, YRDIR2, and YRDIR3 in conjunction with the

following

55

a1:———l———+___1.__ l

YRDIR22 YRDIR32 1011911912

,2=___1_,__1_ 1

YRDIR32 YRDIRlz YRDIR22

,3._._1_.,__1__ 1

(2.17)

 

 

YRDIR12 YRDIR22 YRDIR32
_ 2 _ 2 _ 2
a4 - ———2 a5 — ——-5 a — ———2
YRDIR3 YRDIR2 YRDIRl
where
O' 0' 1 0'
YRDIRl = Y”: YRDIR2 = "’R YRDIR3 = 332 (2.18)

OYP GYP OYP

I I X

For double-reduced steel, YRDIR1=1.0, YRDIR2=0.8405, and YRDIR3=0.8405 were
used.

After the material points in the cylinder reach the plastic region of the stress-strain
curve, the material will begin to harden. The amount of hardening the material
experiences must also be described through a mathematical model. Isotropic hardening,
used in this model, allows the yield locus to retain its center point throughout the analysis,
but the size of the locus will grow with increasing strain hardening. Inputting a table of
plastic strain versus stress into the MARC data deck allowed the program to internally
account for these effects. Table 2.7 contains the measured values for the double reduced
tinplate sheet steel. When this data was plotted in a stress versus plastic strain curve, as in
Figure 2.6, it was observed that the curve was nearly flat. This would indicate that the
strain hardening will not have a significant impact on the results obtained for the model.
Proof of this hypothesis will be left until the model has been completed and the analyses

run.

56

Table 2.7 Stress versus Plastic Strain for Double-Reduced Tinplate Sheet Steel,
Type L.

 

Stress Plastic Strain

 

 

i —

7.7854x104 Psi 0-0

 

 

7.902411104 Psi 6.4000x10'4

 

7.9727x104 Psi 1.3700x10'3

 

8.0079x104 Psi 2. 190011103

 

8.0109x104 Psi 3.10001110-3

 

 

 

 

 

Stress (Psi)
(a) # UI

to
1

 

0 0.5 1 1.5 2 2.5 3
Plastic Strain x 104

Figure 2.6 Stress versus Plastic Strain for Double Reduced Tinplate Sheet Steel.

57
2.4 Solution Parameters and Analysis Options

2.4.1 Solution Method

As previously mentioned, because the problem under consideration is nonlinear,
the Solution path is also nonlinear. In this case a method of tracing this curve was required
in order to obtain solutions. Therefore, the Newton-Raphson method, which is one of
several techniques, was employed to approximate a solution to the resulting nonlinear
algebraic equations. The advantage of the Newton-Raphson method is that it provides a
quadratic rate of convergence to the solution (Zienkiewicz, 1989). In order to understand
the Newton-Raphson method, consider a fictitious load-deflection curve as shown in

Figure 2.7. The desired load is given as PC. In cases where the entire path must be known,

such as in plasticity, the curve may be approximated by loading the structure
incrementally. Within each load step the algorithm ensures the equilibrium by reaching the
solution point in an iterative sequence which is terminated when some convergence
criterion is satisfied.

Accordingly the motion of the forming tools was applied incrementally. Recalling
the definitions of the tool motions in Section 2.2.3, the first stage tool has two separate
translations, one onto the tool and the second a releasing action. In the Newton-Raphson
solution the initial tool movement was divided into 250 equal load steps and the second
was done in 2 unequal steps. The first increment of the releasing action was very small in
order to ensure that the motion started in the correct direction and the second increment
completed the releasing motion. For the second tool, three movements were given:
bringing the tool near the cylinder, the forming operation, and the releasing action. The

initial movement of the second stage tool was performed in one load step, the forming was

58

Load

 

 

 

 

 

 

 

 

Deflection

 

UA UB UC
U1 U2 U1 U2

Figure 2.7 Newton-Raphson solution method.

again divided into 250 equal increments, and the releasing action also performed in two
increments, one small and the other large enough to complete the requirements. Within all
of the load increments a maximum of 20 iterations were allowed, but this was set as upper
bound. After each iteration, the maximum displacement increment was divided by the
maximum displacement change (for the load increment) within the model to calculate a
convergence tolerance. When this tolerance was found to be less than or equal to 0.10 the
increment was said to have satisfied the convergence criterion and the next load increment
was begun. In the event that convergence was not achieved, upon reaching the end of the

20th iteration the program terminated giving a message that the convergence criterion was

59

not met. This leads to some comments about the number of load increments used. 250
increments was found to be a lower limit for convergence reasons. When a smaller
number of increments was used, the load steps were too large for the algorithm to attain a
converged result within the specified number of iterations. Using a greater number of
increments could result in a more accurate solution, with the trade-off being increased
solution time. Therefore, the initial simulations were done with 250 increments, but the

effects of increased load steps on the solution were also investigated.

2.4.2 Analysis Options

The large deformations that are inherent in the process of swaging the cylinder
required a reference frame from which the displacements may be measured. In addition,
the particular measures of strain and stress that are to be used must also allow for the large
deformation effects. MARC offers a set of options to account for this phenomenon, of
which three were used in combination in the simulation model. These were the large
displacement, finite strain plasticity, and updated Lagrange options. The first two cards
account for the large deformations and strains as their names imply. The last option,
updated Lagrange, gives a reference frame for the measurement of deformations and
calculation of the stiffness matrix. Here the mesh is connected to the material throughout
the analysis and is updated at the end of each load increment. By using these three options
in combination the entire spectrum of large deformation were taken into account in the
model. And with that the definition of the model for simulating the swaging process was
completed. The forming process may now be analyzed and an understanding of the

deformations and stress state of the cylinder during forming can be sought.

Chapter 3

Numerical Simulation of Swaging

3.0 Introduction

The model for the simulation of the swaging process has been described. Using
this model preliminary deformations of the cylinder as it progresses through the process
may be observed. The deformations are preliminary in that up to this point two rather
important details have not yet been determined. The element that is to be used is yet to be
finalized, with the remaining choices being elements 10, 116, and 1. In addition to the
element decision, the size of the mesh is also a lingering question to this point. Upon
making a decision about these components, the final deformations and stress state of the
cylinder during and after the necking process may be analyzed. The information revealed

may then be used to ascertain possible causes for the formation of wrinkles in the cylinder

60

61

and ultimately may lead to candidate solutions.

3.1 Simulation of the Swaging Process

Now that a model has been defined, the swaging process can be simulated.
Although, it should be noted that the results shown are preliminary until final selection of
an element and mesh size. Thus, a baseline mesh of 150 elements along the length and 2
elements through the thickness, or a 150 x 2 mesh, was used in concert with element 10 to

produce the initial simulation shown in Figures 3.1 (a) through (i). In Figure 3.1(a) the

th increment of the simulation, before the forming

cylinder and tools are shown at the zero
operation had begun. As the first stage tool progressed forward, the leading edge of the

cylinder was bent down and traveled along the tool ramp toward the “alley” as shown in

Figure 3.1(b), the 50th increment of the simulation. Figure 3.1(c) depicts the leading edge

of the cylinder as it was just being bent back into the “alley”, increment 100, so that the

diameter reduction of the first stage could be achieved. By the 250th increment, the leading

edge had slid down the “alley” until the end of the tool was reached, as illustrated in Fig-

ure 3.1(d). Finally the first stage was completed as the tool was removed in the 252“d
increment of analysis, as seen in Figure 3.1(e).
In Figure 3.1(i) the second stage tool has moved to a position near the cylinder so

as to prevent numerical difficulties in attaining convergence, increment 253. The tool then

pushed the cylinder’s leading edge down until the “alley” was again reached at the 352“d
increment of analysis, as shown in Figure 3.1(g). The advancement of the second tool was

discontinued when the angled portion of the cylinder was parallel to the ramp of the tool

62

 

 

 

 

 

Figure 3.1(a) First Stage of Swaging at the Zeroth Increment.

 

M

 

 

 

 

Figure 3.1(b) First Stage of Swaging at the 50th Increment.

63

 

.0
‘
i 08

 

 

 

 

 

Figure 3.1(c) First Stage of Swaging at the 100th Increment.

 

ii '.
' ‘
[to

 

 

 

 

 

 

 

Figure 3.1(d) First Stage of Swaging at the 250th Increment.

 

“.
ital

 

 

 

 

 

 

Figure 3.1(e) Release of the First Stage Tool.

 

e

e
u
a

a
0
2.1“
:w

 

 

 

1m

 

 

 

Figure 3.1(f) Tool Moved into Position in First Increment of Second Stage.

65

 

£133
igir

 

 

 

8.24

L

”I

 

 

 

Figure 3.1(g) Second Stage of Swaging at the 352nd Increment.

 

: ~

0
:SJM
m

 

 

 

 

#4

L.

”1

 

 

 

Figure 3.1(h) Final Increment of the Second Stage.

66

 

iii?"
l1”

 

 

 

 

 

 

 

Figure 3.1(i) Release of Second Stage Tool.

which occurred at the 503rd increment and is shown in Figure 3.1(h). The forming process

was completed in the 505th increment when the second stage tool was removed from the

cylinder, as illustrated in Figure 3.1(i).

3.2 Final Element Selection

Three elements were still available for use in the modelling of the swaging
process, elements 10,116, and 1. Each of these was used the model previously'discussed.
The results from each of these analyses were then evaluated to determine which provided

the most accurate solution at the best computational cost.

67
3.2.1 Forming Analysis Results Using Element 1

After establishing the baseline with element 10, element 1 was used in the simula-
tion of the forming process to determine if the efficiency which was exhibited in the test
cases still held true. In review, Table 2.1 indicates that element 1 is a two-noded, isopara-
metric, and axisymmetric element which employs Simpson’s rule integration through the
thickness. As such, when substituting element 1 into the model, the number of layers
within the element needed to be defined. For the metal forming analysis to be conducted
11 layers were specified. This was equivalent to using 11 Gauss points through the thick-
ness which was attractive in that an accurate representation of the material effects should
be the result. However, when the analysis was run, the model consistently failed to pro-
vide a solution. Either the solution failed to converge or the stiffness matrix of the model
became singular (or nearly singular). Several attempts were made to overcome these diffi-
culties, but no resolution was found. Therefore the decision was made to abandon further

use of element 1.

3.2.2 Forming Analysis Results Using Element 116

The second element that was used in the model was element 116. Recalling Table
2.1, element 116 is a four-noded, isoparametric, and axisymmetric element employing
reduced integration and hourglass control. In the previous studies of Section 2.1.1, this
element also provided good convergence characteristics within a reasonable amount of
computational time, and the same was hoped of the swaging simulation. Once again,
however, the analysis deteriorated in all cases with element 116 providing the kinematical

model for the cylinder. Inspection of Figure 3.2 reveals the presence of hourglass modes

68

 

INC: I!
808: 0
TIIE:3.9200-01
FREQ:0.0000-o-00

 

am

 

 

 

10M

 

Figure 3.2 Swaging Analysis Using Element 116 with a 150x2 Mesh.

 

69
within the model at the 96th increment of the analysis. Thus it may be concluded that the

stiffness that was added to the element in order to prevent zero-energy modes from
occurring was not sufficient for this particular model. Thus, alternative approaches of
prevention were tested. The first run using element 116 only had two elements through the
thickness. It was believed that this may not provide enough stiffness in the presence of the
plasticity model, nonetheless, increasing the number of through the thickness elements did
not eliminate the hour glass modes. Therefore, a‘second attempt was made using a mixed
model with two elements of type 10 at the leading edge of the cylinder and the remaining
elements being type 116. This again proved to be a fruitless effort. As a result, element

116 was also no longer considered a viable element for the analysis at hand.

3.3 Convergence of Forming Problem

Element 10 was left as the only available element to discretize the model of the
cylinder. \Vrth an element now selected, the mesh size needed to be finalized by
conducting another convergence study similar to that of Section 2.1.1. Again, for the case
of the forming process no tractable analytical solution is available, and experimental
studies of the problem do not provide an adequate description of the internal stress state of
the cylinder. Therefore, to gain confidence in the finite element solution, the model must
be evaluated until a converged solution is found. For the swaging simulation, convergence
must be achieved in two manners.

The mesh used to model the geometry should be refined in both the axial and radial
directions to gain monotonically convergent results. Several reasons exist for the need to

smooth the model in both coordinates. The cylinder initially has regions of plastic stress

70

due to the rolling of the flat sheet into a cylinder at the inner and outer surfaces as
discussed in section 2.2.1. These plastic regions will grow radially inward as the cylinder
is subjected to the forces of the forming dies. In order to capture the material effects of the
growing plastic region, the model should have a large number of Guass points through the
thickness of the cylinder cross section. The convergence of the calculated values in the
radial coordinate would indicate the accuracy of the material model. In addition, the
number of elements in the axial direction were increased to capture the path dependency
involved. This was done in conjunction with a study on the number of load steps that were
necessary for reliable results.

The number of load used to model the forming process will also contribute to the
path dependency of the results obtained from the model. This is because the physical
process of swaging the cylinder is carried out in two smooth, continuous motions.
However, in the finite element model of the process, the tools were given a velocity which
was divided into a number of equal load steps. By this method, the displacement of the
rigid dies is not a continuous function of time, but rather a piecewise linear approximation.
Of course a greater number of time steps taken will increase the smoothness of the tool
displacement which affects the path dependency of the cylinder. Larger time steps will
move the tool a greater distance along the axis of the cylinder and will result in greater
incremental displacements. Therefore it is recommended that a convergence study on the

number of load steps also be conducted.

71
3.3.1 Mesh Refinement Study

A study of the elemental convergence of the first stage was conducted to determine
which mesh would provide for the most accurate results. For this convergence study, the
baseline mesh remained 150 x 2. The elements were then increased through the thickness
of the cylinder up to eight. The number of axial elements were then expanded to 300 and
finally to 400. In the case of 300 elements the through the thickness elements were the
same as for the 150 meshes. For the meshes containing 400 elements, analyses with 2 and
3 elements through thickness were not run because it was already determined that more
elements through the thickness were needed to account for the material variations in that
direction. Upon completion of each analysis, the values of axial and circumferential strain
and stress were collected on the outer surface, mid-plane, and inner surface of the cylinder
at two points within the mesh. The first of these points, labeled ‘Point A’ in Figure 3.3,
was within the first bending region and the second was in the critical region, ‘Point B’ in
Figure 3.3. Also it should be noted that the data was gathered at the 240th increment of the
analysis and not the final 250th increment as specified in the model. The reason for this
was that some of the models exhibited slight penetration through the end of the tool,
which resulted in extraneous axial stresses in the cylinder. Therefore, to make a fair
comparison between the models, the 240th increment was used to avoid the penetration
problem. The strain and stress data resulting from the convergence study can be found in
the Appendix Tables A.1 through A.24. For compactness, the discussion that follows will

refer to that data, but specific data values are not mentioned.

72

 

 

THE : LN
FREQ : 0.000900

 

 

 

 

Figure 3.3 Sampling Locations for Convergence Study of Forming Analysis.

 

73
3.3.2 Discussion of Mesh Convergence Study Data

The analysis of the gathered data began with an evaluation of the strain values,
where very little variation resulted from increasing the number of through-the-thickness
elements. In Figures 3.4 through 3.7, the strain values versus mesh size were plotted for
the critical region and appear to contradict this assessment. While the data for those
meshes with 150 elements and 400 elements along the length were in good agreement
with one another, meshes with 300 elements along the length exhibited large variations in
the strain values with increasing mesh refinement through-the-thickness. This is, however,
misleading due to the scaling used for the generation of these plots. In fact, an
examination of the percentage change in the strain resulting from refining the mesh size
would be more revealing as has been done in Tables 3.1 and 3.2. The observation could
then be made that the changes in strain magnitude were actually less than four percent for
each component on both of the cylinder surfaces. The exception was the axial strain on the
outer surface of the cylinder. Here the values varied by up to 350 percent, however, the
magnitude of strain in this region was an order of magnitude smaller than the other strain
measures. Thus, this strain component was not dominant in the swaging simulation and
therefore did not alter the earlier conclusion. Similar conclusions may be drawn for the
strains in the first bending region. In any case, these small changes in strain were not
surprising because the displacements in the first stage of the forming process were
controlled by the geometry of the tool. The route followed by the cylinder would not be

significantly altered by an increase in the number of elements through the thickness.

74

 

 

 

 

 

 

 

3 0
a
t
3
w
3 -1
'5
O
l
E -2 P — 150 Elements
I; - - 300 Elements
'2 3 ~- ' - 400 Elements
4 l-
_5 A 1 l l l l
2 3 4 5 6 8
Elements Through the Thickness

Figure 3.4 Axial Strain on the Outer Surface in the Critical Region vs. Mesh Size.

 

 

 

 

 

 

 

-0.0615P
a. - - - - ‘9
I \
I 1
-0052 X ‘
8 ' I ‘,
g ’I g- - .. ‘
w A l 4% -1‘\~ k A A
3 m v 1 ‘ If; v
5 ~4 1 7
04.0525- . I
1 ‘ /
7.! X 7’
‘E ‘ I
D
g -o.osa- ‘. ,’
”S —150Elements “ ,
,s_. - - aooElements 1 ,’
Q-0.0535» r-~-400Elements ‘, ,’
\
t”
4
_o.m‘ l J 1 J 4 1
2 3 4 5 6 B
EIemntsThroummeThldmess

Figure 3.5 Circumferential Stress on the Outer Surface in the Critical Region vs. Mesh
Size.

0.044 -

.09
5%

.o

9

d
I

Axial Strain - Inner Surface
0
. 0
§ 2

.0
r

0.037 -

 

Figure 3.6 Axial Strain on the Inner Surface in the Critical Region vs. Mesh Size.

 

75

 

—-— 150 Elements
- - 300 Elements
- - ~ - 400 Elements

 

 

 

 

vol

4 5 6 8
Elements Through the Thickness

 

 

 

 

 

 

-0.0517 -
- r ‘ ‘1
4.1518 " IT— 1
\
I
\
8 4.051% ,' 1 —150Elements
E [I X - - 300 Elements
‘3 .0052- , “ ----4OOE|ernents
g 1
' 4.0521 - ‘1
g 1
a \
3 4.0522) 1‘
S 1
3.‘ -o.0523 - 1
E l
E H- - — -‘ l
\
o -O.m24 " ‘3‘ \
v \ 3 fi-
‘k '\ \ .. — - -+
\ I ' ‘1'- -----
-o.0525L \ ‘ , x'
‘4
-O.m26 ' ' ‘
2 3 4 5 6 8
Elements Throum the “midmess

Figure 3.7 Circumferential Strain on the Inner Surface in the Critical Region vs. Mesh

Size.

Table 3.1 Percentage Change in Axial Strain Values in the Critical Region. (%)

 

Elements Through the Thickness

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Cylinder Axial
Surface Increase 2 3 4 5 6 8
150 to 91.46 267.71 286.70 81.60 34.98 30.70
3 300
5 30010 444.70 -348.90 -21.83
400
1: afi-
150 to 0.70 14.05 -0.95 -140 1.53 0.92
S 300
t:
5 300 to 2.90 3.34 0.60
400

 

 

 

Table 3.2 Percentage Change in Circumferential Strain Values in the Critical Region.

 

 

 

 

 

 

 

 

 

 

 

Elements Through the Thickness
Cylinder Axial
Surface Increase 2 3 4 5 6 8
150 to 0.06 -0.26 1.16 1.22 -2.69 -0.06
3 300
5 300 to -1.03 -1.07 2.59
400
150 to 0.06 -0.23 1.24 1.30 -0.09 -0.06
s 300
t:
E 300 to -1.09 -1.13 -0.06
400

 

 

 

 

 

 

 

 

 

 

 

77
The formation of wrinkles in the cylinder is directly attributable to the build up of

stresses and strains in the critical region. Therefore attention was also focused on the
calculated stress values obtained from the convergence study where it was found that the
trends contrasted those of strains. The change in stress states due to an increase in
elements along the length was found to be minimal, but a change in the number of
elements in the radial coordinate had a distinct impact on the stress values as shown in
Figure 3.8 through 3.11. For the cases of 150 and 300 elements along the length, these
plots show the stresses found on the outer and inner surfaces in the critical region versus
the number of elements through the thickness. The effect of increasing the axial elements
had very little consequence on the stress values obtained from the model as evidenced by
the close agreement between the two curves in each of the plots. In contrast the
multiplication of the number of elements through the thickness seemed to have a dramatic
effect. As discussed in Section 2.2.1, the cylinder enters the forming process with initial
stresses which induce zones of plasticity near the outer and inner surfaces of the cross-
section. During the forming of the cylinder these plastic regions grow radially inward and
the resulting stress gradients were adequately captured only by increasing the number of
elements (thereby increasing the Gauss points) through the thickness. Accordingly, as the
mesh was refined through the thickness, the values of stress began to converge. However,
as can been seen in Figures 3.8 through 3.11, a converged solution was not attained with

the meshes used.

78

I
01
O

1

3

i

 

 

Axial Stress (szi) - Outer Surface
1

 

 

 

 

 

-250 .—
/ — 150 Elements
“300 ” - - 300 Elements
-350 .
4m 1 1 l 1 J I
2 3 4 5 6 8
Elements Through the Thickness

Figure 3.8 Axial Stress on the Outer Surface in the Critical Region vs. Mesh Size.

 

—150Elements
-- 300Elements

2'?

 

 

 

l. 1.
s 1.3

erential Stress (szi) - Outer Surface
1

 

 

 

-250.
E 300.
§
6 _350 _
-400-
2 3 4 5 6 8
Elements Through the Thickness

Figure 3.9 Circumferential Stress on the Outer Surface in the Critical Region vs. Mesh
Size.

79

 

 

 

 

 

 

 

350~

300~
O
0250*
E
3
a)
2
.E 200'
1
6
a
5150-
8
2
a
2100-

50- — 150 Elements
- - 300 Elements
0 1 l I l 1 1
2 3 4 5 6 8
Elements Through the Thickness

Figure 3.10 Axial Stress on the Inner Surface in the Critical Region vs. Mesh Size.

i

 

i

 

 

 

ii

3

 

Circumferential Stress (szi) - Inner Surface

1?

 

o 1 l l l l
2 3 4 5 6

Elements Through the Thickness

@1-

Figure 3.11 Circumferential Stress on Inner Surface in the Critical Region vs. Mesh Size.

3.5

1"
to 01

CPU Time (seconds)
in

0.5

 

80

 

 

 

 

 

x 10
_ , ,+
,+’ ’ ’
/
.. /
I
r: ,
_ , — 150 Elements
/ - - 300 Elements

’ 4t

2 3 4 5 6 8

Elements Through the Thickness

Figure 3.12 Computational Time vs. Mesh Size.

81

The question still remains as to which mesh should be used in the simulation of the
forming simulation. Two final analyses of the convergence data provided the answer. First
of all, was the rather minimal increase in strain and stress values obtained from increasing
the number of elements along the length computationally efficient? This query was
answered by plotting the CPU time (the amount of time the computer took to solve the
finite element equations) versus the mesh size, as was done in Figure 3.12. The CPU times
were obtained by performing each analysis on a Hewlett-Packard Series 9000 model 715
Unix workstation running at 75 MHz with 32 megabytes of RAM and 1 Gigabyte of disk
swap space. In Figure 3.12 the meshes containing 300 elements were observed to require
an average of 1.57x104 CPU seconds (4.4 CPU hours) more CPU time to complete than
those with 150 elements. Therefore, the 150 element meshes appeared to be more

attractive, but final judgement was reserved until a comparison of the results was made.

3.3.3 Comparison of Stress Trends Resulting from Mesh Refinement

Finally, the axial and circumferential stresses in the critical region for three
meshes; 150x8, 300x8, and 150x2, were plotted to see the variations present. By zooming
in on the critical region, the stress gradients were magnified to produce Figures 3.13
through 3.18. In each of the figures, the color bar indicates that compressive stresses were
present on the outer surface, while the inner surface was subjected to tensile stresses.
Comparing the stresses for the 150x8 mesh in Figure 3.13 and 3.14 to those of the 300x8
mesh in Figure 3.15 and 3.16 the only evident improvement was in the lengthwise
gradients, but the gains were minimal. This information coupled with the run time data

lead to the conclusion that 150 element meshes should be used for any subsequent analysis

82
of the forming process. In looking at the results from the 150x2 mesh in Figure 3.17 and

3.18 the trends remained the same, but the contours were more banded than in the
previous plots. This was a result of only two elements through the thickness being used.
As such, the stresses can only be linearly interpolated between four points rather than
sixteen, thus yielding poorer gradations of stress. But the trends still show compressive
stresses on the outer surface and tensile stress on the inner. This appeared to indicate that
the 150x2 mesh would be appropriate for the simulation, provided the conclusions are to
be based on the stress trends and not the absolute values.

Earlier, concern was expressed about the presence of locking in meshes where the
element aspect ratio was greater than 1.0. The 150 x 2 mesh has an aspect ratio of 1.33 and
as a consequence locking should be a concern in the second stage of the forming. In order
to judge the effects of locking an 800 x 2 mesh, aspect ratio = 0.25, was run. The
displacements from that analysis were plotted in Figure 3.19 for increment 303 and in
Figure 3.21 for increment 373. Analogous plots were constructed for the 150x2 mesh and
are shown in Figures 3.20 and 3.22. The results were very nearly the same in both meshes,
thereby indicating that the 150 x 2 mesh was not overly stiff. A possible explanation for
the non-existence of locking could be the complicated state of plasticity in the model, but
an in depth evaluation for this phenomenon was not performed. But, the absence of
locking certainly solidified the argument to use the 150 x 2 mesh in the swaging

simulation.

83

 

 

 

 

...... L

uterine!“ I

 

 

 

 

Figure 3.13 Axial Stress in the Critical Region for the Second Stage of the 150x8 Mesh.

84

 

 

 

 

4M

4m

"M L

“We!“ 1

 

 

 

Figure 3.14 Circumferential Stress in the Critical Region for the Second Stage of the
150x8 Mesh.

85

 

 

 

 

 

per

 

 

 

I‘M“... 1

 

Figure 3.15 Axial Stress in the Critical Region for the Second Stage of the 300x8 Mesh.

86

 

 

 

 

 

" L

 

 

”Genet..- 1

 

Figure 3.16 Circumferential Stress in the Critical Region for the Second Stage of the
300x8 Mesh.

87

 

 

 

 

 

 

 

 

or“
,.
M
V
4.01m“
Id!
tdcerlpd ”as 1

 

Figure 3.17 Axial Stress in the Critical Region for the Second Stage of the 150x2 Mesh.

88

 

 

 

 

 

 

 

1m

 

 

“Met...- 1

 

Figure 3.18 Circumferential Stress in the Critical Region for the Second Stage of the
150x2 Mesh.

89

 

a: ..
....
ee

ii.‘

 

 

 

 

Figure 3.19 Deformations of the 800x2 Mesh at Increment 303.

 

i ..
ee

 

 

 

 

 

Figure 3.20 Deformations of the 150x2 Mesh at Increment 303.

90

 

11::

 

 

 

 

 

Figure 3.21 Deformations for the 800x2 Mesh at Increment 373.

 

-sn

- m
' O

.
1'7...
.
'm
.

 

 

 

 

y:

L

psi

 

 

 

Figure 3.22 Deformations for the 150x2 Mesh at Increment 373.

91
3.3.4 Convergence Study on Number of Load Increments

The last component of the mesh that needed to be finalized was the number of
increments required. Some indications of path independence were expressed above. To
prove this was the case, both stages of the 150x2 mesh were used in an analysis with 500
load steps for each stage of forming instead of the usual 250. Once again, very little
variation in the strains was evident. The stress data in the critical region for the first stage
was collected in Table 3.3 and the second stage in Table 3.4. The data indicated that the
largest change for any of the values was less than 3.0 percent, and therefore it could be
said that the increase in the number of increments used to perform the analysis made no
change in the results. As such, the simulation of the swaging process could confidently be

carried out using the 150x2 mesh and 250 increments.

3.4 Results from Forming Simulation

3.4.1 Strain and Stress States in the Critical Region

Once a mesh was settled upon, results from the simulation were analyzed to
determine possible mechanisms for wrinkling.The results from the final analysis are
shown in Figures 3.23 through 3.26. In each of the plots the cylinder was in the final
increment of the forming process (increment 503) and the data was collected in the critical
region of the cylinder. The first contour plot, Figure 3.23, consists of the axial strains. As
before, the color bar on the left hand side designates the dark red (almost a purple) region
on the outer surface of the cylinder to be in a state of tension which increases in magnitude
to the inner surface that is yellow. Intuitively, one would think the outer surface should be

in compression as a thin-walled structure in bending would exhibit. However, the severe

92

Table 3.3 Stress Values for 250 vs. 500 incs for the First Stage in the Critical Region.

 

 

 

 

 

 

Number Oxx 03“ Cu 0’06 0'90 969
of Ines (outer) (mid) (inner) (outer) (mid) (inner)
250 -288.52 8.38 272.58 -366.38 -64.87 208.42
500 -283.95 8.14 267.11 -361.82 -65.36 204.21
% change 1.58 -2.86 -2.01 1.24 -0.75 -2.02

 

 

 

 

 

 

 

 

Table 3.4 Stress Values for 250 vs. 500 incs for the Second Stage in the Critical

 

 

 

 

 

 

 

 

 

 

Region.
Number (rxx ox, (1,,x 699 699 099
of Incs (outer) (mid) (inner) (outer) (mid) (inner)
250 #26348 :-4.12 240.95 477317 =--—--¥-63. l7IE 175.29
500 -263.56 -4.23 240.22 -276.62 -63.04 174.92
% change -0.03 -2.67 -0.30 0.32 0.21 0.21

 

 

 

 

 

 

 

 

 

state of bending to which the cylinder was being subjected has created membrane forces
which overwhelm the bending forces resulting in tensile axial strains. Figure 3.24 shows
the circumferential strains as increasing in magnitude moving left to right, corresponding
to increased necking. Each value of strain was approximately constant through the
thickness, analogous to a cylinder undergoing a uniform diameter reduction. The axial
stress distribution displayed in Figure 3.25 indicates that the outer surface of the cylinder

is in compression, as shown by the blue region, and the inner surface is in tension, as

93

shown by the red zone. In this case the membrane forces were as dominant as in the case
of strains and so the bending forces were significant enough to create the compressive
stresses on the outer surface. Finally, the stress gradients in the circumferential direction
are given in Figure 3.26 and once again the outer surface was blue, or compressive, and
the inner surface was red indicating tension. Recalling the circumferential strains were
found to be uniformly compressive through the thickness, the presence of tensile

circumferential stresses was initially surprising. These will be explained in the sequel

3.4.2 Assessment of Strain Hardening and Initial Stress Effects on
Strain and Stress States

The effects of strain hardening and the initial stresses were examined by carrying
out additional analyses, once without the work hardening table and the other without the
initial stresses. In Table 2.5 the global maximum and minimum values of strain and stress
are shown for the strain hardening models. From the data tabulated it can be seen that the
hypothesis of Section 2.3.2 was indeed correct, the strain hardening had very little effect
on the outcome of the analysis. Further comparison can be made by contrasting the
contour plots of the model with strain hardening, Figures 3.23 - 3.26, to those without
strain hardening which may be found in the Appendix, Figures A] through A4.

The effects of the initial stresses were also evaluated. In Section 2.2.2 two methods
of calculating the initial stresses were given. At that time it was stated that either method
would not significantly affect the final results obtained from the model. Table 3.6 contains
the global minimum and maximum strain and stress values from two models, the first
being the model run with the initial stresses calculated by the first method (linear

distribution through the thickness) and the second a model analyzed with no initial

94

model are not dependent on either method.

and Stress Values

stresses. From the data it may be concluded that the impact of initial stresses on the strain
and stress values was minimal. More importantly the trends of stress distribution did not
change, which can be seen by comparing the plots in Figures 3.25 and 3.26 to those of
Figures 3.27 and 3.28. As previously stated, the model used does not yield a converged
solution, and any explanations offered for the formation of wrinkles in the cylinder must
be based on the trends of stress. Therefore, despite the fact that an alternate method of

calculating initial stresses was available, the conclusions drawn from the results of the

Table 3.5 Effects of Strain Hardening on the Global Minimum and Maximum Strain

 

 

 

 

 

 

 

 

5xx(m3X/ 899(max/ oxx(rnax/ 099(rnax/
Model . rmn) rmn) mm)
mm)
With Strain 4.686x10‘2 1.210x10'3 3.357x105 2.768x105
Hardening
-1.568x10'2 -8.7O7x10'2 -3.811x105 -3.481x105
Without 4.685x10'2 1.207x10'3 3.39211105 2.862x105
Stain Hard-
ening 1574;110-2 -8.703x10'2 -3.863x105 -3.508x105
0.0 -0.2 1.0 3.3
% Change
0.4 0.5 -1.4 0.8

 

 

 

 

 

 

95

and Stress Values

Table 3.6 Effects of Initial Stresses on the Global Minimum and Maximum Strain

 

 

 

 

 

 

2.2mm 86991:“ “MW stew
Model . mm) mm) mm)
mm)
4.686x10'2 1.210x10'3 ' 3.357x105 2.768x105
With Initial
Stresses -1.568x10’2 -8.707x10’2 -3.811x105 -3.48lx105
4.471x10'2 1.089x10'3 3.463x105 2.867x105
Without Ini-
“al Stresses -1.563x10'2 -8.705x10'2 3923:1105 -3.597x105
4.8 -11.1 3.1 3.5
% change
3.2 0.0 -2.9 -32

 

 

 

 

 

 

96

 

 

 

 

 

 

 

 

£171.43
4M
V
I“ I
post
id Gen. ofTeU w 1

 

Figure 3.23 Axial Strains in the Critical Region at the 503rd Increment.

97

 

 

 

 

 

 

 

4041.02
of”
Y
arch-ea l
per
31 On. of Ted huh 1

 

Figure 3.24 Circumferential Strains in the Critical Region at the 503rd Increment.

98

 

 

 

 

 

 

“To.“ ‘

Idcondfinee 1

 

 

 

 

Figure 3.25 Axial Stresses in the Critical Region at the 503rd Increment.

99

 

Figure 3.26 Circumferential Stresses in the Critical Region at the 503rd Increment.

 

 

 

 

 

 

 

income!“ 1

 

Figure 3.27 Axial Stresses in Critical Region at 503rd Increment, Calculated without
Initial Stresses.

101

 

 

 

 

V
Mstfl'ho-fl I

”We!“ 1

 

 

Figure 3.28 Circumferential Stresses in the Critical Region at the 503rd Increment,
Calculated without Initial Stresses.

 

Chapter 4

Discussion of Wrinkling

4.0 Introduction

An assessment of the forming operation reveals that two basic modes of
deformation are responsible for the diameter reduction of the cylinder. The first is a
uniform diameter reduction and the second is the bending needed. to get the leading edge
of the cylinder into the “alley”. While these two modes are, in fact, coupled due to the
nonlinear nature of the problem, to understand their contribution to the formation of
wrinkles, two arguments will be made. The first is based on the existence of some critical
circumferential strain in the cylinder. The second will reveal the contribution of the

bending energy to the formation of wrinkles.

102

103
4.1 Existence of Critical Strain

For illustration purposes, consider a straight column under an axial load, as shown
in Figure 4.1. From strength of materials it is known that this column has a critical load,

PCR, which, when exceeded, will cause the column to snap into a new configuration upon

the application of some small lateral load. This critical load may be found using Euler’s

equation,

- —, (4.1)

where E is Young’s modulus of the material, I is the moment of inertia of the column, and

L is the length of the column. The moment of inertia for a rectangular column is given as,

I = 1—2bh , (4.2)

where b is the width of the column and h is the thickness. Substituting Equation 4.2 into

4.1 and knowing that the area of the column is given as:

A = bh , (4.3)

Equation 4.1 may be written as:

104

 

1.

 

 

 

Figure 4.1 Euler Buckling Column.

Mxx

 

'80

Figure 4.2 Shell Element Subject to Bending and Resultant Moments.

105

it h A
PCR = (4.4)

12L2 '

 

Finally the critical load may be written in terms of stress, which in turn may be related to a

strain by the following equations,
PCR = oCRA em = £80,. (4.5)

Thus, the critical strain relationship is found by,

em = 12(2):. (3.6)

From Equation 4.6 it may be seen that the critical strain of the column with a constant
length, L, is a function of the thickness, h, squared. So as the thickness of the column is
increased the critical load of the column grows in a quadratic manner.

The issue of wrinkling occurs in a cylinder under external pressure. For this
geometry, the critical strain is proportional to the square of the thickness, h, over the
radius, R, of the cylinder. Thus, if the induced strain is equivalent in the two cylinders, the
thicker cylinder may not exceed its buckling strain, whereas the thin-walled cylinder
under consideration could. From the viewpoint of a cylinder subject to a uniform diameter

reduction the importance of the wall thickness is realized. However, the cylinder was also

106

loaded in bending, which gives rise to a more complicated state.

4.2 Bending Effects

Isolating the applied bending loads, a second approach is proposed to explain the
results obtained for the variation of circumferential stresses as shown in Figure 3.26. The
geometry to consider for this explanation is that of a cylindrical shell element depicted in
Figure 4.2. Imagining that this shell was participating in the forming process, then as the

leading edge of the shell was entering the “alley” it could be considered to be subjected to
a bending moment M H = M , shown in blue. The moment-curvature relations are given
35,

M = D (xxx+v1cee)

, (4.7)
M99 = D (lc66 +vxxx)

where xxx and K99 are the changes in curvature and D is the flexural rigidity given by

3

12(1—v2)

The changes in the circumferential curvature were assumed to be neglible in comparison

to xxx so let, K99 = 0. So K99 was eliminated from both equations in Equation 4.8,

yielding,

107

Mn = D(Kxx+0)

. (4.9)
M69 = D (0+vrcxx)

Examination of the second of Equations 4.9 reveals the existence of a significant resulting

circumferential moment,

M99 = Dvrcn, (3.10)

shown in red in Figure 4.2.

M99 is responsible for the regions of tensile circumferential stress in Figure 3.26,

while in the presence of uniform compressive circumferential strains, Figure 3.24. As the
cylinder was initially bent toward the “alley”, corresponding to Point A in Figure 3.3, it
was being loaded. Then, as it was bent back to go into the “alley”, Point B in Figure 3.3, it
was unloaded. The action of unloading the cylinder allows for the tensile stresses while at
the same time producing compressive strains. This phenomena is demonstrated in Figure
4.3 which shows a history plot of the leading edge node on the inner surface as it travels
through the second stage of the forming process. As can be seen at the zeroth increment
the strains were zero and the stresses took on the value of the initial stresses. The first load
step of the analysis takes both the strains and stresses into the compressive region and by
the one hundredth load step the cylinder has reached the maximum values of strain and
compressive stress. This corresponds to the point at which the leading edge was ready to

enter the “alley”. Subsequent load steps take the leading edge down the “alley”, and the

108

stress strain curve enters the region of tensile stresses and compressive strains.

It has been shown that some critical strain exists within the cylinder. Upon
exceeding that strain, any disturbance, whether due to the process or material
imperfections, will cause the cylinder to buckle. However, the portion of the cylinder in
the “alley” region of the tool does not exhibit the wrinkling phenomena. Thus, it is
hypothesized that the wrinkling is due to the combined effects of diameter reduction and

severe bending. Additionally, the resultant moment, M99, is responsible for the wrinkles

protruding inward rather than outward. While the simple models used up to this point
present possible explanations of the wrinkling problem, they do not, however, allow the
prediction of the onset of wrinkling. They only provide a means to understand the
phenomenon, and, at the same time, the suppositions made allow the possibility to present

some recommendations.

4.4 Introduction of Stiffening Rib

Based on the level of understanding of the mechanics involved in the formation of
wrinkles during the swaging process, a recommendation was proposed as a possible
means of alleviating the wrinkling phenomenon. Because the stability of the cylinder was
governed by a complex relationship between stress, material properties, and geometry,
techniques for increasing the level of stability could be derived from one or more of these

factors. After evaluation of several different solution methods, the addition of a circumfer

109

 

 

and Carpetsueeeutoow)

M1 mm

 

 

 

 

 

 

T
i
/
/
/

 

/

 

/

 

/

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4.0“

mmurwsuungw;

Figure 4.3 Circumferential Stress versus Strains for the First Stage of Forming.

 

110

ential stiffening rib presented itself as the most viable option.

In Section 4.2, the cylinder was shown to have a critical circumferential strain
above which the cylinder could buckle. In addition the severe bending necessary to
accomplish the diameter reduction also contributed to buckling. The combination of these
two mechanisms is believed to give rise to wrinkling. The diameter to which the cylinder
is reduced is an integral feature of the finished product and cannot be altered. Therefore,
the options would be to reduce the compressive/tensile stress gradients or increase the cyl-
inder’s resistance to localized bending associated with the wrinkling mode. Stress reduc-
tion could be achieved by actually thinning the wall thickness of the cylinder in localized
zones. But this would then weaken the resistance of the cylinder to axial loading, for
example in the press fitting of a top, and as such this option did not seem to be acceptable.
Also, a means of thinning the material in a high volume production setting could not be
established. Thus, stress reduction did not appear to be practical. As a result attention was
then focused on increasing the cylinder’s resistance to bending.

Again, several methods were available to increase bending stiffness. An inspection
of the equation which relates the changes in curvature to the bending moments shed light

on these possibilities,

(4.11)

7:
ll
Ell:

Equation 4.11 indicates that with a constant bending moment either the Young’s modulus,
E, or the moment of inertia, I, could be increased to result in smaller changes in curvature.

Because the material used produced other desirable characteristics, the decision was to

lll

attempt to increase the moment of inertia. This concept is reinforced by considering a
beam in bending. In Figure 4.4(a) a bending moment is applied to the wide, flat face of the
beam, whereas in Figure 4.4(b) the bending moment is applied to the thin edge of the
beam. Intuitively, the beam of Figure 4.4(a) should bend more easily than the one in Fig-
ure 4.4(b) because of the increased vertical thickness of the beam in Figure 4.4(b). This is

confirmed by Equation 4.12 which gives the moment of inertia for the beam,

I = 1—2bh . (4.12)

where b is the width of the beam and h is the vertical thickness. Thus, Ia for Figure 4.4(a)
is smaller than 1b for Figure 4.4(b). Substituting Ia and lb into Equation 4.11 individually,

it can be seen that the curvature for the second beam is indeed smaller than that of the first
beam.

In a flat sheet, a stiffening rib can be used to increase the effective moment of iner-
tia (see Figure 4.5). The addition of the central depression in the plate will act as a stiffen-
ing rib in the plate and give an increased resistance to bending. These same principles
could be applied to the cylinder. For the cylinder, the moment of inertia is a function of
both the wall thickness and the (current, or instantaneous) cylinder radius. Increasing the
wall thickness has already been deemed undesirable and a global change to the diameter is
not feasible. So, a local change in the geometry was proposed to increase the moment of
inertia. This would produce the desired increase in bending stiffness, while, at the same

time, possibly enhancing the aesthetics for commercial applications.

112

Figure 4.6 shows a model of the cylinder with the addition of a circumferential
stiffening rib. The placement of the rib was such that upon the completion of the swaging
process the rib would approximately be at the midpoint of the ramped portion into the
“alley”. The axialwidth of the stiffening rib was chosen to be one third of the overall
length of the ramped portion of the cylinder and the depth was arbitrarily chosen to be one
quarter of the width of the rib length. The initial stresses induced by the formation of the
stiffening rib were not calculated as the primary purpose of the following numerical
results was to observe the effect of the stiffening rib on the forming paths and the final
geometry. Figures 4.6 through 4.8 show the initial geometry of the cylinder with the stiff-
ening rib, the formed geometry, and the final geometry upon tool release. The final geom-
etry has a characteristic stairstep shape. However, despite the apparent large change in
geometry, the strain and stress distribution was not significantly affected, as evidenced by
comparing Figures 4.9 through 4.12 to Figures 3.23 through 3.26 which did not include
the stiffener. Nevertheless, the objective behind the stiffening rib was not reduction in
stress and strain level, but an increase in the stability of the critical region.

Introduction of the stiffening rib into the final geometry of the cylinder may be
achieved in several manners. The simplest would be to use rollers or a stamping operation
on the flat sheet before bending into the shape of the cylinder. A consequence of this may
be difficulties in obtaining strong weld seams in the swaged zone of the cylinder. The stiff-
ening rib may also be produced in the cylinder after the welding process has been finished.
One possibility is to use a clamp-like tool with appropriately shaped complimentary roll-
ers on opposite sides. The clamp can be closed, rotated around the circumference of the

cylinder to create the desired shape, and then removed. Finally, the stairstep shape may be

113

realized by using a multi-stage forming process similar to that of this study to gain addi-
tional “steps” in the critical region. For example in the current tooling the amount of diam-
eter reduction produced by the first stage tool could be decreased and the translation of the

second stage truncated. This would require only small modifications to the current tooling.

114

M

_ C\.

 

M
(a)

Figure 4.4 Simple Beam Subjected to a Bending Moment.

M M
\q_\‘\]
a) (b)

(

 

 

 

Figure 4.5 Plate Subjected to a Bending Moment.

115

 

 

 

 

 

 

Figure 4.6 Initial Geometry of Cylinder with Added Stiffening Rib.

 

 

116

 

INC: 503 “I”
SIB: 0 w
“:11“

FREOHMOOOo-N

 

 

 

 

Figure 4.7 Formed Shape of the Cylinder with Added Stiffening Rib.

117

 

iiii
i!”

 

 

 

 

 

 

\ m\\ .
L... 0 1

L

 

 

 

Figure 4.8 Final Shape of the Cylinder with Added Stiffening Rib after Elastic Snapback.

118

 

 

 

 

lab!

taco-monarchs. 1

 

 

 

 

Figure 4.9 Axial Strains in the Critical Region of the Cylinder with Added Stiffening Rib.

119

 

Ii"
it"

 

~‘-'.
‘

 

 

 

 

 

 

-1
w
W
m
4.1m
am
one
a“
v
L_J
lebi
ace-parental» 1

 

Figure 4.10 Circumferential Strains in the Critical Region of the Cylinder with Added
Stiffening Rib.

120

 

 

 

 

 

 

 

 

“tee-potatoes 1

 

Figure 4.11 Axial Stresses in the Critical Region for the Cylinder with Added Stiffening
Rib.

121

 

 

 

 

41M

“Carpet..- 1

 

 

 

Figure 4.12 Circumferential Stresses in the Critical Region for the Cylinder with Added
Stiffening Rib.

Chapter 5

Summary and Future Work

5.1 Summary

The stated objectives of the project were to determine the deformations and the
stress state of the cylinder during the swaging operation. Each of these goals has been
accomplished through the use of a numerical simulation of each step in the manufacturing
process. Simple examples were then used to interpret the results obtained from the analy-
ses. From the explanations offered, a hypothesis of the mechanism responsible for the for-
mation of the wrinkles was presented. Finally, candidate solution methods were proposed,
with the most reasonable being partially studied using in the finite element model.
Throughout this simulation the cost was in man time to develop the model and to subse-
quently run the model on the computer. At no time was any physical simulation per-

formed, thereby yielding a cost reduction in raw materials for test construction of the

122

123

cylinder and the forming tools. In addition, any follow up work or additional proposals
can be performed using the developed model as a basis. Shortcomings have, however,
been demonstrated within the model.

The use of element 10, which is based on a continuum formulation, has proven to
be inherently expensive for the type of analysis being done, as pointed out by Rebolo, et
al. (1990). In the initial convergence studies done to confirm the most appropriate element
for use in the simulation, element 10 required a large number of elements to reach the
exact solutions sought, while the shell elements converge to the closed form solution
quickly. Again, in the convergence studies performed to determine mesh size, conver-
gence was slow and actually never achieved using element 10. Although acceptable solu-
tions were found in each case, the computational cost involved was considerable. The
conclusion must be drawn that an alternative should be found. The model developed does
not predict the formation of wrinkles in the cylinder. This does not diminish the results

obtained up to this point, but only presents the opportunity for further work.

5.2 Future Work

It is proposed to develop a new finite element based on nonlinear shell theory and
a novel formulation that appears to be ideally suited for analysis of sheet metal forming
processes. The theory is based on a displacement field that allows a cubic distribution
through the thickness of the inplane displacement components and satisfies the shear trac-
tion boundary conditions on the top and bottom surfaces of the shell. For cylindrical shells

undergoing axisymmetric deformations, the displacement field takes the form,

124
2 3
z 42 (1W
“x (xx) = u (x) 4.wa (x) '1' Zk51+3—h2[k52—WX(x) -5]

2 3 (5-1)
uy (x,z) = v (x) + 2% (x) + Zzku "' 3411-2 [1‘42 _ ‘15”) — Y%]

uZ (x,z) = w(x)

where ux, try, and uz are the displacements of a point in the shell in the x (axial), y (circum-

ferential), and z (thickness) directions respectively, and,

 

 

1 1‘ 1" 1 1‘ “I"

k =_ 5 _ 5 k =_ 5 + 5
51 2( (N) (1)] 52 2[ (N) (1)]

Q55 Q55 Q55 Q55
t b
1 T4 T4

t b (5.2)
k“ = 2[Q(N) 'Qm) k42 = 2[Q(N) +Q11)
44 44 44 44

 

 

1 b . .
t4 , I; and 1;, 1:5 are the shear tractrons on the top and bottom surfaces of the shell in x

and y directions, respectively and Q 3) , Q :4") and Q5? , Q55?) are the transverse shear

stiffnesses in the first and Nth layer of the shell in the yz and xz planes, respectively. The
shear tractions must be known, and may be due to friction between the tool and the cylin-
der. In the current analysis friction was ignored, so 1251 = k52 = k4, = 1:42 = 0.

The nonlinear shell theory based on the above displacement field allows for trans-
verse shearing effects and a nonlinear distribution through the thickness of deformations
due to bending. Such higher order effects are usually neglected in the analysis of tirin-

walled structures. However, these effects may be important for nonlinear deformations,

125

and the inclusion of these effects allows the satisfaction of shear tractions on the top and
bottom surfaces of the shell, yielding a more accurate physical description of the problem.
This element, once formulated, could be used in the current analysis to provide a very
accurate, efficient, and robust model for thin cylinders undergoing axisymmetric deforma-
tions.

The use of the new element would remove the current barriers to reaching a con-
verged solution, but does not allow for the prediction of wrinkling in the swaging process.
To accomplish this, a method similar to that proposed by Adams (1993) could be
employed. Once a converged solution to the axisymmetric forming problem has been
found the analysis could be rerun, seeking the onset of wrinkling. This could done by
checking the solution after each stage of the forming process, to determine if circumferen-
tial buckling has occurred. Tracing the determinant of the characteristic system of equa-
tions for cylindrical shells to determine when (or if) it goes to zero would indicate that the
cylinder has buckled or wrinkled.

The characteristic system of equations could be formulated from the non-
axisymmetric form of the governing equations of the nonlinear shell theory discussed
above. The circumferential variations of the deformations, in this case, could be
approximated using a Fourier series method. In this method, the problem is analyzed for
several independent harmonics which are superimposed to give a final solution. Because
of the independence of the harmonics only those harmonics which influence the buckling
of the cylinder need be analyzed. Use of this method is possible because the
circumferential direction (the rolling direction) is also a principal material direction, and
the cylinder geometry is axisymmetric. Therefore, each of the five degrees of freedom

may be represented by a Fourier series in y,

126

ny
uinPl. (x) cos?

Ms
.MN
I

u(x.y) = 2 17:,(x)COS'-'R-y- =

n=1

3
II
N.
~
ll
-

0° m 2
v(x, y) = 2 17;!(x)sin'% = 2 z PEP, (x) sin-{1R3
n = 1 n = Ii = I
co m 3
w(x, y) = 2 uTn(x)c0s’% = 2 2 EN, (x) cos? (5.3)
n = 1 n = 11': I

_ — 22 ._. — . '2’
‘11,“: 1’) - 2 Wxn(x)005 R WxinPl (x) cos R

n=I

11MB 3M5

ny
"Pl. (x) cos 7?-

MNEMe

:5

W,(x. y) = 2 \7; (x) 603% =

n=1

3
N
-
ll
1...

where m is the number of Fourier harmonics, P,- are the linear Lagrange interpolation
functions, and N,- are the quadratic Lagrange interpolation functions. In the above

equations the DOF are approximated in the circumferential direction by the product of a
trigonometric function and an amplitude term, with the axial variation of the amplitudes

approximated by the finite element method. Also of note is the approximation w“. In this

case, an interdependent interpolation scheme is employed to allow for greater accuracy
without an increase in computational cost. The third DOF in the element approximation of

Wu would be eliminated before assembly by imposing a constraint on the variation of

transverse shear strain (see Averill, 1994 and Tessler, 1983)

Due to the form of the displacement approximations in Equation (5.3), integration
in the thickness and circumferential directions can be performed analytically, so that the
model is reduced to one dimension (along the axis of the cylinder). The final model is

given as,

127

16.1: = p” (5.4)

J j

where n is the Fourier harmonic, K the global stiffness matrix, d the Fourier series
amplitude, and F the force matrix. The characteristic system of equations is obtained for
each harmonic to measure the stability of the current configuration. When the determinant
goes to zero for a given harmonic, buckling is predicted to occur, with a buckling pattern
wavelength related to n.

The development of the axisymmetric, nonlinear shell element is currently being con-
ducted by other researchers. Upon completion of the development, the element can be
implemented into the current model for solution of the axisymmetric forming problem.
The solutions from this analysis must then be coupled with the proposed wrinkling analy-

sis formulation to complete the solution.

APPENDIX

128

A.l Mesh Refinement Study Data

In Section 3.3.1, an elemental convergence study was conducted on the first stage

of the forming process. In this study a series of meshes were analyzed using the developed

model. From those analyses data was collected on the outer, mid-plane, and inner surfaces

of the cylinder at locations within the first bending region and the critical region. Tables

A.1 through A.24 contain all of that data collected.

Table A.l Axial Strain on the Outer Surface in the First Bending Region. (x10'3)

 

 

 

 

 

 

 

 

 

 

 

Elements 2 3 4 5 6 8
150 14.3964 15.2937 15.2429 15.2084 15.2017 15.2091
300 14.2752 15.6635 15.5345 15.6663 15.6150 15.7637
400 ----- ----- 15.6079 15.7318 15.8140 ' -----

 

 

 

Table A.2 Axial Strain on the Mid-Plane Surface in the First Bending Region. (x10'3)

 

 

 

 

 

 

 

 

 

 

 

Mesh 2 3 4 5 6 8
150 0.838835 1.233395 0.997778 0.98309 0.953442 0.923100=H
300 0.645666 1 .45242 1.04027 1 .00798 0.990998 0.884100
400 ----- mu 1 . 12246 1 . 158305 1 .080070 -----

 

 

 

 

 

129

Table A.3 Axial Strain on the Inner Surface in the First Bending Region. (x10'3)

 

 

 

 

 

Mesh 2 3 4 5 6 8
150 43.2124 43.1500 43.5015 43.4648 43.4725 43.5070
300 43.4868 43.0982 43.7462 43.9155 43.9076 44.1870
400 ---------- 43.6394 43.6624 43.8760 ------

 

 

 

 

 

 

 

 

Table A.4 Circumferential Strain on the Outer Surface in the First Bending Region.

 

 

 

 

 

(x10'3)

Mesh 2 3 4 5 6 8
150 -6.98879 -6.93845 -6.91135 -6.96230 -6.97252 -6.64920
300 -6.9025 8 -6.8201 1 -6.70422 -6.66099 -6.68490 -6.63974
400 ---------- -6.85317 -6.85937 -6.81014 -----

 

 

 

 

 

 

 

 

Table A.5 Circumferential Strain on the Mid-Plane Surface in the First Bending

Region. (x10'3)

 

 

 

 

 

 

 

 

 

 

 

Mesh 2 3 4 5 6 8
150 -6.92393 -6.87485 -6.84736 -6.899825 -6.90983 -6.91335
300 -6.83654 -6.75465 -6.63764 -6.59537 -6.6l930 -6.57391
400 ---------- -6.78859 -6.79609 -6.74601 ---—-

 

 

 

 

130

Table A.6 Circumferential Strain on the Inner Surface in the First Bending Region.

 

 

 

 

 

(x103)

Mesh 2 3 4 5 6 8
150 6.89687 6.84650 -6.82324 6.87584 -6.88715 6.89121
300 6.80847 6.72467 6.68394 6.56921 6.59467 6.54976
400 ----- 6.76480 6.77228 6.72332 -----

 

 

 

 

 

 

 

 

Table A.7 Axial Stress on the Outer Surface in the First Bending Region. (szi)

 

 

 

 

 

 

 

 

 

 

 

 

Mesh 2 3 4 5 6 8
150 240.268 187.504 158.886 138.823 126.935 111.188
300 240.424 191.156 162.397 144.685 132.495 115.812
400 ---------- 157.724 140.057 128.731 -----

 

Table A.8 Axial Stress on the Mid-Plane Surface in the First Bending Region. (szi)

 

 

 

 

 

 

 

 

 

 

 

 

Mesh 2 3 4 5 6 8
150 ' 0.17758 -4.8125 W 8.37254 -10.2651 -11.9106
300 1.14157 -4.7234 -65.3091 -8.6686 -10.8890 -13.7535
400 ---------- 45.8308 -6.54295 -8.965.84 ------

 

 

 

131

Table A.9 Axial Stress on the Inner Surface at in the First Bending Region. (szi)

 

 

 

 

 

 

 

 

 

 

 

 

Mesh 2 3 4 5 6 8

T50 -268.912 -208.256 476.320 453.901 440.543 422.824
300 -271.538 -213.056 480.558 459.968 446.372 427.406
400 ---------- 471.748 452.187 439.352 -----

 

 

 

Table A.10 Circumferential Stress on the Outer Surface in the First Bending Region.

 

 

 

 

 

 

 

 

 

 

 

 

 

(szi)

Mesh 2 3 4 5 6 8
150 155.628 106.460 787539—385953— 47.0817f 31.7309
300 156.503 109.640 81.7260 64.5423 52.8537 362061
400 ----- 75.4760 58.2574 47.2932 -----

 

 

 

Table A.11 Circumferential Stress on the Mid-Plane Surface in the First Bending

 

 

 

 

 

Region. (szi)
Mesh 2 3 4 5 6 8
== 1: 3;:
150 -55.8569 -56.7963 -62.6269 -6l4201 -64.9001 -657393
300 -58.4541 -58.0932 -65.3584 -63.0065 -67.8407 -69.l656
400 ---------- 65.8639 63.9364 68.2242 --..-

 

 

 

 

 

 

 

 

 

132

Table A.12 Circumferential Stress on the Inner Surface in the First Bending Region.

 

 

 

 

 

 

 

 

 

 

 

 

(szi)

Mesh 2 3 4 5 6 8
150 -256.441 — 198.253 —165.400 - 145.595 -131795 -114095
300 -258.999 -199.263 -167.479 - 147.434 -132.657 -116.286
400 ----- ----- - 172.804 7151.758 - 138.334 ----

 

 

Table A.13 Axial Strains on the Outer Surface in the Critical Region. (x10'3)

 

 

 

 

 

 

 

 

 

 

 

Mesh 2 3 4 5 6 8
150 -4. 18674 -1.04515 -0.278518 0.101525 0.336483 0.644162
300 -.35759 1.75286 0.14918 0.551638 0.517526 0.929506
400 ----- ----- -0.333748 0.122887 0.424808 -----

 

 

 

Table A.14 Axial Strains on the Mid-Plane Surface in the Critical Region. (x10'3)

 

 

 

 

 

 

 

 

 

 

 

 

 

Mesh 2 3 4 5 6 8
150 17.8305 1W T9696 :- 18.9899 19.0050 18.9928=
300 18.2751 23.7753 18.8816 18.8507 19.3794 19.2492
400 ----- ----- 19.2907 19.3346 19.3919 ----

 

 

 

133

Table A.15 Axial Strains on the Inner Surface in the Critical Region. (x10‘3)

 

 

 

 

 

 

 

 

 

 

 

 

Mesh 2 3 4 5 6 8
150 38.8382 38.3533 37.7740 37.6701 37.6560 37.5003
300 39.1136 44.6235 37.4199 37.1498 38.2427 37.8484
400 ---------- 38.5379 38.4350 38.4716 -----

 

Table A.16 Circumferential Strains on the Outer Surface in the Critical Region.

 

 

 

 

 

 

 

 

 

 

 

 

(x103)

Mesh 2 3 4 5 6 8
150 62.3648 -52.3186 -52.—3061 -52.3227 62.3256 62.3274
300 62.3355 62.4560 61.7063 61.6914 63.7710 62.3589
400 ----- 62.2466 62.2544 62.4124 -----

 

Table A.17 Hoop Strains on the Mid-Plane Surface in the Critical Region. (x10'3)

 

 

 

 

 

 

 

 

 

 

 

 

Mesh 2 3 4 5 6 8
150 -52.4591 -52.4024 -52.3916 -52.4067 -52.4104 -52.4121
300 -52.4281 -52.5357 -51.7063 -51.7540 -52.4607 -52.4428
400 ---------- -52.3256 -52.3312 -52.4950 -----

 

 

 

 

134

Table A.18 Hoop Strains on the Inner Surface in the Critical Region. (x10'3)

 

 

 

 

 

 

Mesh 2 3 4 5 6 8
150 -52.4937 -52.4353 -52.4211 -52.4368 -52.4393 -52.4409
300 -52.4615 -52.5579 -51.7774 -51.7640 -52.4881 -52.4710
400 ----- ----- -52.3482 -52.3537 -52.5213 -----

 

 

 

 

 

 

 

 

 

Table A.19 Axial Stresses on the Outer Surface in the Critical Region. (szi)

 

 

 

 

 

 

Mesh 2 3 4 5 6 8
150; -360.213 -246.750 -197.342 -171.537 -153.706 -127.704
300 -354.848 -232.841 —208.305 - 177.613 -155.114 -130.090
400 ----- ----- -216.718 -186.542 - 164.428 -----

 

 

 

 

 

 

 

Table A.20 Axial Stresses on the Mid-Plane Surface in the Critical Region. (szi)

 

 

 

 

 

 

 

 

Mesh 2 3 4 5 6 8
I = r I
150 12.007 4.2907 2.18321 1.0688 -0.00944 -2.0748
300 12.5069 4.06294 5.9629 4.2178 1.6014 0.0673
400 ---------- 6.5923 4.3753 3.7378 -----

 

 

 

 

 

 

 

 

135

Table A.2] Axial Stresses on the Inner Surface in the Critical Region. (szi)

 

 

 

 

 

 

 

 

 

 

 

 

Mesh 2 3 4 5 6 8
150 335.940 226.770 180.082 157.769 140.193 117.602
300 331.314 266.730 197.465 168.970 144.993 122.361
400 ---------- 201.974 174.301 154.099 -----

 

Table A.22 Circumferential Stresses on the Outer Surface in the Critical Region.

 

 

 

 

 

 

 

 

 

 

 

 

(szi)

Mesh 2 3 4 5 6 8
150 -394. 100 -284.650 -233.855 -204.445 - 185.684 - 160.038
300 -397.440 -282.371 -235.476 -204.863 -189.751 - 162.522
400 ----- ----- -238.453 -209. 154 - 188.076 -----

 

Table A.23 Circumferential Stresses on the Mid-Plane Surface in the Critical Region.

 

 

 

 

 

 

(szi)
Mesh 2 3 4 5 6 8
150‘ 68.99% 63.8572 68.5703 67.9680 -70.7456 20.9200
300 69.7367 20.7880 67.0310 66.2861 21.3535 22.1101

 

400

 

 

 

 

 

 

 

 

 

 

 

136

Table A.24 Circumferential Stresses on the Inner Surface in the Critical Region.

 

 

 

 

 

 

(szi)

Mesh 2 3 4 5 6 8
150 253.175 155.333 109.250 83.8303 66.7538 44.0962
300 247.643 204.532 119.618 91.9736 65.4983 43.4282
400 ---------- 123.298 96.5859 76.1358 -----

 

 

 

 

 

 

 

 

137
A.2 Analysis Results without Strain Hardening

The effects of strain hardening on the analysis results were evaluated in Section
3.4.2. It was determined at that time the effects were minimal at most. This may be con-
formed by comparing the plots that follow in Figure A.1 - A.4 to those of Figures 3.23 -
3.26. Doing so will show that while the magnitudes of the strain and stress values have
small incremental changes, the trends remain constant. Therefore, strain hardening effects

are categorized as minimal.

138

 

 

 

 

 

you

 

 

Incend‘l’ehl“ I

 

Figure A.1 Axial Strains in the Critical Region Without Work Hardening.

139

 

 

 

 

 

 

4‘1.“

4“

 

 

 

 

‘1 ado-garment. 1

Figure A.2 Circumferential Strains in the Critical Region Without Work Hardening.

 

 

 

 

 

 

 

I‘M“... 1

 

Figure A.3 Axial Stresses in the Critical Region Without Work Hardening.

141

 

 

 

 

 

1.51

“001.01“ 1

 

 

 

Figure A.4 Circumferential Stresses in the Critical Region Without Work Hardening.

REFERENCES

142

Adams, G. G., (1993),“Elastic Wrinkling of a Tensioned Circular Plate Using von Karman
Plate Theory”, Journal of Applied Mechanics, 60. pp. 520 - 525

Averill, R. C. and Reddy, J. N., (1990), “Behaviour of Plate Elements Based on the First
Order Shear Deformation Theory”, Engineering Computing, vol. 7, pp. 57-74

Averill, R. C., (1994) “Static and Dynamic Response of Moderately Thick Laminated
Beams with Damage”, Composites Engineering, vol. 4, pp. 381-395

Butman, J ., (1991),Car Wars, Grafton Book, London, p 11

Chan, S. L., (1993), “A Non-Linear Numerical Method for Accurate Determination of
Limit and Bifurcation Points”, International Journal for Numerical Methods in
Engineering, vol. 36, pp. 2779 - 2790

Dym, C. L., (1974), Introduction to the Theory of Shells, Pergamon Press, Oxford, p.
26.

Eterovic, A. L. and Bathe, K. J ., (1991), “On the Treatment of Inequality Constraints
Arising form Contact Conditions in Finite Element Analysis”, Computers and
Structures,vol.40, pp. 203-209

Green, C. M., (1956), Eli Whitney and the Birth of American Technology, Little,
Brown and Company, Boston, p. 120

Kalpakjian, S., (1991), “Manufacturing Processes for Engineering Materials”, Second
Edition, Addison-Wesley Publishing Company, Reading, MA, pp. 461-462

Keck, P., Wilhelm, M., and Lange, K., (1991), “Application of the Finite Element Method
to the Simulation of Sheet Forming Processes: Comparison of Calculations and
Experiments”, International Journal for Numerical Methods in Engineering, vol. 30,
pp. 1415-1430

Lee, J. K., Cho, U. Y., and Hambrecht, J., (1991), “Recent Advances in Sheet Metal
Forming Analysis”, Advances in Finite Deformation Problems in Materials
Processing and Structures, ASME 1991, AMD vol. 125

MARC User’s Guide Volume A: User Information, (1994), MARC Analysis Research
Corporation, Palo Alto, CA

Mendelson, A., (1968), Plasticity: Theory and Application, Robert E. Krieger
Publishing Company, Malabar, FL.

Mentat H Command Reference, (1994), MARC Analysis Research Corporation, Palo
Alto, CA

143

Neale, K.W. and Tugcu, P., (1990), “A Numerical Analysis of Wrinkle Formation
Tendencies in Sheet Metals”, International Journal for Numerical Methods in
Engineering, vol. 30, pp. 1595-1608

Padovan, J., (1974), “Quasi-Analytical Finite Element Procedures for Axisymmetric
Anisotropic Shells and Solids”, Computers and Structures, vol. 4, pp. 467 - 483

Parr, G., (1958), Man, Metals and Modern Magic, American Society for Metals,
Cleveland and The Iowa State University Press, Ames, Iowa

Rebolo, N ., Nagetaal, J. C., and Hibbitt, H. D., (1990), “Finite Element Analysis of Sheet
Forming Processes”, International Journal for Numerical Methods in Engineering,
vol. 30. PP. 1739 - 1758

Reddy, J. N., ( 1984), Energy and Variational Methods in Applied Mechanics, John
“Wiley and Sons, New York

Riccobno, R., (1992), “A Comparison Between Two Different Models for Metal Forming
FEM Analysis”, Computers and Structures, vol. 44, pp. 429-433

Sedaghat, M. and Herrmann, L. R., (1983), “A Nonlinear, Semi-Analytical Finite Element
Analysis for Nearly Axisymmetric Solids”, Computers and Structures, vol. 17, no. 3,
pp. 389 - 401

Simo, J. C., Wriggers, P., and Taylor, R. L., (1985), “A Perturbed Lagrangian Formulation
for the Finite Element Solution of Contact Problems”, Computer Methods in Applied
Mechanics and Engineering, vol. 50, pp. 163 - 180

Tessler, A. and Hughes, T. J. R., (1983), “An Improved Treatment of Transverse Shear in
the Mindlin-Type Four-N ode Quadrilateral Element”, Computer Methods in Applied
Mechanics and Engineering, vol. 39, pp. 311-335

Toh, C. H., (1989), “Prediction of the Forming Limit Curves of Sheet Materials Using the
Rigid-Plastic Finite Element Method”, International Journal of Machine Tools and
Manufacture, vol. 29, pp. 333-343

Wood, M. D., (1994), “The Effect of Microstructure and Texture on the Mechanical
Properties of Double-Reduced Tinplate Steel”, Master’s Thesis, Michigan State
University, East Lansing, Michigan

Zhang, L. C. and Yu T. X., (1988), “The Plastic Wrinkling of an Annular Plate Under
Uniform Tension on Its Inner Surface”, International Journal of Solid Structures, vol
24, no. 5, pp. 497 - 503

Zienkiewicz, O. C. and Taylor, R. L., (1989), The Finite Element Method, McGraw-Hill,
New York