.«MN Vﬂﬁm‘mﬂw . V . .V V . . V , .. . V ... . . , .. .V . V I “1.14.55.13- V‘

VV 932‘.» MW: . . V V V . . V . a . ﬁfﬁs§ VV .r. V

V ﬁr. , . . _ . . . _ V . V . ”mwﬂaamﬁ.
V . . V . I V N I. V Vf.\.«..

is“

:5 . . . V .V . , . V am;
. . . tV:-m
, _

. .
maﬁa. V .
H! . V : V3. . V .
. .mmém V , V V an“
, . . c . g.

2:... .

1.. Va.
Jainnmr:
1%....
1.13%.. «v :.
. . w A. a V
r... :m Vilngmvx

t...
payout
V. ,1

ud’u) x

. I):

V _ . . , . . . E. Ev? . V.
vs.“ ”mun. x , , z . V V . VV V V . V ., V . Ream. w§.?mv 5%“?
5‘. 1 1m. . . . V . . . V V V .3. .wm.%u1 ﬁwMWmWWWTZﬁ? VVoV
. :n . \ (“3.1V V V. v V V . V V , . . , . «W “"4;
V129, 3.? a V . £9 V . . ,. . . . . 3.. V V V Emu,
V . Einawwuaawpﬂwﬁ , V . . V , . V .. . V, .V V .umww . ﬁn%mﬁé x;
M1. , . V . . V , V . . x... .mMmug V

.r

. f...»

t. V
. 5:
1?

. Asa...
u. .L. «pvmx;

.‘l- i.

. : .é. ..

IV: 12:.

t

5
. pV
.V 5;. ﬂyﬂthr.
of. . . )V. .5.
. . 4.

V20“
9% ; 435$
7. it if}. u.
inn“? $.4th it
. Jun... \
I
V . 3.13 . .mnvnu
‘2.ng 3&5...
93527 . »

:5: V

, .. .329
mug?

-r t.)
)an a...
r. !V ,

 

 

u .3 .7 .
>f|\
neiv.

 

. ‘ ,.. {.5»...P
V1. .VV‘v. I. .u z}. _

 

 

THESIS

l
.2 cm

This is to certify that the

thesis entitled

A CROSS SECTION ANALYSIS FINITE ELEMENT CODE FOR TUBE
HYDROFORMING

presented by

Y abo Guan

has been accepted towards fulﬁllment

of the requirements for
Masters Mechanical Engineering
Jemee in

yo.

 

G.\@Wwo&/\

Major professor

Date ’2/” /2000

 

0-7639 MS U is an Afﬁrmative Action/Equal Opportunity Institution

 

 

LIBRARY

Mlchlgan State
Universlty

 

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before .date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

11m m.“

A CROSS SECTION ANALYSIS FINITE ELEMENT CODE
FOR TUBE HYDROFORMING
By

Yabo Guan

A THESIS
Submitted to
Michigan Sate University
In partial fulﬁllment of the requirements
For the degree of
MASTER OF SCIENCE
Department of Mechanical Engineering

2000

ABSTRACT

A CROSS SECTION ANALYSIS FINITE ELEMENT CODE
FOR TUBE HYDROFORMING

By

Yabo Guan

Hydroforming of metal components is widely used in several industries. While
numerous variations of hydroforming exist, the basic principle remains the same: utilize
ﬂuid pressure to form a component. In tube hydroforming, interrml pressure is applied to
form a tubular metal blank into a structural component having a closed cross section.

The thesis presents a "cross-section" analysis method for calculating tube
stretching and bending followed by hydroforming and expansion to ﬁll the die in the
hydroforrning process. The material is assumed to be elastic-plastic and to satisfy the
plasticity model that takes into account rate-independent, work hardening and normal
anisotropy. Based on the virtual work principle and constitutive equation, the equilibrium
equation can be derived. Numerical solution obtained ﬁ'om a Total-Lagrangian
formulation of a general shell theory is employed. The boundary friction condition is
introduced into the formulation in the form of penalty ﬁinction, which imposes the
constraints directly into the tangent stiﬂ‘ness matrix. The Newton-Raphson algorithm is

used to solve the nonlinear equation.

ACKNOWLEDGMENTS

The author wishes to express his appreciation to the following people and institutions for
making the completion of his Master of Science possible.

To my major professor, Dr. Farhang Pourboghrat, for his help, patience and guidance,
without which this project would not have been possible.

To Ford Company and MSU Manufacturing Research Consortium, for their ﬁnancial
support of my master studies.

To Haiyan Miao, my wife and closest friend, for her support and encouragement, when I

need it most, and for help in preparing this manuscript.

TABLE OF CONTENTS

ACKNOWLEDGMENTS .......................................................................... iii
LIST OF CONTENTS .............................................................................. iv
LIST OF FIGURES ................................................................................. vi
INTRODUCTION .................................................................................. 1
CHAPTER 1
LITERATURE REIVEW ........................................................................... 4
CHAPTER 2
THIN SHELL THEORY ............................................................................ 8
Principal Stretch (Total Lagrangian) ....................................................... 8
Principal Curvature and Stretch (Updated Lagrangian) ................................ 11
CHAPTER 3
KINEMATICS ....................................................................................... 15
Straight Segment ............................................................................ 16
Circular Arc Segment ...................................................................... l7
Constraints .................................................................................. 20
CHAPTER 4
CONSTITUTIVE EQUATION ................................................................... 21
CHAPTER 5
CONTACT ALGORITHM ........................................................................ 23
Numerical Procedures for Contact ....................................................... 23
Detection of Contact ....................................................................... 25
Contact Tolerance ................................................................. 25
Normal Vector and Tangential Vector .......................................... 26
Condition for Contact Violation ................................................ 28
The Projection Algorithm ................................................................. 29
The Iterative Method .............................................................. 30
The Directive Method ............................................................ 30
Implementation of Constraint ............................................................ 3]
Friction Model .............................................................................. 32
Separation ................................................................................... 34
CHAPTER 6
NUMERICAL SIMULATION .................................................................... 35
Principle of Virtual Work ................................................................. 35
Newton-Raphson Algorithm .............................................................. 36

iv

Three Types of Modeling ................................................................. 37

Tube Bending Modeling ......................................................... 37
Pressure Loading Modeling ...................................................... 39
Frictional Contact Modeling ..................................................... 40
CHAPTER 7
PROGRAM STRUCTURE ....................................................................... 44
Input File .................................................................................... 44
Geometric, Numerical and Material Input ..................................... 44
Loading Input ...................................................................... 48
Program Flow Chart ........................................................................ 48
Output File .................................................................................. 51
CHAPTER 8
EXAMPLES ......................................................................................... 54
Bending Examples ......................................................................... 54
Free Expansion Case ....................................................................... 56
Contact Examples .......................................................................... 56
CHAPTER 9
CONCLUSIONS .................................................................................... 67
APPENDIX A ....................................................................................... 68
APPENDIX B ....................................................................................... 73
REFERENCES ...................................................................................... 75

LIST OF FIGURES

Figure 1.] Representative hydroforrned components and assemblies ......................... 2
Figure 1.2 Principal hydroforming tooling diagram ............................................. 2
Figure 2.] Shell element at original time to and current time 1
(Total Lagrangian) ...................................................................... 10
Figure 2.2 Shell mid-surface at reference time 0! and current time t
(Total Lagrangian) ..................................................................... 10
Figure 3.1 A rectangular tube and the assumed global coordinate system .................. 14
Figure 3.2 Geometry of a straight segment ...................................................... 18
Figure 3.3 Geometry of a circular segment ...................................................... 18
Figure 5.1 Contact tolerance ....................................................................... 25

Figure 5.2(a) The computational mesh

(b) The normal vector and tangent vector ............................................ 26

Figure 5.3 The associated die segment with contact node ..................................... 28
Figure 5.4 The projection algorithm .............................................................. 30
Figure 5.5 Coulomb ﬁiction model ............................................................... 32
Figure 6.1 A schematic of tool-workpiece check ............................................... 40
Figure 6.2 Flow chart for contact iteration ...................................................... 42
Figure 7.1 Flow chart of the program ............................................................ 49
Figure 8.1 Stretching/bending of a rectangular section ........................................ 57
Figure 8.2 Stretching/bending of a circular section ............................................. 57
Figure 8.3 Bending a tube to collapse ............................................................ 58
Figure 8.4 Strain and displacement result compared with ABAQUS ........................ 59

vi

Figure 8.5a A round aluminum tube was expanded into a square die

with sharp—comers ..................................................................... 60
Figure 8.5b Final shape from ABAQUS ......................................................... 60
Figure 8.63 Hoop strain distribution compared with experiment ............................. 61
Figure 8.6b Intermediate shapes of the tube during hydroforming ........................... 62
Figure 8.7a Circumferential strain distribution after bending ................................. 63
Figure 8.7a Circumferential strain distribution after bending ................................. 64

vii

INTRODUCTION

Hydroforming offers a way to cut material and manufacturing costs while
improving product performance in a variety of applications. Advantages include material
savings, making shapes that cannot be made by conventional manufacturing methods,
and the ability to shape dissimilar tailor-welded materials. It reduces tooling, part, and
labor costs, while signiﬁcantly improving product performance (Anon, 1997).

The process currently is either being used or being considered for making a wide
range of complex cylindrical, ﬂat and tubular components. Cylindrical parts include gas
cylinders, washing-machine drums, and cooker cavities. Flat components include auto
body panels, fuel tanks, and aerospace parts.

Tube hydroforming is receiving the greatest attention, especially in the auto
industry (V ari-Form Inc, 1998), because existing multi-piece, stamped/welded assemblies
in auto body and frame structures potentially can be replaced with less expensive
hydroformed parts that are lighter, stronger and more precise (Figure 1.1).

In tube hydroforming, internal pressure is applied to form a tubular metal blank
into a structural component having a closed cross section. As shown in Figure 1.2, a
straight or prebent tube is placed in an open die, the die is closed, the tube ends are
sealed, and the tube is ﬁlled with a hydraulic ﬂuid. Internal hydraulic pressure forces the
tube to conform to the shape of the die cavities, producing a part having different cross-

section shapes.

BOI

 

Figure 1.1 --- Representative hydroformed components and assemblies
From left to right: instrument panel support beam, engine-cradle assembly.
Photo courtesy Van-Form, Woodstock, Ontario, Canada.

.Punch plate centering

   
 
 

_' Top tool
Top die Base plate
Horizontal cylinder

_ Pivotable
Cylinder bracket

Sealing punch
Bottom die

Bottom tool
Base plate

Table plate ASE component

Figure 1.2 --- Principal hydroforming tooling diagram
Courtesy Schafer Hydroforming, Wilnsdorf, Germany.

The thesis presents a "cross-section" analysis method for calculating tube
stretching and bending followed by hydroforming and expansion to fill the die in the
hydroforming process.

The present method employs a number of simplifying assumptions. The principal
geometrical assumption is that the cross-section of the extrusion can be of any shape
consisting of straight-line segments and circular arcs. All segments making up the cross-
section are assumed to be relatively thin and of constant thickness. The deformation of
the member is assumed not to vary along the length. Hence, the analysis can be
considered to be a "section" analysis.

The material is assumed to be elastic-plastic and to satisfy the plasticity model
that takes into account rate-independent, work hardening and normal anisotropy
(Pourboghrat et a], 2000).

At the tool-tube interface, the frictional stress is assumed, on the basis of ordinary
Coulomb friction, to be proportional to the contact pressure whenever relative sliding
occurs.

By using the virtual work principle and constitutive equation, the equilibrium
equation can be derived. Numerical solution obtained from a Total-Lagrangian
formulation of a general shell theory is employed. The boundary friction condition is
introduced into the formulation in the form of penalty function, which imposes the
constraints directly into the tangent stiffness matrix. The Newton-Raphson algorithm is

used to solve the nonlinear equation.

 

 

h}
kn
an

en,

prc

Bile!

desi‘:

CHAPTER 1

LITERATURE REVIEWS

1.1 Hydroforming Concept

Hydroforming has been one of the fundamental sheet metal forming processes for
quite a long time, having been developed at least since pre-World War H, when it was
applied to the German aircraft industry (Zhang, 1997). About 30 years ago, the first
hydroformed parts were fabricated by expanding a straight or prebent tube to make
manifold elements and components of similar geometry for sanitary use. Today, many
hydroformed parts are already in mass-production in the automotive industry. Well-
known hydroformed automotive applications include exhaust manifolds, exhaust pipes
and subframes, such as radiator enclosures, space frames, dash assemblies, frame rail,
engine cradles, etc.(Dohman et al, 1997)

According to its features, hydroforming can be classified as the following
processes (Zhang, 1999 and Siegert et al, 2000):

1. Flat sheet hydroforming or hydroforming deep drawing

2. Tube hydroforming or hydrobulging of tubular components

3. Hydromechanical deep drawing process

4. New variation of hydroforming including the integral hydrobulge forming

(IHBF) of sheet shell products and viscous pressure forming (VPF) processes

In this thesis, we focus on tube hydroforming that has become an economic

alternative to various stamping processes with the availability of advanced machine

design and controls.

 

IC

fe

Tube hydroforming offers several advantages as compared to conventional
manufacturing via stamping and welding (Ahmetoglu et a1, 2000). These include:

(a) part consolidation (stamped and resistance welded two or more pieces of a

box section can be manufactured in one operation from a hollow component)

(b) weight reduction through more efficient section design and tailoring of the

wall thickness

(c) improved structural strength and stiffness

((1) lower tooling cost due to fewer parts

(e) fewer secondary operation (no welding of sections required and holes may be

punched during hydroforming)

(f) tight dimensional tolerances and low springback

(g) reduced scrap

While numerous variations of hydroforming exist, the basic principle remains the
same: utilize ﬂuid pressure to form a component. Both traditional high-pressure
hydroforming and pressure sequence hydroforming are presented in Longhouse (1998)
and Morphy (1997).

If hydroforming technology is to be applied economically, it is essential to have
knowledge of the avoidance of failure cases as well as of the behavior of the tube in the
tool under the compressive stress and forces that are exerted by the machine (i.e., axial
feed). Several failure modes, buckling, wrinkling and bursting are shown in Dohmann

(1996).

1.2

dra‘

£155

inc

6X2

W3

me

thi

CO

an:

ph

of

of

 

 

dc]

fOr

(19

W11

fOn

 

1.2 Numerical Simulation

Kobayashi (1989) has worked on several sheet forming problems including
drawing, bulging, stretch forming and bending, solving the problems using a ﬁnite
element method based on incremental strain theory. In his work, the material was
assumed as rigid-plastic and anisotropy and work hardening characteristics were
included. An elastic-plastic FEM formulation of stretch forming for a punch and a die of
arbitrary shape was introduced by Wang and Bundiansky (1978), which is a classic
example veriﬁed by ABAQUS. A rigid-viscoplastic material model FEM formulation
was presented by Germain, Chung and Wagoner (1989).

All the numerical analyses mentioned above were obtained using simple
membrane theory which does not include the variation of stresses through the sheet
thickness. This approximation is justiﬁed when the thickness of sheet material is small
compared to the curvature of the sheet surface, the radius of either the punch or the die,
and when in-plane stretch dominates bending. However, in order to simulate bending
phenomena such as the early deformation in stretch forming, spring-back etc., inclusion
of the bending stress is necessary. A total Lagrangian formulation based on an extension
of a thin shell theory which takes into account large membrane and bending strain was
derived by Wang and Tang (1988). An updated Lagrangian and plane-strain bending
formulation based on nonlinear thin-shell theory was employed by Choudhry and Lee
(1994). A general purpose ﬁnite element program for analysis of forming processes,
which was implemented in MARC, was introduced by Nagtegaal and Rebelo (1988).

Built on the successful application of ﬁnite element method on sheet metal

forming, explicit ﬁnite element method is being applied to hydroforming of tubular parts

by Ni (1994), bulging process by Hu et a1. (1997), and a sheet metal operation involving
punch stretching, drawing and hydroforming by Hsu and Chu (1995). In Noh and Yang
(1998) work, a general formulation for hydroforming of arbitrarily shaped boxes was
expressed.
1.3 Friction and Contact Problem

Friction and contact between sheet and die play an important role in hydrofornring
processes. Classical Coulomb’s law is capable of describing only friction effects between
effectively rigid bodies and gross sliding of one body relative to another, However if
Coulomb’s law is applied pointwise in contact problem, then the contact stress developed
normal to the contact surface is ill deﬁned. In Oden and Pires (1983), nonlocal and
nonlinear friction laws and variational principles for contact problems were introduced.
A method of discretizing the die boundary condition considered for the analysis of metal

forming processes with arbitrary shaped dies was employed by Oh (1982).

 

IN)

I(:

the

Wt.

the

In}

CHAPTER 2
THIN SHELL THEORY

The deformation of the mid-surface of an element in the plane containing the meridian
curve will be considered. The deformation of the element will be obtained from the
theory of thin shells and Figure 2.1-2.2 show the mid-surface of an element of the sheet at
the initial time to, reference time °t and at the current time t as it bends and stretches.
Derivation of principal curvatures, stretches of a shell element undergoing an
axisymmetric deformation was discussed in (Pourboghrat, F. et a1, 2000) using both total
and updated Lagrangian formulations. These formulations are used in the program.

2.1 Principal Stretch (Total Lagrangian Formulation)

Figure 2.1 shows a ﬂat element at the initial time t0 and at the current time
t(=°t+At) as it bends and stretches. By referring to Figure 2.1, the following relationship
between the current and the initial position vectors and the total displacement vector for
the mid-surface of the element can be written:

x = X + U (2-1)
where the initial position vector X and the total displacement vector U are, deﬁned in
the global coordinate system,

X—X‘ 11—" (22
~—O’ -—w .)

In Eq.(2.2), u and w are the total displacements of the mid-surface. By substituting from

Eq.(2.2) into (2.1) , the position vector x becomes

[X1 +u]
x: (2.3)

W

To calculate the principal stretch in the radial direction, we will consider the unit tangent

vector, (i , of the current mid-surface, i.e.

dx dx dx 1+“
= ~ ~d5 ~_1_=[ “ii (2.4)

where A? = ds/dS . Here, 3 and S are travel length of the current and initial surfaces from
center, respectively, and the comma stands for the differentiation with respect to X,.

Then, by calculating the magnitude of vector ii from Eq.(2.4) and setting it equal to
unity, A? can be calculated to be:

"a" = a- a = 1.0 (2.5)

NI—

13 = [(1+ 14.x, )2 + wij (2.6)

 

NI

 

 

 

 

0
X2
A d8 time t
U

” w

z

u
F 1 § X1
dS
time to

Figure 2.1 Shell element at original time to and current time t(Tota1 Lagrangian)

 

 

E
ii 6 ,
ds trme t
4""-~....
Aw Au "~ -.
.0.
92 '° w
. u, ..... :7"
A .............
dS time ° t
81 ...............
r

 

Figure 2.2 Shell mid-surface at reference time °t and current time «Updated Lagrangian)

IO

2.2 Principal Curvature and Stretch (Updated Lagrangian)

After bending and stretching, the principal mid-surface curvature, k1 , of a shell

element at the current conﬁguration, t(=°t+At), can be calculated from the known

information about the element at the reference conﬁguration (time °t, see Figure 2.2), as

follows:

r=R+uA+wIQ

(2.7)

where R is the reference conﬁguration at time °1: and u and w are incremental

displacements deﬁned in Figure 2.2. In Eq.(2.7), the unit tangent vector A and the unit

principal normal vector If] to the mid-surface of the reference are deﬁned as

r>)
ll
mIQJ
%rz
ll
'70

l>)

.S

10:

W |

I

where K1 is the curvature. Therefore,
10. A = o

and

11

(2.8)

(2.9)

(2.10)

(2.11)

Is

To calculate the current curvature, k], the unit tangent vector ‘ = and the

I D

 

 

 

 

 

 

 

 

n
unit principal normal vector of the mid-surface of the shell " = 4 are calculated.
n
Using Eq.(2.7), the tangent vector 0 is found as:
a r ,. A . ,.
9=E=f5 :5,+u~5 A+uA'S+w'S 1Y+w1§l5 (2.12)

By substituting from Eqs.(2.8), (2.9) and (2.11) into (2.12), and re-arranging, the

following expression results:
a=—”—= S =(l+u.s—K1w)A+(w.S+K,u)1il=cA+d19 (2.13)

The principal incremental stretch of the mid-surface in the radial direction calculated

from the magnitude of the base vector a in Eq.(2.13) is

’11:

 

 

 

 

a

= [9-9 =ch +212 =[(l+u.s —1<,w)2 +(w5 +K,u)2]i (2.14)

The current length of the mid-surface of the shell in the radial direction, ds, is calculated
from the reference length, dS, and A; as follows:
('15 = Ads 2.15)

The unit principal normal vector of the surface of the current shell, ii , is

73 =-———:———"- 2.16)

12

which, from Eqs.(2.l3), (2.14) and (2.16), shows that pic? =0. The current principal
curvature of the shell, k1, can now be found as:

~r'i =-r ~ri =——a-ii (2.17)
where a is given by Eq.(2.12) and ti can be derived from Eq.(2.l6):
~ ~ .s

as 2, (d, 21+ c, 10- K,c21— [WIND—Au -(-d A+ch7)
—~= * ~ ~ 2 ~ ~ ~ (2.18)
.s S ,1]

 

l3>

In Eq.(2.18), xi” is assumed to vanish within an element and the above expression

simpliﬁes as :

(at, + ch) 21+ (—c,, + [(4)10
'5 = - A (2.19)

By substituting from Eqs.(2.13), (2.14) and (2.19) into Eq.(2.17), the current centerline

 

13>

curvature of the shell, k], can be found:
Cd's -dc.s + K1212
= a:

It should be noted that it is possible to recover Eq.(2.6), for total Lagrangian formulation,

k, (2.20)

 

From Eq.(2.14). This is done by setting the reference curvature K1 equal to zero for a ﬂat
sheet, replacing S with X; and treating the incremental displacement u and w as total

displacements in Eqs.(2.13) and (2.14). Similarly, the current curvature of the shell, k1,

can be calculated from the initial conﬁguration by setting K1 equal to zero, replacing A,
with if and S with X1 and treating u and w as total displacements in Eq.(2.20)

_[_w'xlu-X1x1 + (1+u.X. )M’Xﬂll ]

k, 3 2.21
(2:) ( ’

 

l3

CHAPTER 3
KINEMATICS

Based on the general thin shell formulations in Chapter 2, the formulations for
ﬁnite strain and small rotation have been derived for the straight line segment and are
segment that are used in the program. (Derivation in details can be seen in Appendix A).

The formulation in Corona and Pourboghrat (1998) focuses on a closed structural
member, such as a thin-walled hollow tube, that is subjected to combined tension and
bending loadings. The cross-section of the tube can be of any closed or open shape

consisting of straight-line segments and circular arcs (see Figure 3.1).

 

Arc Element
f \/

 

Y { d—Straight Element
Bending Surface

3.

 

 

 

 

L J
////vZ////

Figure 3.1. A rectangular tube and the assumed global coordinate system.

For the purpose of modeling, all segments making up the cross-section are
assumed to be thin and of constant thickness. The length of the straight segments should
be greater than ten times the thickness, while the radius of the circular arcs should be

greater than eight times the thickness. The deformation of the structural member is

14

assumed not to vary along its length (in X-direction). Hence, the analysis can be done for

a section along the axial dimension, as shown in Figure 3.1.

3.1 Straight Segment

The rectangular cross-section shown in Figure 3.1 lies in the plane deﬁned by the
Y and Z-axes, while the X—axis is in the axial (length) direction of the member. This ﬁgure
indicates that the curvature K2 of a tube bent about the Y-axis will match that of the
bending surface. A simultaneous bending about the Z-axis is also possible in order to
impart an additional curvature K, to the tube. In this formulation, the initial geometry of
the cross-section is speciﬁed by ﬁrst deﬁning the Y-Z coordinates of a set of nodes and
then deﬁning the segments. Two nodes must bind each segment.
For the straight segment in Figure 3.2, the local coordinates of the segment are 5 along
the segment and z in the through-thickness direction. The initial segment geometry is
deﬁned by the location of the points = 0 given by Y0 and Z0 , the orientation angle 6, the
length 10 and the thickness t of the segment. The position 2:0 is located at the mid-surface
of the segment, indicated by the dashed line. The displacement components of the mid-

surface are u and w along the original 3 and 2 directions respectively. The rotation of a

line is given by B
B = —— = —w. (3.1)

The axial strain component of a point is given by:

8,, = £3 + K25 — Kyn (3.2)

where 81’ is the axial strain (axial feed) along the length of the tube (X-direction) and

15

f = Z0 +(s+u)sin0—wcos€+z(-cost9+ﬁsin 6) (3.3a)
17 = Y0 +(s+u)cost9+wsin0+ z(sin6+ﬂcos6) (3.3b)

The strain components, along .9 is given by
83 = a? + K32 (3.4)

where £3 is the membrane strain component given by
0 1 2 1 2
a, = u: + 2“" + 2w" (3.5)

where the local curvature is given by

W

K, = ——‘”2—, (3.6)
(1+ wJ )3

3.2 Circular Arc Segments
For the circular arc segment shown in Figure 3.3, the local coordinates of the

segment are the angle ¢ and the through-thickness z. The initial geometry of the
segment is speciﬁed by 1) the coordinates of the center of the arc, Y, and Zc, 2) the
angle (00 , which locates the line where (I) = O , 3) the radius R, 4) the angle spanned by the
arc A¢ , 5) and the thickness of the segment t. As in the straight segment, the mid-surface

is the origin for z and is indicated by the dashed line. The rotation of a line, ,8 , can be

 

given as
dw
“—25 u—w
.a
= = 3.7
16 R R ( )

The axial strain of a point is given by Eq.(3.2), but for an are

16

f = Z( +(R+w)sin$+ucos¢7+ z(sin(i)_+ﬂcos$)
77 = Y( +(R+w)cos(5—usin¢7+z(cos$—,Bsin(3)
where 5 = ¢o +¢ . The tangential strain is given by
8, =62+K,z

where

 

 

 

2 2
u +W u-w w+u
6:): .¢ +_1_ .¢ +1 .¢
R 2 R 2 R

Cd -dc +2R
and IQ: " A; 11/

 

where c=l+%(u" —w), d=-I-l€(w" +14) and xi, =\/c2 +612

17

(3.8a)

(3.8b)

(3.9)

(3.10)

(3.11)

 

 

 

 

 

 

VZ

Figure3.2 Geometry of a straight segment

 

 

 

 

 

 

 

Figure3.3 Geometry of a circular arc segment

18

3.3 Constraints
Since the cross—section is made of several independent segments which have their
own local coordinate systems and variables, it is necessary to enforce compatibility of
deformations at junction of two or more segments. Two constraint equations are used to
ensure compatibility of displacements and one to ensure compatibility of rotation
between two segments. Therefore, at a junction where M segments come together, 3(M-
1) constraint equations need to be enforced. At a junction between two straight segments
having orientation angle 6, and 62 and components u,,wi and u2,w2 , the two
displacement compatibility conditions can be written in terms of the displacement
components in the Y and Z directions as follows:
u2 00502 +w2 sintS?2 —ul costSl1 -wl sin 6, =0 (3.12)
u2 sine2 —w2 c036?2 —ul sint9l +wl c036l =0 (3.13)
The rotation constraint requires that the angle between segments at a junction
remain unchanged. For the current example, if the rotations of the two members at the

junction are ,6, and ﬂ, , compatibility condition is

,6, — ,6, = O (3.14)

19

CHAPTER 4
CONSTITUTIVE EQUATION

The elastic-plastic model which has been implemented in the program has been
developed by Pourboghrat et a1 (2000). It is based on the isotropic hardening model. The

uniaxial stress-plastic strain curve of the material is given by

5 = K(E+Eo )” (4.1)
where a" is the effective stress and E is the plastic strain. Parameters K, N and 50 are
material constants. The elastic strain increment is related to the stress increment through
the equations of linear, isotropic elasticity with Young’s modulus E and Poisson’s ratio.

The yield function allows for anisotropic yielding of the material. The yield function is

given by

2+ 2+R — 2
= 0x 0': (0.x 0', ) — 0': = 0 (4'2)
1+ R

 

f

where 0', is the yield stress and R is an anisotropy parameter.

During loading, Hooke's law is used to calculate stress below the elastic limit; i.e.,
5" 5 0y , where 0'y is the initial yield stress of the sheet obtained from a uniaxial tensile
test. Beyond the elastic limit; i.e., E > 0,, the co—rotational time derivative of stress

(Jaumann stress rate) is calculated, for a given strain rate, from an elasto-plastic

constitutive equation (Becker, 1992):

 

a: L— .. : :D (4.3)

 

 

20

V
Here 0' and D(= D‘+D” ) are the Jaumann rate of stress and strain rate tensors,

respectively, 0' is the stress tensor, 5 is the material flow strength, h(= 83/35) is the

plastic hardening parameter, L is the fourth order elastic tensor and

an

 

Ii: (Mr/3 g/"dqb/a 9'"), where “am 9'“ = Jim/do : d¢/Bg , is the second order tensor

representing the unit normal to the ﬂow potential surface. The plastic strain rate,

associated with Eq.(4.3) , is calculated from the following expression:

V
. P:L: -P:0’
5:” ” ” ‘ " (4.4)

 

The fourth order elastic tensor L(= LU“) used in this work is the standard tensor for the

isotropic elasticity, which has only two independent components.

21

CHAPTER 5

CONTACT ALGORITHM

The simulation of hydroforming processes requires the ability to model the contact

phenomena. The analysis of contact behavior is complex because of the requirement to

accurately track the motion of multiple geometric bodies, and the motion due to the

interaction of these bodies after the contact. The numerical objective is to detect the

motion of the bodies, apply a constraint to avoid penetration, and apply appropriate

boundary condition to simulate the frictional behavior (Marc’s Theory and user

information).

5.1 Numerical Procedures for Contact

Two widely used procedures are available, namely the Lagrange multiplier

method and penalty method, which will be introduced here brieﬂy. Both, the Lagrange

multiplier method and penalty method operate on variational formulation of problem

under consideration (Bathe, 1982).

For purpose of demonstration, consider the total potential energy l'lp of a discrete

system for steady-state analysis,

11, =%UTKU—UTR

The constraint equations can be expressed as
CU = O
In the Lagrange multiplier method we amend the r.h.s. of (5.1) to obtain

11; =—;-UTKU-UTR +1 TCU

22

(5.1)

(5.2)

(5.3)

where ,1 contains as many Lagrange multipliers as there are constraint equations, and

invoke that 5 fl; = 0 , which gives

T U R
K C = (5.4)
C 0 1. 0
The second of Eq.(5.4) is Eq.(5.2), the equation of constraint. Eq.(5.4) are solved for both

U and A. The l 1 may be interpreted as force of constraint.

In the penalty method we also amend the r.h.s. of (5.1) but without introducing an
additional variable. Now we use

n; = éUTKU-UTR+%tTat (5.5)

in which t = CU , and t = 0 implies satisfaction of the constraint. The usual potential IIp
of the system can be augmented by a penalty function ét'rat , where a is a diagonal

matrix of penalty numbers a, . The condition 5 I]; = 0 now yields

[K + cTac J{U} = {R} (5.6)
If 0. =0, the constraints are ignored. As it grows, U changes in such a way that the
constraint equations are more nearly satisﬁed.

It is difﬁcult to assess which method, the Lagrange multiplier method and penalty
method, is better in the ﬁnite element implementation. However, it can be mentioned that
the Lagrange multiplier method introduces a large number of independent variables.
Also, the diagonal terms of the stiffness matrix corresponding to the 2. always become
zero and special attention is required during the assembly of the stiffness matrix so that

an equation corresponding to any A does not become the ﬁrst one in the stiffness

23

equations after applying the boundary conditions (Kobayashi et al, 1989). In comparison
with Lagrange multiplier, penalty functions have the advantages of introducing no new
variables and requiring no special care for the order of unknowns for successful equation
solving. Their implementation can be as easy as assigning a high modulus for an element
already in the program. Thus, we use penalty functions to implement the constraints in
contact problem.

Effectively, the penalty procedure constrains the motion by applying a penalty to
the amount of penetration that occurs in the hydroforming processes. The penalty
approach can be considered as analogous to a nonlinear spring between the two bodies.
Using the penalty approach, some penetration occurs with the amount being determined
by the penalty constant or function. The choice of the penalty value can also have a
detrimental effect on the numerical stability of the global solution procedure.

5.2 Detection of Contact

During the contact process, each potential contact node is ﬁrst checked to see
whether it is near a contact segment. The contact segments are either edges of the other
2D deformable bodies, faces of 3D deformable bodies, or segments from rigid bodies
that stand for the dies in tube hydroforming process. The motion of the nodes is checked
to see whether it has penetrated a surface by determining whether it has crossed a
segment. The determination of when contact occurs and the calculation of the normal
vector are critical to the numerical simulation. In this part, the normal vector of potential
contact node has been pre-deﬁned and a so-called cross-product algorithm is used to

determine which segment is associated with the potential contact node.

24

5.2.1 Contact Tolerance

During the contact process, it is unlikely that a node exactly contact the surface.
For this reason, a contact tolerance is associated with each surface (Figure 5.1). If a node
is within the contact tolerance, it is considered to be in contact with the segment. The

contact tolerance is calculated by the program as 5% of the thickness of segment.

 

2 x Tolerance

 

///////f/////////}

Figure 5.1 Contact tolerance
5.2.2 Normal Vector and Tangential Vector
In order to track the motion of each potential contact node, it is important to
deﬁne the normal vector and tangential vector of the node (Sheet-S).
The normal and tangential vectors of the nodes of the mesh (Figure 5.2(a)) can be

calculated as follows: The element i has nodes i and k and length Li , The element j has
nodes k and j and length L j . The unit normal vector, 11 , at the contact node k, as shown
in Figure 5.2(b), is evaluated by averaging the normals of elements i and j, Ni and N j ,

as follows:

Ni +Nj
n:| (5.7)

|N1+lel

By rewriting the normal vector, n, in the following component form:

2
"bl—hi ’1}; ZN=,/z§+1 (5.8)

25

 

 

 

(b)

Figure 5.2 (a) The computational mesh and

(b)The normal vector and tangent vector

26

where Zx which is interpreted as the slope at node k, can be obtained from Eq.(5.7), in

terms of the coordinates, (x,z), of the nodes i, k and j, as follows:

_ Lj(zk —zi)+Li(zj —zk)
" Lj(xk —x,)+Li(xj -xk)

 

(5.9)

where Li and L J. are respectively lengths of elements i and j. As seen in Eqs(5.8) and

(5.9), the normal vector, n, at the contact node k, can be expressed in terms of the
neighboring nodal coordinates.

The tangent unit vector, t, at the node k, whose direction is normal to the normal
vector, 11, and which is oriented opposite to the sliding direction, is expressed using the

orthogonality as follows:

r—i 1 (510)
—ZN zx ’

where ZN and Z, are deﬁned in Eq.(5.8) and Eq.(5.9), respectively.

Actually, as we can see in Appendix B, when the Fourier series expansion is used
to approximate the nodal coordinates, it is easy to exactly deﬁne the normal of a node by
differentiating the Fourier series without any numerical difﬁculty.

5.2.3 Condition for Contact Violation

As shown in Figure 5.3, when a node is within the contact tolerance, it is
considered to be in contact with the segment. The location of the closest point on the die
corresponding to the contact node needs to be determined by using the cross-product
algorithm. B“, B, and BM are die nodes, A), is one of the contact nodes and n is its

normal vector.

27

If (A,,B,_l xn)-(A,B, xn)>0 (5.11a)

 

and (A131 xn)o(AkB xn) < 0 (5.11b)

1+1

Then, we deﬁne segment Bi Bmto be the segment associated with contact node Ak.

 
  

 

Sheet Segment

/

/Die Segment

 

Bi+l

Figure 5.3 The associated die segment with contact node

5.3 The Projection Algorithm
A nodal position produced by the trial solution may penetrate the die. By using

the cross-product algorithm, the closest point on the die corresponding to the node can be

28

found. The nodal coordinates are then modiﬁed by a projection scheme such that the
node just touches the die surface. There are two ways to bring the penetrated node back
to the die surface.
5.3.1 The Iterative Method

The position of a penetrated node is changed to locate it at the tool (die) surface
where the normal vector for that node intersects it. After calculating the normal vector,
u = (nvnz) , using Eq.(5.8), the tool surface point to be projected (t-nx,t-nz) can be
iteratively found by solving for the length parameter, t, from the nonlinear equation:

S(t-n,‘)-—t-nz =0 (5.12)

where S(x) is the tool surface equation.
5.3.2 The Directive Method
As shown in Figure 5.4, PQ is assumed to be the tool segment associated with the
penetrated node A, point B is the intersection point between normal vector and PQ , and

O is the original point. Based on the following vector equation, the coordinate of point B

could be calculated:
5P+P§=DX+TB (5.13a)
ﬁrm—£2, , 7A3=—n,n (5.13b)
IIPQII

where 771,172 are scalar. Once I], , 772 are solved from Eq.(5.13), the coordinate of B could

be determined.

29

  
   

 

 

2 / Sheet Segment
Q
‘/ Die Segment

P 711
0%

Figure 5.4 The projection Algorithm
5.4 Implementation of Constraint

A node located on the sheet at time t +At , is constrained to move in the tangent direction
deﬁned by the trial solution, Au'. The tangent direction incorporates the Z, in Eq.(5.9).
The constraint on 6 u = (8 u, ,5 u) , for contacting nodes is then:

5 u -n = 0 (5.14)
which by substituting in Eq.(5.8) we will have

Su Z = 2x ~6u x (5.15)
Actually, Eq.(5.15) is the alternative constraint equation, Eq.(5.2), in the contact

problem.

30

5.5 Friction Model
Friction is a complex physical phenomena that involves the characteristics of the surface
such as surface roughness, temperature, normal stress, and relative velocity. The
numerical modeling of the friction has been simpliﬁed to two idealistic models.
The most popular friction model is the Coulomb Friction model. This model is used for
most applications. The Coulomb model (Marc’s Theory and user information) is:
oﬁS—uon-t (5-16)

where o n is the normal stress

a f, is the tangential(friction) stress

11 is the friction coefﬁcient

t is the tangential vector in the direction of the relative velocity

V

f

t = l— , vr is the relative sliding velocity.
vl’

 

The Coulomb model is also often written with respect to forces
f,S-1.rfn-t (5.17)
where f , is the tangential force, fn is the normal reaction.

For a given normal stress, the friction stress has a step function behavior based upon the

value of vr (see Figure 5.5).

31

‘ ftOI'O'ﬁ

Slip

+— Stick

 

 

 

T Slip

Figure 5.5 Coulomb Friction Model

This discontinuity in the value of of, can result in numerical difﬁculties so a
modiﬁed Coulomb friction model (Oh, 1982) is implemented:

f, =—mk£— =-3mktan"[Av’] (5.18)
lAv: 7t u0

 

 

where Av, is slipping velocity, m is friction factor, k local ﬂow stress in shear and uo is

very small positive number compared to Av, .

When the Coulomb model is used with the stress based model, the integration point
stresses are used to calculate the normal stress component of the contact node. The

tangential stress is then evaluated and a consistent nodal force is calculated.

32

5.6 Separation

After a node comes into contact with a surface, it is possible for it to separate in a
subsequent iteration or increment. Mathematically, a node should separate when the

reaction force between the node and surface becomes tensile or positive.

When contact occurs, a reaction force associated with the node in contact balance
the internal stress of the element adjacent to this node. When separation occurs, this
reaction force behaves as a residual force (as the force on a ﬁee node should be zero).

This requires that the internal stresses in the deformable body be redistributed.

33

CHAPTER 6
NUMERICAL SIMULATION

In this chapter, based on the virtual work principle, the equilibrium equation can be
derived. The numerical solution obtained from a Total Lagrangian formulation of a
general thin shell theory is employed. In contrast to traditional ﬁnite element method, the
nodal displacement approximate function based on Fourier series is used. The boundary
friction condition is introduced into the formulation in the form of penalty function,
which imposes the constraints directly into the tangent stiffness matrix. The Newton-
Raphson algorithm is used to solve the nonlinear equilibrium equation.

6.1 Principle of Virtual Work

An equivalent approach to express the equilibrium of the body is to use the principle of
virtual work. This principle states that the equilibrium of the body requires that for any
compatible, small virtual displacements imposed onto the body, the total internal virtual
work is equal to total external virtual work. In the Lagrangian incremental analysis
approach we express the equilibrium of the body at time t+At using the principle of

virtual work as:

J‘HArTij‘sHA’eU r+ArdV _—_ J'HArfB‘i‘iHArdV + J'HArﬁSaliS r+ArdS (61)

I
r+Arv r+ArV ”Ar;

Where l.h.s is the total internal virtual work and r.h.s is total external virtual work. The

””1”. are the Cartesian components of the Cauchy stress tensor, the ”A, e”. are the

Cartesian components of an infinitesimal strain tensor, the ”‘37." and "mfs are the

components of the externally applied body and surface force vectors, respectively, and

ﬁn, is the ith component of the virtual displacement vector.

34

6.2 Newton-Raphson Algorithm
Substituting the element coordinate and displacement interpolation into the principle of

virtual work expression Eq.(6.1), the equilibrium equation can be written as
HA! R - HAIF = O (6.2)
where the vector ”‘5’ R lists the externally applied nodal point forces in the conﬁguration

at time t+Ar , and the vector ”A’F lists the nodal point forces that correspond to the
element stresses in this conﬁguration. For the nonlinear problems, the basic approach for
an incremental step-by-step solution is to assume that the solution for the discrete time t
is known, and that the solution for the discrete time t+At is required, where At is a

suitably chosen time increment. Hence, considering (6.2) at timer + At we have

f(U') = 0 (6.3)
where f(U') = ”A' R(U') - ”A’ F(U’) (6.4)
Assume that in the iterative solution we have evaluated ”A“ U('"’ ; then a Taylor series

expansion gives

f(U’) =r('*~ U““’) + [35-] (U‘-'*~ U“’”) (6.5)

3U

 

MN UU-l)
where higher-order terms are neglected. Substituting from (6.4) into (6.5) and using (6.3)
we obtain

a_R_aF

O: I+AIR_ HAIFU-l) + __
[6U 3U

] (u‘ _I+Al Uu-u) (6.6)
0+6! UU-l)
when we assume that the externally applied loads are independent of the displacements,

which means 113‘ = O , (6.6) becomes
MAI ”(I-I)

35

3F . . ._
__ U _!+Atu(r-l) : t+AtR_ HAtFt l 6.7
[an]...u..-.,( ) ( >

We now define

AU“) 2U. _ t+Atu(i-l) (6.8)
and 3: ="A’ Km” (6.9)
MN UU-I)

where "“'K“"’ is the tangent stiffness matrix in iteration H, which we assume to be
nonsingular. Thus (6.7) can be written as

HA: K(i-1)AU(i) = (+21! R _ r+At F(i-l) (6.10)
Since (6.5) represents only a Taylor series approximation, the displacement increment
correction is used to obtain the next displacement approximation

HA! u“) = HAtuﬁ-l) +AU“) (6.11)
The relation in (6.9) and (6.10) constitute the Newton-Raphson solution of (6.2).

6.3 Three Types of Modeling

Three different models, i.e. for bending, pressure loading and frictional contact,

are derived respectively.

6.3.1 Tube Bending Modeling

Based on the kinematics and constitutive equation discussed in Chapter 3 and Chapter 4,
the principle of virtual work is used to satisfy equilibrium and takes the following form:

i L, (0,165; +6; 55; )dA‘ +22ij =raef (6,12)
,.

i=l

36

where I is the number of segments in the cross-section, J is the number of constraint

equations, xi]. are Lagrange multipliers, C j are the constraint equations(from Eqs.(3.12-

3.14), T is the tension (compression) applied, dA‘ = dsidz, for the straight segment or
dAi = Ridoidz, for arcs.
The displacement component w, ,u, for each straight segment are approximated using the

following Fourier series expressions:

             

               

       

w = a0 + ”261' cos— (6.13a)
n=1 1‘0 It
. N‘
u, = y; + 2y; (6.13b)
n=l "=1
while for an arc, they are given by
w =a' + N26? cos— ‘+ 2 ,8; mm
0 Sin— (6.143)
"=14: M
u = yo +:y; cosl— IZH'Z‘S‘ sinnA— (6.14b)

=1
After substituting from Eqs.(6.133—b) and (6.14a-b) into the principle of virtual work
Eq.(6.12) , a nonlinear expression of the following form will result

r(c,:c,,xZ ,T) = o (6.15a)
where c={a3,a;,ﬂ,j,y3,y;,6j,ef,lj} (6.15b)
Eq.(6.15) should be solved for c for given values of dKY ,dKZ and dT.

Since Eq.(6.15) is highly nonlinear, it will be numerically solved by the Newton-Raphson

method. The Newton-Raphson iterative method used to solve c will look like this:

37

(6.16)

In!

11511491441: F. — F
where K is the stiffness matrix, (Se is the incremental c and R is the force residual.

F Ex, is externally applied nodal forces, F ,m is the nodal forces that correspond to the

element stress in this conﬁguration ( Internal forces).

6.3.2 Pressure Loading Modeling
Pressure loading is modeled as an external force to expand the tube. The external work

done by a pressure p applied to the inside of the tube is equal to

3(Au‘)

C

 

I . . . I .
Mil 48441044; =21 44'

6921.4; (6.17)

where 6(Au‘) is the virtual displacement increment having two components 6(A w‘)

and 6 (A v’) , 73' is the unit outward normal to the surface AL.

The principle of virtual work for bending and pressure loading modeling becomes

$191581”; 58 )dA'+Z/1 5C =1’T58 +2jpn 6(Au )dA‘ (6.18)

P
[=1 A;

The variation of the virtual external work due to internal pressure loading is

FE; = nil] ITEM (6.19)

Due to the follower forces effect (Hibbitt, 1979), the load stiffness matrix is

 

' 66 .,-82A‘ ,.
1‘31:le Barz"c‘)(u‘)+ (u)

 

i=1 A, dc n 829 p
I 301 W803”
6.20
=12}. a? a? dA; ( )

b

P

38

Eq. (6.19) and (6.20) will appear on the right-hand and left-hand side of the Newton-

Raphson expression (Eq. (6.16)) respectively, namely,

[15+ K3,][d€]=[13]= Ff,“ + F5, — Fm (6.21)

When we implemented above formulations into the code, we found that the load stiffness

matrix Kg, has nearly no effects on the solution which leaded us to simplify the

formulation. Actually, if we approximated ﬁ‘dmu‘) = 5(Aw‘) , the whole formulation

32 (Aw‘)

82

 

above can be simpliﬁed. The proof can be seen in Appendix B. Since = O , there

C

will be no stiffness matrix due to pressure (K g, ) going to the left side of the Newton-

Raphson expression (Eq.(6.21)). Both methods are implemented in the code, the results

are almost the same.

6.3.3 Frictional Contact Modeling

The most challenging task when developing a numerical code for metal forming
processes is to model frictional contact. To model the tooling-workpiece frictional
contact correctly, the following two conditions must be continually checked during each
equilibrium iteration:

1. Penetration of the workpiece nodes into the tooling.

2. Nodal contact forces becoming tensile at the contact boundary (separation).

Once the penetration of the workpiece nodes into the tooling has been detected, the

penetration nodes must be returned to the tooling surface and constrained to stay on the

39

tooling surface for the remainder of the equilibrium iterations. The nodes, which are
returned to the tooling surface, are constrained to move only tangent to the tooling
surface and only condition 2 stated above can cause a contacting node to be detached
from the tooling surface. Figure 6.1 shows an example of a typical contact check during

the Newton-Raphson equilibrium iteration.

Tooling Surface

Workpiece (

 

P1 ‘

 

 

 

'—._I—I-I-._._._.—.-.—I-C_.-O—I

.

:
.. I

(a) (b) (C) ((1)
Figure 6.1 A schematic of tool-workpiece contact check
(a) shows the tube and the tooling, (b) shows initial penetration of some of the nodes,
(c) shows how those nodes are returned to the tooling surface and ﬁnally (d) shows how

the equilibrium shape is obtained aﬁer several iterations.

The external work done by the frictional contact is added into the virtual work principle

equation (6.12) as following:

40

I J I l
2L5(‘7;65;+0;.651.}1Ai +2216c, = was +2 Pﬁ,5(Au')dA; +2 J‘TiMAui )dA; (6.22)
i=1

l=l [=1 I=l A;

where 1,. is the traction on the surface of the tube and 6 (Au‘ )is the virtual incremental

displacement of the contacting nodes.

In order to improve convergence, a special algorithm is introduced (see Figure 6.2). For
each trial set of contacting and non-contacting nodes, equilibrium iteration is performed.
After equilibrium is satisﬁed, the nodes are reexamined for non-penetration condition.
The contact set is then updated by releasing or projecting certain nodes and another
equilibrium iteration is initiated.

In each contact iteration, the trial displacements are ﬁrst updated according to the
Newton-Raphson procedure and the non-penetration contact condition is then applied to
these trial values by projecting the contact nodes to the tool surface along the normal
vector. The modiﬁed trial solutions are used for Newton-Raphson iteration. Within this
force equilibrium iteration, the internal force is calculated. The signs of the sheet normal
force at contact points are checked so that the nodes having non-compressive force are

released.

41

 

TIME Iteration

 

 

 

 

 

 

 

CONTACT Iteration

 

 

 

 

 

 

 

N-R Iteration

 

 

 

 

 

 

 

 

 

 

New Trial

 

Converge

 

 

 

 

 

N-R loop

 

 

 

 

 

penetration

 

 

CONTACT loop

 

 

 

 

 

 

 

.__,

 

 

adjust time

+

 

 

 

 

 

 

 

 

 

 

TIME loop end

 

 

 

Figure 6.2 Flow chart for contact iteration

42

CHAPTER 7

PROGRAM STRUCTURE
The program includes input ﬁle, main program and output ﬁles.
7.1 Input File
The input is divided into two parts. The ﬁrst part contains all the geometric, numerical
and material information for a run. It is written in an input ﬁle that is read by the program
prior to speciﬁcation of the loading. The second part contains the loading history. This
input is provided interactively by the user. All input is free-formatted.
7.1.1 Geometric, Numerical and Material Input

1. Title of this case
TITLE

 

 

 

 

TITLE: Array of 72 characters used to identify a case.

2. Number of nodes
NODE

 

 

 

 

NODES: Number of nodes in the cross-section. Each segment in the cross-section is
bounded by two nodes. Maximum: 20.

3. Node coordinates. Input this record N ODES times before proceeding to 4.

 

 

 

YN.ZN

 

YN: Y global coordinate of the location of this node.
ZN: Z global coordinate of the location of this node.

4. Number of segments.
NSEG. MU. DELTAV

 

 

 

 

43

NSEG: Total number of segments in the cross-section. The segment may be straight lines
or circular arcs. Maximum: 40.

MU: The frictional factor.

DELTAV: The relative sliding velocity

5. Type of segment. Input NSEG item {5,6} pairs before proceeding to 7.

 

STYPE

 

 

 

STYPE: Enter str if the segment being described is straight or are if it is a circular are.
6. Segment parameters.

For straight segments input the following record:

 

 

 

INODE. JNODE. T. NCOEFF. IS. IZ

 

INODE: First node of the segment (3:0).

JNODE: Second node of the segment (szlo ).

T: Thickness of the segment (I).

NCOEFF: Number of terms in the Fourier series expansions for the displacement
components(N). Maximum: 10. Recommended: 4 to 6.

IS: Number of Gauss integration points along the segment. Maximum: 48.
Recommended: 8 to 12.

ll: Number of Gauss integration points along thickness of the segment. Maximum:
48. Recommended: 3 to 7.

For arc segments input the following record:

 

 

 

INODE. JNODE. T. NCOEFF. IS. IZ. YC. ZC. IAO

 

INODE: First node of the segment (¢=0).

JNODE: Second node of the segment (¢=A¢).INODE and JNODE must be chosen
so that the coordinate ¢ has a counterclockwise direction.

T: Thickness of the segment (t).

NCOEFF: Number of terms in the Fourier series expansions for the displacement
components(N). Maximum: 10. Recommended: 4 to 6.

IS: Number of Gauss integration points along the segment. Maximum: 48.
Recommended: 8 to 12.

ll: Number of Gauss integration points along thickness of the segment. Maximum:
48. Recommended: 3 to 7.

YC: Y-coordinate of the center of the arc(Yc).

ZC: Z-coordinate of the center of the arc(Zc).

IAO: Input a if A¢S 7trad., input 0 if A¢>7£rad.

7. Number of point constraints

 

NPC

 

 

 

NPC: Number of point constraints. At least three constraints are needed to prevent rigid
body motion. Also input the constraints needed to maintain contact with the mandrel or
symmetry conditions. All constraints are enforced at nodes.

8. Point constraint information. Input this record NPC times before proceeding to 9.

 

NSPC. IVARPC

 

 

 

NSPC: Segment number for which the point constraint will be enforced.
o If NSPC is positive, the constraint is enforced at the node with s = O or ¢ = O.

o If NSPC is negative, the constraint is enforced at the node with s = lo or (D =A¢.

IVARPC: Variable which will be constrained.

45

0 To constrain the rotation 13, input 1.
0 To constrain the displacement component u, input 2.
0 To constrain the displacement component w, input 3.

9. Iteration parameters:

 

 

 

TOL. ITEMAX. JEVAL

 

TOL: Convergence tolerance of Newton-Raphson iteration. Recommended:0.0001.
ITEMAX: Maximum allowable number of iterations.

JEVAL: Multiple of the iteration count at which the Jacobian is evaluated. If JEVAL =
1 the Jacobian is evaluated every iteration, If JEVAL = 2 it is evaluated every other
iteration, etc.

10. Material properties.
E. V. SIGMY. R, XK. XN. EPSO

 

 

 

 

E: Young’s modulus.

V: Poisson’s ratio.

SIGMY: Initial yield stress.
R: Anisotropy parameter.
XK: Material constant.

XN: Hardening exponent.
ESPO: Initial plastic strain.

1 1. Input I Output options.
IWR. IZSYMM

 

 

 

 

IWR:

- IWR = 1: print coefﬁcients and Lagrange multipliers.

46

0 IWR = 0: do not print coefﬁcients and Lagrange multipliers.
IZSYMM: This variable indicates whether the cross-section deﬁned by the input above
is a whole section or on-half of a symmetric section. The value of IZSYMM affects input
/ output operations only.

0 IZSYMM = 1: the cross-section is symmetric about the Z axis.

0 IZSYMM = 0: no symmetry about the Z axis exists.
7.1.2 Loading Input
The input of the loading parameters is interactive. Loading is prescribed in an
incremental manner. The loading record below may be repeated as many times as
necessary.

Loading parameters

NINC. CZINC. CYINC. TINC. PINC

 

 

 

 

NINC: Number of increments. Input 0 to stop execution.
CZINC: Increment size of the variable [(2. Note that the sign of CZINC must be

decided by the location of the mandrel. For example, CZINC must be negative for
the case shown in Fi g.4.1.

CYINC: Increment size of the variable X, . This will usually be zero.

TINC: Tension increment.
PIN C: Internal pressure increment.

7.2 Program Flow Chart
The program includes a main program and 24 subroutines.

Subroutine input: input geometric and material parameters.
Subroutine store: store integration, trigonometric and second derivative information.

Subroutine newrap: solve system of nonlinear equations using Newton-Raphson method.

47

Subroutine build: evaluate the system of equation and the jacobian.

Subroutine constr: add constraint contribution to system of equation and jacobian matrix.
Subroutine decomp: ﬁnd the inverse matrix using Gauss-Jordan method.

Subroutine disp: evaluate displacements and their derivatives with respect to s at
integration points along a segment.

Subroutine dstrn and dstrnC: derivatives of displacement components and their
derivatives with respect to the coefﬁcients of the series expansions at one point along s
for straight and are segment respectively.

Subroutine d2curvature: evaluate the second derivatives of curvature.

Subroutine incm: increment of all variables.

Subroutine solve: solve the nonlinear equations.

Subroutine stndsp: ﬁnd strains at the integration points along the local coordinate.
Subroutine stsstn: ﬁnd stress increments and force residual at integration points.
Subroutine cnstutv: constitutive equation.

Subroutine contact: deal with contact problem.

Subroutine coloumb: use Coloumb friction model.

Subroutine dieinfo: calculate die geometric information.

Subroutine forcecheck: check nodal force at each contact point.

Subroutine getcoord: calculate the coordinates of each integration point.

Subroutine nodalf: calculate nodal force at each contact point.

Subroutine setpenalty: set penalty function.

Subroutine violation: check each integration point to see if it violate contact.

Subroutine bringback: bring each contact point to die surface.

48

Following is the flow chart of the program:

 

 

SIOI‘C

 

 

input

 

 

 

 

 

 

 

 

 

dieinfo

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

build

 

 

 

 

 

 

 

 

COl’lStl'

 

 

 

 

 

 

 

 

 

 

 

newrap

 

 

 

 

decomp

 

 

 

 

 

solve

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

contact

 

 

 

 

 

 

 

 

 

incm

 

 

 

 

 

 

 

 

 

 

 

 

 

stndsp disp
ern or dstrnC

stsstn __J cnstutv
violation bringback
setpenalty
.oloumb
forcecheck nodalf

 

 

 

 

 

 

Figure 7.1 Flow chart of the program

 

7.3 Output File
The program has three output ﬁles: betges.out, betgesmk, betges.coeff. The contents of
each ﬁle are described below:
betges.out: This ﬁle is formatted and contains the output of the program in the form of
tables. The ﬁrst part of this ﬁle is written prior to the speciﬁcation of the loading and
contains an augmented echo print of the input. It is useful to look at this part of the ﬁle to
ensure that the input ﬁle was read correctly. The second part includes a list of the
curvature, axial strain, bending moments, tension and internal pressure at each loading
step. If the variable IWR is set to 1 then the ﬁle will also contain the coefﬁcients of the
trigonometric series expansions for each segment. If the variable IZSYMM is set to 1
then the value of T and M y are reported for the full section, while M2 is set to zero.
betges.mk: This ﬁle contains one line per loading step containing the value of the
following parameters:
KZ,K,,£f,M,,MZ,T,P
betges.coeff: This ﬁle contains information necessary to make plot of the deformed
cross-section of the member at any loading step. It also contains information to plot the
shape of the cross-section prior to loading. The ﬁle is organized as follows:
Line 1: (Format: i2)

Number of segments.

For each segment, the ﬁle contains the next set of three records:

50

 

Line 2.1: (i3, e135, a4, 3el3.5)
Number of terms in this segment, length of the segment, type of segment, location
of the segment in Y, location of the segment in Z, radius of the segment. For straight

segments, the location of the segment is partially given by the coordinates Y0 and Z0 as

shown in Figure 3.2. The radius is assigned a value of one. For arc segments, the length
of the segment is in angular measure, the location of the segment is partially given by the

coordinates of the center Yc and Zr as shown in Figure 3.3.

Line 2.2: (4e13.5)

Coordinates of the nodes bounding this segment in the following order: Y-
coordinate of ﬁrst node, Z—coordinate of ﬁrst node, Y-coordinate of second node, Z-
coordinate of second node.

Line 2.3: (2el3.5)
Sine and cosine of the angleB deﬁned in Figure 3.2 for straight segments. For

arcs,sin 19 = cos ¢0 and cost? = —-sin ¢o , where ¢ois deﬁned in Figure 3.3.

This concludes the part of this ﬁle which contains the original, undeformed description of
the cross-section. The second part of the ﬁle contains, for each loading step, additional
information necessary to describe the deformed shape of the cross-section. For each
loading step, the ﬁle contains the following:
Line 3.1: (i2)

Step number.
Line 3.2: (i2)

Segment number.

51

 

For each segment in the cross-section the output ﬁle contains the following set of four
records:
Line 3.3.1: (6(e22.14))
a0 ,( an , n = 1, number of terms).
Line 3.3.2: (6(e22.14))
(ﬂu , n = 1, number of terms).
Line 3.3.3: (6(e22.14))

70 .( 7,. . n = 1, number of terms).

Line 3.3.4: (6(e22.14))

( 5n , n = 1, number of terms).

52

CHAPTER 8
EXAMPLES
There are several examples shown in this chapter, including stretch-bending, expansion
under internal pressure and contact example.
8.1 Bending Examples
Example 1 : Stretch-bending a tube with rectangular cross-section
This example concerns a tube with rectangular cross-section. The section is 1.97 in. wide

and 1.18 in. high (outside dimensions). The wall thickness is 0.07 in. See Figure 8.1. The
material properties in plasticity implementation 0' = K (E + 8,, )N are: E = 10 x 106 psi, v:
0.33, a, = 5,500 psi, R=0.7, K = 64 x 103psi, so = 0.0 and N = 0.25 .

The input ﬁle used to run this case is shown below:

’ Rectangular section ’
4

0.950,0.555
0.950,-O.555
-0.950,-O.555

-0.950,0.555
4,0,0

str
1,2,0.070,4,8,3
str
2,3,0.070,4,8,3
str
3,4,0.070,4,8,3
str
4,1,0.070,4,8,33
1,2

1,3

—1,3
0.0001,20,1
l.e7,0.33,5500.,O.7,64000,.25,0.0
0,0

53

 

The loading history consisted of the following
1. Apply tension to 4200 lbs.

2. Bend to a radius of 20 in.

The interactive input records were as follows:
21, 0., 0., 200., 0.

100, -0.0005, 0., 0.

Example 2: Stretching I bending a tube with circular cross-section
This example concerns a tube with circular cross-section. The diameter is 0.6092 in. The

wall thickness is 0.0356 in. See Figure 8.2.
The material properties(aluminum-alloy 0' = K (6 + £0)” ) are E = 10.295 x 106 psi, v =
0.3, 0', =19,618psi, R = 1.0, K: 83,634 psi and N = 0.3593, 80 = 0.01658.

The ﬁnal state is x = 0.064 in“ and T = 1000 lb.

The input ﬁle used to run this case is shown below:

’Circular case’
2 .

0.0001,20,1
10.295e6,0.3,l96l8,l.0,83634,0.3593,0.01658
0,1

54

Example 3: Bending a tube with circular cross-section to collapse

Figure 8.3 shows example of bending a tube to collapse. The material and dimension are
similar to the one in example 2.

8.2 Free Expansion Case

Example 4: Formulated, implemented and tested the internal pressure loading and
expansion of the tube.

Figures 8.4 shows examples of expanding symmetric tube with internal pressure. The
material and dimension are similar to the one in example 2. Final internal pressure is
lMMpﬁ.

Same example was run in ABAQUS by using 3D shell element. Maximum hoop strain
are 16% and 18%, axial strain are —8% and —9% and the change of diameter are 0.183in
and 0.2m, in this 2D-Hydroforming code and ABAQUS respectively. It shows that
results are in good agreement.

8.3 Contact Examples

Example 5: A round aluminum tube expanded into a square die with sharp-corners.
Final internal pressure is 7,000 psi, material is similar to aluminum-alloy in example 2.
Final shape is seen in Figure 8.5a. Compared with the results of ABAQUS(see Figure
8.5b). The ﬁnal shape is in good agreement. Twenty-six S4 shell elements were used in

ABAQUS.

55

Example 6: Hoop strain distribution and intermediate shapes in a specialized

contact example (Compared with GM experiment data)
The material properties (aluminum—alloy 0' = K (8 + £0)” ) are E = 10.3 x 106 psi, v: 0.3,
0’, = 20,000 psi, R = 1.0, K = 63657 psi and N = 0.25, 80 = 0. 0.

Figure 8.6a shows hoop strain distribution and Figure 8.6b shows the intermediate

shapes. Figure 8.7a-b are the experiment data of strain distribution.

56

Curvature = 0.05 I in, Tension = 4200 lb

 

I Y T Y I I Y Y I

’ Original shape

 

 

Final Shape

0.2 .

-0.4 »

 

 

 

0.5 Q .

 

I J 1 1 l l 1

-0.8 —0.5 -0.4 -0.2 0 0.2 0.4 0.5 0.8

Figure 8.1 Stretching / bending of a rectangular section

Curvature = 0.02lln, Tension = 1cm":

 

 

 

 

 

 

06 V V ’1th V WK; V Y . .
’ Orlgrnal shape
.57
0.4 '- I/
y
0.2 - .
l
f
0 h k
\r
0.2 '- I:\\. l! \
\~ 7% Final Shape
‘ . /
\. .2
0.6 1 a \M— r J/ r r
-0.6 -04 -0 2 0 0.2 0.4 0.5

Figure 8.2 Stretching / bending of a circular section

57

 

 

 

 

 

 

 

Figure 8.3 Bending a tube to collapse

58

mDO<m< £3, 2:388 :32 EoEooﬁamE use cﬁbm Wm Bswi

:3 :83 w

as- em- “68: seem 3.2

as: as “seem 86: .52
as $2

 

 

8:05 0956 u 305.02% .855 good u 36am 3
in coca 0 8.58.5 a

 

wmmzod new .mmmmd M Z .md u > in 36mm u v— ia @202 n m> in camomdd n m “we? 05 3:232

md

A

 

1 mo-

 

 

 

 

an 08.m n 2333 3525

59

lntemal pressure = HID psi

 

 

W‘
0.. 3"},
./ ,/"
f. /’
: I/
o
1 - . -
.//
//
I/JI
/
,/
l/
0.5 " / "‘
//

 

 

 

-1.6 -1.4 -1.2 -1 08 05 -04 -02 0

Figure 8.5a A round aluminum tube was expanded into a square die with sharp-comers

 

 

 

 

Figure 8.5b Final shape from ABAQUS

60

E05596 :33 33:88 59536.6 58; noon 3.x aswi

sees a: _ .o u 365.25. seeps at: u 33. 65 co 8:3:
a: coed .l. 8:395
one .2212 .31, 2.88%" x .2 38?? .8 8:87:

 

 

 

Do
I

1l
L/
7

 

 

 

 

652w:

6&3

65%:

used—

$02

$5.3

semi

$93

on:

 

some
see:
sane
.32
.32
see:
.3:
s3:

Bu:

3
gala
”321%:

seems.

653

$0.3

$5.3

65?»:

655.2

$02

awe—.3

 

men:

236568 8.3::

653-0?

$3-0m

sew—-0:

osNTa

$2-0

§E-m

oxemTM

 

~585me

 

8836

61

 

8:2: :2 .o u 3222.: .855 E: u peas 65 .6 see: a
begin»:

a: Sod pee as 83 88% :2. 83“ u 2:28: welsﬁ
due .deuz

.89 .882? M .2: 38?? .2: 8:02": :2: 8: sees:

 

z?» + eka H b ”—2002 Etoﬁz

wEE:£o::>E misc 2.3. of Mo 3926 23:08:85 nod enmﬁ

 

 

 

 

 

 

 

62

ugens JOUHN

:oﬁwuﬁEoo

 

c222

wcﬁcon “can couanEmG Egm EaceﬂEzoto a ﬁw 2sz

cap

our

2;

oww

8P

cm

005% 60:30..

on ow ON

 

q

_ _ —

 

 

1111

 

 

 

 

 

 

 

 

 

3
352+ M
352:? H
- ow
- 8

 

 

 

% ‘UIWIS

63

1s .Iougw

UIBJ

388$ $053802 5 833.5% 58; Eﬁéesea £5 2&E

 

00.320500
com F

 

00.002

O3

3500 60:000..

 

 

 

 

 

 

 

 

Oop 8P ONF OOF OO 8 cc ON O
k ( K\

M

m 3:..0D 1

m .5500 Sun In]

$00.00: no. .

 

 

 

 

 

OF

mp

ON

mN

% ‘UIWIS

CHAPTER 9
CONCLUSIONS

Based on the formulations presented in the thesis, a cross section analysis ﬁnite element
code had been developed. This code is very efficient in calculating tube stretching and
bending followed by hydroforming and expansion to ﬁll the die in the hydroforming
process. The code has been tested for several cases. The results from this 2D code were in
good agreement with both commercial code ABAQUS and the limited experiment data
available to us.

In the future, the failure model and instability (bifurcation, post-buckling etc) analysis
will be implemented into the code. It is also planned to extend this code for the 3D

analysis tube hydroforming process.

65

 

APPENDIX A
According to the total Lagrangian formulation discussed in Chapter 2 (Thin shell theory),
we derived the membrane strain, rotation of normal vector and the current principal
centerline curvature used in Chapter 3 (Kinematics) and compared with the formulation
developed by Corona and Pourboghrat (1998) and Brush and Almroth (1975).
1. Membrane strain

For straight segment, the curvature K1=O, then from Eq.(2.l4) ,
,1] =[(l+u.s)2 +w‘25]é =[l+u‘25 +2115 +w.25]% (1)
Let’s approximate the square root as:

(1+x)% = l+%x—-~ forx<<1, Eq.(l) becomes:

1 1 2
A] =l+u.s +5113; +§W§ (2)
Since, A, = ds/dS and the engineering deﬁnition e = ﬂ = ﬂ --1 = 2., —1
dS dS
We have the membrane strain:
1 l
e z u's +5185 +5“): (3)

which is used in Chapter 3 as Eq.(3.5) and Eq.(l.6) in Brush and Almroth (1975).
For circular arc segment, the undeformed shell is non-ﬂat and has an initial curvature K r,

then rearranging Eq.(2.l4) and approximating the square root , we obtain
__ 1 2 2 2 _ _ 2 2 2
A] ~l+—2-[u.s +2145 +K1w 2Klw 2K1wu5 +w,s +Klu +2K1uw‘s] ,

then membrane strain is

66

 

e zémi +w‘25)+%K12(u2 +W2)+u.s +K1(uw.s ’W“.s _W) (4)

Consider the circular arc segment shown in Figure 3.3 Since Rd¢ = dS , then

30 _ 1 30) 30)

BS R W) ‘35, Eq.(4) becomes:

e z $19204; + w;) +%K12 (u2 + w2 ) + K114, + Kl (Kluw, — K,wu., - w)
=-;-Kl2[(u.2, +wi,)+(u2 +w2)+ 2(uw, —wu.¢)]+ K1(u', —w)

: u,,—w +1 u+w.' 2+1 w—u, 2 (5)
R 2 R 2 R

As shown in Figure 3.3, if the positive w direction is the opposite direction of the normal

 

 

 

to the mid-surface vector, ﬁ , then Eq.(S) becomes:

2 2
e: u',+w +1 u—w, +1 u',+w (6)
R 2 R 2 R

which is used in Chapter 3 as Eq.(3.10) and Eq.(4.3) in Brush and Almroth (1975).

 

 

 

Based on the following assumptions, we could derive the normal vector rotation and
current principal centerline curvature for both straight segment and circular arc segment.
Assumption 1: The rotation angle B is very small,

Assumption 2: Bending deformation is dominant and stretching is negligible, namely,

u" =0, 14.55 ==Oor/1, fig-=1.

67

2. Rotation of normal vector

We assumed the angle between the current normal vector ii and the S (arc length) is on,

the angle between the current normal vectorn and normal vector 19 at the reference is B.

Then, we have B +32E =a , and Pr A = cosa = -sinB (7)

Substituting Eq.(2.16) into Eq.(7), we obtain:

 

 

 

 

 

w + K w + K
SinB = i = .S I“ = .5 I“ (8)
I11 11. Mm, -K,w)2 +(w,_. +K1u)2
For straight line segment (K 1:0), according to assumptionl, Eq.(8) becomes
. W s
srnB z = ‘ (9)

According to assumption 2 and wi <<1, Eq.(9) becomes: B = w‘s , for the coordinates -
shown in Figure 3.2, we have the normal vector rotation for straight segment

B = -W,s (10)
which is used in Chapter 3 as Eq.(3.1).

For circular arc segment ( K 1 ¢ 0) , according to assumption 1 and 2, Eq.(8) becomes:

B = K,u +w.s Since d5 = Rd¢ and K1 =% , we obtain

 

 

_ u + w", (11)
_ R
For the coordinates shown in Figure 3.3, we have
B-"-”‘ an
- R

which is used in Chapter 3 as Eq.(3.7) and Eq.(4.4) in Brush and Almroth (1975).

68

 

3. The current principal centerline curvature

For straight segment (K 1:0), the centerline curvature at the current conﬁguration, k, will

 

 

be obtained from Eq.(2.20):
, = Cd'sédc's : (1+u.5 )WZ: — w‘Su‘SS (13)
According to assumption 2(u_S = uss = O ), 2, = W , the Eq.(l3) becomes:
k. = $7- (14)
(1+ Wis )3

which is used in Chapter 3 as Eq.(3.6).

For circular arc segment, it seems not easy to get the simpliﬁed curvature from
Eq.(2.20).We begin from Eq.(12),which is based on the small rotation, however, for ﬁnite
rotations, we use

u—w

I
R (15)

 

sinB =
According to Flugge (1973, pp.362), we have the change of curvature:

K, =% (16)

By differentiating Eq.(lS) , we obtain:

_“¢

- Rcosﬁ

“W” (17)

where

 

 

cosﬂ=‘/l-sin2,l3=‘ll—(uuRw‘i)2 (18)

Substituting from Eqs.(l7),(l8) into Eq.(l6), we have

69

 

x = R (19)
_ u—w‘c, 2
RJ1 (———R )

which is Eq.(l l) in Corona and Pourboghrat (1998) for small rotation case.

 

 

70

 

APPENDIX B

Here we try to prove the equivalence of two formulations in the pressure modeling, ie:

F;,=p ZIMd/l ~—p2[ ig‘g’ja;

(1)

From the Chapter 3 kinematics, for the point on mid—surface (circular segment), we have

4"" =26 +(R+w‘)sin$+v‘cos¢7

77‘ = Yc +(R+w‘)cos$—visin¢7

where (3 = ¢0 + at . Replacing the w, u with w‘ , v’ respectively from Eqs.(3.8a-b).

Deﬁne ii" is the unit outward normal to the surface A; , then we have

 

 

(_a__§' _3___77'
.,- 8¢ 3¢
(_ _8_§' 877_'_
E)¢ 8(1)
where
3¢_ 29¢ srn¢ (R+w)cos¢ B¢COS¢+V srn¢
L3).=_é;c osJ—(R+w‘)sin¢7-%%lsin5-Vi005$
-9:er _ aw' av :-
"( a(1),a¢)"_{(a—¢—--v‘2) +[(R+w‘)+a¢]2“} R+w
Wealsohave
amu"): 6(An") am?)
a? 39 ’ a?
where

71

(2.a)

(2.b)

(3)

(4.a)

(4.b)

(4.c)

(5)

 

A? = Aw‘ sin e) + Av' 003;? (6.a)
An‘ = Aw‘ cosE—Av‘ sin 5 (6.b)

Combining (3), (4.a-c),(5) and (6.a-b), we have

 

 

 

 

 

 

 

ﬁ,a(Au)=_ l . (R+w,)3(Aw)+v,-8(Av)+av 3(Aw)_8w 8(Av)
a c R + w‘ 8 c a c a,» a c 59¢ a c
3(Aw‘) ‘
,_ _ 7 E
ac ( ) 7
For the straight segment, same expression can be obtained. Thus, we have proven the (1). E
. 32 (Aw‘) . . . . .
Since 32 c = 0 , there Will be no stiffness matrix due to pressure gorng to left srde

of Newton-Raphson expression.

 

72

REFERENCE

Ahmetoglu M. and Altan T., 2000, Tube hydroforming: state-of-the-art and future
trends, J Mater Process Tech 98: (1) 25-33.

Anon, 1997, Hydroforming technology, Advanced Material Process, 151(5) ,50-53

Bathe, K-J, 1982, Finite Element Procedures in Engineering Analysis. Prentice-Hall,
New Jersey.

Becker, RC, 1992, An efﬁcient semi-implicit integration algorithm for anisotropic ﬂow
potentials

Brush, DD. and Almroth, 1975, Buckling of Bars, Plates, and shells, McGraw-Hill.

Choudhry, S. and Lee, J.K., 1994, Dynamic plane-strain ﬁnite element simulation of
industrial sheet metal forming processes, Int. J. Mech. Sci. Vol.36, No. 3, pp.189-207.

Corona, E. and Pourboghrat, F., 1998, Bending-Stretch Forming of Extrusions- BET GES
Computer Program Manual, Alcoa internal report, March.

Dohmann, F. and Hartl, Ch., 1997, Tube hydroforming: research and practical
application, J. Mater Process Tech 71 174-186.

Dohmann, F. and Hartl, Ch., 1996, Hydroforming: a method to manufacture light-weight
parts, J Mater Process Tech 60 669-676.

Flugge, W., 1973. Stresses in Shells, 2nd Edition. Springer-Verlag, New York.
Germain,Y.,Chung, K. and Wagoner, R.H.,l989,Int. J. Mech. Sci. Vol. 31, N01, pp. 1-24.

Hibbitt, H. D.,l979, "Some follower forces and load stiffness", Int. J. num. Mech. Engng,
14,937-941.

Hsu, Tze-Chi and Chu, Chan-Hung, 1995, A ﬁnite element analysis of sheet metal
forming process, J Mater Process Tech 54 , 70-75.

Hu, P.,Liu, Y.Q., Li, Y.X. and Lian, J ., 1997, Rigid viscoplastic ﬁnite element of the gas-
pressure constrained bulging of superlastic circular sheets into cone disk shape dies, Int.
J. Mech. Sci. Vol. 39, No.4, pp. 487-496.

Kobayashi, S., Oh, 8.1 and Altan, T., 1989, Metal Forming and the Finite Element
Method, Oxford University Press, New York.

Longhouse, B.,l998, Applied pressure sequence hydroforming, 3rd Annual automotive
tube fabricating.

73

 

Marc’s Theory and user information
Morphy, G., 1997,Hydroforming high strength steel tube for automotive structural
applications using expansion, SAE International Congress & Exposition.

Nagtegaal, Joop C. and Rebelo, N.,l988, On the development of a general purpose ﬁnite
element program for analysis of forming processed, Int. J .Num. Methods Eng.,25, 113-
131.

Ni, C.M., 1994, Stamping and hydroforming process simulation with a 3D ﬁnite element
code, SAE Paper #940753, SAE SP-1021

Noh, TS. and Yang, D.Y., 1998, A general formulation for hydroforming of arbitrarily
shaped boxes and its application to hydroforming of an elliptic-circular box, Journal of
Manufacturing Science and Engineering, August, Vol. 120 481-488.

Oden, J .T. and Pires, EB, 1983, Nonlocal and nonlinear friction laws and variational
principles for contact problems in elasticity, ASME J. Appl. Mechs., Vol. 50,67-76.

Oh, S.I., 1982, Finite element analysis of metal forming processes with arbitrarily shaped
dies, Int. J. Mech. Sci. Vol.24, No.8, pp. 479-493.

Pourboghrat, F., Karabin, M.E., Becker, R.C.and Chung, Kwansoo, 2000, "A hybrid
membrane/shell method for calculating springback of anisotropic sheet metals
undergoing axisymmetric loading", International Journal of Plasticity 16 677-700.
Sheet-S, Ohio State University

Siegert K., Haussermann M. and Losch B., et al, 2000, Recent developments in
hydroforming technology, J Mater Process Tech 98: (2) 251-258.

The ABCs of Hydroforming, Vari-Form Inc.
Wang, NM. and Tang, SC, 1988, Analysis of bending effects in sheet forming
operations, Int. J. Num. Methods Eng., 25 253-267.

Wang, NM. and Budiansky, B., 1978, Analysis of sheet metal stamping by a ﬁnite
element method, ASME J. Appl. Mechs.,Vol. 45 73-82.

Zhang, SH, 1999,Developments in hydroforming, J Mater Process Tech 91:(1-3) 236-
244.

74

IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

ylllluallwljﬂlumlﬂw(will