A new element for analyzing large deformation of thin Naghdi shell model. Part II: Plastic

Applied Mathematical Modelling 35 (2011) 2650–2668

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

A new element for analyzing large deformation of thin Naghdi shell model. Part II: Plastic Aazam Ghassemi ⇑, Alireza Shahidi, Mahmoud Farzin Department of Mechanical Engineering, Najafabad Branch, Islamic Azad University, Isfahan, Iran

a r t i c l e

i n f o

Article history: Received 17 August 2009 Received in revised form 12 October 2010 Accepted 15 November 2010 Available online 3 December 2010 Keywords: Large deformation Thin Cosserat shell Plastic Constrained director

a b s t r a c t In this paper a new element is developed that is based on Cosserat theory. In the finite element implementation of Cosserat theory shear locking can occur, especially for very thin shells. In the present investigation the director vector is constrained to remain perpendicular to the mid surface during deformation. It will be shown that this constraint yields accurate results in very large deformation of thin shells also the rate of convergency is very good. For plastic formulation, the model introduced by Simo is used and it has been reduced for constrained director vector and the consistent elasto-plastic tangent moduli is extracted for finite element solution. This model includes both kinematic and isotropic hardening. For numerical investigations an isoparametric nine node element is employed then by linearization of the principle of virtual work, material and geometric stiffness matrices are extracted. The validity and the accuracy of the proposed element is illustrated by the numerical examples and the results are compared with those available in the literature. Ó 2010 Elsevier Inc. All rights reserved.

1. Introduction Large deformation analysis of shells is usually studied using two different approaches: – three-dimensional theory; – direct method or Cosserat theory in which a director is assigned to a non-euclidean plane [1–4]. Similarly in numerical analysis of large deformation of shells and plates two methods have also been employed. For three-dimensional theory, a three-dimensional modified element was developed by Ahamd [5]. Other researchers like Hughes and Liu [6] and Hughes and Carnoy [7] developed this element. This element is recorded in standard finite element text books like, Bathe [8], Hughes [9] and Blytschko [10]. Another approach for this theory was developed, namely, co-rotational method. Numerical formulation of this approach was firstly presented by Wempner [11] whose article introduces co-rotational finite element in nonlinear analysis of shells. The work of Argyris [12] can also be mentioned along this approach. This approach is suitable for large deformation with small strains. Some examples of this approach are works of Parish [13], Beuchter et al. [14], Sansour et al. [15], Peng et al [16], Jiang et al. [17] and Liu et al. [18]. Another theory is direct method. This method is one of the best theories for modeling large deformation of shells. This theory was presented by Cosserat for first time and further elaborated upon by a number of authors such as Naghdi [19],

⇑ Corresponding author. E-mail addresses: [email protected] (A. Ghassemi), [email protected] (A. Shahidi), [email protected] (M. Farzin). 0307-904X/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2010.11.029

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668

2651

Antman [20]. The basic assumption of this theory is that the mid surface of the shell is regarded as an inextensible one-director surface. Typically this approach yields an exact analytical definition of the initial geometry of the shell and the representation of the stress and strain state in curvilinear coordinates and stresses are entirely in term of stress resultants and stress couples. Green and Naghdi [21] derived a general form of constitutive equation for an elastic perfectly plastic material with attention to thermo dynamical constraints. When dealing with stress resultants, it is of course important to be able to identify a yield surface which remarks the limiting values of the stress resultants. Ilyushin [22] derived an exact form of the yield surface for a linear elastic, perfectly plastic isotropic material which obeys Vonmisses yield criterion. Crisfeild [23] improved this criterion by adding a pseudo hardening effect due to progression of yielding across the shell’s thickness. Simo [24–26] extended Ilyushin criterion for an isotropic and kinematic hardening materials also he developed finite element formulation of this theory. Since the Ilyushin criterion is a multifunction and it’s surface has corner, several researchers such as Mohammed and Skallerud [27] modified Ilyushin ’s criterion to one function. In the above works, shear deformations in the direction of the thickness are taken into account. In these analyses, as the thickness approaches zero, their numerical analyses mostly experience shear locking. According to the well known Kirichhoff’s hypothesis, straight lines perpendicular to the mid-surface remain perpendicular to the deformed mid-surface. This hypothesis yields satisfactory results only when the thickness approaches zero and the deformation is not large. This hypothesis can lead to numerical difficulties, if used for large deformations. However it will be shown in this paper that by employing Cosserat’s surface and constraining the director vector to remain perpendicular to mid surface during deformation, very good results can be obtained for large deformation of thin plates and shells without any locking. This constraint is in fact a limiting analysis of the Cosserat theory in which Kirichhoff’s hypothesis is enforced and hence the shear strains in the direction of the shell’s thickness are ignored. For plasticity solution, the model extended by Simo [26], is used and it is modified for a constrained director surface also the consistent elasto-plastic tangent moduli is extracted for this modified surface. Using principle of virtual work and linearization process stiffness matrices are extracted. For numerical solution a 9 node isoparametric element has been used. The outline of this paper is as follows. In Section 2 the theory is explained. In this section the algorithm of return mapping for plastic solution is extracted and the elasto-plastic tangent moduli is derived for numerical solution. In Section 3 finite element scheme for solution a constrained Cosserat shell is developed. In this section by using virtual work, the geometric and material stiffness matrices are derived through a linearization process. In Section 4 several numerical examples are presented and the results are compared with literature. Finally conclusions are drawn in Section 5. 2. Theory In this section the stress and strain vectors have been illustrated and then elasto-plastic constitutive equations are presented. The return mapping algorithm is explained and finally elasto-plastic tangent moduli is derived for presented model. 2.1. Kinamatic relations Fig. 1 shows geometry of a three-dimensional shell with a mid surface (M). On the mid surface the convective coordinate system h1, h2 is considered which has the base vectors a1, a2 and a3 which is orthogonal to a1 and a2. The position vector of any point with respect to O is [28]:

R ¼ rðh1 ; h2 Þ þ h3 a3 :

ð1Þ

Fig. 2 shows the mid surface of an arbitrary shell in equilibrium states before and after deformation (t = 0, t respectively). In this figure x, y, z represent reference Cartesian coordinate system and h1, h2 are the convective coordinate system. a3

h/2

S

θ3 M

−h/2 R r

O Fig. 1. Geometry of a three-dimensional shell.

2652

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668 0 a =0 d 3 t a =t d 3

z

0a

ta

2

θ1

θ

2

0a

θ

2

ta

2

1

θ

1

1

E3 E2

y

E1

x Fig. 2. Equilibrium state of a quadrilateral plate at times zero and t.

The base vectors of convective coordinate system, in initial configuration are denoted by 0ai. Similarly, tai denotes base vectors of convective coordinate system at time t. It should be noted that the director vector is constrained to be perpendicular to the mid surface at each time, so ta3 = td. The Position vector of a material point, which is a function of h1 and h2, is:

2 t

t

xðh1 ; h2 Þ

3

6 7 rðh1 ; h2 Þ ¼ 4 t yðh1 ; h2 Þ 5: t

1

ð2Þ

2

zðh ; h Þ

The base vectors can be written as:

)

t

t

@ r aa ¼ t r;a ¼ @h a

t

a3 ¼ t a3 ¼ kt aa11 t aa22 k

t

t

ð3Þ

:

Components of the first and the second fundamental tensors of the surface and also components of membrane and bending strains are written as [1]: t

aab ¼ t x;a  t x;b ;

t

bab ¼ t a3;a :t ab ¼ t a3 :t aa;b ;  1 t t x;a :t x;b  0 x;a :0 x;b ; 0 eab ¼ 2 1 t t ffi ½ x;ab :t x;1  t x;2 ; jab ¼ t a3 :t x;ab ¼ pffiffiffiffi ta t 0

ð4Þ

qq ¼ t jab  0 jab :

In the above relations, lower left subscripts in the strain components denote reference configuration. Let computational strain vector, 0te, according to t 0





e ¼ t0 e11 ; t0 e22 ; 2t0 e12 ; t0 q11 ; t0 q22 ; 2t0 q12 T :

Such that

ð5Þ

t 0

e contains elastic and plastic parts and it can be decomposed into:

e ¼ t0 e_ e þ t0 e_ p :

t _ 0

In this formulation

ð6Þ

t p 0

e is the plastic part of the strain and it will be explained as follow.

2.2. Stress resultants and stress couples In Cosserat theory, membrane and bending stresses are defined in terms of stress resultants in the direction of thickness [1]. Fig. 3 shows the effective Cauchy stresses at a material point of the deformed configuration and also the effective symmetric Piola stresses corresponding to the Cauchy stresses. In the above figure, ttnab, ttmab are membrane stresses and bending moments per unit length in the deformed configuration, respectively. The invariant forms of these stresses are:

2653

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668 t 0 0

a2 0

t 0

m

n12 t 0

a1

t

n

11

21

t t

a2 t

m 21

t t

a1 t t

t 0

m 22

t t

n12

n11

m 22

Fig. 3. Effective Cauchy and Piola stresses at a material point of the body.

9 t t ab t aa  t ab ; > t nn ¼ t n = t t at : t q ¼ t q aa ; > ; ab t t t t aa  ab : tm ¼ tm

ð7Þ

Similarly 0tnab, 0tmab are membrane stresses and bending moments per unit length in the reference configuration, respectively. The invariant forms of these stresses are: t t ab 0 aa  0 ab ; 0n ¼ 0n t t a b0 aa  0 ab : 0m ¼ 0m

) ð8Þ

:

These stresses are related to each other according to the following relations:

)

abab ¼ Jtt nab : : abab ¼ Jtt mab ;

ð9Þ

t

where J ¼ dd0SS is the transformation Jacobean, namely, the ratio between the element surface after and before deformation. For a convective coordinate the Jacobean term can be written as:

J ¼ detðFÞ where Fji ¼ t ai  0 aj :

ð10Þ

In this formulation F is the deformation gradient tensor. For simplicity computational stress vector, 0tr, has been defined according to: t 0

r¼

D

11 t 22 t 12 t 11 t 22 t 12 t 0n ; 0n ; 0n ; 0m ; 0m ; 0m

E

ð11Þ

:

2.3. Constitutive equations In the previous sections the stress and strain vectors were determined. Because the stress components are in terms of stress resultants and stress couples, the constitutive equations must be formulated directly according them. The generalized Ilyushin–Shapiro elasto-plastic model is entirely in terms of stress resultants and stress couples. It should be noted that this yield surface is multifunction. 2.3.1. Elasto-plastic constitutive equations For an elasto-plastic material which the elastic region is linear we can write [24]: t 0

where





r_ ¼ C t0 e_  t0 e_ p ;

ð12Þ

t p 0

e is the plastic part of the strain and: 

C¼

Cm

0

0

Cb

 ð13Þ

; 2

Eh 6 6 Cm ¼ ð1  m2 Þ 4

ð0 a11 Þ2 sym

mð0 a11 Þð0 a22 Þ þ ð1  mÞð0 a12 Þ2

ð0 a11 Þð0 a12 Þ

ð0 a11 Þ2

ð0 a22 Þð0 a12 Þ 1þm 0 12 2 ð a Þ 2

þ 12 m ð0 a11 Þð0 a22 Þ

3 7 7; 5

2

Cb ¼

h Cm : 12

To define the plastic part of the strain, the yield function must be defined. In general for an elasto-plastic material that its yield surface is one function, the yield function can be defined as follow:

2654

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668

f = f(0tn, 0tm, Ps, U0) where Ps, s = 1, 2, . . . n characterizing the hardening response and U0 is the character of initial configuration. If the associative flow rule is considered then: 0 t _p s _s _ @f _ s t t 0 e ¼ c @ t r and a general hardening law as P ¼ ch ð0 n; 0 m; p ; U Þ. 0

In the above equations c_ P 0 is the plastic consistency parameter; a function satisfying Kuhn–Tucker complementary conditions:

c_ P 0; f 6 0; c_ f ¼ c_ f_ ¼ 0:

ð14Þ

The above condition must be satisfied through plasticity model. According to Simo [26] the set of Ps, s = 1, 2, . . . n is supplemented by a conjugate set of internal variables as, s = 1, 2, . . . n through the transformation: P = rH(a) = Da, where H(a) is hardening potential and for simplicity, it will be assumed that it is strictly quadratic (HðaÞ ¼ 12 aT DaÞ such that D 2 Rn  Rn is constant. According to Simo [26], the hardening relation can be written as:

  Da k0 where D ¼ 0 ¼  P Da k

  p

P¼

and D ¼

2 0 H I6 ðI6 is the unique matrixÞ: 3

ð15Þ

The constants k0, k0 , H0 are yield parameters. k0 is the uniaxial yield stress, k0 is the linear isotropic hardening moduli and H0 is  2 R6 are associative with isotropic and kinematic hardening of the the kinematic hardening moduli. Variables a 2 R and a yield surface, respectively. For a yield surface with multiple functions, the elastic domain can be defined as follow:

n

o

Xr ¼ ðr; PÞ 2 R6  R6 jfl ðr; PÞ < 0

for all l 2 f1; 2; . . . ; mg; n o 6 P @ Xr ¼ ðr; PÞ 2 R  R jfl ðr; PÞ ¼ 0 for all l 2 f1; 2; . . . ; mg;

ð16Þ ð17Þ

where oXr is the boundary of the yield surface. The functions fl(r, P) are smooth functions which are assumed to define independent constraints at any (r, P) 2 oXr and may intersect in a nonsmooth fashion. The closure Xr [ oXr is assumed to be a closed convex set. If the associated flow rule is used, according to Koiter rule, the plastic components of strains and hardening characters can be written as below:

e_ p ¼

q X

c_ l

l¼1

a_ ¼

q X

@fl ; @r

ð18Þ

c_ l @ P fl ðr; P; pÞ;

ð19Þ

l¼1

where q = {l 2 mjfl = 0} (active surfaces). 2.3.2. Generalized Ilyushin–Shapiro elastoplastic model Let the back stress ‘‘P’’, then the yield function is determined as:

fl ðr þ P; pÞ ¼ ul ðr þ PÞ 

j2 ðpÞ 6 0; l 2 f1; 2g: j20

ð20Þ

In the above formulation:

ul ðr þ PÞ :¼ ðr þ PÞT Al ðr þ PÞ;

ð21Þ

Such that:

2 Al ¼ 4

1 n20

signðlÞ pffiffi 2 3n0 m0

where signðlÞ :¼

signðlÞ pffiffi 2 3n0 m0

p



1 m20

p

þ1; 1;

if if

p

p

3 5;

ð22Þ 2

3

1 l ¼ 1 and p :¼ 4 11  2 0 5. 2 1 0 l¼2

0

0

3

In above relations m0 and n0 are the yield parameters associated with membrane and bending response respectively. These yield parameters are typically related to the uniaxial yield parameter j0 through the relations: 2 n0 = hj0 and m0 ¼ h4 k0 where h is the shell thickness. Also

jðpÞ ¼ j0 þ j0 p; which defines the radius of the yield surface.

ð23Þ

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668

2655

If the associated flow rule is used, according to Koiter rule, the plastic components of strains can be written as below:

e_ p ¼ c_

X @fl : @r l2f1;2g

ð24Þ

So from Eq. (20) we have:

X

e_ p ¼

c_ l 2Al ðr þ PÞ:

ð25Þ

l2f1;2g

The parameters p and P according to the generalized Ilyushin–Shapiro yield function are determined as below:

a_ ¼

m X

c_ l @ p fl ðr; P; pÞ ! a_ ¼

l¼1

a_ ¼

m X

X

c_ l

2j0 jðpÞ

l2f1;2g

c_ l @ p fl ðr; P; pÞ ! a_ ¼

X

j20

;

ð26Þ

c_ l 2Al ðr þ PÞ:

ð27Þ

j0 2 a and P ¼ D0 a ¼  H0 I6 a ðI6 is the unique matrixÞ: 3 j0

ð28Þ

l¼1

l2f1;2g

Then we have:

p ¼ Da ¼ 

For solving the illustrated yield function and finding stresses and plastic strains an iteration step must be done. This iteration process is discussed in return mapping algorithm. 2.3.3. Return mapping algorithm The algorithm has the standard geometric interpretation of a closest-point-projection, in the energy norm, of a trial state onto the elastic domain. Because the illustrated surface has corner, the active surface must be defined during return to the yield surface. 2.3.3.1. Discerete algorithmic problem. Let current mid-surface of the shell is known. Consider a time disceretization of the  n are known. Let DUn+1, Dtn+1 be a given increment interval [0, T]  R of interest. We assume that the variables en ; epn ; an ; a  n must be updated to enþ1 ; epnþ1 ; in the displacement and time on the interval t 2 [tn, tn+1]. Then the variables en ; epn ; an ; a anþ1 ; a nþ1 at tn+1 2 [tn, T]. To this end, application of an implicit, backward Euler difference scheme leads to the following nonlinear coupled system:

epnþ1 ¼ epn þ

q X l

q X l

l¼1

l¼1

a nþ1 ¼ a n þ

cnþ1 @ r fl ðr; P; pÞnþ1 ¼ epn þ

cnþ1 2Al ðr þ PÞnþ1 ;

q X l

cnþ1 @ P fl ðr; P; pÞnþ1 ;

l¼1

rnþ1 ¼ Cðenþ1  epnþ1 Þ; Pnþ1

2  nþ1 ¼  H0 a  nþ1 ; ¼ D a 3

ð29Þ

0

anþ1 ¼ an þ

q X l

cnþ1 @ p fl ðr; P; pÞnþ1 ;

l¼1

p ¼ Danþ1 ¼  l

j0 a ; j0 nþ1 l

where we have set cnþ1 ¼ Dtc_ nþ1 . Because the use of convective coordinates, it does not need to objective rates. It is considerable that the discrete counterpart of the Kuhn–Tucker loading/unloading takes the following form:

clnþ1 P 0; f l ðr; P; pÞnþ1 6 0 and clnþ1 fl ðr; P; pÞnþ1 ¼ 0 ðno sum on lÞ l 2 f1; 2g:

ð30Þ

Convexity of the yield surface is guaranteed by a positive definite Al(l 2 {1, 2}) and consequently the solution is unique. The trial state in the interval t 2 [tn, tn+1] can be determined as below: At the first stage set:

2656

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668 trial

 epnþ1 :¼ epn ; a trial atrial nþ1 :¼ an ; nþ1 ¼ an ; trial

eenþ1 :¼ enþ1  epn ; 



p rtrial nþ1 :¼ C enþ1  en ;

ð31Þ

2 0 trial j0 trial Ptrial ptrial nþ1 :¼  H anþ1 ; nþ1 ¼  0 anþ1 ; 3 j  trial trial  trial fl;nþ1 :¼ fl rnþ1 ; Pnþ1 ; p : As the notion of Simo [26], the yield function, fl ðr; P; pÞnþ1 , can be expressed in terms of the consistency parameters c1nþ1 ; c2nþ1 . By this, the return mapping reduces to solution of the following nonlinear system: fl ðrnþ1 þ Pnþ1 ; pnþ1 Þ ¼ f l ðc1nþ1 ; c2nþ1 Þ ¼ 0l 2 f1; 2g as follow: At the first stage it is obvious that:

X

p rtrial nþ1 ¼ rnþ1 þ CDenþ1 ¼ rnþ1 þ C

clnþ1 2Al ðr þ PÞnþ1 ;

ð32Þ

l2f1;2g

2 0 2 0 X l  Ptrial c 2Al ðr þ PÞnþ1 : nþ1 ¼ Pnþ1 þ H Danþ1 ¼ Pnþ1 þ H 3 3 l2f1;2g nþ1

ð33Þ

Define parameter g as:

gnþ1 ¼ rnþ1 þ Pnþ1 ¼ where rnþ1 ¼



n m





gn gm

ð34Þ

; nþ1

and Pnþ1 ¼



nþ1

Pn Pm



trial trial and gtrial nþ1 ¼ rnþ1 þ Pnþ1 .

nþ1

By this definition, from Eqs. (32)–(34) we have:

0 trial nþ1

g

¼ @I6 þ

X 4 clnþ1

3

1 H Al þ 2CAl Ag 0

ð35Þ

nþ1 :

From Eqs. (13) and (22)1:

2 CAl ¼ 4

1 n20

signðlÞ pffiffi 2 3n0 m0

Cn p

signðlÞ pffiffi 2 3n0 m0

1 m20

Cm p

Cn p

Cm p

3 5:

ð36Þ

For simplicity consider orthogonal and diagonal matrix Q as:

2

1

1 0

3

6 7 Q ¼ 4 1 1 0 5: pffiffiffi 0 0 2 Matrices p and C can be rewritten as:

p ¼ Q Kp Q T

and C ¼ Q KC Q T :

ð37Þ

where Kp and KC can be computed as below:

23

3 2 E 0 0 1þm 2 6 7 6 Kp :¼ 4 0 12 0 5 and KC :¼ 4 0 0 0 0 3

0 E 1m

0

0

3

7 0 5:

ð38Þ

E 2ð1þmÞ

Such that p and C are commute; i.e., pC ¼ Cp Let nn :¼ QTgn and nm :¼ QTgm, so from Eq. (35) we have:



C1

C2

C3

C4



nn nm

"

 ¼ nþ1

ntrial n ntrial m

# ; nþ1

ð39Þ

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668

where:

  2 2 0 ¼ I3 þ cpnþ1 w1 ; H K þ hK K p p C n20 3   2 2 0 H Kp þ hKC Kp ¼ cmnþ1 w2 ; :¼ ðc1  c2 Þnþ1 pffiffiffi 3n0 m0 3 " # 3 2 2 0 h 1 2 H Kp þ :¼ ðc  c Þnþ1 pffiffiffi K Kp ¼ cmnþ1 w3 ; 12 C 3n0 m0 3 " # 3 2 2 0 h 1 2 H Kp þ :¼ I3 þ ðc þ c Þnþ1 2 K Kp ¼ I3 þ cpnþ1 w4 ; 12 C m0 3

2657

C1 ðc1 ; c2 Þnþ1 :¼ I3 þ ðc1 þ c2 Þnþ1 C2 ðc1 ; c2 Þnþ1 C3 ðc1 ; c2 Þnþ1 C4 ðc1 ; c2 Þnþ1

ð40Þ

where cpnþ1 ¼ ðc1 þ c2 Þnþ1 and cmnþ1 ¼ ðc1  c2 Þnþ1 From Eq. (39) we have:



nn nm



 ¼

nþ1

C1

C2

C3

C4

1 "

ntrial n ntrial m

# ð41Þ

: nþ1

With these results, first part of Eq. (20) can be rewritten as:



ul;nþ1 ¼



nn

 Xl

nm

nm

nþ1

2



nn

where Xl :¼ 4

nþ1

1 n20

signðlÞ pffiffi 2 3n0 m 0

Kp

signðlÞ pffiffi 2 3n0 m 0

1 m20

Kp

Kp

Kp

3 5:

ð42Þ

By replacing Eq. (41) into (42) we have:

"

ul;nþ1 :¼

ntrial n

#T 

ntrial m

nþ1

By defining Hðc1 ; c2 Þ :¼

"

ul;nþ1 :¼

ntrial n



#T

C1

C2

C3

C4

C2 C4

C1 C3

"

H Xl H nþ1

 Xl

1

T

ntrial m

T

C1

C2

C3

C4

1 "

ntrial n

#

ntrial m

ð43Þ

: nþ1

we have:

ntrial n ntrial m

# ð44Þ

; nþ1

where the above equation is only function of c1nþ1 and c2nþ1 . Similarly, the second term of the yield function can be written in terms of c1nþ1 and c2nþ1 that is discussed as follow. According to the given yield function we have:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

^ l ðc1 ; c2 Þnþ1 jðpÞnþ1 ¼ j0 ul ðr þ PÞnþ1 so j^ ðc1 ; c2 Þnþ1 ¼ j0 u

ð45Þ

And from Eq. (29)5 we have:

anþ1 ¼ an þ

X

cl

 2j0 jðpÞ

j20

l2f1;2g

so pnþ1 ¼ pn þ

X

clnþ1

l2f1;2g

2jðpÞ

j0

:

ð46Þ

So from Eq. (45) we can write:

X

^ðc1 ; c2 Þnþ1 ¼ pn þ p

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

^ ðc1 ; c2 Þnþ1 ¼ j0 þ j0 p ^ l ðc1 ; c2 Þnþ1 and j ^ðc1 ; c2 Þnþ1 : clnþ1 2 u

ð47Þ

l2f1;2g

By Eqs. (44) and (47) the yield function can be written completely as a function of c1nþ1 and c2nþ1 :

j^ ðc1 ; c2 Þnþ1 1 2 ^f l ðc1 ; c2 Þ ^ ¼ 0; nþ1 ¼ ul ðc ; c Þnþ1  2 j0

l ¼ 1; 2:

ð48Þ

Because ^f l;nþ1 monotonically decrease with clnþ1 , for increasing hardening laws, Eq. (48) has a unique solution clnþ1 P 0; clnþ1 P 0. For determination of c1nþ1 and c2nþ1 the Eq. (48) must be solved. These equations are nonlinear in terms of c1nþ1 and c2nþ1 . @^f For solving Eq. (48) and for finding c1nþ1 and c2nþ1 Newton Raphson algorithm is used. So the term @ clb b 2 f1; 2g must be computed. By differentiating of Eq. (46) with respect to cb we have:

@ ul @ cb

" ¼2

ntrial n ntrial m

where Hðc1 ; c2 Þ :¼



#T HT Xl nþ1

C1 C3

C2 C4

1

" # @H ntrial n ; @ cb ntrial m nþ1 . Computation of

ð49Þ @H @ cb

has been discussed in Appendix A.

2658

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668 ^ c1 ;c2 Þ @ kð

nþ1 Also for Newton Raphson iteration to find c1nþ1 and c2nþ1 ; must be determined that is computed at follow. @ cl If active surface is single (l = 1 or l = 2) rename k = {ljfl,n+1 > 0} then

1

2

k nþ1

^ðc ; c Þnþ1 ¼ pn þ 2c p

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ^ k ðc1 ; c2 Þnþ1 and u

0 1 ^ k;nþ1 @u qffiffiffiffiffiffiffiffiffiffiffiffiffi ^ @j @ ck 0@ ^ k;nþ1 A: ¼ j 2 pffiffiffiffiffiffiffiffiffiffiffiffiffi þ 2 u @ ck ^ k;nþ1 u

ð50Þ

And if two surfaces are active, then

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

X l

^ðc1 ; c2 Þnþ1 ¼ pn þ p

^ l ðc1 ; c2 Þnþ1 : cnþ1 2 u

ð51Þ

l¼1;2

So

0 1 ^ l;nþ1 @u qffiffiffiffiffiffiffiffiffiffiffiffiffi X l ^ @j @ cb 0@ ^ b;nþ1 A: ¼j cnþ1 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ 2 u @ cb ^ l;nþ1 u l¼1;2

ð52Þ

By this for one active surface:

^ l 2j ^ @j ^ @^f l @ u ¼  2 @ ck @ ck j0 @ ck (

Dck ¼ fk;k Dcv ¼ 0

@fk;k @ ck

l ¼ k;

 ;

ð53Þ

ð54Þ

v – k;

And for two active surfaces we have:

^ l 2j ^ @j ^ @^f l @ u ¼  2 ; @ cb @ cb j0 @ cb (

Dc1 Dc2

)

 ¼

f1k f2k



2

ð55Þ

@f1

1

@f1 @ c2

@f2 @ c1

@f2 @ c2

4 @c

31 5 :

ð56Þ

Now, the algorithm of return mapping can be written as below: p trial trial trial 1. Compute rtrial nþ1 ¼ Cðenþ1  en Þ; pnþ1 ¼ pn ; Pnþ1 ¼ Pnþ1 and fl;nþ1 for l = 1, 2 from Eq. (20). trial 2. If fltrial ;nþ1 6 0 (for l = 1, 2) then we are in elastic phase and set ð. . . Þnþ1 ¼ ð   Þnþ1 and Exit else go to step 3.

3. This step is start of Newton Raphson iteration for finding c1nþ1 and c2nþ1 . At the first set k = 0 and k = {ljfl,n+1 > 0} 4. cknþ1;k ¼ 0 and Dcknþ1 ¼ 0. @f 5. Compute fk,n+1 from Eq. (20) and @k;nþ1 ca from Eqs. (53) or (55), according its condition.  1 Dc 6. Compute from Eqs. (54) or (56),according its condition. Dc2 l l nþ1 ¼ cnþ1;k þ Dcl where l = 1,2. 7. Let c lnþ1 < 0 for l = 1, 2 then reset k ¼ fljc lnþ1 > 0g and go to step 4 else clnþ1;k ¼ c lnþ1 and set k = k + 1. 8. If c 1 1 9. Check convergency, if (Dc + Dc ) 6 tolerance exit, else go to step 5. By this the algorithm of the return mapping is completed and the parameters c1nþ1 and c2nþ1 are determined. 2.3.4. Elastoplastic tangent moduli nþ1 For linearizing the weak form of equilibrium equations the expression ddrenþ1 or elasto-plastic tangent moduli, is needed. This moduli is computed for an isotropic and kinematic hardening rule and is given at follow. The process of extracting elasto-plastic tangent moduli is discussed in Appendix B.

X X @fa;nþ1 drnþ1 ¼ Br  Br z1 y Br ; denþ1 @ r ab b b2k a2k where k = {ljfl,n+1 = 0} 1 B1 r ¼C þ

X a2k

ð57Þ

l = 1, 2. And

canþ1

@ 2 fa;nþ1 : @ r2

ð58Þ

2659

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668

Also: T

yk ¼

T

T

@fk @fk 0 X a @ 2 fa;nþ1  D cnþ1 @r @r @ r2 a2k

1 where B1 þ canþ1 p ¼ D

and zak ¼ yk Be

@fa;nþ1 @fk 0 @fa;nþ1 @fk @fa;nþ1 þ D þ Bp ; @r @r @r @p @p

ð59Þ

@ 2 fa;nþ1 . @p2

3. Finite element implementation In this section the numerical solution is discussed. For the numerical solution the principle of virtual work is used to obtain the weak form of the governing differential equations and material and geometric stiffness matrices are derived through a linearization process. 3.1. Variatonal form of the virtual work After defining stress and strain components, by using the principle of virtual work, at time t we have:

Z

tS

ðtt nab dt0 eab þ tt jab dabjab Þdt S ¼ t Rext :

ð60Þ

Or it can be written as:

Z

0S

ðt0 nabdt0 nab þ t0 mab dt0 jab Þd0 S ¼ t Rext ;

ð61Þ

where Rext is virtual work of the external forces and can be written in terms of boundary tractions according to the following relation:

Z

 :dt U þ t m   dt dÞdt S þ ðt n

tS

Z

 :dt U þ t m:d  t dÞd@ t S ¼ t Rext : ðt n

ð62Þ

@t S

 ; tm  are distributed force and moment vectors at time t per unit area of the deformed surface, respecIn the above relation, t n t t  tively, and n; m are distributed moment and force vectors at time t per unit length of boundary, @ t S, respectively. Computational stress and strain vectors, 0tr and 0te, are defined as: t 0 t 0





e ¼ t0 e11 ; t0 e22 ; 2t0 e 12 ; t0 q11 ; t0 q22 ; 2t0 q12 T ; r¼

D

11 t 22 t 12 t 11 t 22 t 12 t 0n ; 0n ; 0n ; 0m ; 0m ; 0m

E

ð63Þ ð64Þ

:

So the formulation (61) can be written as:

Z

0S

t 0

r T dt0 e d0 S ¼ t Rext :

ð65Þ

At follow, the variational form of the strain vector is computed. An arbitrary element with liner boundaries in Cartesian coordinates can be mapped into a standard isoparametric 9 nodes element. For simplicity natural coordinates f, g are considered to be convective coordinate by the simple boundary equations of: ha = ±1 (see Fig. 4). 0

x ¼ N 1 0 x1 þ N 2 0 x2 þ    þ N 9 0 x9 ;

ð66Þ

0

y ¼ N1 0 y1 þ N2 0 y2 þ    þ N9 0 y9 ;

ð67Þ

where Ni is ith shape function of the isoparametric 9 nodes element.

η 4

y0

7

3 6

ξ

9

8

2 1

5

x0 Fig. 4. Isoparametric nine node element.

2660

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668

So the base vectors in reference configuration are:

2 0

N1;1 0 x1 þ N2;1 0 x2 þ    þ N9;1 0 x9

3

6 7 a1 ¼ 4 N1;1 0 y1 þ N2;1 0 y2 þ    þ N9;1 0 y9 5 and

2 0

N1;2 0 x1 þ N2;2 0 x2 þ    þ N9;2 0 x9

3

6 7 a2 ¼ 4 N1;2 0 y1 þ N2;2 0 y2 þ    þ N 9;2 0 y9 5:

0

ð68Þ

0

Thus the first fundamental form of the surface in reference configuration is:

0

"



aab ¼

0

a1  0 a1

0

a1  0 a2

0

a2  0 a1

0

a2  0 a2

# ð69Þ

:

Let U be the displacement vector, then the position of any point can be written as: t

rðh1 ; h2 Þ ¼ 0 rðh1 ; h2 Þ þ t Uðh1 ; h2 Þ;

ð70Þ

2

t

3 u 6 7 U¼4v 5 w

ð71Þ

In-plane displacements are interpolated as follow:

u ¼ N1 u1 þ N2 u2 þ N3 u3 þ    þ N9 u9 ; v ¼ N1 v1 þ N2 v2 þ N3 v3 þ    þ N9 v9 ;

ð72Þ

where Ni is ith shape function of the isoparametric 9 nodes element. For out of plane displacements the Hermitian shape functions are employed for the 4 corners as follow:

w¼

4 X

N1i wi þ N2i

i¼1

! @w @w @2w ; þ N3i þ N4i @ gi @fi @fi @ gi

ð73Þ

where N 1i to N 4i are the Hermitian shape functions. For example,

N11 ¼ H01ðfÞH01ðgÞ; N21 ¼ H01ðfÞH11ðgÞ H01ðfÞ ¼ 1=4ð2  3f þ f3 Þ;

ð74Þ

N31 ¼ H11ðfÞH01ðgÞ H11ðfÞ ¼ 1=4ð1  f  f2 þ f3 Þ; N41 ¼ H11ðfÞH11ðgÞ: Let’s assume t

b U ¼ N U:

ð75Þ

In the above relation, N is the shape function matrix and can be written as follow:

2

3 Nu 0 0 6 7 N ¼ 4 0 Nv 0 5; 0 0 Nw

ð76Þ

where

Nu ¼ ½N1

N2

N3

N4

N5

N6

N7

N9 ;

Nv ¼ ½N1 N2 N3 N4 N 5 N6 N7 N9 ; h 01 11 01 11 Nw ¼ H01 H01 H11 H11 ... f1 H g1 f1 H g1 f1 H g1 f1 H g1

ð77Þ ...

02 H01 f4 Hg4

12 H01 f4 H g 4

02 H11 f4 H g 4

i 12 H11 f4 H g 4 :

ð78Þ ð79Þ

By this definition from Eq. (4) we have:

b ¼ Eab d U; b dt eab ¼ ðrT;a N;b þ rT;b N;a Þd U  qab 1 b ¼ Kab d U; b dt jab ¼ pffiffiffiffiffi Q ab  pffiffiffiffiffiffiffi A d U ta 2 t a3

ð80Þ ð81Þ

where:

Q ab ¼ t

t

T r;1  t r;2 N;ab þ ðt r;2  t r;ab ÞT N;1 þ ðt r;ab  t r;1 ÞT N;2 t

t

qab ¼ r;ab : r;1  r;2

ð82Þ ð83Þ

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668

2661

The symbol ‘‘T’’ is the transpose symbol. Also from Eq. (3)2 we have:

dt d ¼ dt a3 ¼



1 1 b ¼ Yd U: b pffiffiffiffiffi bT12  pffiffiffiffiffiffiffi AT ðr;1  r;2 ÞT d U ta 2 t a3

ð84Þ

Such that

2

t

0

6t 1 1 t b12 ¼ 6 4 z;1 N;2  z;2 N;1 t

y;2 N1;1  t y;1 N1;2

t

z;2 N2;1  t z;1 N2;2

t

0

t

x;1 N2;2  t x;2 N2;1

y;1 N3;2  t y;2 N3;1

3

7 x;2 N3;1  t x;1 N3;2 7 5

ð85Þ

0

And

A ¼ 2t a11 t rT;2 :N;2 þ 2t a22 t r T;1 :N;1  2t a12 t r T;2 :N;1  2t a12 t r T;1 :N;2 ;

ð86Þ

where ‘‘a’’ is determinant of aab. So by substituting (80), (81) and (84) in Eq. (65) we have:

Z 0S

t ab 0 n Ea b

Z  ab þ t0 m Kab d0 S ¼

t

 T N þ tm  T Y dt S þ n

@t S

Z

t  T N þ t mYÞd@  ðt n S;

ð87Þ

@t S

where N, Eab, Kab and Y are determined through Eqs. (80), (81) and (84) respectively. 3.2. Local cartesian system For simplicity, the curvilinear convective coordinate system can be mapped to a local Cartesian system [25]. Let’s define a local Cartesian system {xa, x3} with base vectors {t n1, t n2, t n3} by means of the orthogonal transformation:

Kt ¼ ðE3  t nÞI3 þ ½E3  t n þ tr

1 ðE3  t nÞ  ðE3  t nÞ; 1 þ E3  t n

ð88Þ

t r

where t n ¼ t a3 ¼ kt r;1;1 t r;2;2 k is the normal to the mid surface and E1, E2 and E3 are the base vectors of reference Cartesian coordinate:

2 3 2 3 2 3 1 0 0 6 7 6 7 6 7 E1 ¼ 4 0 5 E2 ¼ 4 1 5 E3 ¼ 4 0 5: 0 0 1

ð89Þ

Also [E3  tn] is skew-symmetric tensor corresponding to E3  tn vector. Observe that Kt maps E3 ? tn = KtE3 without drilling and tr,a.t n = 0 and tna = KtEa such that tna.tnb = dab Also it can be seen that at time t:

@xa t t ¼ n a : r;a @ha

and

t

r;a ¼

@xa t na ; @ha

ð90Þ

where ha is a curvilinear system. So in the local Cartesian system 0aab = dab 3.3. Geometry and material stiffness matrices In this section stiffness matrices are extracted by linearization the virtual work, Eq. (87). It is obvious that Eq. (87) at time b and should be linearized for the numerical analyses. t, is nonlinear in term of U For linearization, the Newton Raphson method is employed as follow: b k is known where ‘‘k’’ is iteration number, then It is assumed that U

b kþ1 ¼ U b kþ1 : b k þ DU U

ð91Þ

Rename:

Z 0S

t ab 0 n Eab

ab

þ t0 m Kab d0 S 

Z tS

t T tm N

 T Ddt S  þn

Z @t S

t T tn N

b  T Dd@ t S ¼ fð UÞ: þ tt m

ð92Þ

Then

b kþ1 Þ ¼ fð U b kÞ þ fð U

 b  @fð UÞ b ¼0)  DU b bk @U U

b ¼ DU

 !1 b  @fð UÞ b k ÞÞ: ðfð U  b bk @U U

ð93Þ

2662

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668

The stiffness matrix

  @fðb U Þ  k is computed as follows: U b @b U



 b  @fð UÞ  ¼ KM þ KG ; b bk @U

ð94Þ

U

where KM can be computed as below:

Z

KM ¼

@ t0 r 0 Bd S ; b @U

0S

ð95Þ

where B is defined as:

3 E11 6E 7 6 22 7 7 6 6 E12 7 7 B¼6 6 K 7: 6 11 7 7 6 4 K22 5 2

ð96Þ

K12 R R For elastic deformation, KM ¼ 0 S BT CBd0 S and for plastic deformation, KM ¼ 0 S BT Cep Bd0 S, where Cep is determined from Eq. (57). R And KG ¼ 0 S t0 r @B d0 S, this term is computed as below: U @b

Z

KG ¼

ab

0S

E ab ¼

ab

ðt0 n E ab þ t0 m Kab Þd0 S;

ð97Þ

 1 T N;a N;b þ NT;b N;a ; 2 0 1 0

ð98Þ 1

t

0

1

0

1

t

1 1 1 B 3 qab C B qab C B B C C Kab ¼ @ qffiffiffiffiffiffiffiffiffiffiAAT :A  @ qffiffiffiffiffiffiffiffiffiffiAA  @ qffiffiffiffiffiffiffiffiffiffiAAT :Q ab  @ qffiffiffiffiffiffiffiffiffiffiAQ Tab :A þ pffiffiffiffiffiffiffiffi Qab ; t aÞ 5 3 3 3 ð t t t t 2 ð aÞ 2 ð aÞ 4 ð aÞ 2 ð aÞ

ð99Þ

A ¼ 4NT;2 t rT;2 t r;1 N;1 þ 4NT;1 t rT;1 t r;2 N;2  2NT;1 t rT;2 t r;1 N;2  2NT;1 t rT;2 t r;2 N;1  2NT;2 t rT;1 t r;1 N;2  2NT;2 t rT;1 t r;2 N;1 ;

ð100Þ

Qab ¼ BT12 N;ab þ BT2ab N;1 þ BTab1 N;2 ;

ð101Þ

2

Bab1

0 6t 1 t ¼6 z N  z;1 N1;ab 4 ;ab ;1 t y;1 N1;ab  t y;ab N1;1 2

0

6t 1 1 t B2ab ¼ 6 4 z;2 N;ab  z;ab N;2 1 1 t t y;ab N;2  y;2 N;ab

t

z;1 N2;ab  t z;ab N2;1

t

0

t

t

x;ab N2;1

t

z;ab N2;2  t z;2 N2;ab

t

0

t

t

x;2 N2;ab





t

t

x;1 N2;ab

x;ab N2;2

y;ab N3;1  t y;1 N3;ab

3

7 x;1 N3;ab  t x;ab N3;1 7 5;

ð102Þ

0 y;2 N3;ab  t y;ab N3;2

3

7 x;ab N3;2  t x;2 N3;ab 7 5:

ð103Þ

0

Rename:

Fk ¼

Z 0S

  ðnab Eab þ mab Kab Þd0 S  bk

ð104Þ

U

and

Fext ¼

Z tS

 T N þ tt m  T DÞdt S þ ðtt n

  T N þ tm  T DÞd@ t S Þ : ðtt n t bk @t S U

Z

ð105Þ

So from Eq. (93):

b ¼ ðKM þ KG Þ1 ðFext  Fk Þ: DU Therefore the algorithm of solution can be summarized as below: 1. Consider interpolation matrix, N, from Eq. (76). b ¼ 0 and D U b ¼ 0. 2. For n = 0, let U

ð106Þ

2663

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668

3. Compute KM and KG and Fext and Fk. b from Eq. (106). 4. Compute D U b < tolerance exit else, let U bn ¼ U b and n = n + 1 and go to 3. b n1 þ D U 5. check for convergence, if norm ðD UÞ 4. Numerical examples In this section the presented method is tested with some numerical examples. The most advantage of this method is that the convergency rate is very high and shear locking is eliminated. 4.1. Pinched cylinder with end rigid diaphragms A short cylinder with pinching vertical force at the middle section, and two rigid diaphragms at the ends, is studied. The geometric and material property of cylinder is as below: L = 600 Radius = 300, thickness = 3, E = 3000,m = 0.3, k0 = 24.3, k0 = 300 and H0 = 0. Because of symmetry only one octane of cylinder is modeled. The results are derived for a 30  30 mesh. It can be seen (form Fig. 5) that the results are very close to the results presented in [29, 30]. 4.2. A simply supported plate under uniform lateral load In this example the deformation of a rectangular plate, under uniform lateral load is studied. The material and geometric properties are as below:

a ¼ b ¼ 407 mm thickness ¼ 7:6 mm E ¼ 2:11  105 N=mm2 ; k0 ¼ 250N=mm2

0

and k ¼ 50

Because of symmetry, only one quarter of plate was modeled. The results are found for a 20  20 mesh. The results are shown in Fig. 6 and have been compared with literature. 4.3. A simply supported trapezoidal plate under uniform lateral load In this example the deformation of a simply supported trapezoidal plate under uniform lateral load is studied. The geometry and material property of plate are as below:

a ¼ 1m; ;E ¼ 2  105 Mpa;

k0 ¼ 250Mpa;

h ¼ 0:01; a

0

k ¼ 1000:

This problem experience very large elasto plastic deformation. The results are shown in Fig. 7. In this figure W0 denotes out of plane displacement of the center of the plate. 4.4. A simply supported skew plate under uniform lateral load In this example the deformation of a simply supported skew plate under uniform lateral load is studied. The geometry and material property of the plate are as below:

a ¼ 1m;

b ¼ 1 m; E ¼ 2  105 Mpa;

k0 ¼ 250Mpa;

h ¼ 0:01; a

0

k ¼ 1000:

9000 8000 7000 Rigid diaphragm support

6000

F

5000

presented by[12]

4000

present mthod

3000

Simo[5]

2000

Brank[13]

1000 0 0

100

200

300

400

Rigid diaphragm support

W0 Fig. 5. Vertical deflection at center of the pinched cylinder.

2664

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668

a

b 12

Collapse load (0.509) )

10 8 6 4

elastic present method

2

Mohammed[9] 0 0

1

2

3

4

normalized central displacement(w/h) Fig. 6. Vertical displacement at the center of a simply supported plate under uniform lateral load.

0.025

Fa^4/(Dh)10^(-8)

0.02 0.015

isotropic hardening (b/a=0.2) elastic (b/a=0.2)

0.01

isotropic hardening (b/a=0.6) elastic (b/a=0.6)

0.005 0 0

20

40

60

80

100

w0/h Fig. 7. Vertical displacement at the center of a simply supported trapezoidal plate under uniform lateral load.

The problem is solved for a = 30°, a = 60° and a = 45° and the results are obtained for elastic and plastic deformation and they are shown in Fig. 8. In this figure W0 denotes out of plane displacement of the center of the plate. This problem also experience very large deformation. 5. Conclusion A new non-linear method based on the Cosserat theory, with constrained director, has been presented for large elastoplastic deformation with isotropic and kinematic hardening. The most advantage of this new method is that it eliminates the shear locking problem during the thin shell analyses and the convergency rate is very good. The material and geometric

2665

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668

a

0.12

α

Fa^4/(Dh)10^(-8)

0.1

a elastic(Alfa=30deg)

0.08

plastic(Alfa=30deg)

0.06

elastic(Alfa=60deg) plastic(Alfa=60deg)

0.04

elastic(Alfa=45deg) 0.02

plastic(Alfa=45deg)

0 0

20

40

60

80

100

w0/h Fig. 8. Vertical displacement at the center of a simply supported skew plate under uniform lateral load.

stiffness matrices, for finite element solution, have been derived through linearization of virtual work equation. For numerical solution, a nine node isoparametric element was implemented. Consistent elasto plastic tangent moduli is derived for FE solution. The method is computationally efficient and the numerical results exhibited very good agreement with the known values in the literature. Appendix A According relation (41) rename:

 H¼

C1

C2

C3

C4

"

1 ¼

# H2 : H4

H1 H3

ðA:1Þ

By using Eq. (40) we can write:

H1 ði; iÞ ¼ H2 ði; iÞ ¼ H3 ði; iÞ ¼ H4 ði; iÞ ¼

1 þ cpnþ1 w4 ði; iÞ gði; iÞ cmnþ1 w2 ði; iÞ gði; iÞ cmnþ1 w3 ði; iÞ

ðA:2Þ

;

ðA:3Þ

;

; gði; iÞ 1 þ cpnþ1 w1 ði; iÞ gði; iÞ

ðA:4Þ ðA:5Þ

;

Where gði; iÞ ¼ ð1 þ cpnþ1 w1 ði; iÞÞð1 þ cpnþ1 w4 ði; iÞÞ  c2mnþ1 w2 ði; iÞw3 ði; iÞ:

ðA:6Þ

So

2 @H @ cpnþ1

@H1

6 @cp ¼ 4 nþ1 3 @H

@ cp

nþ1

@H2 @ cp

nþ1

@H4 @ cp

3 7 5;

ðA:7Þ

nþ1

where 2

2

2

@H1 ði; iÞ w1 ði; iÞ  cmnþ1 w2 ði; iÞw3 ði; iÞw4 ði; iÞ  2cpnþ1 w1 ði; iÞw4 ði; iÞ  cpnþ1 w1 ði; iÞðw4 ði; iÞÞ ¼ ; @ cpnþ1 ðgði; iÞÞ2

ðA:8Þ

@H2 ði; iÞ cmnþ1 w2 ði; iÞw4 ði; iÞ þ cmnþ1 w1 ði; iÞw2 ði; iÞ þ 2cmnþ1 cpnþ1 w1 ði; iÞw2 ði; iÞw4 ði; iÞ ¼ ; @ cpnþ1 ðgði; iÞÞ2

ðA:9Þ

@H3 ði; iÞ cmnþ1 w3 ði; iÞw4 ði; iÞ þ cmnþ1 w1 ði; iÞw3 ði; iÞ þ 2cmnþ1 cpnþ1 w1 ði; iÞw3 ði; iÞw4 ði; iÞ ¼ ; @ cpnþ1 ðgði; iÞÞ2

ðA:10Þ

2

2

2

@H4 ði; iÞ w4 ði; iÞ  cmnþ1 w1 ði; iÞw2 ði; iÞw3 ði; iÞ  2cpnþ1 w1 ði; iÞw4 ði; iÞ  cpnþ1 w4 ði; iÞðw1 ði; iÞÞ ¼ @ cpnþ1 ðgði; iÞÞ2

ðA:11Þ

2666

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668

And

@H1 ði; iÞ 2cmnþ1 w2 ði; iÞw3 ði; iÞð1 þ cpnþ1 w4 ði; iÞÞ ¼ ; @ cmnþ1 ðgði; iÞÞ2

ðA:12Þ

2 2 2 @H2 ði; iÞ w2 ði; iÞ  cmnþ1 w3 ði; iÞðw2 ði; iÞÞ  cpnþ1 w1 ði; iÞw2 ði; iÞw4 ði; iÞ  cpnþ1 w2 ði; iÞðw1 ði; iÞ þ w4 ði; iÞÞ @H3 ði; iÞ ¼ @ cmnþ1 @ cmnþ1 ðgði; iÞÞ2

¼

w3 ði; iÞ  c2mnþ1 w2 ði; iÞðw3 ði; iÞÞ2  c2pnþ1 w1 ði; iÞw3 ði; iÞw4 ði; iÞ  cpnþ1 w3 ði; iÞðw1 ði; iÞ þ w4 ði; iÞÞ ðgði; iÞÞ2

;

@H4 ði; iÞ 2cmnþ1 w2 ði; iÞw3 ði; iÞð1 þ cpnþ1 w1 ði; iÞÞ ¼ @ cmnþ1 ðgði; iÞÞ2

ðA:13Þ

ðA:14Þ

And then we can write:

@H @H @H ¼ þ @ c1 @ cpnþ1 @ cmnþ1

@H @H @H ¼  @ c2 @ cpnþ1 @ cmnþ1

and

ðA:14Þ

Appendix B In this appendix we derive explicit expression for elasto-plastic tangent moduli. From Eq. (12) we have:

drnþ1 ¼ Cðdenþ1  depnþ1 Þ:

ðB:1Þ

Also from Eq. (29)1 we can write:

depnþ1 ¼

X

dcanþ1

a2k

! @fa;nþ1 @ 2 fa;nþ1 ; þ canþ1 d r @r @ r2

ðB:2Þ

where k = {ljfl,n+1 = 0} l = 1, 2 (active surfaces). So from (B.1) and (B.2) 1

C drnþ1 ¼ denþ1 

X a2k

So

C1 þ

X

canþ1

a2k

! 2 @fa;nþ1 a @ fa;nþ1 dcnþ1 þ cnþ1 dr : @r @ r2 a

! X @ 2 fa;nþ1 @fa;nþ1 drnþ1 ¼ denþ1  dcanþ1 : 2 @r @r a2k

ðB:3Þ

ðB:4Þ

By renaming: 1 B1 r ¼C þ

X

canþ1

a2k

drnþ1 ¼ Br denþ1 

@ 2 fa;nþ1 ; @ r2 X a2k

dcanþ1

ðB:5Þ !

@fa;nþ1 : @r

ðB:6Þ

For an active surface fk = 0 so dfk,n+1 = 0 So

 T  T @fk @fk @fk dr þ dP þ dp ¼ 0; @r @p @P Also; da ¼

X a2k

dcanþ1

@fa;nþ1 @ 2 fa;nþ1 dp: þ canþ1 @p @p2

Also dp ¼ Dda so  D1 dp ¼

X a2k

So 1

 D

ðB:7Þ

dcanþ1

@fa;nþ1 @ 2 fa;nþ1 dp: þ canþ1 @p @p2

! X @ 2 fa;nþ1 @fa;nþ1 dp ¼ þ cnþ1 dcanþ1 : @p2 @p a2k a

ðB:8Þ

ðB:9Þ

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668

2667

By renaming

@ 2 fa;nþ1 ; @p2 X @fa;nþ1 dp ¼ Bp dcanþ1 @p a2k

1 B1 þ canþ1 p ¼ D

ðB:10Þ ðB:11Þ

presented yield function, from Eq. (29)1, (29)2 and (29)4, it is obvious that

dP ¼ D0 dep :

ðB:12Þ

By replacing (B.6), (B.11), (B.12) in (B.7) we have:

! !   T T X X @fk @fk  0 p  @fk a @fa;nþ1 a @fa;nþ1 Br denþ1  Br dcnþ1 D de þ dcnþ1 Bp ¼ 0: þ @r @ r2 @r @p @p a2k a2k

ðB:13Þ

By replacing dep from (B.2) in the above equation we have:

! !!  T  T 2 X X @fk @fk @fk 0 a @fa;nþ1 a @ fa;nþ1 a @fa;nþ1 ¼ 0: þ dr þ D dcnþ1 þ cnþ1 dr dcnþ1 Bp @r @r @r @ r2 @p @p a2k a2k

ðB:14Þ

So

! !  T T T X X @fk @fk 0 X a @ 2 fa;nþ1 @fk @fk 0 a @fa;nþ1 a @fa;nþ1 ¼ 0: dr þ þ  D cnþ1 D dcnþ1 dcnþ1 Bp @r @r @ r2 @r @r @p @p a2k a2k a2k

ðB:15Þ

By renaming: T

yk ¼

T

@fk @fk 0 X a @ 2 fa;nþ1  D cnþ1 : @r @r @ r2 a2k

ðB:16Þ

We have:

yk dr þ

X a2k

T

@fk 0 @fa;nþ1 @fk @fa;nþ1 dcanþ1  D  Bp @r @r @p @p

! ¼ 0:

ðB:17Þ

And by replacing dr from Eq. (B.6) we have:

yk Br denþ1  yk Br

X a2k

@fa;nþ1 X a @fk 0 @fa;nþ1 @fk @fa;nþ1 dcnþ1 þ dcnþ1  D  Bp @r @r @r @p @p a2k T

a

! ¼ 0:

ðB:18Þ

So

yk Br denþ1 þ

X

T

a

dcnþ1

a2k

@fa;nþ1 @fk 0 @fa;nþ1 @fk @fa;nþ1 yk Br  D  Bp @r @r @r @p @p

! ¼ 0:

ðB:19Þ

By renaming T

zak ¼ yk Br

@fa;nþ1 @fk 0 @fa;nþ1 @fk @fa;nþ1 þ D þ Bp : @r @r @r @p @p

ðB:20Þ

So

yk Br denþ1 

X

dcanþ1 zak ¼ 0;

ðB:21Þ

a2k

when only one surface is active then:

dcanþ1 ¼ z1 aa ya Br denþ1

no sum on a:

If both of surfaces are active then:

dcknþ1 ¼

X

z1 ka yk Br denþ1

ðB:22Þ

a2k

And then:

X X @fa;nþ1 drnþ1 ¼ Br  Br zab yb Br : denþ1 @r b2k a2k

ðB:23Þ

2668

A. Ghassemi et al. / Applied Mathematical Modelling 35 (2011) 2650–2668

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]

P.M. Naghdi, The theory of shells and plates, in: C. Truesdell, Handbuch der Physik, Via/2, 1972. J.L. Ericksen, C. Truesdell, Exact theory of stress and strain in rods and shells, Arch. Ration. Mech. Anal. 1 (1959) 295–323. J.L. Sandres, Nonlinear theories for thin shells, Arch. Ration. Mech. Anal. 21 (1962) 21–36. A.E. Green, P.M. Naghdi, A general theory of a Cosserat surface, Arch. Ration. Mech. Anal. 20 (1965) 287–308. S. Ahmad, B.M. Irons, O.C. Zienkiewicz, Analysis of thick and thin shell structures by curved finite elements, Int. J. Numer. Methods Eng. 2 (1970). 619– 451. T. J Haghes, W.K. Liu, Nonlinear finite element analysis of shells: Part I Three-dimensional shells, Comput. Methods Appl. Mech. Eng. 26 (1981) 331– 362. J. R Hughes, E. Carnoy, Nonlinear finite element analysis of shell formulation accounting for large membrane strains, Comput. Methods Appl. Mech. Eng. 39 (1983) 69–82. K.J. Bathe, Finite Element Procedures, Prentice-Hall, Englewood Cliffs NJ, 1996. T.J. R Hughes, The Finite Element Method, Prentice-Hall, Englewood Cliffs NJ, 1987. T. Belytschko, W.K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, John Wiley & Sons LTD, 1987. G. Wempner, Finite elements, finite rotations and small strains of flexible shells, Int. J. Solids Struct. 5 (1969) 117–153. J. Argyris, An excursion into large rotations, Comput. Methods Appl. Mech. Eng. 32 (1982) 85–155. H. Parisch, An investigation of a finite rotation four node assumed strain shell element, Int. J. Numer. Methods Eng. 31 (1991) 127–150. N. Buechter, E. Ramm, Shell theory versus degeneration comparison in large rotation finite element analysis, Int. J. Numer. Methods Eng. 34 (1992) 39– 59. C. Sansour, H. Bufler, An exact finite rotation shell theory; its mixed variational formulation and its finite element implementation, Int. J. Numer. Methods Eng. 34 (1992) 73–115. X. Peng, M.A. Crisfield, A simple four nodded co-rotational formulation for shells using the constant stress/constant moment triangle, Int. J. Numer. Methods Eng. 35 (1992) 1829–1847. L. Jiang, M.W. Chernuk, A simple four-noded co-rotational shell element for arbitrarily large rotations, Comput. Struct. 53 (1994) 1123–1132. C.S. Liu, H.K. Hong, Using comparison theorem to compare co-ratational stress rates in the model for perfect elastoplasticity, Int. J. Solids Struct. 38 (2001) 2969–2987. P.M. Naghdi, Foundations of elastic shell theory, Prog. Solid Mech. 4 (1963) 1–90. S.S. Antman, Ordinary differential equations of nonlinear elasticity; Part I: Foundations of the theory of non-linearrly elastic rods and shells, Arch. Ration. Mech. Anal. 61 (1976) 307–351. A.E. Green, P.M. Naghdi, Theory of an elastic–plastic Cosserat surface, Int. J. Solids Struc. 4 (1968) 907–927. A.A. Ilyushin, Plasticity, Gostekhizdat, Moscow, 1948 (in Russian). M.A. Crisfield, Non-linear Finite Element Analysis of Solids and Structures, vol. 2, John Wiley and Sons, New York, 1997. J.C. Simo, D. Fox, On a stress resultant geometrically exact shell model. Part I: Formulation and optimal parametrization, Comput. Methods Appl. Mech. Eng. 72 (1982) 267–304. J.C. Simo, D. Fox, M.S. Rifai, On a stress resultant geometrically exact shell model. Part II: The linear theory: computational aspects, Comput. Methods Appl. Mech. Eng. 73 (1989) 53–92. J.C. Simo, J. G Kennedy, On a stress resultant geometrically exact shell Model. Part V: Nonlinear plasticity formulation and integration algorithms, Comput. Methods Appl. Mech. Eng. 96 (1992) 133–171. K.A. Mohammed, B. Skallerud, Simplified stress resultants plasticity on a geometrically nonlinear constant stress shell element, Comput. Struct. 79 (2001) 1723–1734. A.E. Green, W. Zerna, Theoretical Elasticity, Oxford University Press, 1968. K.D. Kim, G.R. Lomboy, A co-rotational quasi-conforming 4-node resultant shell element for large deformation elasto-plastic analysis, Comput. Methods Appl. Mech. Eng. 195 (2006) 6502–6522. B. Brank, D. Peric, F.B. Damjanic, On large deformation of thin elasto-plastic shells: implementation model for quadrilateral shell element, Int. J. Numer. Methods Eng.. 40 (1997) 689–726.