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

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 chIRC/DateDuost-p. 15

FOURIER SERIES BASED IMPLICIT FEM PROGRAM FOR AXISYMMETRIC
ANALYSIS OF TUBE HYDROFORMING
By

Manish Sharma

A THESIS

Submitted to
Michigan Sate University
in partial fiilfillment of the requirements
for the degree of

MASTER OF SCIENCE
Department of Mechanical Engineering

2003

ABSTRACT

F OURIER SERIES BASED IMPLICIT FEM PROGRAM FOR AXISYMMETRIC
ANALYSIS OF TUBE HYDROFORMING

By
Manish Sharma
Tube Hydroforming has become a preferred alternative to stamping and
joining processes in industry for manufacturing closed cross section structural
components. This process requires proper combination of part design, material selection,
application of internal pressure and axial feeding. By using finite element analysis, these

process parameters can be detemlined efficiently in the very early stage of production.

The thesis presents a Fourier series based axisymmetric analysis of Tube
Hydroforming. The material is assumed to elastic-plastic and to satisfy the plasticity
model that takes into account rate independent work hardening and normal anisotropy.
Numerical solution obtained from Total-Lagrangian formulation of general shell theory is
employed. The axial stroke is incorporated using Lagrange multipliers. Contact violation
and boundary fn'ction condition are introduced into the formulation based on penalty
functions, which impose constraints directly into the tangent stiffness matrix. The Force
Limiting Curves based on shear instability and on experimental results are used as
fracture criteria. The Axisymmetn’c Hydroforming program developed is verified against
results from Abaqus and experimental measurements for isotropic and anisotropic

(transverse) materials.

ACKNOWLEDGMENTS

The author wishes to express his appreciation to the following people 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 Guan Yabo for his help in extending this code and also in preparing this manuscript.

iii

TABLE OF CONTENTS

LIST OF FIGURES ................................................................................. vi
INTRODUCTION .................................................................................. 1
CHAPTER 1
LITERATURE REVIEW ........................................................................... 4
CHAPTER 2
THIN SHELL THEORY ............................................................................ 9
Principal Stretch (Total Lagrangian) ....................................................... 9
Principal Curvature and Stretch (Updated Lagrangian) ................................ 12
CHAPTER 3
KINEMATICS ....................................................................................... 16
Kinematics Assumptions ................................................................... 16
Principal Strains ............................................................................. 16
Constraints .................................................................................. 18
CHAPTER 4
CONSTITUTIVE EQUATION ................................................................... 20
CHAPTER 5
CONTACT ALGORITHM ........................................................................ 22
Numerical Procedures for Contact ....................................................... 22
Detection of Contact ....................................................................... 24
Contact Tolerance ................................................................. 24
Normal Vector and Tangential Vector .......................................... 25
Condition for Contact Violation ................................................ 27
The Projection Algorithm ................................................................. 28
The Iterative Method .............................................................. 28
The Directive Method ............................................................ 29
Implementation of Constraint ............................................................ 3O
Friction Model .............................................................................. 3O
Separation ................................................................................... 32
CHAPTER 6
NUMERICAL SIMULATION .................................................................... 33
Principle of Virtual Work ................................................................. 33
Newton-Raphson Algorithm .............................................................. 34
Four Types of Modeling ................................................................... 35
Pressure Loading Modeling ...................................................... 35
Axial Loading Modeling ........................................................... 38

iv

Displacement Stroke Modeling .................................................. 38

Frictional Contact Modeling ..................................................... 39

Fracture Criteria ........................................................................... 43
CHAPTER 7

PROGRAM STRUCTURE ....................................................................... 45

Input File .................................................................................... 45

Geometric, Numerical and Material Input ..................................... 45

Loading Input ...................................................................... 49

Program Flow Chart ........................................................................ 51

Output File .................................................................................. 52
CHAPTER 8

EXAMPLES ......................................................................................... 54

Aluminium Tube Example ................................................................ 54

Steel Tube Example ....................................................................... 57
CHAPTER 9

CONCLUSIONS .................................................................................... 73

APPENDIX A ....................................................................................... 75

APPENDD( B ....................................................................................... 77

REFERENCES ...................................................................................... 8O

LIST OF FIGURES

Figure 0.1 Representative hydroforrned components and assemblies ......................... 2
Figure 0.2 Principal hydroforming tooling diagram ............................................. 2
Figure 1.1 Hydroforming Process Sequence ...................................................... 6

Figure 2.1 Shell element at original time t0 and current time t

(Total Lagrangian) ...................................................................... 11
Figure 2.2 Shell mid-surface at reference time °t and current time t

(Total Lagrangian) ..................................................................... l 1
Figure 3.1 Representative Meridian for axisymmetric case ................................... 16
Figure 3.2 Geometry of a straight segment ...................................................... 17
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 .............................................................. 29
Figure 5.5 Coulomb friction model ............................................................... 31
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 ............................................................ 51
Figure 8.1 Current Radius vs. Tube Length .................................................... 60
Figure 8.2 Hoop Strain vs. Tube Length (Aluminium) ......................................... 61
Figure 8.3 Axial Strain vs. Tube Length (Aluminium) ........................................ 62
Figure 8.4 Thickness vs. Pressure ................................................................. 63

vi

Figure 8.5 Hoop strain distribution compared with Abaqus ................................... 64
Figure 8.6 Axial strain distribution compared with Abaqus ................................... 65

Figure 8.7 Hoop strain comparison for different yield functions and

Experimental Measurements ......................................................... 66
Figure 8.8 Current Radius vs. Tube Length ....................................................... 67
Figure 8.9 Displacement Stroke vs. Internal Pressure History ................................. 68
Figure 8.10 Hoop Strain vs. Tube Length (Steel) ................................................. 69
Figure 8.11 Axial Strain vs. Tube Length (Steel) ................................................ 70
Figure 8.12 Thickness vs. Pressure at maximum bulge point ................................... 71
Figure 8.13 Major vs. Minor Strain ................................................................ 72

vii

INTRODUCTION

Hydroforming offers a way to cut material and manufacturing costs while
improving product performance in a variety of applications. Advantages include weight
reduction due to improved part design, part consolidation where a single component
replaces an assembly, reduced tooling cost as a result of part consolidation, and improved
structural strength and stiffness of the hydroformed component. It thus reduces tooling,
part, and labor costs, while significantly improving product performance (Anon, 1997).

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

Tube hydroforming is receiving the greatest attention, especially in the auto
industry (V ari-Fonn 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 0.1).

In tube hydroforrning, internal pressure is applied to form a tubular metal blank
into a structural component having a closed cross section. As shown in Figure 0.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 filled with a hydraulic fluid. Internal hydraulic pressure forces the
tube to conform to the shape of the die cavities, producing a part having different cross-

section shapes.

 

 

 

 

Figure 0.1 — Representative hydroformed components and assemblies
From left to right: instrument panel support beam, engine-cradle assembly.
Photo courtesy Vari-Fomr, Woodstock, Ontario, Canada.

 

» , , . ' ' Base plate

   

Horizontal cvlinder
. Pivotable
Cylinder bracket

”3°.“ , _ ~ 1
Sailingpunch ._ ' ........-
Bottomdie ' I I

 

 

‘ g " . . I 3 Bottom tool
*‘ ' ' ’ Base plate

 

Table plate ASE component

Figure 0.2 ~- Pn'ncipal hydroforming tooling diagram
Courtesy Schafer Hydroforming, Wilnsdorf, Germany.

The thesis presents an axisymmetric analysis method for calculating strains and
expansion as the tube is subject to simultaneous internal pressure and axial displacement

stroke. Only free forming stage of tube hydrofomring is treated.

The present method employs a number of simplifying assumptions. The principal
geometrical assumption is that the representative meridian of the tube is initially straight.
All segments making up the meridian are assumed to be relatively thin and of constant
thickness. The deformation of the member is assumed not to vary along the cross-section
of the tube. Hence, the analysis can be considered to be an axisymmetric 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 fi'om a Total-Lagrangian
formulation of a general shell theory is employed. Axial stroke is incorporated into the
equilibrium equation using Lagrange Multiplier technique. The boundary fiiction
condition is introduced into the formulation in the form of penalty firnction, which
imposes the constraints directly into the tangent stiffness matrix. The Newton-Raphson
algorithm is used to solve the nonlinear equation. Force Limiting Curve equations based
on the onset of shear instability and on experimental measurements are incorporated as

the fiacture criteria.

CHAPTER 1
LITERATURE REVIEW

1.1 Hydroforming Concept

Hydroforming has been one of the fiindamental sheet metal forming processes for
quite a long time, having been developed at least since pre-World War II, 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 subfi'ames, such as radiator enclosures, space fi'ames, 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):

g—a

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

Tube hydroforming offers several advantages as compared to conventional
manufacturing via stamping and welding (Ahmetoglu et al, 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)
(t) tight dimensional tolerances and low springback
(g) reduced scrap
While numerous variations of hydroforming exist, the basic principle remains the
same: utilize fluid pressure to form a component. Both traditional high-pressure
hydroforming and pressure sequence hydroforming are presented in Longhouse (1998)
and Morphy (1997). Typical Hydroforming sequence is shown below (Fig. 1.1), where
internal pressure and axial feed are applied simultaneously to improve the material

shaping abilities.

1) Place the blank tube into tooling 2) Seal the ends and fill with
fluid

 

4) Open the tooling, remove the part 3) Increase pressure, move
punches
Fig. 1.1 Hydroforming process sequence (Leitlofl', F. 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) and Asnafi (1999).

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 finite
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 verified 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 justified 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 finite 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 finite element method on sheet metal
forming, explicit finite 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 hydroforming
processes. Classical Coulomb’s law is capable of describing only fiiction 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 defined. 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).

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 fi'om 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 and stretches of a shell element undergoing an
axisymmetric deformation was discussed in (Pourboghrat, F. et al, 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 flat element at the initial time to 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, defined in
the global coordinate system,

#le mu
~ 0 ~ 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

§:[Xl +21] (2.3)

W

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

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

d" dde_ d5 1 _[1+“.x.]_1_ (24)

where if = 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 I.

Then, by calculating the magnitude of vector 5 from Eq.(2.4) and setting it equal to

unity, if can be calculated to be:

"a” = (a a = 1.0 (2.5)

NI—

1.? = [(1+ u,x, )2 + Win] (2.6)

The principal stretch in the circumferential direction, 1.3, is calculated fiom the change in

the radial position of the element; i.e.

0_ Xr+u _ _u_
2.2—K X] )]—1+XI (2.7)

 

10

Nl

 

 

 

 

 

 

0
X2
A ‘15 time t
x
‘ U
~ w
z
X i
r u 1 X]
T r >
dS
time to

Figure 2.1 Shell element at original time to and current time «Total Lagrangian)

 

z-
I? 6 .
2 ds time t
4' ..........
Aw Au --------
Z "
(32 w
A? ~ 3! r
A A M“ W
dS time ”t
e, ..............
R ............
“’ r

 

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

ll

2.2 Principal Curvature and Stretch (Updated Lagrangian)

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

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

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

follows:

r=R+uA+wN

(2.8)

where R is the reference configuration at time °t and u and w are incremental

displacements defined in Figure 2.2. In Eq.(2.8), the unit tangent vector A and the unit

principal normal vector IV to the mid-surface of the reference are defined as

1%)
ll
Q)
all”
I

(hi)

,8

A7:

’i l

l

where K; is the curvature. Therefore,

117-21:

0

and

12

(2.9)

(2.10)

(2.11)

(2.12)

In:

To calculate the current curvature, In, the unit tangent vector “ = and the

l

..’__s

unit principal normal vector of the mid-surface of the shell “ are calculated.

n
Using Eq.(2.8), the tangent vector a is found as:

ar A A A A
az—izr :R +uSA+uA +wSN+wN (2.13)
~ as ~,s ~,. ' ~ ~,s ' ~ ~,s

By substituting fi'om Eqs.(2.9), (2.10) and (2.12) into (2.13), and rearranging, the
following expression results:

6r .. . . .
az-a—gzrs =(1+u’s-K,w)/~1+(w’S+Klu)N=c/:1+dN (2.14)

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

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

Ni—

i‘ Z ”a" = J“ = 1’62 +d2 =10 +21. —K.w)2 +(w,. +K.u)2]

(2.15)

The current length of the mid-surface of the shell in the radial direction, ds, is calculated

from the reference length, dS, and 2.1 as follows:
(IS = AldS (216)

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

a=—:——~ (2.17)

13

which, from Eqs.(2.l4), (2.15) and (2.17), shows that fr-[r = 0. The current principal

curvature of the shell, k1, can now be found as:

k,=—&-fi :—r 4‘1 =——a-ii (2.18)

where a is given by Eq.(2.13) and f1 S can be derived from Eq.(2.17):

Q.

a ,1, (d, 21+ c, A?— K,c2r— Kid/V) — 2,, -(—d 21+ c1?)
: -——-:- : ~ ~ ~ 2 ~ ~ ~ (2. 19)
.s S A,

 

13>

In Eq.(2.19), 31.3 is assumed to vanish within an element and the above expression
simplifies as :

(d, + K1c);1+(—c, + K,d)1i/
= — ~ ~ (2.20)
.s ,11

By substituting from Eqs.(2.l4), (2.15) and (2.20) into Eq.(2.18), the current centerline

 

1:)

curvature of the shell, k], can be found:

_ eds “dcs +[(1/112

k] 213

 

(2.21)

It should be noted that it is possible to recover Eq.(2.7), for total Lagrangian formulation,
From Eq.(2. 15). This is done by setting the reference curvature K1 equal to zero for a flat
sheet, replacing S with X 1 and treating the incremental displacement u and w as total

displacements in Eqs.(2.l4) and (2.15). Similarly, the current curvature of the shell, k1,

can be calculated from the initial configuration by setting K1 equal to zero, replacing A,
with 211" and S with X1 and treating u and w as total displacements in Eq.(2.21)

{—W.Xr uerXr + (1 + uvxr )w,Xl.\’l ]
k] z (103 (2.22)

14

 

The current principal stretch increment of the mid-surface of the shell in the

circumferential direction, 1;, is found from the geometry (in Fig. 2.2) as:

,1, = — (2.23)

The principal curvature of the shell in the circumferential direction, k2 , is found from
Fig. 2.1 (Flugge, 1973), as follows:

12 42-5"”) blitz—“"3 (2.24)

 

 

where w and s were defined previously and R is defined as (see Fig. 2.1):
R = X , + u (2.25)
Since it is known that as = 21.?Xm , Eq. (2.24) becomes

raw rawaX, raw w},

k, = ——— — ———————

Ras RaX, as :_,1,°R5)—rf:-2,°R

 

 

(2.26)

Finally, by substituting from Eq. (2.25) into Eq. (2.26), we get the following expression
for the principal curvature in the circumferential direction:

W

 

k, = —l (2.27)

.(Xr + 101?]

15

CHAPTER 3

KINEMATICS

Based on the general thin shell formulations in Chapter 2, the formulations for
finite strain and small rotation have been derived for the straight-line segment that are

used in the program. (Derivation in details can be seen in Appendix A).

 

¢
v

Y K
...,. :> .........

v2 vz

 

 

 

 

 

Fig 3.1 Representative Meridian for axisymmetric case
3.1 Kinematics Assumptions

For the purpose of modeling, all segments making up the meridian are assumed to
be thin and of constant thickness. The length of the straight segments should be greater
than ten times the thickness. The deformation of the structural member is assumed not to
vary along its cross-section for a given point on the meridian. Hence, the analysis can be
done for a straight meridian parallel to the axis of the tube.
3.2 Principal Strains

The representative meridian shown in Figure 3.1 lies in the plane defined by the Y
and Z-axes. In this formulation, the initial geometry of the meridian is specified by first
defining the Y-Z coordinates of a set of nodes and then defining the segments. Two nodes

must bind each segment.

16

 

 

 

 

 

 

 

VZ

Figure3.2 Geometry of a straight 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

defined by the location of the point 3 = 0 given by Y0 and Z 0 , the orientation angle 0, the

length 10 and the thickness t of the segment. The position F0 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 s and 2 directions respectively. The rotation of a

line is given by B

.6: ——= -w,. (3.1)

The coordinates of a point on a segment undergoing displacements u and w fi'om the
original configuration are (see Fig. 3.2) :
5: Zo +(s+u)sin9—wcos6+z(-cos€+,Bsin 6) (3.2a)

77=Yo +(s+u)cos€+wsint9+z(sint9+,BcosQ) (32b)

17

The axial strain component ex along the length of the tube is given by
8x = 83 + sz

where 63 is the membrane strain component given by
0 1 1

a = u + -—u 5 + —w

2 ’ 2 '

where the local curvature is given by

W

3

(1+w’5)2

The strain component in the circumferential direction is given by
0
a, = a, + 192
where a? is given by

0 _ 5-(2O +ssin6—zcost9) _ usin6—wcost9+zflsin6

 

gr
R R

and the local curvature K, is given by

3.3 Constraints

(3.3)

(3.4)

(3.5)

(3.6)

(3.7)

(3.8)

Since the whole length of the tube is made up 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-l) constraint equations need to be enforced. At a junction

between two straight segments having orientation angle 6, and 62 and components

18

u1 ,WI and u,,w2 , the two displacement compatibility conditions can be written in terms
of the displacement components in the Y and Z directions as follows:

u2 c0362 + w2 sin 62 — ul c0319l — w1 sin (91 = O (3.9)

u2 sin 02 —w2 c05192 —ul sin 01 +1121 0050, =0 (3.10)

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 fl, , compatibility condition is

flz-flr:0 (3-11)

19

CHAPTER 4
CON STITUTIVE EQUATION

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

uniaxial stress-plastic strain curve of the material is given by
6=K(E+§o)” (4.1)
where 6' 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 firnction allows for anisotropic yielding of the material. The yield firnction is
given by

2 2 _ 2
:0”: +03+R(0x 0'3) _0_::0 (42)
1+R

 

f

where a, 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.,

6" 3 0y , where 0y is the initial yield stress of the sheet obtained from a uniaxial tensile
test. Beyond the elastic limit; i.e., E > cry, 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):

— -1

<1
IN
3":
e'c
It“

 

0' = L— ° Z :1) (4.3)

+
“B
2's
is

 

 

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

respectively, or is the stress tensor, 6 is the material flow strength, h(= 66/85) is the

~

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

~

 

 

 

), where Harri/6g" = Pip/6g : art/6g , is the second order tensor

1%: 6¢/6 g/“ags/ag

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

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

1Q<

R: :D—P:

 

E:

lth

(4.4)

P:0'

I "c

:L:P

The fourth order elastic tensor L(= LW) 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 briefly. 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 I]1D of a discrete

system for steady-state analysis,

11, 2%UTKU-UTR (5.1)

The constraint equations can be expressed as
CU = 0 (5.2)

In the Lagrange multiplier method we amend the r.h.s. of (5. 1) to obtain

11; = é—UTKU—UTR + xTCU (5.3)

22

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

invoke that 5 I]; = 0, which gives

T U R
K C = (5.4)
C 0 2. 0
The second of Eq.(5.4) is Eq.(5.2), the equation of constraint. Eq.(5.4) are solved for both
U and l. The A ,. 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
1 T T 1 T
l];=—U KU—U R+—t at (5.5)
2 2
in which t 2 CU , and t = 0 implies satisfaction of the constraint. The usual potential I],
of the system can be augmented by a penalty function étTat , where a is a diagonal

matrix of penalty numbers a, . The condition 6 II; = 0 now yields

[K + CTaC Kn} = {R} (5.6)
If a = O, the constraints are ignored. As a grows, U changes in such a way that the
constraint equations are more nearly satisfied.

It is difficult to assess which method, the Lagrange multiplier method and penalty
method, is better in the finite 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 3. 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 first 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 successfirl equation
solving. Their implementation can be as easy as assigning a high modulus for an element
already in the program. Thus, we use penalty firnctions 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 first 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 fi'om 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-defined 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). Ifa 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

 

/////////////////IL

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
define 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 1 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, N 1 and N j ,

as follows:

Ni+Nj

n = W (5.7)

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

F143}; zN=,/z: +1 (5.8)

ZN

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 fiom 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)
x 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, 11, 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:

t—1 l 510
"zfiz, (')

where ZN and Zx are defined 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 define the normal of a node by
differentiating the Fourier series without any numerical difficulty.

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, Ak is one of the contact nodes and n is its

normal vector.

27

If (AkBH xn)-(AkBi xn) > 0 (5.11a)

and (Akai x n)-(AkBH, xn) < o (5.1 lb)

Then, we define segment Bi Bmto be the segment associated with contact node Ak.

n
Akf Sheet Segment

/Die Segment

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
found. The nodal coordinates are then modified 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,
11 = (nx,nz), 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, fiom the nonlinear equation:

S(t-nx)—t-nz=O (5.12)

where S(x) is the tool surface equation.

28

5.3.2 The Directive Method

As shown in Figure 5.4, PO 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:

'0—15+1_>1§=5}1'+E (5.13a)
1713' = r), If? , 7413 = —q,n (5.13b)
IP

 

 

where 711,772 are scalar. Once I], , 772 are solved from Eq.(S. 13), the coordinate of B could

be determined.

 

 

Figure 5.4 The projection Algorithm

29

5.4 Implementation of Constraint
A node located on the sheet at time t +At , is constrained to move in the tangent
direction defined by the trial solution, A u’. The tangent direction incorporates the Zx in
Eq.(5.9). The constraint on 5 u = (5 11K ,6 uz), for contacting nodes is then:
6 u - n = 0 (5.14)
which by substituting in Eq.(5.8) we will have
5u z = Zx ~6u , (5.15)
Actually, Eq.(5.15) is the alternative constraint equation, Eq.(5.2), in the contact
problem.
5.5 Friction Model
Friction is a complex physical phenomenon 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 simplified to two idealistic models.
The most popular fiiction model is the Coulomb Friction model. This model is used for
most applications. The Coulomb model (Marc's Theory and user information) is:
afiS-an-t (5.16)
where o n is the normal stress
a f, is the tangential(friction) stress
11 is the friction coefficient

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

v . . . . .
t = l—il- , vr rs the relatrve slrdrng velocrty.
v!’

The Coulomb model is also often written with respect to forces

30

1't s—ufn-t (5.17)
where f, is the tangential force, fn is the normal reaction.

For a given normal stress, the fiiction stress has a step firnction behavior based upon the

value of vr (see Figure 5.5).

f, or c f,
Slip T

 

 

 

T Slip

Figure 5.5 Coulomb Friction Model

This discontinuity in the value of of, can result in numerical difficulties so a

modified Coulomb fiiction model (Oh, 1982) is implemented:

 

 

A A
f3 =—mk v, e_—2-mktan“[ 1),] (5.18)
IAV, It 110

 

where Av: is slipping velocity, m is friction factor, k local flow stress in shear and u0 is

very small positive number compared to Avs.

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.

31

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 free node should be zero).

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

32

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 finite 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'HAITUJHAteU ”de : J‘HNflBa‘i t+ArdV+ INN-fisalfnmcis (61)

r+Aty HAIV t+Ats
Where l.h.s is the total internal virtual work and r.h.s is total external virtual work. The

”wry. are the Cartesian components of the Cauchy stress tensor, the 1+ A, ea. are the

Cartesian components of an infinitesimal strain tensor, the ”mfg and ”A’fi" are the

1

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

611,. is the ith component of the virtual displacement vector.

33

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
”A’ R - "A’F = O (6.2)
where the vector ”A’ R lists the externally applied nodal point forces in the configuration

at time t+ At, and the vector "A’ F lists the nodal point forces that correspond to the
element stresses in this configuration. 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 timet + 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‘ U0"); then a Taylor series
expansion gives

f(Ue):f(r+ArU(tl))+[§_] (Uo_r+ArU(i—1)) (65)

t+AtU(i-l)

8U

 

where higher-order terms are neglected. Substituting from (6.4) into (6.5) and using (6.3)

we obtain

0: t+AtR_ ”INFO-1) +[§R;_ 6F] (U.-!+AIU(F1)) (66)

HA! U(i—l)

auBfi

 

when we assume that the externally applied loads are independent of the displacements,

which means ER 2 O , (6.6) becomes
6U t+AtU(i-1)

34

fl (Us-HAgUa—n): r+AtR_ t+AtFi—1 (67)
6U ”mun—1)

We now define

AUU) : U‘ _ t+AtU(i-l) (68)
and a—F =“A’KU‘” (6.9)
6U 1+A1U(1—1)

 

where "’3’ K0”) is the tangent stiffness matrix in iteration i-l, which we assume to be
nonsingular. Thus (6.7) can be written as
”A’ K""’AU"’ = 1+4: R — ”’1‘ F0") (6.10)
Since (6.5) represents only a Taylor series approximation, the displacement
increment correction is used to obtain the next displacement approximation
”A’U“) = ”A’UU") +AU") (6.11)
The relation in (6. 9) and (6. 10) constitute the Newton-Raphson solution of (6.2).
6.3 Four Types of Modeling
Four different models, i.e. for pressure loading, axial force, displacement stroke
and fiictional contact, are derived respectively.
6.3.1 Pressure Loading Modeling
Based on kinematics and constitutive equation discussed in Chapter 3 and Chapter
4, the principle of virtual work is used to satisfy the equilibrium equation and takes the

following form for pressure loading:

21“, [(a; as; +6; 55; )R‘dL’ +itj 6C}. = 5m; (6. 12)
1'21

i=1 L‘

35

where I is the number of segments along the tube length, I is the number of constraint

equations, ,1]. are Lagrange multipliers, C J. are the constraint equations(from Eqs.(3.9-
3.11), 5W3: is virtual work due to external force (pressure), dL’ is the incremental length
for the straight segment and Ri is the radius of the tube at current segment.

The displacement component w, ,u, for each straight segment are approximated
using the following Fourier series expressions:

mas, mzs,

 

 

NI
_ x i
w, —a0+ E ancos

 

NI
[0 + 2,6; sin 10 (6- 13a)
n2] i ":1 i
i N‘ r‘ "”3: N! i - "7131'
11,270 +27" COS [0 +25" Stu—15— (6.13b)
n:l r' n=1 1‘

After substituting from Eqs.(6. 13a-b) into the principle of virtual work Eq.(6.12) ,

a nonlinear expression of the following form will result
f(c, P) = 0 (6. 14a)

where

c={a5,a;, $.7s,r;,6;,1,} (6.14b)

Eq. (6.14) should be solved for c for given values of P.
Since Eq. (6.14) 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:
[’Sfldfi- [6]: Fg- Fm (6.15)

where K is the stiffness matrix, 6c is the incremental c and R is the force residual.

36

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

~

the element stress in this configuration ( Internal forces).

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

W;,=:Ipn 6(Au~)RdL‘ =ij fi‘a—L—g‘z)ac RdL" (6.16)

i=1 Ll, i=1 Ll

~

where 6(Au’) is the virtual displacement increment having two components 6(A 145')

and 6 (A v‘) , ft" is the unit outward normal to the segment Lip .

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

17;: pzjn" ———)a(aA:i‘R dL" (6.17)

i=lL‘

~

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

 

’ am") 6(Au") .. 62(Au‘) . .
KP = '— R'dL’ 6.18
E“ p: [ 60 ac +n 62c p ( )

N

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

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

[K+K§u][dc]:[1~i]=EEu+FEfl,— Fm (6.19)

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

stiffness matrix K g, has nearly no effects on the solution that led us to simplify the

E

37

formulation. Actually, if we approximated ii‘ 6(Au’) z 6(Aw‘) , the whole formulation

above can be simplified. The proof can be seen in Appendix B. Since 2 O , there

 

62 (Aw')
2

C

~

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

Raphson expression (Eq. (6. 19)). Both methods are implemented in the code; the results

are almost the same.

6.3.2 Axial Loading

Application of axial force with internal pressure is critical to improve material
shaping abilities for deeper draw configurations and higher expansion in localized
regions. Axial compressive feed F is modeled as an external force at both the ends of the

tube as follows:

i=nse ,k= 1 i: use ,kzl
5W5; = (—1)"Fi‘6(Au‘) g =-—( 1)‘r Ft a_(____§uc )6 g (6.20)

~ 1:11:20 0 ~ i=l,k=0

A:

where S(Aui ) is the virtual displacement increment having two components (S(Awi )

and 5(A v‘) , 5‘ is the unit tangent vector along the segment L; and nseg is total number

of segments. Since the force is compressive at both ends, when i =1(i.e. first segment),
k = O and when i = nseg (i.e. last segment), k =1.

6.3.3 Displacement Stroke Modeling

The displacement stroke at both ends of the tube is included in the equilibrium equation
using Lagrange Multiplier technique. The stroke is specified in the u direction (along the

length of the tube).

38

For the left end ofthe tube(i =1):

ll

 

_. *
s=0,i=l — ll

             

 

_ *
s O,i—'l — u

.11." $1.]

:[73 +271.

NI
1 l
=>ro+Zrn=u*

nzl

Ni
=>73 +27}. -u*=0 (6.21)

n=l

For the right end of the tube (i = nseg ):

 

   

 

         

u “in... z _u ,.
:> [73 +217; "2:9. :1 561%] WW8 : _u*
375% + 2(1),)”, =—u*
2’ 70 Q + N21- 1)"7,',“‘3 +u*_ — (6.22)

Eqn. (6.21) and Eqn. (6.22) are incorporated in the equilibrium equation using

two Lagrange Multipliers.

6.3.4 Frictional Contact Modeling

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

39

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

    
     

\

Nodes on !
The TM

 

 

 

I
l
I
I l
i Po ‘ : P1
1 :
! I
! i
! I
! i
y
3mm ”in (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 finally (d) shows how the equilibrium

shape is obtained after several iterations.

40

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

principle equation (6.12) as following:

J+2

$1,101.66; +01,58j,)R’d1’ +2166, 2
i=1 1:]

tweak—J + Pi] fit a(Aui)R"d[jp +212IT,6(Aui)RidL'} (6'23)

(—1)"Fi"5(Au")
~ 66.

 

r=l,k:0 i=1 L1

where r, 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. Afier equilibrium is satisfied, 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 first 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 modified 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
.___fi

Adjust time

 

 

 

 

 

 

 

 

 

 

 

¢

 

 

 

 

 

TIME loop end

 

 

 

Figure 6.2 Flow chart for contact iteration

42

6.4 Fracture Criteria

Force Limiting Curve equations based on the onset of shear instability and on
experimental measurements are incorporated in the program as fracture criteria:
6.4.1 Shear Instability:

Force Limiting Diagram based on the onset of local shear instability, using

Hosford yield criterion and critical shear stress under plane strain (Lee & Kim, 1989), is

“1'

" k2 + p)V~ (1 + p)“/5»(a;,) (6.24)

given by:
(na- "—1 )/

1111.1 1 “11.11.

“)Hawg’ [I + R(1— at)“1 |1 - a

l

1....
,6

 

 

 

(1+1ar +R|1-a

 

 

where
R = Anisotropy parameter
n 2 Strain Hardening Exponent
a = Stress ratio

p 2 Strain ratio

,6 = 1 + R X"

a 2 Order of Yield Function (=2 for Yield Function used)
a: = Maximum principal strain at shear instability

8].}, = Maximum principal strain fi'om plane strain test

The effects of thickness, surface roughness and strain gradient through the

thickness are not accounted for in the equation above.

43

6.4.2 Experimental Measurements:
The Force Limiting Curve derived from experimental measurements for Hot-dip

galvanized DP6OO (Asnafi, 1999) is used for Steel in the second verification problem.

44

 

CHAPTER 7
PROGRAM STRUCTURE

This chapter describes the input file, control flow of the main program and the
output files.
7.1 Input File

The input is divided into two parts. The first part contains all the geometric,
numerical and material information for a run. It is written in an input file that is read by
the program prior to specification 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

 

NODES. RADIUS

 

 

 

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

RADIUS: Tube radius

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

 

YN, ZN

 

 

 

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

45

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

4. Number of segments.

 

NSEG. MU. DELTAV

 

 

 

N SEG: 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.

6. Segment parameters.

For straight segments input the following record:

 

INODE. JN ODE. T. NCOEFF. IS. IZ. NC

 

 

 

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

JN ODE: Second node of the segment (s=lQ).

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.

46

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

NC: Number of contact nodes on segment (t).

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 IfNSPC is positive, the constraint is enforced at the node with s = 0.

o IfNSPC is negative, the constraint is enforced at the node with s = 10.
IVARPC: Variable which will be constrained.

0 To constrain the rotation ,6, input 1.

0 To constrain the displacement component it, 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.00000001
47

ITEMAX: Maximum allowable number of iterations.
JEVAL: Multiple of the iteration count at which the Jacobian is evaluated. IfJEVAL =
l 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.

11. Hardening Model.
I_choice, A, B, C

 

 

 

 

I_choice :
o I_choice=1: Ludwig-Holloman Model
0 I_choice=2: Voce’s Model

A, B, C : Material constants for Voce’s model

48

11. Output option.
IWR

 

 

 

 

o IWR = 1: print coefficients and Lagrange multipliers.

o IWR = 0: do not print coefficients and Lagrange multipliers.

7.1.2 Loading Input
Loading is prescribed in an incremental manner.

Loading parameters

 

NINC. TINIT. TINC. PINC

 

 

 

NINC: Number of increments / Load-steps
TlNIT: Initial Tension
TINC: Tension increment.
PlNC: Internal pressure increment.
Displacement stroke increment is calculated by “Ioad_curve.f’ based on the loading

curve programmed.

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.
Subroutine build: evaluate the system of equation and the Jacobian.

Subroutine constr: add constraint contribution to system of equation and Jacobian matrix.
49

Su
Sul

F01]

Subroutine decomp: find 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 coefficients of the series expansions at one point along 3
for straight and are segment respectively.
Subroutine d2cu rvature: evaluate the second derivatives of curvature.
Subroutine incm: increment of all variables.
Subroutine solve: solve the nonlinear equations.
Subroutine stndsp: find strains at the integration points along the local coordinate.
Subroutine stsstrr: find 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 firnction.
Subroutine violation: check each integration point to see if it violates contact.
Subroutine bringback: bring each contact point to die surface.
Subroutine load_curve.f: calculates stroke increments based on current pressure.

Following is the flow chart of the program:

50

 

 

input

 

 

 

store

 

 

 

 

 

 

 

 

newrap

 

 

 

 

 

 

 

 

 

dieinfo

 

 

 

 

 

 

 

 

 

 

 

 

 

build

 

 

 

 

 

 

 

 

 

 

constr

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

decomp

 

 

 

solve

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

contact

 

 

 

 

incm

 

 

 

 

 

 

 

 

 

Figure 7.1 Flow chart of the program

 

 

 

 

 

 

 

 

 

 

Load curve
stndsp disp
dstrn
stsstn cnstutv
violation bringback
setpenalty
coloumb
forcecheck nodalf

 

 

 

 

 

 

 

 

of
111:
L11
C01

off

7.3 Output File
The program has two output files: betges.out and betges.ooeff. The contents of each file
are described below:
betges.out: This file is formatted and contains the output of the program in the form of
tables. The first part of this file is written prior to the specification of the loading and
contains an augmented echo print of the input. It is useful to look at this part of the file to
ensure that the input file was read correctly. The second part includes a list of the
curvature, axial strain, tension and internal pressure at each loading step. Ifthe variable
IWR is set to 1 then the file will also contain the coefficients of the trigonometric series
expansions for each segment.
betges.coeff: This file 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 meridian prior to loading. The file is organized as follows:
Line 1: (Format: i2)

Number of segments.
For each segment, the file contains the next set of three records:
Line 2.1: (i3, e135, a4, 3e13.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. The location of

the segment is partially given by the coordinates Y0 and Z 0 as shown in Figure 3.2.

Line 2.2: (4e13.5)
Coordinates of the nodes bounding this segment in the following order: Y-coordinate

of first node, Z-coordinate of first node, Y-coordinate of second node, Z-coordinate

52

$1

Ln

111

For.
offi)

Lme

Lme3

Lhe}

Lme33

of second node.
Line 2.3: (2e13.5)

Sine and cosine of the angle6 defined in Figure 3.2 for straight segments. This
concludes the part of this file which contains the original, undeformed description of the
cross-section. The second part of the file contains, for each loading step, additional
information necessary to describe the deformed shape of the meridian. For each loading
step, the file contains the following:

Line 3.1: (i2)

Step number.
Line 3.2: (i2)

Segment number.

For each segment along the length of the tube, the output file 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))
( ,6" , n = 1, number of terms).
Line 3.3.3: (6(e22. 14))
70 ,( y" , n = 1, number of terms).
Line 3.3.4: (6(e22. 14))

(6”, n = 1, number of terms).

53

 

J Jib LII N

'O :--_- to g.

A

CHAPTER 8
EXAMPLES
The AXHD (Axisymmetric Hydroforming) program has been validated using the
following two examples:
8.1 Bulging of an Aluminium Tube
This example concerns bulging of a straight Aluminium tube under pressure and
stroke loading. The tube is 8.0” long with an outside diameter of 2.0” and a wall
thickness of 0.049’ ’. The following are the material properties used for Aluminium 6061-
T4 material in plasticity (Ludwig Holloman) model:
E = 10300000 psi, v = 0.3, 0,: 18730 psi, R = 0.8179, K=69183 psi, N=0.2646 and 80:0
The input file to run this case is shown below (the definition of input parameters
at each line is given in Chapter 7):

Title: 'Aluminium Straight tube subject to 2030 psi pressure and 0.1181 in. stroke’
5,1,000

4.0,1.000
3.0,1.000
2.0,1.000
1.25,].000
0.0,1.000
4,0.1,0.01,0

Stt’
1,2,0.049,5,10,5,41
str
2,3,0.049,5,10,5,41
str
3,4,0.049,5,10,5,41
str
4,5,0.049,5,10,5,41
3

1,1

1,3

-4,2

0.000000001, 150, 1
10.3e6,0.3,18730,0.8179,69183,0.2646,0

54

1,0,0,0

This tube was subjected to a pressure loading of 2030 psi and 0.1181” axial
displacement stroke. The deformed shape and strain results obtained from AXHD were
compared with the results from ABAQUS and experimental measurements.

To study the effect of Aluminium’s anisotropy on the strain distribution along the length
of the tube, the ABAQUS code was run with the following three yield firnctions:

0 Von Mises (Isotropic)

o Hills (with transverse anisotropy)

o Barlat (with planar anisotropy)

Due to the inability of ABAQUS/ Standard to account for the current thickness and
mid-surface offset for shell elements in contact analysis, the same simulation (with S4R
shell element and Hill’s yield function) was run with ABAQUS/Explicit and
ABAQUS/ Standard. The results were compared and it was found that the true strain
distribution along the length of the tube is offset by a finite magnitude (within an
acceptable error limit) when the mid-surface (instead of the external surface) contacts the
die.

Figure 8.1 shows the intermediate shapes of the deformed tube predicted by AXHD
program at various pressure levels. Similar to the experiment the tube was pressurized to
the maximum pressure level of 2030 psi while axially being compressed to a maximum
stroke of 0.1 181” (3.0 mm). The predicted deformed shape of the tube at the maximum

pressure of 2030 psi is compared against the actual one. The radius of the predicted

55

deformed tube is 1.17” (29.72 mm), while that of the actual tube was measured to be

1.25” (31.75 mm). This corresponds to an underestimation of 6.4%.

Corresponding distribution of hoop and axial strains along the length of the tube are
shown in figures 8.2 and 8.3. The maximum hoop and axial strains occur in the middle of
the bulged tube and they are 17.54% and —7.4%, respectively. Figure 8.4 shows the
thinning of the tube at the maximum bulge point (middle of the tube) as a function of the
applied pressure. It is interesting to note that thinning accelerates exponentially beyond
the pressure level of 1000 psi. The thinning strain at the pressure level of 2030 psi is

10%.

Figures 8.5 and 8.6 show a comparison of hoop and axial strain distributions
predicted with AXHD code and ABAQUS/Explicit code using axisymmetric element
(SAXl) and four-node shell element with reduced integration (S4R). The hoop strain
distribution predicted by the AXHD code matches those of ABAQUS very well, despite
the fact that a S4R shell element takes into account the shear deformation. Similarly, the
axial strain distribution predicted by the AXHD code again matches those of ABAQUS
very well. The two spikes (at 4:15") shown in figure 8.6 for the ABAQUS/Explicit code

are attributed to the dynamic effect.

Figure 8.7 shows a comparison of the measured and predicted hoop strain
distribution, as a function of material models used. For comparison, in the numerical
simulations isotropic and two types of anisotropic material models were used. Hill’s yield
criterion was reduced to Hosford yield criterion (Eqn. 4.2), for transverse anisotropy, for
both SAX 1 and S4R shell elements and then compared with AXHD results. Barlat’s 1996

56

yield function was also used to account for planar anisotropy. As expected, the isotropic
material model (von Mises) had the worst performance and significantly underestimated
the magnitude of the maximum hoop strain (14.9%). Those models (Hill’s SAXl, S4R
and AXHD) accounting for transverse anisotropy showed some improvement but still
underestimated the measured hoop strain (17.4%). Although not perfect, the best strain
predictions were obtained with Barlat’s yield firnction accounting for planar anisotropy
(20.4%). The coefficients used in Barlat’s yield function for the aluminum 6061-T4 tube
were obtained from reference [Pourboghrat, F. et a1, 2002]. The maximum hoop strain
measured experimentally was about 25%. The high error percentage between the
experimental and FEA results can be attributed to the high anisotropy of Aluminium.
Also, the maximum attainable pressure in the used hydroforming equipment is 30,000 psi
while the working range in Aluminium’s case is 2000 psi. So, even an error of 0.5% (150
psi) in the pressure reading from the equipment will cause a high change in the strain and
deformation results. This could account to the discrepancy in the simulation and
experimental results too, since the experimental values obtained could well be for

pressure beyond 2000 psi.

8.2 Steel Tube

This example is concerned with the bulging of a straight steel tube under pressure
and displacement stroke. The tube is 8.66” long with an outside diameter of 2.36” and a
wall thickness of 0.05787”. The material properties used for hot-dip galvanized DP600
(HG/2140) in plasticity (Ludwig Holloman) model are:

E = 31465000 psi, v = 0.3, 0,, = 351040 psi, R = 1.0, K=143260 psi, N=0.182 and 80:0

57

Simula
(HO)
meaSUru
figlues

18811115 fc

The input file to run this case is shown below:

Title: 'DP600 isotropic material under 5220 psi pressure and 0.24 in. displacement stroke’
7,1.146

4.33,].146
3.608,].146
2.8867,].146

2. 165,]. 146
1.4433,].146
0.722,].146

0.0,]. 146
6,0.0,0.00],0

str
1,2,0.05787,6,]2,7,41
str
2,3,0.05787,6,12,7,4]
Stt'
3,4,0.05787,6,12,7,4]
str
4,5,0.05787,6,12,7,41
str
5,6,0.05787,6,12,7,41
str
6,7,0.05787,6,12,7,41
3

1,1

1,3

-6,2

0.00000000001, 150, 1
31.465e6,0.3,51040,1.0,143260,0.182,0
1,0,0,0

0,0

The stroke vs. pressure loading history and experimental results for this
simulation are obtained from [Asnafi, N. et al, 1999]. The Forming Limit Diagrams
(FLD), based on the onset of shear instability [Lee & Kim et al., 1989] and experimental
measurements [Asnafi, N. et al, 1999] are used as fiacture criteria. As shown in the
figures described below, the results from AXHD are very close to the experimental

results for this isotropic material.

58

an
sho
mee
(FL:
of 52

middl

Figure 8.8 shows the intermediate shapes of a bulged steel tube at various
pressure levels. Those deformed shapes were predicted by AXHD program using the
displacement stroke-intema] pressure history shown in figure 8.9, obtained fi'om
analytical equations given by [Asnafi, N. et a1, 1999] in order to maintain a proportional
axial to hoop strain ratio of ——0.5. The final displacement stroke was 0.24” (6.] mm) at the

maximum pressure level of 5220 psi in the tube bulging process.

Figures 8.10 through 8.12 show the evolution of hoop strain, axial strain and the
maximum thickness reduction as a function of the internal pressure. The maximum
thickness strain at the maximum pressure level of 5220 psi is 15.57%. The maximum
hoop and axial strains predicted by AXHD program were 28% and -10.9%, respectively.
These predicted values compared very well against measured maximum failure hoop and
axial strains of 29% and —10% reported by [Asnafi, N. et al, 1999]. Finally, figure 8.13
shows the predicted major and minor strains by AXHD program against experimentally
measured [Asnafi, N. et al, 1999] and numerically predicted Forming Limit Curves
(FLC) based on shear instability [Lee & Kim et al., 1989], at the maximum pressure level
of 5220 psi. These results show that the steel tube is about to the fail (burst) in the

middle of the tube where maximum bulge occurs.

59

3 .5“.
59.6.. 3:» .555»: .5653 855m 52 .5551“. 585.65» .555“:
.855. 3:... 65 .6 59.3 .55..

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

3 N m... F no o no- F. 3- N- on.
o , >
No
.5 an}
to .5 was"...
a .5 95...":
m $-88 5555:?
w . an Sing 11
w o o
w an 821". A
n 1515111111,. Fungus. .52
N. no an 82".“. 11
m. an 828 A
an 80an A
V Ilfifulrll: .11111111
3
ma

598.. .a> matey. E250

60

mod

cod

. no.0

_ «.0

name doou

_
f
A
l «to
L :6

L 3.0

2.0

In.» .FIIII- r 1111

a.» .2".

50:01. Etc—Io; £82.35 38.50.22 35:79.. Jam-£03... EEC—«cu

€2.65. 33. 65 .6 593..
n a P o F.

 

 

 

8 8o Fun. 4

Ba 83......”— A

 

 

 

 

8 83am 4

ed 25 Fun. 1 5

 

in ooouua A

 

 

 

.5 ans.
.5 23....
.5 «3...»...
45-38 555532

 

 

 

 

 

 

.8 fl
. omow n. A

 

59.6.. .u> £85 goo:

61

“19118 I?!“

no .2".
52.3 .655"... .5255 .555"... 685.65... .555":

 

 

 

 

 

 

 

32.65. 33 2.36 :55:
N F o F- N- a- v.
\111\\
.5 our.
.5 2.3"...
.5 Scan:
”5.58 5555:?

 

4 a
8.0. L.
883. 1
So- Sam"; 4
08 W6 A
8.? 1
8» Va 4
3o.
3.?
8 o- 7d ooiua 11
mod-
3.?
ed 08 Fun. A

 

 

5.0

 
 

_

59.6.. .2 5an .52

 

 

 

undquan _u_x< .an.

 

62

comm

55:0.— .QEE Io.— ..GuoEuE costs» 22 _u_::.lou.uuocxu_:._. .055":

ooow

323.513.. 38...:

v.» .2".
«585 fiancee—«e d 6 6352...... 2250!

.35 2:32.“.

oom _. 000 v

_

com

 

 

.5 on...
.5 £3»...
.5 63.6...
45.58 555532

 

 

 

 

2:32“. .2, :52. 02:9 E:E_xaE .3 3052......

 

 

mead

__ £36

_ 33°

33.0

nomad

m Rod

(sorrow) sseuxauu

mwmod

63

1‘ ‘l‘c:

arms d00H

m.» .9“.
.5305 2.5.823... 056553.365?
59.3 3:... .555"! £82.35 32.5922 355%.. 635.02... fiEEu:

 

 

 

 

 

 

 

 

av u: c a
a. a ~ P. .. scan 33 p F- m- n- 4-
O _ . a .. , a . . . .. . . . . . Z . I
85
3o
86
3°
to . Sonar.
.5 83......
9.0 .5 93...»:
p-38 5555:?
Ed
253% 1 /
dune A
2...
5 6. x< a < <1
No

.2523 .u> 9.39: 585 goo:

“19118 I?!“

a.» .5".

E802“. 055.283: org—5:55.205?
50:3 =5.an 0825.0 35.50 Es. .2555. .3352... .325»:

.855. 59.3 33-
3. 3 m... 3-

93341

5316334
2.0.0.. A

3.33131 3

000.0-

 

 

. mood.

63o-

63o.

39o- ,/ \5
a _

Rod-

9.0.0-

mood- 1 - a r

A6.00.0

6 3.6 .3662 .2, 9.32. 526 .52

 

 

.5 on...
.5 33.6
.5 63.9.:
45-38 5555::

 

 

65

is 1100”

.3 .5”.
EuEefi =30 0085.52.2- §<m 5:038:25 03:02. «55:5 :30 3.360.100
50:3 3:... 355%.. c.2350 32.5955. _:_u_:_no0 £35.02... .555":
330:: 50:3 3:...

0.0 0N oé 0.0 0.0- v. T 0N-

   

.5 on
.5 23"....
.5 25...":
5.53 555532

 

 

 

  

 

£235-59 5...: A

 

cBiEEvm: 51.5.: A

 

52mm: iEwémmaxm A
”No 3334‘ E2: 3:05am 0:: 6:28:3— 0_o_>
22220 .8 3303500 58.0 :00...

womodH

     
  

66

a.» .2“.
55:0.— 33. 35:79.. .bwquED ESE—lot Jaw-3.3:... EEC—no“

 

 

 

 

 

 

 

 

 

 

.35.... 59.3 3.5
on N 3 F 9° o .3- F- n. F- N-
an 88": A .5 2.» n J
.5 23 u on
a: 89.": A .5 2:3... u :
.SEQ... 82:
a: 83": A
a: 88".“. A
a: 88": A
an 8%": A

 

 

  

 

 

 

E has»...

 

 

 

 

 

 

 

 

:60 £0 05 :_ one. 03::an :5 :0 cans» 02.3.39

 

 

«a o.
O O
(”noun 8mm: wanna

P

N
F

V.—

8..

w...

67

oooo

 

ooom

a.» .2“.

593.. can... Etc—non. £30520 33.5%.. 6:05.02... 33.5.3

83.

:3. 2:32..

coon

 

5:. 2.0.8..

ooom

ooo _.

 

 

.5 5... u S
.5 23 u 2.
.5 253. u 3
Gina... 2.2.:

 

 

a mod-

1 mod

0

dsu

o'
axons wawaoel

I m...o

(«noun

1
a!
o

. mud

 

r 9.0

68

2.» .5...
:55... 83 .555 u 9. .3855: .555 u .u 685.55 3...... n o.
.855. 3.... 2.. .8 5:5...

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.. n. v-
.8 88"“.
.8 Senna.
.8 88":
.8 83am
.8 83".".
.5 S... u c..
.5 «an: n on
.5 .23... n o.
.8 888. side... 88:
.8 omwmnn.

 

520.. a> Exam :8...

mod-

mod

3
1s doou

36

N6

mud

nd

69

 

3... a...
50:3 3:... .0551... c825... 3558.0 635.95... .255»...
.855. 33 o... .o 58....
..

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.. m N o .- N. m- _ .1 m-
j... . . . i ,.. . . . . _ NF.°I
.8 888 r >
3-
.8 88"..
8...-
.8 88a“. .5 S... n 7.
.5 N8.“ a on ‘ 8...- m
.5 3.8... u o.
3380... 82.:
.8 83.8. 3...-
.8 88.8
f... _ .w... N0.0:
.8 88".. .. \,
f1... 4.»...
.8 88 n. . m.
. v .. i n 4 o m
m.
8...

58.... 2. 58.... .52

up... .5“.
2......» 8.5.. a goo... 8.2.5.83... 685.55 2.2.3...
58.... 8.5 35...... 52:5: 355.... £85.25 355-...
an... 2:32..

 

 

 

 

 

 

88 88 8o... 88 88 8.: o
. . . . . . o

5...
8...
8...

.5 S.» n 0.. .

.5 «93 u on 3 o

..5 83.98.850.555“; 8%»... .5 .28... u a
8380:. 8......

8...
8...
B...

2:32... .2. «35.05...

71

(swam) ssaunouu

2... .9".
598.. 3:... 355%.. £82.35 3...:.uo0 £35.05... 3.55%.

 

 

 

 

 

 

 

.8... 2.. .o 5.... 5.85.5... 2.. 9.2.. «.52. 5.5.8.5 .a. 5....» 35....
.8... 3.....- 8...... 8...... 2.3.. 8.....-
90 O
O
8... o
o .5 3... u 9..
.5 «8.. a on
2. o .5 .88... u...
33.56.... 8...:
0
.t...
m .
m 0
E .
m. m o
u In!
8...
no
8...

:3...m .052 .m> .53....

 

3.6.

72

CHAPTER 9
CONCLUSIONS

An Axisymmetn'c Hydroforming finite element analysis program has been
developed based on the formulations presented in the thesis. Fourier series interpolation
fimctions, which reduce the size of global stiffness matrix and the number of variables
drastically, are employed for approximating the displacements. The program is very
efficient in predicting deformations for free-forming stage of tube hydroforming under
simultaneous action of internal pressure and displacement stroke. Failure model based on
FLD is incorporated in the code. The results from the program are compared with the
results from commercial code (ABAQUS) and experimental measurements and are found
to be in good agreement.

In the future, this code is to be extended for 3D analysis of tube hydroforming
process. It is also planned to develop an ability to optimize the loading curve (pressure

vs. stroke) for improving the material shaping abilities.

73

APPENDICES

74

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 K 1 =0, then fi'om Eq.(2. 15) ,
2. =[(1+u,.)2 +w§1i =U+u§ +221. +w§1i (1)
Let's approximate the square root as:

(1+x)’g z1+%x—~- forx<<1, Eq. (1)becomes:

 

A, z1+u§ +%u:. +%wfs (2)
Since, A, = ds/dS and the engineering definition e = (is — d5 = ES— —1 = x1, —1
d8 dS
We have the membrane strain:
ezu +—1-u2 +—1-w2 (3)
,S 2 ,S 2 ,S

which is used in Chapter 3 as Eq. (3.4) and Eq. (1.6) in Brush and Almroth (1975).
2. Rotation of normal vector

We assumed the angle between the current normal vector {rand the S (arc length)

is a, the angle between the current normal vectorii and normal vector Ill at the reference

is B.

75

Then, we have B +% 2a , and 23/} = cosa = -sinB (4)

Substituting Eq.(2.17) into Eq.(4), we obtain:

 

 

 

 

 

, d w S + K ]it w S + K 111
sm|3 = — = ' = ' (5)
,1, 1. J(1+ us — K,w)2 + (w, + K,u)2
For straight line segment (K 1 =0), according to assumption], Eq.(S) becomes
. w s
SlnB zB = (6)

JG + 11.3.)2 + w:
According to assumption 2 and W; <<1, Eq. (6) becomes: B = W’s , for the coordinates
shown in Figure 3.2, we have the normal vector rotation for straight segment

B = -W.s (7)

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

3. The current principal centerline curvature

For straight segment (K 1 =0), the centerline curvature at the current

configuration, k, will be obtained fi'om Eq. (2.21):

cd_s —dc.s (1+u's)w,ss —w_su,ss

 

 

3 [113 A: (3)
According to assumption 2( as = ass 2 O ), 1, z 1/1 + w; , the Eq.(8) becomes:
ks = _:V_§§_3 (9)
(l + w; )3

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

76

APPENDIX B

Here we try to prove the equivalence of two formulations in the pressure modeling, i.e.:

421542—- —)‘RdL‘fl pzja 5(A—w— RdL’

':1 LI ‘:1 LI

~

From the Chapter 3 kinematics, for the point on mid-surface, we have
5’ = Z0 +(s+u‘)sin6—w’ 0030
77’ = YO +(s+u‘)cos€+w‘ sin6l

Define n‘ is the unit outward normal to the meridian segmentL‘p ,then we have

 

 

 

.55 _
fit- _ as as
I_as 66 6877’

   

     

where

ai = [1+ 2)sin 0 -— 9”.ng
6s 6s 6s

6_77_ : [I +Qjcos6+flsin0
6s 6s 6s

ca: a—”-)= +a”]2 {93) ms]
696 6¢ 61s [65

We also have

 

5(Aui): a(A'Ii) 6(A4")
6e 6e ’ 6c

77

(1)

(2.a)

(2b)

(3)

(4.a)

(4.b)

(4.0)

(5)

 

where
Ar,“ 2 Au' sin6—Aw' cost9 (6.a)
An' = Au’ c056 + Aw’ sin 6 (6b)

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

if 6(Au’) _ _ 6(Aw‘) + 6wi 6Au' _ 611" 6Awi
6c 66 6s 66 6s 6c

 

 

~ _ 6(Aw')

62‘ (7)

 

Thus, we have proven Eqn. (1).

62 (Aw')

62

 

Since 2 0 , there will be no stiffness matrix due to pressure going to lefi side

C

~

of Newton-Raphson expression.

78

 

 

REFERENCES

 

79

 

REFERENCES

Asnafi N. and Skogsgardh A., 1999, Theoretical and experimental analysis of stroke-
controlled tube hydroforming, Materials Science and Engineering A279 (2000) 95-110.

Asnafi N., 1999, Analytical Modelling of Tube Hydroforming, Thin-Walled Structures
34 (1999) 295-330.

Lee D. N. and Kim Y. K., 1989, Forming Limit Diagrams for Stainless Steel clad
Aluminium Sandwich sheets, Intemal Report, Seoul National University.

Bressan J. D. and Williams J. A., Int. J. Mech. Sci. 25 (1983) 155.

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 efficient semi-implicit integration algorithm for anisotropic flow
potentials

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

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

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

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

Dohmann, F. and Hart], 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.,1989,Int. J. Mech. Sci. Vol. 31, N01, pp. 1-24.

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

80

 

Hsu, Tze-Chi and Chu, Chan-Hung, 1995, A finite 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 finite 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, 8., Oh, 8.1 and Altan, T., 1989, Metal Forming and the Finite Element
Method, Oxford University Press, New York.

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

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.,1988, On the development of a general purpose finite

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 finite 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, IT. 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, 8.1., 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., Hausserrnann 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.

81

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

82

IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

1111111111111111111311111|