An elastic-plastic finite element analysis on hydrostatic bulging of rectangular diaphragms by using layered degenerated shell elements

Int. J. Mech. Sci. Vol. 32, No. 1, pp. 49--64, 1990 Printed in Great Britain.

0020-7403/90 $3.00+.00 © 1990 Pergamon lh'-~s pie

A N ELASTIC-PLASTIC FINITE E L E M E N T ANALYSIS O N HYDROSTATIC BULGING OF RECTANGULAR DIAPHRAGMS BY U S I N G LAYERED D E G E N E R A T E D SHELL E L E M E N T S H. B. SmM and D. Y. YANG Department of Production Engineering, Korea Advanced Institute of Science and Technology, Seoul, Korea

(Received 25 April 1989; and in revised form 12 September 1989) Abstract--An incremental formulation incorporating the effect of shape change is derived for the non-steady large deformation problem of a rate-independent isotropic elastic-plastic shell. The deformation during a step is studied using the natural convected coordinate system. To account for the spread of the plastic zone in the thickness direction, a layered model is used. Change of thickness is considered by using the incompressibility condition. The corresponding finite element equations are found in order to analyse the sheet metal deformation in which the effects of shape change and bending are significant. Based on the derived equations, the FEM program is developed to analyse the hydrostatic bulging process for rectangular diaphragms. The computational results are then compared with the available experimental results.

NOTATION B L,

BNL matrix relating Lagrangian strain with nodal displacement for linear part and nonlinear part, D

E~, E ~ e~, e ~

F G~p, G~ hk .K

L N R

S,S Sa#

T,T t,i 6 Ul

W X

C~p6~ E

H

H 0• ~J

v~,~ X,x X~,x k

respectively rate of deformation tensor covariant and contravariant base vectors at time t o, respectively covariant and contravariant base vectors at time t o + At, respectively nodal force vector due to stress divergence term covariant and contravariant metric tensors at time to covariant and contravariant metric tensors at time to + At shape function in the shell mid-surface stiffness matrix spatial gradient of velocity normal vector at time t o nodal force vector due to external loading 2nd Piola-Kirchhoff stress and its rate, respectively contravariant components of deviatoric stress tensor 1st Piola-Kirchhoff stress and its rate, respectively surface traction vector and its rate, respectively column vector for Cartesian components of displacement of nodal points components of displacement in Cartesian coordinate components of displacement in convected coordinate spin tensor column vector for position vector of a material point column vector for position vector of nodal points rotations of normal vector about V] and V~ from t o to t o + At, respectively constitutive tensor related with Jaumann rate of stress tensor Lagrangian strain tensor covariant and contravariant components of Lagrangian strain tensor linear part of Lagrangian strain e~p nonlinear part of Lagrangian strain e,p shape function matrix for coordinate . shape function matrix for displacement constitutive tensor related with Truesdell rate of stress tensor natural convected coordinates Kirchhoff stress tensor and its rate Jaumann rate of Kirchhoff stress tensor Truesdell rate of Kirchhoff stress tensor unit normal vector to the shell mid-surface at nodal point k at time to and time t o + A t , respectively position vector of a material point at time t o and time t o + At, respectively position vector of nodal point k at time to and time to + At, respectively 49

50

H. B, SHIMand D. Y. YANG INTRODUCTION

Two major approaches are usually employed in analysing the large deformation of metals by the finite element method; the rigid-plastic FEM using a rigid-plastic material model in which the elastic deformation is ignored and the large strain elastic-plastic FEM. In the rigid-plastic FEM, it is not necessary to check the yield condition during the computation procedure and a larger step of deformation can be taken and thus the computing time can be greatly reduced. These advantages have made the rigid-plastic FEM more favored for analysing metal forming processes recently. Despite the numerous advantages of the rigid-plastic FEM, it suffers a severe drawback in that it fails to predict the stress history whenever elastic loading or plastic unloading are encountered. Accurate computations of the complete stress and deformation history are useful in predicting the final mechanical state of the product and for warning against possible defects in it and also it is not appropriate for problems that have significant effect of bending and shape change. In analysing the non-steady deformation of a solid numerically, continuously changing modes of deformation should be approximated by a finite number of deformations and a step-by-step procedure should be employed accordingly. A plastically deforming solid has path-dependent properties and each step size should be sufficiently small. In sheet metal deformation the displacement for each step is large even though the strain increment during a step is very small. In cases of a large displacement occurring during a step, the effect of shape change must be considered properly. Consideration of the effect of shape change implies two things; the effect of rigid body rotation of a material element in computing strains and the force equilibrium at the deformed geometry should be considered. In elastic-plastic FEM, general incremental formulations to include the effect of shape change are available. Hibbit et al. [1] introduced the complete large strain finite element formulation by using a total Lagrangian formulation in which the reference state is the original undeformed configuration. An alternative procedure of using the current state under consideration as the reference state, which is known as the updated Lagrangian formulation, was developed by McMeeking and Rice [2]. Wifi [3] analysed axi-symmetric stretching and deep-drawing of sheet metal using solid elements. Wang and Budianski [4] suggested a membrane shell formulation for sheet metal forming based on the plane stress assumption and analysed the process of axisymmetric stretching using hemispherical punch. But the formulation does not describe shell behavior properly, since the variation of stress and strain in thickness direction was neglected. Wang and Tang [5] developed a shell bending element and applied it to axisymmetric and plane strain problems, since it is based on the principal curvature. Iseki et al." [6, 7] proposed an elastic plastic FEM formulation for hydrostatic bulging of sheet metal and carried out the analysis for various shapes of diaphragms. Nakamachi et al. [8, 9] analysed stretching and deep drawing of sheets using nonlinear membrane shell theory [10], in which the effect of bending is considered using the principal curvature and the shape of die is used to find the principal curvature of diaphragm. Kim and Yang [11, 12] analysed hydrostatic bulging of sheet on the basis of convected coordinates by using the general rigid-plastic membrane model. Of all the elements developed for plate and shell structures, the degenerated shell element family has been widely considered to be the most efficient, since Ahmad et al. [ 13] developed the element. Despite the encouraging results obtained by this type of element, some weaknesses were noted. In particular the element resulted in too stiff models for application to moderately thin elements. Furthermore, the element behavior was worsened by the reduction of the thickness of elements. This is the result of what is known as shear and/or membrane locking behaviors. An improvement was obtained by reducing the order of integration by the so-called reduced or selective integration method [14]. Ofiate and Zienkiewicz [ 15] analysed axisymmetric stretching using an axisymmetric shell element and analysed stretching of rectangular diaphragm using a 9 node degenerated shell element based on Cartesian coordinates. Huang and Hinton [16] suggested elements based on the use of assumed covariant transverse shear strains referred to the natural coordinate system. The foregoing formulations are based on Cartesian coordinate system with coordinate transformation of local to global Cartesian coordinate systems or global to local coordinate

Elastic-plastic FE analysis

51

systems and much computation time is required. The possibility of locking behavior still remains to some extent. Park et al. E17] showed that the use of assumed natural coordinate strains improved transverse shear locking phenomena without using reduced integration in the shell problem. He expressed natural coordinate strains in the form of Cartesian strains by tensor transformation. Jang and Pinsky [18] showed that the shear and/or membrane locking phenomena were eliminated by using the assumed natural strains which are obtained by tensor transformation from the global Cartesian coordinates. In the present paper, an incremental formulation for finite element analysis of large deformation of an elastic-plastic shell employing the natural convected coordinate system which excludes the coordinate transformation procedure for strain components and constitutive equations is derived and associated finite element equations are found accordingly. For integration in the thickness direction and to account for the spread of the plastic zone through thickness direction, the element is divided into layers parallel t o t h e middle surface of the shell. The stress state of each layer is assumed to be a plane stress state. Elimination of the transverse shear stress by an assumption improves the locking behavior. In order to check the validity of the present formulation the computational results are compared with existing experimental results I-7] of the hydrostatic bulging process for rectangular diaphragm. THEORY Updated Lagrangian form of field equation for elastic-plastic solids A deformed metal surface is considered in three-dimensional space (Fig. 1). In analysing the non-steady state deformation by a step-by-step procedure, let us consider the deformation during one step from time t o to time t o + At. In Fig. 1, 01 and 02 are taken as surface convected coordinates and the 03-axis is taken to be the direction normal to the sheet surface. Let G~a and g~p be the metric tensors of the undeformed and deformed configurations and let G~a and g~B be their respective inverses. Base vectors in the undeformed configuration are denoted by E, and their reciprocals by E ~. Similarly, the base vectors in the deformed body are denoted by e~ and their reciprocals by e ~ dX E~ = d0~

dx e~ = ~30,~-

(1)

G~p = E~. Eo

g~p = e~'ea

(2)

E ~= G~PE~

e~= g~e~.

(3)

The displacement vector u from the undeformed configuration is (4)

u = u~E~ = u~E ~ uiEi, =

X3,x~, 03

x

03 X2 It2

FIG. 1. Coordinatesystemsin the deformationprocess.

52

H.B. SHIMand D. ¥. YANG

where u'=G'Pup. It is implied that the Greek indices are referred to the convected coordinate system and Latin indices are referred to the rectangular Cartesian coordinate system. The Lagrangian strain tensor e in the convected coordinate system is then given by [ 19] e = e~pE~E# = ~#E~E#

(5)

e,# = ½(9,#- G,#) = ½(u,,# + u#,~ + uy, uy.#),

(6)

where comma denotes covariant differentiation with respect to the undeformed metric. Decomposing the Lagrangian stress tensor into a linear part and a nonlinear part, then we have ~# = e,# + ~/,#,

(7)

where e,# = ½(u,.# + u#.,), q,# = ½u~,u~,#. The updated Lagrangian equation for large deformations is given as follows:

= f s~ t~'°+a') ru, dS - fvO ~'# ae ~#d V.

(8)

Detailed derivation for equation (8) is given in Appendix A. If the constitutive tensor is taken for the 2nd Piola-Kirchhoff stress increment and Lagrangian strain, we have

AS ~# = L~#~'~e~,

(9)

where the detailed derivation of L "ar~ is given in Appendix B. Combining equations (8) and (9), a resultant updated Lagrangian equation is obtained for the elastic-plastic solid:

= f so tl'° +a° fu,dS - fvo r'# ,~e,#d V.

(10)

As the components of the 2nd Piola-Kirchhoff stress tensor are equal to the components of the Kirchhoff stress tensor in the convected coordinate system [ 11 ], the stress integration procedure is very simple as compared with other coordinate systems. Once the approximations to the displacement increments have been calculated using equation (10), the stresses corresponding to time t o + At are calculated using equation (9)

(S,#)(to +ao = (S,#),o + AS,#

(1 1)

(C#),o +a0 = (S,#),o + ~o.

(12)

The working material can be assumed to be incompressible for a plastically deforming solid. Then the Kirchhoff stress tensor x is equal to the Cauchy stress tensor a.

Degenerated shell element formulation The element behavior is based on the assumption that lines which are originally normal to the midsurface of the shell element remain straight during the element deformations and transverse normal stress is neglected. In addition to these common assumptions for shell elements, the transverse shear stress is assumed to be absent. To take into account the change of the material property through the thickness direction, the elements are divided into several layers. In other words, a plane stress state is assumed at each layer.

Elastic-plastic FE analysis

53

X3,u~

. i FIG. 2. Nodal degrees of freedom of degenerated shell element.

•

Layer j

FIG. 3, Finite element model of layered shell structure.

The coordinates of a material point at time t o and t o + At are given, respectively, in the following form [20] 03 X = k= ~ l hkXk + - 2 k~Z akhk Fk" (13) x =

hkXk +

akhkV~,

k=l

k

(14)

1

where a k is the thickness at node k. The displacement u during At is given by (15)

u=x-X.

Substituting equations (13) and (14) into equation (15), u

hk uk +

:

k=l

akhk(v~ -- vk),

(16)

2k=l

where Vk and v.k are unit vectors normal to the shell mid-surface in direction 0 a at the kth nodal point at t o and to + At, respectively. Let fie and fig be the rotations of the normal vector about the vectors Vk and V~ from the configuration at time to to the configuration at time t o + At. Then the following approximation is made for small incremental angles ffl and fl~ v k. - V .k = -

v k llk dr.

2~t

k k Vlfl2.

(17)

54

H.B. SHIMand D. Y. YANG

Substituting equation (17) into equation (16), the incremental displacements for an internal element are obtained in terms of incremental displacements and rotations of a nodal point, 03 ,

u:

k:~ l h~uk+Tk~=~ akhk[-- Vkzflkl+ V]flk2]"

(18)

Finite element equations X~ are the orthogonal Cartesian coordinates of a material point at time to and u~ are the orthogonal Cartesian components of the displacement vector of the material point during At as shown in Fig. 4. Since the foregoing expressions are based on an arbitrary convected coordinate system and the natural coordinate is also a convected coordinate system [21], the natural coordinates (01, 02) in Fig. 4 can be taken as the convected coordinates without loss of generality. In the natural convected coordinate system, the Lagrangian strain is then given by

ox, G¢

ou,

2 \ O 0 " O0~ + ~0~ ' ~30p ~ ~30" ~30e ]"

(19)

The Cartesian coordinates and displacements of a material point are expressed in the following matrix form X=

fxlt

X 2 =HX

(20)

X3

u=

u2

(21)

=Hfi,

U3

where H is the shape function for the coordinates and H is the shape function for the displacements, .~ and ~ are nodal point coordinates and displacement vectors and are expressed as follows

2~ = {xl, x~, xL xL xL x~ . . . . . x;, xL x ~ 7 ~T = {U~, U 1, b/~, fl~, ill, U1,2 U2,2 //23, f12, f12 . . . , U~, U~, U~, fl•, fl~}T.

The linear part of the Lagrangian strain e can be expressed as follows: e ~ = ~ \ ~ 0 ~ ~0~ + ~0~ ~ j .

(22)

02

~o 3

FIG. 4. Natural convectedsystem.

Elastic-plastic FE analysis Then

{el1 I"L'0} e~21

55

(23)

~BL~ ~

BL~, BL2 and BL3 are given by

xT 0HT 0H BL~ = BL2 = ~T

'001

'001

~H T OH 002 002

(24)

B L 3 = I ( x T 0H T J H F~ T 0H r 0 H ) ~01 ~0 2 002 ~ • 2\ For the nonlinear part of the Lagrangian strain

~-

1

= t BNL' -

= BNLU,

(25)

where

~H BNL1-~01,

0H BNL2=~'

Since in previous work [15-18] the rectangular Cartesian coordinates were taken as a reference frame, the stress and strain components in the natural basis must be transformed to the Cartesian components. Since these transformations should be taken at each iteration, the procedure requires much computation time. Computation time can then be reduced greatly and the finite element coding procedure can be made simpler. Through a Newton-Raphson procedure, the matrix expression of equation (10) is reduced to the following equation K" All (0 = R - F tl -

1),

(26)

where F " - 1~is the stress divergence term at time (to + At) and (i- 1)th iteration and R is the externally applied load vector at time (t o + At).

with 0 < fl ~< 1 denoting a deceleration coefficient. APPLICATION TO HYDROSTATIC BULGING In order to check the validity of the present formulation, computations are carried out for hydrostatic bulging of a rectangular diaphragm. The computational results are compared with existing experimental results. Iseki et al. [7] made the experimental study of the hydrostatic bulging of elliptic and rectangular diaphragms. Computational procedure

For the finite element analysis 8-node serendipity elements are used which are divided into 5 membrane layers in thickness direction. In the Newton-Raphson iteration procedures, the iteration is terminated when the ratio of correction displacement vector norm to current displacement vector norm becomes less than 5 × 10 -6. As the transverse shear stress is neglected in the formulation, the stiffness matrix is singular at the first step. To avoid this singularity, it is assumed that a very small amount of stress is applied to the sheet initially. The slight change of initial configuration due to the small initial stress does not affect the

56

H.B. SHIMand D. Y. YANG

final state of the sheet appreciably. This assumption makes the problem free of shear locking. In order to check the efficiency of the present formulation, the solutions for hydrostatic bulging of square diaphragms are compared with those by conventional degenerated shell formulation based on rectangular Cartesian coordinates. Figure 5 shows the mesh system for the test and the test material is annealed aluminum whose material properties are listed in the following section. Computational results are compared only for the elastic range. To reach 1 MPa of hydrostatic pressure 100 incremental steps are used and computations are performed on SUN-3/110 micro computer with a floating-point accelerator. Table 1 shows the comparison between present formulation and conventional degenerated shell formulation. At 1 MPa, the pole height differs only by about 1% and the present, formulation needs only about 27% of the computation time that is spent for the conventional degenerated shell formulation. The change of thickness which is important in shell bending problems is considered by using volume constancy. The thickness is assumed not to change during the deformation step and the thickness at the nodal point is updated at the end of the deformation step by considering the Jacobian determinant of area at the nodal point.

Comparison with the existiny results The material constants of the test diaphragm (annealed aluminum) are as follows I-7] Young's modulus Poisson's ratio work hardening initial thickness

E = 7000 kgf/mm 2 v =0.314 # = 18(0.000769 + gp)O.29 kgf/mm 2 t o = 0.314 mm,

and the dimensions of the dies are shown in Table 2. Figure 6 shows the mesh systems used in the analysis. About 200 steps are required to reach the maximum pressure and the whole computation time for the given model in Fig. 6 is about 6000 CP seconds on a CRAY-2S computer.

TABLE 1. COMPARISON OF POLE HEIGHT AND COMPUTATION TIME BETWEEN PRESENT FORMULATIONAND CARTESIAN-BASEDFORMULATIONAT 1 M P a

Present formulation Cartesian-based formulation

Pole height (mm)

Computationtime (CP rain)

0.55172 0.55816

52.8 195.2

"//////,

"/

// /

/

I

/

i

/"

200 mm

1

FIG. 5. Finite element mesh system for test.

Elastic-plastic FE analysis TABLE 2.

57

DIMENSIONS OF DIES IN THE CALCULATION

Aspect ratio

(b/a)

Major axis (2a) Minor axis (2b) (mm) (mm) 56 65 79

1.0

0.74 0.51

56 48 40

Y b/a= 071,

b

clamped

clamped

u i : l')z=0

0

UZ= fli=O

a

X

FIG. 6. Typical finite element meshes for computation.

2.~ blz = OSl

,~

).2)

b/z:07~

1.8

b/i

z X:

= 1.0

.iP o.

eL 1,2

0.6

oeb/i=ON • : b l a = 0.51

0

0

0.2

O.Z~

1

0.6

I

0

0.8

NORMALIZED POLAR HEIGHT, hc/a

FIG. 7. Variation of polar height with respect to hydrostatic pressure for rectangular diaphragm of various aspect ratios: comparison among the present analysis and the experiment.

Figure 7 shows the relationship between pressure and polar height of the diaphragm for three aspect ratios. The present analysis is in good agreement with the existing experimental results. The proposed element does not exhibit any shear locking behavior. Figure 8 shows the strain distributions along the axis of a square diaphragm for the polar thickness strain utc = -0.1. Figures 9-12 show the strain distributions along the major and minor axes of a diaphragm for aspect ratio = 0.74. Figures 9 and 10 are for the polar thickness strain utc = -0.1. Figures 11 and 12 are for the polar thickness strain et~ = -0.2. The strain distributions are generally in good agreement with the experiment except at the boundary of the die. As the die dimension used in calculation is slightly adjusted from the actual die dimension considering the die corner radius as shown in Table 3, the distance between the actual measuring point and die edge is slightly different from that in the

58

H . B . SHIM and D. Y. YANG 0.1

l

I

I

I

b/a= tO, Etc= -0.1 o

005

0

o

~

o

- Present analysis Experiment ( 7 ]

-

-o.os

a a

Ex Ey

o

E~

~

o

o

o

o

-01

-0.15 [ 0

I

l

1

l

0.2

0 t~

0.6

08

CURRENT

10

COORDINATE ( Y / a )

FIG. 8. Strain distribution along Y axis for square diaphragm when polar thickness strain t',c = - O . I.

01 b/a = 0 7t, , Etc: - 01 005

: Present analysis Experiment [ 71

z

-

-oo5

-

tx

:

Ex

o

:

Ey

o :Et -0.1

-015

I

I

I

I

02

0~

06

08

CURRENT

FIG. 9. Strain distribution along X axis for aspect ratio •~;tc = - 0 . 1 .

0.1

10

COORDINATE ( X / a )

I

i

b a 07 o,

b/a=0.74

I

I

o

°

o

when polar thickness strain

a

°

o

0.05

0

z

-0.0S

--

Present analysis Experiment [71 •

Ex

a o

Ey Et

-0.1

-0.I~

I 02

I 04

I 0.6

] 0.8

1D

CURRENT COORDINATE ( Y / b )

FIG. 10. Strain distribution along Y axis for aspect ratio e,~ = 0.1.

b/a =0.74

when polar thickness strain

Elastic-plastic FE analysis 0.2

r

I

59 I

I

b / a = 0. % . G¢=- 0.2

~ ~ , ~ , _ ,

0.1

0.0

z

- -

Present

•

•

•

analysis

D

Experiment

~ -0.1

o

Ex E

o

Et

[7]

y

~

-02 t

I

[

I

02

0l,

06

08

-0.3

CURRENT

1.0

X - cOORDINATE ( X / a )

FIG. 1I. Strain distribution along X axis for aspect ratio b/a = 0.74 when polar thickness strain e,,¢= -0.2.

0.2

1

I

h/a=0%,

n

I

I

G~=-0.2

g

g

,

°

o

a

o

°

o

0.1

z

0.0

- -

: Present

Experiment

m -0.1

• o

Ex ~y

o

I~t

analysis

[71

-0.2

-0.3 0

I I I I (]2 0.t, Q6 0.8 CURRENT Y - COORDINATE { Y / b )

1D

FIG. 12. Strain distribution along Y axis for aspect ratio b/a=0.74 when polar thickness strain ~,, =

-

0.2.

TABLE 3. D I M E N S I O N S OF DIES IN THE EXPERIMENT

Aspect ratio (b/a) 1.0 0.74 0.51

Major axis (2a) Minor axis (2b) (mm) (mm) 56 66 80.6

56.2 49 41.4

calculation. The a s s u m p t i o n of a s h a r p edge at the die b o u n d a r y causes the e r r o r in c o m p u t a t i o n . As a whole the present analysis shows g o o d a g r e e m e n t with the experiment. F i g u r e 13 shows the relation between the p o l a r thickness strain a n d h y d r o s t a t i c pressure for two aspect ratios. As the pressure increases, the present solution is in better a g r e e m e n t with the experiment. F i g u r e s 14, 15 a n d 16 show the d i s t r i b u t i o n of local effective strain in the d i a p h r a g m s for three aspect ratios at the given Values of pressure. The a m o u n t of d e f o r m a t i o n is largest at the center of pole for b/a = 1 a n d it decreases with distance a w a y from the center. In the case

60

H.B. SHIM and D. Y. YANG 2.~

b/a = 0.51 1.8

O2~ E E

z Z

==

o.

.

1.2

o

0_

°°

"'~/~ 7/

P,esent anatysis Experiment171

-0.1

~ 0.6 I

I

I

I

I

02 0 t. 0.6 POLARTHICKNESS STRAIN, -Etc

I

0

FIG. 13. Variation of polar thickness strain with respect to hydrostatic pressure.

\

0 /ii 0

FIG, 14. Effective strain distribution for aspect ratio P = 1.4 M Pa.

b/a= 1.0 when hydrostatic pressure

FZG. 15. Effective strain distribution for aspect ratio P = 1.4 MPa.

b/a=0.74 when hydrostatic pressure

of b/a = 0.51 the amount of deformation near the center of the longer side edge is as large as the value at the sheet center as shown in [l 1]. The comparison of the present solution with the reported experiments shows that there is good agreement between the present theory and the experiment.

Elastic-plastic FE analysis

61

FIG. 16. Effective strain distribution for aspect ratio b/a=0.51 when hydrostatic pressure P = 1.4 MPa.

CONCLUSION

An incremental formulation incorporating the effect of shape change has been derived for the non-steady large deformation problem of a rate-independent isotropic elastic-plastic shell. The formulation of the strain-displacement relation is simplified by using the natural convected coordinate system. The assumption of plane stress at each membrane layer is applied to the multi-layered shell element and locking behavior is eliminated. The proposed scheme allows the thickness to change properly. The developed shell formulation has been successfully applied to the analysis of hydrostatic bulging of non-circular diaphragms, i.e. rectangular diaphragms of different aspect ratios. The proposed formulation is found to be far more efficient than the conventional degenerated shell formulation based on rectangular Cartesian coordinates. The relation between the bulging pressure and deformation has been investigated for three aspect ratios. The present solution has been shown to be in good agreement with existing experimental results. It has been thus shown that the proposed shell formulation considering large deformations is applicable to the analysis of sheet metal processes.

REFERENCES 1. H. D. HIBBIT, P. V. MARCALand J. R. RICE, Finite element formulation for problems of large strain and large displacements. Int. J. Solids Struct. 6, 1069 (1970). 2. R. M. MCMEEKING and J. R. RICE, Finite element formulations for problems of elastic-plastic deformations. Int. J. Solids Struct, 11, 601 (1975). 3. A. S. WIFI, An incremental complete solution ofthe stretch forming and deep drawing of a circular blank using a hemispherical punch. Int. J. Mech. Sci. 24, 23 (1976). 4. N. M. WANG and B. BUDIANSKI,Analysis of sheet metal stamping by a finite element method. J. appl. Mech. 45, 73 (1978). 5. N. M. WANG and S. C. TANG, Analysis of bending effects in sheet forming operations. Proc. of the NUMIFORM '86 Conference, Gothenberg, Sweden, p. 71 (1986). 6. H. ISEK1,T. JIMMA and T. MUROTA,Finite element method of analysis of the hydrostatic bulging of a sheet metal Part 1. Bull. J S M E 17, 1240 (1974). 7. H. ISEKI, T. MUROTAand T. JIMMA, Finite element method in the analysis of the hydrostatic bulging of a sheet metal Part 2. Bull. J S M E 20, 285 (1977). 8. S. TAKEZANO, E. NAKAMACHIand T. YAMAGUCHI,Elasto/viscoplastic analysis of thin circular plates under large strains and large deformations. J. appl. Mech. 47, 741 (1980). 9. E. NAKAMACHI, Finite element modelling of the punch press forming of thin elastic-plastic plates. Proc. NUMIFORM '86, p. 333 (1986). 10. B. BUDIANSKY, Notes on nonlinear shell theory. J. appl. Mech. 35, 393 (1968). 11. Y.J. KIM and D. Y. YANG, A rigid-plastic finite element formulation considering the effect of geometric change and its application to hydrostatic bulging. Int. J. Mech. Sci. 27, 453 (1985). 12. D.Y. YANG and Y. J. KIM, Analysis of hydrostatic bulging of anisotropic rectangular diaphragms by the rigidplastic finite element method. J. Engng Ind. 109, 148 (1986). 13. S. AHMAD, B. M. IRONS and O. C. ZIENKIEWICZ,Analysis of thick and thin shell structures by curved finite elements. Int. J. Num. Meth. Engn 9 2, 419 (1970). 14. D. BRAISSOULIS,Machine locking of degenerated thin shell elements. Int. J. Num. Engn 9 26, 1749 (1988). 15. E. O~ATE and O. C. ZIENKIEWICZ,A viscous shell formulation for the analysis of thin sheet metal forming. Int. J. Mech. Sci. 25, 305 (1983).

62

H.B. SHIM and D. Y. YANG

16. H.C. HUANG and E. HINTON, Elastic-plastic and geometrically nonlinear analysis of plates and shells using a new nine node element. In Finite Element Methods for Nonlinear Problems (edited by BERGAN et al.), p. 283. Springer, Berlin (1986). 17. K.C. PARK, G. M. STANLEY and H. CABINESS,A family ofC ° shell element based on generalized Hrennikoff's method and assumed natural-coordinate strains. In Finite Element MethodsJbr Nonlinear Problems (edited by BERGAN et al.), p. 265. Springer, Berlin (1986). 18. J. JANG and P. M. PINSKY, An assumed covariant strain based 9-noded shell element. Int. J. Num. Meth. Engnf.I 24, 2389 (1987). 19. A. E. GREEN and W. ZERNA, Theoretical Elasticity, 2nd edn. Oxford University Press, Oxford (1968). 20. K. J. BATHE, Finite Element Procedures in Engineering Analysis. Prentice-Hall, Englewood Cliffs, NJ (1982). 21. W.J. CHUNG. Y. J. KIM and D. Y. YANG, Rigid-plastic finite element analysis of hydrostatic bulging of elliptic diaphragms using Hill's new yield criterion. Int. J. Mech. Sci. 31, 193 (1989). 22. L MALVERN,Introduction to the Mechanics of a Continuous Medium. Prentice Hall, Englewood Cliffs, NJ (1969). 23. K. WASHIZU, Variational Methods in Elasticity and Plasticity, 3rd edn. Pergamon Press, Oxford (1982). 24. J. W. HUTCmNSON, Finite strain analysis of elastic-plastic solids and structures. In Numerical Solution oJ Nonlinear Structural Problems (edited by R. F. HARTUNG), p. 17. ASME (1973). APPENDIX A U P D A T E D LAGRANGIAN F O R M U L A T I O N IN GENERAL COORDINATES If the current configuration (t -- t 0) is taken as a reference configuration, the 1st Piola-Kirchhoffstress rate tensor T, and the 2nd Piola-Kirchhoff stress rate tensor S are given by [22] T='~-L.'t

(AI)

;S= - L-~ +.~-~- L"r,

(A2)

where L is the spatial gradient of the velocity and P=~+x.W-W.~

(A3)

f T r = f _ x . LT--L..~ '

(A4)

From equations (A2) and (A4), the 2nd Piola-Kirchhoff stress rate tensor is equal to the Truesdell rate of the Kirchhoffstress tensor when the current configuration (t = to) is taken as a reference configuration. From equations (A3) and (A4), the Jaumann rate of the Kirchhoff stress tensor is related to the Truesdell rate of the Kirchhoff stress tensor '~J= fTr +'~' D + D-.L

(A5)

In the absence of body and inertia forces, the rate form of the field equation to be satisfied in solid mechanics is as follows [23]: V.T=0

in V°

(A6)

N-T=i

on S O.

(A7)

Applying the divergence theorem, the following weak form of equations (A6) and (A7), i.e. rate form of the principle of virtual work, can be obtained.

fv ':6LTdV= I i.c3vdS. o

(AS)

,)sO

Substituting equations (A 1) and (A2) into equation (AS) and using the symmetry of S and z, the following expression is obtained

fvO('+z'LT):c~LVdV=fvO':'DdVq'-fvo

L''f:t~LTdV

f~o i' 6v dS.

(A9)

Multiplying both sides of equation (A9) by (At) 2, and approximating u = v At, e = D At, AS = SAt and At = iAt, then

fvoAS"6e~adV+fror~au,~fukadV=fs?At,fu, dS. The right-hand side of equation (A10) can be replaced by the divergence theorem

oAt,'u'dS=f f o't'+At""'dS-fsot"'dS ,,,

j vo

(A10)

Elastic-plastic FE analysis

63

where t (r°÷~) is the traction vector at time t o + A t which is measured at time t o

fv° AS'# 6e~#d V + .:v "ff~'uk~Suk~dV =fs tl'°+A"t~u~dS--I °t

(AI2)

C#Se~#dV'

JV o

APPENDIX

B

TENSOR FORM OF CONSTITUTIVE RELATIONSHIP Deviatoric stress and the J2 invariant are generalized as

s~#= C#--~g~#gT,:~

(BI)

J 2 = l,qayg#aSa#S'/a'

(B2)

In the absence of plastic deformation, the material is assumed to be hyperelastic, the Jaumann stress rate is related with the rate of deformation tensor by [24] (V#)J =

C~D~a

(B3)

C'#~a=2# g ~ g ' q ~ + ~ - O * P g ~6 -2,uF

(B4) S

with F=I

ifsY~Dra>~0

and

F=0

ifs~aD~.a<0

or

J2=(J2)max J2<(J2)max

2 2/' h''~ s=~0 ~1+~),

and

where p is Lame's constant and h' is the slope of equivalent stress-plastic strain curve. Combining equations (A5), (B3) and (B2), the rate-constitutive equation for Truesdell's rate of the Kirchhoff stress is given as (V#) rr =

S~ =

L~P~aD~

(B5)

L~aT~= L'~P = C~#76 _ ½(g~.:z#a+ g#~C6 + g~:,~#.~+ g#6C~).

APPENDIX CHANGE

(B6)

C

OF THICKNESS

AT N O D E S

Jacobian determinant of volume Jv for deformed configuration is

t3x Jv

3x2

3x3

~301 301

001

~x~

~X 3

~02

GX2

~02 ~02 =%'(% ×ez).

(c1)

cx~ dx2 ~x3 ~03 ~ 0 3 ~ 0 3 Covariant base vector for the direction of thickness direction t e I × e2

e3

2 lel × %1

,

(C2)

where t is the thickness at a node. From the volume constancy of small element at node

d V=fffJv

d01 d02 dOS = constant.

(C3)

Therefore, the Jacobian determinant of volume must be constant. From equations (C1), (C2) and (C3) the Jacobian determinant of volume is expressed as: to

t

J v = 2-IE ~ xE21 =21% x%],

(C4)

64

H.B. SHIM and D. Y. YANG

where E 1, E 2 denotes covariant base vectors for the undeformed configuration and e 1, e 2 denote covariant base vectors for the deformed configuration, t and t o denote the thickness of the node in the deformed and the undeformed configurations, respectively. Then, the thickness of node at the deformed configuration is given as follows IE1 ×E21

t = t ° lel xe21 "

(C5)