Revisited Simo algorithm for the plane stress state

Applied Mathematics and Computation 237 (2014) 730–751

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Revisited Simo algorithm for the plane stress state Sanda Cleja-T ß igoiu a,⇑, Nadia Elena Stoicutßa a,b a b

University of Bucharest, Faculty of Mathematics and Computer Science, Academiei 14, 010014 Bucharest, Romania University of Petrosßani, Department of Mathematics-Informatics, Institutului 20, 332006 Petrosani, Romania

a r t i c l e

i n f o

Keywords: Elasto-plastic model Plane stress Variational formulation Finite element method Trial state Radial Return Mapping Algorithm Numerical simulation

a b s t r a c t In this paper a new numerical approach of the elasto-plastic problem with mixed hardening in-plane stress state is proposed within the classical constitutive framework of small deformation. In the formalized problem, the non-zero normal component of the strain (namely, the component in the direction perpendicular to the plane of stress) is considered to be compatible with the plane stress. We eliminate the rate of strains which are developed normal to the plane and solve an appropriate plane problem in the strain setting. The integration algorithm for solving the elasto-plastic problem with mixed hardening realizes the coupling of Radial Return Algorithm, namely an update algorithm, with the finite element method applied to the discretized equilibrium balance equation. The solution of the pseudo-elastic nonlinear problem is then solved by Newton’s method. The numerical algorithm (which is an incremental one) consists of the Radial Return Mapping Algorithm, originally proposed by Simo and Hughes (1998) adapted to the stress plane problem and coupled with a pseudo-elastic problem. We rebuild the algorithmic formula giving rise to the plastic factor in terms of the trial stress state and the algorithmic elasto-plastic tangent moduli which is requested to evaluate the Jacobian that is necessary to determine the displacement. The proposed algorithm has been named the revisited Simo algorithm and is applied to exemplify the mode of integration of the bi-dimensional problem. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction This paper deals with a new approach of the elasto-plastic problem with mixed hardening in the plane stress state, within the classical constitutive framework of small deformation. To solve the boundary and initial value problem, we use an iterative solution for the discretized equilibrium balance equation, a pseudo-elastic problem, which is coupled with the update algorithm required to solve the rate-type constitutive model. Numerous references are devoted to elasto-plasticity with small strains with respect to the constitutive description, functional and numerical methods to solve the initial and boundary value problems. For the history and evolution of concepts that occur in the constitutive equations, as well as different theoretical formulations of elasto-plasticity, the reader is referred to Hill [11], Martin [19], Kachanov [7], Simo and Hughes [25] and Khan and Huang [8]. The evolution equations for tensorial hardening variables, which describe the mixed hardening, were proposed by Prager [22,23], Armstrong and Frederick [1], Chaboche [5]. Numerical tests are realized for elasto-plastic materials with kinematic hardening of Armstrong-Frederick type and scalar hardening function of Chaboche type and are compared with the experimental results ⇑ Corresponding author. E-mail addresses: [email protected] (S. Cleja-T ß igoiu), [email protected] (N.E. Stoicutßa). http://dx.doi.org/10.1016/j.amc.2014.03.126 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

731

obtained in laboratory by Broggiato et al. [4], using FEM code. The variation in time of the irreversible variables has been alternatively defined through the principle of the maximum plastic dissipation. This principle has been attributed by Hill [11] to Mises and considered, for instance, by Johnson [14], Lubliner [18] and Simo and Hughes [25]. This plays a basic role in the mathematical development of the problem in plasticity, see for instance Han and Reddy [10] and Duvaut and Lions [17]. The variational form of governing equations in plasticity in a discrete formalism was performed by Simo and Hughes [25] in Chapter 4, by assuming the coupling between the total free energy functional, the Lagrangian functional associated with the plastic dissipation, together with the Lagrange multiplier interpreted as plastic factor. In Simo and Hughes [25], the authors pay attention to the so-called Radial Return Mapping Algorithm, which is considered to be the central problem in computational plasticity, and they mention that the method was first proposed by Wilkins [26] and Krieg and Key [16]. The trial stress is a basic concept of the Radial Return Mapping Algorithm and is computed from the elastic type constitutive equation for the frozen plastic behaviour at the previously reached value. Within the J 2 plasticity theory the Radial Return Mapping Algorithm makes possible a closed-form solution of the problem, related to the rate-type description of constitutive models in elasto-plasticity. For an in-plane stress state, the algorithm can not be applied since the plane stress condition is violated. The procedure proposed by Simo and Hughes [25] follows Simo and Taylor [24] and performs the return mapping directly in the constrained plane stress space. Simo and Hughes [25] mentioned that the strain components e; ee and ep do not appear explicitly in their formulation since ee33 ¼  1m m ðee11 þ ee22 Þ, where m is the Poisson ratio, and ep33 ¼ ep11  ep22 . Thus e33 ¼ ee33 þ ep33 . The Radial Return Mapping Algorithm in its various versions and extensions is based on the implicit backward Euler method, applied to the system of differential equations which describe the models. The algorithms and their accuracy and stability have been studied for instance by Hughes and Taylor [13] and Ortiz and Popov [21]. The finite element method (FEM) was generally used to solve variational equalities. The main ideas related to the FEM, which have been developed by Bathe [2], Belytschko et al. [3], Fish and Belytschko [9], Hughes [12], Johnson [15] and applied to various elastic problems, say for instance the isoparametric FEM, will be applied in this paper. In Section 2, we formulate the differential type constitutive equations, which describe the behaviour of isotropic elastoplastic material with isotropic and kinematic hardening. In the formalized rate-type constitutive model for in-plane stress, the non-zero normal component of the strain e33 (namely, the component in the direction perpendicular to the plane of stress) has to be consistent with the requirement to have a plane stress state within the considered framework. We apply the procedure developed by Cleja-T ß igoiu [6] within the finite elasto-plasticity to eliminate the rate of strain which develops normal to the plane. We emphasize the peculiar aspects related to the generalized plane strain: (a) the matrix of elastic coefficients for isotropic material is replaced by the constitutive matrix of the elastic moduli E ð2Þ , a fourth-order tensor associated with bi-dimensional symmetric tensors, Sym2 ; (b) the normality of the flow rule to the yield surface is lost for the in-plane stress; (c) the new expression of the plastic factor is provided in terms of the plane strain rate. In Section 3 the initial and boundary value problem, Problem P, is formulated starting from the elastic type constitutive equation and the equilibrium equation to be satisfied by the stress, together with the rate-type constitutive description of the elasto-plastic material. To solve numerically Problem P, we consider the discretized weak formulation associated with the equilibrium equations, namely a pseudo-elastic problem, in order to find the stress rnþ1 and the displacement vector unþ1 at time tnþ1 . In Section 4, for a given incremental strain history, the update algorithm is built to solve the differential-type system which describes the in-plane stress state. We use here the Radial Return Algorithm proposed by Simo and Hughes [25], in a revised form adapted to the in-plane stress state. We rebuild the algorithmic expression of the plastic factor in terms of the trial stress state and derive the new algorithmic elasto-plastic tangent moduli for the isotropic-kinematic hardening model. Next, in Section 5 the FEM is used to solve the pseudo-elastic problem, which is reduced to a non-linear system for the displacement vector unþ1 at time t nþ1 . The displacement unþ1 is then solved by Newton’s method, while to evaluate the Jacobian of the system the previously calculated elasto-plastic moduli are requested. To exemplify the numerical integration of the bi-dimensional problem and prove the efficiency of the numerical algorithm proposed in Section 6, we consider two examples in Section 7. In the numerical simulations a trapezoidal plate, with and without a hole, has been deformed for an in-plane stress generated by periodic forces applied on the horizontal edge of the plate. The numerical results of the simulations for the in-plane stress were performed using the material parameters given by Broggiato et al. [4] and compared to the solution that corresponds to Simo and Taylor [24], the so-called Simo solution. Further we shall use the following notations: R is the set of real numbers, and R60 ¼ fx 2 Rjx 6 0g; þ Lin; Lin are the sets of second-order tensors and the corresponding elements with positive determinant, respectively; Sym  Lin the set of symmetric tensors; Sym2 the set of symmetric plane tensors; r0 ¼ dev r ¼ r  13 ðtrrÞI3 the deviatoric part of the stress tensor; E the tensor of elastic moduli (the fourth order tensor); ka ; la the Lamé elastic constants; E the Young modulus; m the Poisson ratio; fei g a Cartesian basis in R3 ; I2 ¼ dij ei  ej the second-order identity tensor in R2 ; I3 the identity tensor in R3 ;  1 I4 ¼ 2 dik djl þ dij dkl ei  ej  ek  el the fourth-order identity tensor;

732

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

r_ ¼ @ðrÞ=@t the time derivative of r; ruðx; tÞ ¼ ui;j ei  ej the displacement gradient; eðx; tÞ ¼ 12

P3



i;j¼1

infinitesimal strain tensor at time t 2 I; @ C F the partial derivative with respect to the tensorial field C of the function F . qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P3 1 2 kxk ¼ i;j¼1 xij the norm of x 2 Lin; hzi ¼ 2 ðz þ jzjÞ; 8z 2 R the ramp function;

@ui @xj

 @u þ @xj ei  ej the i

H is the Heaviside function, defined by HðzÞ ¼ 0; if z < 0 HðzÞ ¼ 1; if z P 0 a  b and A  B vector and tensor scalar products, respectively; A  B the tensorial product which is defined as a bilinear mapping which associates the element ðA  BÞC ¼ AðB  CÞ to any C 2 Lin. 2. Framework of constitutive elasto-plastic model in a plane stress and a generalized in-plane strain states We consider the body B to be a plate X  ½0; h, with X  R2 , and h the thickness of the plate, considered to be much smaller than a certain characteristic length of the plate. Let e1 ; e2 be the in plane axes and e3 be the unit vector of the normal to the plane. We denote by x ¼ ðx1 ; x2 Þ a material point in X and t 2 I ¼ ½0; T  Rþ . We assume that X is an open and bounded set with a smooth boundary @ X and its closure is denoted by X ¼ X [ @ X. We introduce the following notation: u : X  I ! R2 ; uðx; tÞ  ui ðx; tÞei is the displacement vector, e : X  I ! Sym is the infinitesimal strain tensor, r : X  I ! Sym is the Cauchy stress tensor, a : X  I ! Sym is the kinematic hardening variable, k : X  I ! R is the isotropic hardening variable. We restrict ourselves to a plane stress state, namely

r¼

2 X

rij ei  ej ; r 2 Sym2 ; i:e r0  rð2Þ

ð1Þ

i;j¼1

since

r13 ¼ r23 ¼ r33 ¼ 0. We use the following representation for an in-plane symmetric tensor r

r ¼ r11 e1  e1 þ r22 e2  e2 þ r12 e1  e2 þ r21 e2  e1 ; with r12 ¼ r21

ð2Þ

and we denote r by ðr11 ; r22 ; r12 ; r21 Þ in order to avoid the representation ðr11 ; r22 ; 2r12 Þ. The kinematic hardening variable a is defined by

a¼

3 X

aij ei  ej ; a 2 Sym

ð3Þ

i;j¼1

with a13 ¼ a23 ¼ 0 and is a zero trace tensor, i.e. tr a ¼ 0. The generalized strain state is characterized by a strain tensor e that has a non-zero component in the direction perpendicular to the shell, denoted by e33 , and added to the in-plane strain eð2Þ :

e ¼ eð2Þ þ e33 e3  e3 ; eð2Þ ¼

2 X

eij ei  ej ; eð2Þ 2 Sym2 :

ð4Þ

i;j¼1

The normal strain, e33 , has to be consistent with the requirement that a plane stress state occurs within the considered framework. The elasto-plastic model with mixed hardening is described by the following relations: The rate of elastic type constitutive equation is given by

r_ ¼ E e_ el ; e_ el ¼ e_  e_ p

ð5Þ

with the linear isotropic elastic constitutive equation, E e ¼ ka tre I3 þ 2la e , as a particular case of (5). The rate of the plastic strain tensor is described by associated flow rule el

el

el

e_ p ¼ k@ r0 F ðr0 ; a; kÞ

ð6Þ

The variation of the kinematic hardening variable a is given by the Prager law [22,23]

a_ ¼ C e_ p

ð7Þ

with tra ¼ 0 and C is kinematic hardening modulus. The variation of the isotropic hardening variable k is given by

k_ ¼

pffiffiffiffiffiffiffiffiffiffiffiffi e_ p  e_ p ;

ð8Þ

where the irreversible properties of the material are described in terms of the yield function F : DF  Sym  Sym  R ! R60 , dependent on the stress and hardening variables

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

F ðr0 ; a; kÞ ¼ Jðr0  aÞ  FðkÞ:

733

ð9Þ

Here Jðr0  aÞ ¼ kr0  ak is the second invariant of the deviator ðr0  aÞ and, for an in-plane stress state this reads as:

Jðr0  aÞ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  0        r11  a11 2 þ r022  a22 2 þ r033  a33 2 þ 2 r012  a12 2 :

ð10Þ

Function FðkÞ given by (9) is a strictly increasing and, in the linear case, it is defined by

rffiffiffi 2 ½Hk þ rY ; FðkÞ ¼ 3

ð11Þ

where H > 0 is the hardening parameter and rY represents the initial yield condition. The plastic factor, k, is a function of the state of material, which is defined by the Kuhn–Tucker condition

F ðr0 ; a; kÞ 6 0;

k P 0;

kF ðr0 ; a; kÞ ¼ 0

ð12Þ

and the consistency condition

kF_ ðr0 ; a; kÞ ¼ 0:

ð13Þ

Remark. Further we consider a 2 Sym without any mention about its trace, thus the formulae derived below allow for symmetric notations. _ given by (8), becomes: The rate of the isotropic hardening variable, k,

k_ ¼ k:

ð14Þ

In the following, if we used the tensorial representations for the plane stress state and plane kinematic variable



r0 ¼ r0ð2Þ þ r033 e3  e3 a ¼ að2Þ þ a33 e3  e3

)

r0  a ¼ r0ð2Þ  að2Þ þ ðr033  a33 Þ e3  e3

ð15Þ

and we denote by

8 0 ^ ¼ að2Þ ; > < s ¼ rð2Þ ; a 0 r33 ¼ trs; trs ¼ s11 þ s22 ; > : a33 ¼ tra^ ; tra^ ¼ a^ 11 þ a^ 22

ð16Þ

then we have

r0  a ¼ ðs  a^ Þ  ðtrðs  a^ ÞÞe3  e3 :

ð17Þ

In this case, the norm of the difference ðr  aÞ is given by 0

kr0  ak ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ^ k2 þ ðtr ðs  a ^ ÞÞ2 : ks  a

The passage from matrix

0

2 3

B1 B P¼B 3 @ 0 0

r  rð2Þ  Sym2 to s ¼ r0ð2Þ can be characterized by an invertible mapping P. given by the following

 13 0

0

1

0

0C C C; 1 0A

0

0 1

2 3

ð18Þ

0

s ¼ P r;

r ¼ rð2Þ :

ð19Þ

Hence, the yield function with respect to the variables previously introduced becomes

~ ðs; a ^ ; kÞ ¼ F ðr0 ; a; kÞ  F

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffi 2 ^ k2 þ ðtrðs  a ^ ÞÞ2  ½Hk þ rY : ks  a 3

ð20Þ

If the yield function is written in terms of s, then one can say that

~ ðs; a ~ ðs; a ^ ; kÞ and @ r0 F ðr0 ; a; kÞ ¼ @ trs F ^ ; kÞ; prð2Þ ð@ r0 F ðr0 ; a; kÞÞ ¼ @ s F 33

ð21Þ

  where @ r033 F ðr0 ; a; kÞ ¼  @ r011 F ðr0 ; a; kÞ þ @ r022 F ðr0 ; a; kÞ . Consequently, whenever the yield function F ðr0 ; a; kÞ and its par~ ðs; a ^ ; kÞ and its corresponding derivatives, tial derivatives appear, these can be replaced by the new yield function F respectively.

734

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

Remark. We restrict ourselves to the linear elastic constitutive equation. If we introduce the standard notations



E¼



la 3ka þ 2la k   ; m¼  a  2 ka þ la ka þ la

then we can introduce the new representations in terms of E and

ð22Þ

m instead of ka and la .

Proposition 1. 1. For a given history of the total deformation t ! eðtÞ, the rate of plane stress is calculated by the constitutive equation





r_ ¼ E ð2Þ e_ ð2Þ  kðrÞ @ s F~ ðs; a^ ; kÞ ;

ð23Þ

where the original plastic factor k is replaced by kðrÞ , i.e k ¼ kðrÞ , in order to emphasize that the plastic factor is dependent on the plane components of the stress and strain tensors only. The tensor E ð2Þ : Sym2 !Sym2 ,

0 E ð2Þ

E 1m2

B mE B 1m2 ¼B B 0 @

mE 1m2 E 1m2

0

0

E 1þm

0

0

0

0

0

1

C 0 C C 0 C A

ð24Þ

E 1þm

characterizes the elastic isotropic property of the material in terms of E and m where I4 denotes the identity, i.e. E ð2Þ I4 ¼ E ð2Þ . 2. The history of total deformation has to be compatible with a generalized plane strain, with e33 dependent on eð2Þ only, i.e:

e_ 33 ¼ 

m 1m

ðe_ 11 þ e_ 22 Þ 

1  2m ðrÞ ~: k @ trs F 1m

ð25Þ

~ ðs; a ^ ; kÞ ¼ 0 as 3. The plastic factor kðrÞ : X  I ! R is calculated on the yield surface F

kðrÞ ¼

hbi HðF Þ; h

ð26Þ

where

b¼

   E  ~ þ m @ trs F ~ e_ 11 þ E ~ þ m @ trs F ~ e_ 22 þ 2E ð@ s F ~ Þe_ 12 ; @ s11 F @ s22 F 1þm 1þm 1 þ m 12 1m 1m

rffiffiffi 

 2 2 E E Eð1  2mÞ 2 ~ ~ H: @ trs F þ h¼ þ C @sF þ þCþ 1þm 1þm 1  m2 3

ð27Þ



ð28Þ

Here, the hardening parameter h is supposed to be strictly positive and H is the Heaviside function.

Proof. 1. The components of the rate of the stress tensor can be derived from (5), together with (6), and are given by the following system:

8 r_ 11 ¼ ka ðe_ 11 þ e_ 22 þ e_ 33 Þ þ 2la e_ 11  2la k @ r011 F ðr0 ; a; kÞ; > > > > < r_ 22 ¼ ka ðe_ 11 þ e_ 22 þ e_ 33 Þ þ 2l e_ 22  2l k @ r0 F ðr0 ; a; kÞ; a a 22 > r_ 33 ¼ ka ðe_ 11 þ e_ 22 þ e_ 33 Þ þ 2la e_ 33  2la k @ r033 F ðr0 ; a; kÞ; > > > :_ r12 ¼ 2la e_ 12  2la k @ r012 F ðr0 ; a; kÞ: 2. In (29),

ð29Þ

r33 ¼ 0, and therefore we can obtain the derivative of e33

e_ 33 ¼ 

  ka 2l a ðe_ 11 þ e_ 22 Þ  k @ r011 F ðr0 ; a; kÞ þ @ r022 F ðr0 ; a; kÞ : ka þ 2la ka þ 2la

ð30Þ

We can introduce representation (25) in terms of E and m for (30) if the notations (22) are accounted for. In this case, by eliminating the component e33 from system (29) via (30), equations (29) become

8   > r_ 11 ¼ 1Em2 ðe_ 11 þ me_ 22 Þ  kðrÞ 1Em2 @ s11 F~ þ m@ s22 F~ ; > > <   r_ 22 ¼ 1Em2 ðe_ 22 þ me_ 11 Þ  kðrÞ 1Em2 @ s22 F~ þ m@ s11 F~ ;   > > > ~ : : r_ 12 ¼ E e_ 12  kðrÞ @ s12 F 1þm

ð31Þ

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

735

3. We consider the consistency condition F_ ðr0 ; a; kÞ ¼ 0 on the yield surface F ðr0 ; a; kÞ ¼ 0 for kðrÞ > 0

@ r0 F  r_ 0 þ @ a F  a_ þ @ k F  k_ ¼ 0

ð32Þ

together with (7), (14) and (23). These lead to

#  2 2 2 2  rffiffiffi E 2 H þ C @ r011 F þ @ r022 F þ @ r033 F þ 2 @ r012 F þ 1þm 3  E  @ r011 F e_ 11 þ @ r022 F e_ 22 þ @ r033 F e_ 33 þ 2ð@ r012 F Þe_ 12 ¼ 0 þ 1þm

 kðrÞ

"

ð33Þ

since @ r0 F  r_ ¼ @ r0 F  r_ and @ a F ¼ @ r0 F . Substituting e_ 33 from (25) and @ r033 F via (21), relation (33) yields: 0

rffiffiffi # 

  E E Eð1  2mÞ ~ 2 þ @ s F ~ 2 þ 2 @ s F ~ 2 þ ~ 2 þ 2H @ þ C @ s11 F þ C þ F trs 22 12 1þm 1þm 1  m2 3     E m E m 2E ~þ ~ e_ 11 þ ~þ ~ e_ 22 þ ~ Þe_ 12 ¼ 0: @ s11 F @ trs F @ s22 F @ trs F ð@ s F þ 1þm 1þm 1 þ m 12 1m 1m

 kðrÞ

"

ð34Þ

We finally obtain:

rffiffiffi # 

   2 2 E E Eð1  2 m Þ @ trs F ~ þ ~ þ 2H þ E ~ þ m @ trs F ~ e_ 11 k þ C @ s F þCþ @ s11 F 2 1þm 1þm 1m 3 1þm 1m  E  m 2E ~þ ~ e_ 22 þ ~ Þe_ 12 ¼ 0: þ @ s22 F @ trs F ð@ s F 1þm 1 þ m 12 1m ðrÞ

"

ð35Þ

The expressions of the complementary plastic factor b and the hardening parameter h given by (27) and (28) follow immediately from (35). h

^ , is introduced, then the variation of the plastic strain Proposition 2. If a new variable, s ¼ s  a

e_ p ¼ e_ pð2Þ þ e_ p33 where

8 p ðrÞ ^ sa ffi; > < e_ ð2Þ ¼ k pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ^ k2 þðtrðsa ^ ÞÞ2 ksa

^Þ trðsa > ffi : e_ p33 ¼ kðrÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2

ð36Þ

^ k þðtrðsa ^ ÞÞ ksa

and the variation of the kinematic hardening variable

a_ ¼ a^_ þ a_ 33 where

8 ^ sa > ^_ ¼ C kðrÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi;
^Þ trðsa > ffi : a_ 33 ¼ C kðrÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2

ð37Þ

^ k þðtrðsa ^ ÞÞ ksa

as well as the yield function (20), can be equivalently written in terms of the variable, s. Proof. Eqs. (36) and (37) are obtained as a direct consequence of the constitutive relations (6), (7), together with (20) and (21). h

3. The boundary and initial value problem The boundary conditions are formulated on the boundary @ X ¼ C of the body, C ¼ Cu [ Cr , where Cu \ Cr ¼ ;. We can formulate the following problem: Problem P. Given b : X ! R2 the body forces, f and g functions defined on parts of @ X, and having the values in R2 , find u; r; e; ep ; a; k defined on X  ½0; TÞ and satisfying the differential type constitutive equations

8 p s ; e_ ¼ kðrÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > ksk2 þðtrðsÞÞ2 > ð2Þ > > > trðsÞ > > e_ p33 ¼ kðrÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; > > ksk2 þðtrðsÞÞ2 > > > > s > > a^_ ¼ kðrÞ C pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; > > ksk2 þðtrðsÞÞ2 > > < trðsÞ ðrÞ a_ 33 ¼ k C pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; ksk2 þðtrðsÞÞ2 > > >_ > > k ¼ kðrÞ ; > > > > ^ Þ; > kðrÞ ¼ hbi HðF > h > > qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffi > > > ^ ðs; kÞ ¼ ksk2 þ ðtrðsÞÞ2  2½Hk þ rY ; > >F 3 > > : s ¼ s  a^ ; s ¼ Pr;

ð38Þ

736

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

together with (27) and (28). Here the Cauchy stress tensor, r, is given by the following elastic type constitutive equation



r ¼ E ð2Þ eð2Þ  epð2Þ



ð39Þ

and satisfies the boundary and initial value problem:

8 div ½r þ b ¼ 0 > > > > > > < rn ¼ f u¼g > > > > uð0; xÞ ¼ u0 ðxÞ; > > : p e ð0; xÞ ¼ ep0 ðxÞ;

in X  I; in Cr  I; in Cu  I;

ð40Þ

rð0; xÞ ¼ r0 ðxÞ; eð0; xÞ ¼ e0 ðxÞ; að0; xÞ ¼ a0 ðxÞ; kð0; xÞ ¼ k0 ðxÞ;

where n is the outward normal field. The set of kinematically admissible velocity field at every time t is denoted by:

W ad ¼



wi 2 L2 ðXÞ;

@wi 2 L2 ðXÞ; wi jCu ¼ g i ; i; j ¼ 1; 2  H1 ðXÞ  H1 ðXÞ; @xj

ð41Þ

where H1 ðXÞ is the Sobolev space. To solve Problem P at every time t 2 ½0; T, for any x 2 X, we use an iterative solution of the boundary and initial value problem for the discretized equilibrium balance equation coupled with an update algorithm. We discretize the weak formulation of problem (40) together with (39), i.e. we replace the continuous time interval ½0; T  by the sequence of discrete times t0 ; t1 ; . . . ; tN with t nþ1 ¼ t n þ Dt, n = 0, . . . , N, and formulate the pseudo-elastic problem P enþ1 at time tnþ1 : Problem P enþ1 . Find the displacement field unþ1 , which satisfies the discretized weak formulation:

Z

X

  E ð2Þ ½eð2Þ nþ1  ½epð2Þ nþ1  eðwnþ1 Þdx ¼ Lðwnþ1 Þ

ð42Þ

with Lðwnþ1 Þ given by

Lðwnþ1 Þ ¼

Z

f  wnþ1 dCr þ

Cr

Z

rnþ1 n  gdCu þ

Cu

Z

b  wnþ1 dx;

X

8wnþ1 2 W had ;

ð43Þ

where W had is a finite-dimensional approximation to W ad . At each iteration step n þ 1, one solves a so-called pseudo-elastic problem P enþ1 , where wnþ1 is the test function associated with the ðn þ 1Þ iteration step. 4. The Radial Return Algorithm ^ ; k is derived by using the constitutive equation For a given incremental strain history, an update algorithm for s; epð2Þ ; a (38). Next, we consider the weak formulation associated with equations (40) to find the stress rnþ1 and the displacement vector unþ1 at t nþ1 . System (38) is solved based on the Radial Return Algorithm proposed by Simo and Hughes [25], applied at any time t nþ1 and adapted here to the in-plane stress. The elastic type constitutive equation (39) at t nþ1 leads to

rnþ1 ¼ E ð2Þ ð½eð2Þ nþ1  ½epð2Þ nþ1 Þ:

ð44Þ

In order to numerically solve the differential-type equation (38), we apply the backward Euler method using the initial condition at t n . Thus

h

i snþ1 ; epð2Þ þ ½k nþ1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 2 ksnþ1 k þ ðtrðsnþ1 ÞÞ2  p    trðsnþ1 Þ ; e33 nþ1 ¼ ep33 n  ½k nþ1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ksnþ1 k2 þ ðtrðsnþ1 ÞÞ snþ1 ; a^ nþ1 ¼ a^ n þ ½k nþ1 C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ksnþ1 k þ ðtrðsnþ1 ÞÞ2 trðsnþ1 Þ ½a33 nþ1 ¼ ½a33 n  ½k nþ1 C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; ksnþ1 k2 þ ðtrðsnþ1 ÞÞ2

epð2Þ

i

nþ1

¼

h

knþ1 ¼ kn þ ½k nþ1 ;

h i where ½k nþ1 ¼ kðrÞ

nþ1

Dt ¼ ½kr nþ1 ðt nþ1  t n Þ and

snþ1 ¼ snþ1  a^ nþ1 .

ð45Þ

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

737

System (45) at t nþ1 is associated with the discrete Kuhn–Tucker conditions

^ ðsnþ1 ; knþ1 Þ 6 0; F

½k nþ1 P 0;

^ ðsnþ1 ; knþ1 Þ ¼ 0; ½k nþ1 F

ð46Þ

where the yield function at t nþ1 is defined by

^ ðsnþ1 ; knþ1 Þ ¼ F

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffi 2 2 ½Hknþ1 þ rY : ksnþ1 k2 þ ðtrðsnþ1 ÞÞ  3

ð47Þ

Elastic Predictor. At this step, the discrete problem is reformulated in terms of the trial state. The trial state at t nþ1 corresponds to the constitutive elastic type relations (44), written for the frozen plastic state at tn , namely





rtrnþ1 ¼ E ð2Þ ½eð2Þ nþ1  ½epð2Þ n ; tr

½epð2Þ nþ1 ¼ ½epð2Þ n ;

ð48Þ

a^ trnþ1 ¼ a^ n ; ktrnþ1 ¼ kn ; ½ep33 trnþ1 ¼ tr½epð2Þ n ; ½a33 trnþ1 ¼ tr½a^ n :

ð49Þ

^ tr is given by The expression of the trial yield function F nþ1

^ tr ¼ F ^ F nþ1



tr tr nþ1 ; knþ1

s



rffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    2 2 tr 2 tr  snþ1 þ trðsnþ1 Þ  ½H ktrnþ1 þ rY : ¼ 3

ð50Þ

Using (19), the deviator of the trial stress s is defined (via the formula) by

  strnþ1 ¼ PE ð2Þ ½eð2Þ nþ1  ½epð2Þ n :

ð51Þ

^ tr ¼ 0, then the trial state is elastic. If the trial state is located inside or on the trial yield surface defined by F nþ1 ^ tr > 0, then ½k can not be Plastic Corrector. If the trial state is located outside of the surface of plasticity, namely F nþ1 nþ1 zero and the trial state is taken as the initial condition for the problem of plastic corrector, using formulae (44) and (45). The consistency condition at tnþ1 is restored by the radial return algorithm. The tensor snþ1 and the hardening kinematic variable a^ nþ1 are given by

snþ1

snþ1 ¼ strnþ1  ½k nþ1 PE ð2Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; ksnþ1 k2 þ ðtrðsnþ1 ÞÞ2

ð52Þ

snþ1 a^ nþ1 ¼ a^ trnþ1 þ ½k nþ1 C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ksnþ1 k þ ðtrðsnþ1 ÞÞ2 since to pass from

r0ð2Þ  að2Þ to s  a^ , it is necessary to use the invertible mapping P given by (19). Consequently, we obtain

snþ1 ; snþ1 ¼ strnþ1  ½k nþ1 ðPE ð2Þ þ CI4 Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ksnþ1 k þ ðtrðsnþ1 ÞÞ2 where

ð53Þ

strnþ1 ¼ strnþ1  a^ trnþ1 . From (53) we also obtain

snþ1

½k nþ1 ðPE ð2Þ þ CI4 Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ strnþ1  snþ1 : ksnþ1 k2 þ ðtrðsnþ1 ÞÞ2

ð54Þ

Proposition 3.

^ nþ1 and strnþ1  a ^ nþ1 are collinear and they have the following direction 1. Tensors snþ1  a

^ nþ1 Þ  trðsnþ1  a ^ nþ1 Þe3  e3 ðsnþ1  a ^ nþ1 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi n 2 ^ nþ1 Þk þ ðtrðsnþ1  a ^ nþ1 ÞÞ2 kðsnþ1  a

or

s  trðsnþ1 Þe3  e3 ^ nþ1 ¼ qnþ1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : n ksnþ1 k2 þ ðtrðsnþ1 ÞÞ2

ð55Þ

2. The projections of the appropriate stress and the trial stress on the plane and on the normal to the plane are given by

snþ1

str

nþ1 ; nnþ1  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 kstrnþ1 k2 þ ðtrðstrnþ1 ÞÞ ksnþ1 k2 þ ðtrðsnþ1 ÞÞ2

trðstrnþ1 Þ trðsnþ1 Þ ; n?nþ1  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 kstrnþ1 k2 þ ðtrðstrnþ1 ÞÞ ksnþ1 k2 þ ðtrðsnþ1 ÞÞ2 ^ nþ1 ¼ nnþ1 þ n? where the following relation holds n nþ1 e3  e3 .

ð56Þ

738

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

Proof. As a consequence of formulae (15) and the collinearity of the deviators

str  trðstrnþ1 Þe3  e3 s  trðsnþ1 Þe3  e3 ^ nþ1 ¼ qnþ1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  qnþ1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : n 2 kstrnþ1 k2 þ ðtrðstrnþ1 ÞÞ ksnþ1 k2 þ ðtrðsnþ1 ÞÞ2

0 r0nþ1  anþ1 and rtrnþ1  atrnþ1 , we have

ð57Þ



Proposition 4. The following loading/unloading condition is satisfied ^ tr 6 0, then ½k  1. If F nþ1 nþ1 ¼ 0.

^ tr > 0, then ½k  ^ 2. If F nþ1 nþ1 – 0 and the plastic factor ½k nþ1 is calculated, on the yield surface F ðsnþ1 ; knþ1 Þ ¼ 0, by

½k nþ1

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffi tr 2 kstrnþ1 k2 þ ðtrðstrnþ1 ÞÞ  23½H knþ1 þ rY  qffiffi ¼ : 2 E H þ 3ð1 3 mÞ þ C

ð58Þ

Proof. To determine the expression of the plastic factor ½k nþ1 , we impose the following yield condition given by (47), ^ ðsnþ1 ; knþ1 Þ ¼ 0. In order to derive the expression of the plastic multiplier ½k  , we apply the trace to relation (54) and, F nþ1 accounting for (50)1, this leads to

     1 ½k nþ1 tr PE ð2Þ þ CI4 strnþ1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ tr strnþ1  tr ðsnþ1 Þ: 2 2 tr tr ksnþ1 k þ ðtrðsnþ1 ÞÞ

ð59Þ

If we eliminate trðsnþ1 Þ from (56)2, then we obtain

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    2 kstrnþ1 k2 þ ðtrðstrnþ1 ÞÞ  ksnþ1 k2 þ ðtrðsnþ1 ÞÞ2 : ½k nþ1 tr PE ð2Þ þ CI4 strnþ1 ¼ trðstrnþ1 Þ

ð60Þ

^ ðsnþ1 ; knþ1 Þ ¼ 0, this yields If we take into account relation (47) for F

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffi 2 tr ½H ðknþ1 þ ½k nþ1 Þ þ rY  ksnþ1 k2 þ ðtrðsnþ1 ÞÞ2 ¼ 3

ð61Þ

Finally, the relation (60) becomes

(rffiffiffi   ) qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffi tr PE ð2Þ þ CI4 strnþ1 2 2 tr 2   ¼ kstrnþ1 k2 þ ðtrðstrnþ1 ÞÞ  Hþ ½H knþ1 þ rY  ½k nþ1 3 3 tr strnþ1

ð62Þ

  if and only if tr strnþ1 is not vanishing. To evaluate the trace of the coefficient from (62), we use formulae (19) and (23)

   tr PE ð2Þ þ CI4 strnþ1 ¼



 E þ C trðstrnþ1 Þ 3ð1  mÞ

ð63Þ

Consequently, relation (62) becomes (58). h Algorithm 1. tr ^ tr 6 0, then ½k  ^ nþ1 ¼ a ^ trnþ1 ; ½ep nþ1 ¼ ½ep trnþ1 ; a If F nþ1 ¼ 0 and snþ1 ¼ snþ1 ; nþ1 ^ tr > 0, then ½k  If F cannot be zero and is calculated by (58) and nþ1 nþ1

tr

knþ1 ¼ knþ1 .

str

nþ1 snþ1 ¼ strnþ1  ½k nþ1 PE ð2Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 kstrnþ1 k2 þ ðtrðstrnþ1 ÞÞ

strnþ1 a^ nþ1 ¼ a^ trnþ1 þ ½k nþ1 C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 kstrnþ1 k2 þ ðtrðstrnþ1 ÞÞ trðstrnþ1 Þ ½a33 nþ1 ¼ ½a33 trnþ1  ½k nþ1 C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 kstrnþ1 k2 þ ðtrðstrnþ1 ÞÞ h i h itr strnþ1 ; epð2Þ ¼ epð2Þ þ ½k nþ1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi nþ1 nþ1 2 kstrnþ1 k2 þ ðtrðstrnþ1 ÞÞ  p    trðstrnþ1 Þ ; e33 nþ1 ¼ ep33 trnþ1  ½k nþ1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 kstrnþ1 k2 þ ðtrðstrnþ1 ÞÞ tr

knþ1 ¼ knþ1 þ ½k nþ1 ;

rnþ1 ¼ P1 snþ1

ð64Þ

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

739

The Radial Return Algorithm is completely defined if we show how the associated elasto-plastic moduli can be determined. These elasto-plastic moduli will be used into the discrete weak formulation (42). Algorithm 2. We determine the associated elasto-plastic tangent moduli from the algorithmic constitutive equations. Proposition 5. In the case of a plane stress state associated with the generalized plane strain, the algorithmic elasto-plastic tangent modulus can be expressed in the following modified form

( E ep nþ1

¼

^ tr 6 0; if F nþ1   ^ tr > 0; E ð2Þ ½H1 nþ1  ½H2 nþ1 nnþ1  E ð2Þ P nnþ1 þ n?nþ1 E ð2Þ P Ið2Þ if F nþ1

E ð2Þ

ð65Þ

where

½k nþ1 ffi PE ð2Þ ; ½H1 nþ1 ¼ I4  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 str  þ ðtrðstr ÞÞ2 nþ1 nþ1 0 ½H2 nþ1

1

ð66Þ

½k nþ1 C  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi AE ð2Þ :   2 2 2 E tr tr s  þ ðtr ðs ÞÞ H þ 3ð1mÞ þ C 3 nþ1 nþ1

B ¼ @qffiffi

1

Proof. By differentiating the algorithmic expression rnþ1 ¼ P1 snþ1 with respect to the total strain at tnþ1 , the elasto-plastic tangent moduli E ep are calculated using the following formula

@ rnþ1 @snþ1 ½E ep nþ1 ¼    P1   ; @ eð2Þ nþ1 @ eð2Þ nþ1

ð67Þ

which is referred to as the algorithmic tangent modulus, a notion firstly introduced  by Simo and Taylor [24]. We differentiate the algorithmic constitutive equation (52)1 with respect to eð2Þ nþ1 , in case when the plastic factor is non-vanishing, namely

0 @

eð2Þ

0

11

@str strnþ1 @ B B CC  ¼  nþ1    @½k nþ1 PE ð2Þ @qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi AA:   2 2 @ eð2Þ nþ1 @ eð2Þ nþ1 tr tr s  þ ðtr ðs ÞÞ nþ1 nþ1

@s  nþ1  nþ1

ð68Þ

In relation (68), we observe the presence of the elastic moduli E ð2Þ via the formula

@



@

eð2Þ

ðstrnþ1 Þ ¼ PE ð2Þ :



ð69Þ

nþ1

We calculate the derivative of the plastic factor, starting from the derived formula (58), as follows

0

1 0 qffiffih i 1 tr 2

þ r H k Y nþ1 ½ k  3 1 @ C C nþ1   B   B ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffiA ¼ qffiffi @qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi @1  q A:    2 @ eð2Þ nþ1 @ eð2Þ nþ1 2 E tr 2 tr str 2 þ ðtrðstr ÞÞ2  þ C H þ s þ ðtr ð s ÞÞ 3 3ð1mÞ nþ1 nþ1 nþ1 nþ1 @

ð70Þ

Note that, for any X 2 Sym2 , the following identities hold

@ @

 

@

eð2Þ @

eð2Þ



 nþ1





strnþ1  Ið2Þ  X ¼ PE ð2Þ X  Ið2Þ ¼ E ð2Þ PIð2Þ  X; ð71Þ

   str 2  X ¼ 2str  PE ð2Þ X ¼ 2E ð2Þ Pstr  X: nþ1 nþ1 nþ1

nþ1

Here Ið2Þ is the identity matrix in Symð2Þ . Thus relation (70) becomes

0

1

qffiffih i tr   2 H knþ1 þ rY E ð2Þ Pstr þ tr ðstr Þ E ð2Þ PIð2Þ ½k nþ1 3 B C nþ1 nþ1   @qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi q ffiffi ¼ A

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 :     @ eð2Þ nþ1 2 E str 2 þ ðtrðstr ÞÞ2 H þ 3ð1 str 2 þ ðtr ðstr ÞÞ2 3 mÞ þ C nþ1 nþ1 @



nþ1

Using relations (69) and (72), Eq. (67) can be expressed as

nþ1

ð72Þ

740

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

@

@s  nþ1 

eð2Þ

nþ1

qffiffih i tr   2 H knþ1 þ rY E ð2Þ Pstrnþ1 þ tr ðstrnþ1 Þ E ð2Þ PIð2Þ ½k nþ1 3 tr ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PE ð2Þ snþ1  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 PE ð2Þ PE ð2Þ  qffiffi ¼ PE ð2Þ  q    2 E str 2 þ ðtr ðstr ÞÞ2 H þ 3ð1 str 2 þ ðtr ðstr ÞÞ2 mÞ þ C 3 nþ1 nþ1 nþ1

nþ1

ð73Þ or, after some calculations are performed, the algorithmic elasto-plastic tangent modulus yields

0

1

½k nþ1 B C ffi PE ð2Þ A ¼ E ð2Þ @I4  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 2 tr tr s  þ ðtr ðs ÞÞ nþ1 nþ1 qffiffih i tr 2 H knþ1 þ rY   3 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi E ð2Þ nnþ1  E ð2Þ P nnþ1 þ E ð2Þ nnþ1  n?nþ1 E ð2Þ P Ið2Þ ;  qffiffi   2 2 E str  þ ðtr ðstr ÞÞ2 H þ 3ð1 3 mÞ þ C nþ1 nþ1

½E ep nþ1

ð74Þ

where nnþ1 is a deviatoric tensor, while n? nþ1 is a scalar field. Using again formula (58), relation (74) can be written as

0

1

0

1

½k nþ1 ½k nþ1 1 B C B C ffi q q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffiffi PE ¼ E ð2Þ @I4  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   A @ AE ð2Þ nnþ1 ð2Þ     2 2 E tr 2 tr str 2 þ ðtr ðstr ÞÞ2  þ C H þ s þ ðtr ð s ÞÞ 3 3ð1mÞ nþ1 nþ1 nþ1 nþ1    E ð2Þ P nnþ1 þ n?nþ1 E ð2Þ P Ið2Þ ;

½E ep nþ1

ð75Þ

We next introduce two fourth-order tensors, ½H1 nþ1 and ½H2 nþ1 , respectively, via formulae (66). h The algorithmic tangent modulus given by relation (75), together with the fourth-order tensors (66), reads as (65)2. If we ^ tr > 0. Othertake into account that formula (75) has been derived under the assumption that ½k nþ1 is non-vanishing, i.e. F nþ1 ^ tr 6 0 then ½k  wise, if F vanishes and (65) holds as a consequence of formulae (67)–(69). 1 nþ1 nþ1 Remark. Formulae (65) can be reduced to the formulae that correspond to a plane stress state, which can be found in Simo and Hughes [25], but associated with the plane strain, namely when the strain in the direction perpendicular to the plane is neglected.

5. Finite element method and non-linear equations for unþ1 The FEM is now utilized to solve the discretized weak formulation (42), i.e. the pseudo-elastic Problem P en . Applying the FEM, Problem Pen for unþ1 is reduced to solving an associate non-linear elastic-type problem. Finally, the non-linear equation is solved using a Newton–Raphson type algorithm. Proposition 6. The discrete version of the weak formulation, namely problem Pen reads as

^ ðu ^ nþ1 Þ ¼ F int ðu ^ nþ1 Þ  F ext G nþ1 ¼ 0; where F

int

ð76Þ

is the matrix of internal forces and

F ext nþ1

is the vector of external forces.

^ nþ1 should be the solution of the following Proof. We discretize the weak formulation (42) and the displacement field u equations: nel X 

^ enþ1 w

T

Z

¼

^ enþ1 w

T

e¼1

þ

Z

1

e¼1 nel X 

1

2 X nel X 

1

h

1

Z

1

Z

1

^ enþ1 w

T

þ



i¼1 e¼1

^ enþ1 w

1

T

Z





i



i ^e J e ðn; gÞ dndg ðNðn; gÞÞT Nðn; gÞb nþ1

1

h

i T ðNðn; 2ði  1Þ þ 1ÞÞ Nðn; 2ði  1Þ þ 1Þ^f enþ1 ðJ 1 Þei dn

h

i T ðNð2ði  1Þ þ 1; gÞÞ Nð2ði  1Þ þ 1; gÞ^f enþ1 ðJ 2 Þei dg

1

Z



r u^enþ1 ðn; gÞ Je ðn; gÞ dndg

h

1

i¼1 e¼1 2 X nel X

T

ðBe ðn; gÞÞ

1

1

 ^ 2 W had , where ðJ 1 Þe ¼ ðJ 1 Þe1 for any w

ðJ 1 Þe2

T

 ; ðJ 2 Þe ¼ ðJ 2 Þe1

ðJ 2 Þe2

T

with

ð77Þ

741

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

ðJ 1 Þe1 ðJ 2 Þe1

¼ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2  2ffi ðx1 Þe3  ðx1 Þe4 þ ðx2 Þe3  ðx2 Þe4 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   2ffi e e 2 ðx1 Þ2  ðx1 Þ3 þ ðx2 Þe2  ðx2 Þe3 2

;

ðJ 1 Þe2

¼

;

ðJ 2 Þe2

¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2  2ffi ðx1 Þe1  ðx1 Þe2 þ ðx2 Þe1  ðx2 Þe2 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   2ffi e e 2 ðx1 Þ1  ðx1 Þ4 þ ðx2 Þe1  ðx2 Þe4 2

; ð78Þ :

In (77), the shape function Nðn; gÞ and strain–displacement matrices Be ðn; gÞ are given by

Nðn; gÞ ¼



N1 ðn; gÞ

0

N2 ðn; gÞ

0

N3 ðn; gÞ

0

N4 ðn; gÞ

0

0

N1 ðn; gÞ

0

N2 ðn; gÞ

0

N3 ðn; gÞ

0

N 4 ðn; gÞ

 ð79Þ

with

1 ð1  nÞð1  gÞ; 4 1 N3 ðn; gÞ ¼ ð1 þ nÞð1 þ gÞ; 4

1 ð1 þ nÞð1  gÞ; 4 1 N4 ðn; gÞ ¼ ð1  nÞð1 þ gÞ 4

N1 ðn; gÞ ¼

N2 ðn; gÞ ¼

ð80Þ

and

2

Be11 ðn; gÞ

0

6 0 Be21 ðn; gÞ 6 Be ðn; gÞ ¼ 6 e 4 B21 ðn; gÞ B11 ðn; gÞe Be12 ðn; gÞ B11 ðn; gÞe

Be12 ðn; gÞ

0

0

Be22 ðn; gÞ

Be22 ðn; Be22 ðn;

gÞ gÞ

Be12 ðn; Be21 ðn;

Be13 ðn; gÞ

0

0

Be23 ðn; gÞ

Be23 ðn; Be32 ðn;

gÞ gÞ

gÞ gÞ

Be13 ðn; Be31 ðn;

Be14 ðn; gÞ

gÞ gÞ

0 Be24 ðn; Be42 ðn;

0

3

Be24 ðn; gÞ 7 7 7: gÞ Be14 ðn; gÞ 5

ð81Þ

gÞ Be41 ðn; gÞ

For the sake of the simplicity, we introduce the following notations

8  i  R 1 R 1 h e  T  e > ^ enþ1 ¼ 1 ^ nþ1 ðn; gÞ J e ðn; gÞ dndg; f int u ðB ðn; gÞÞ r u > 1 > > > > i  e R1 R1 h > ^e J e ðn; gÞ dndg; > < f1ext nþ1 ¼ 1 ðNðn; gÞÞT Nðn; gÞb nþ1 1 i  ext  e R1 h T > ^f e ðJ Þe dn; > f ¼ ð Nðn; 2 ð i  1 Þ þ 1Þ Þ Nðn; 2 ð i  1 Þ þ 1Þ > 1 2 nþ1 i 1 i nþ1 > > > i >    R h > : f ext e ¼ 1 ðNð2ði  1Þ þ 1; gÞÞT Nð2ði  1Þ þ 1; gÞ^f e ðJ Þe dg; 3 nþ1 2 i 1 i nþ1

ð82Þ

Hence Eq. (77) becomes nel n X 

^ enþ1 w

T 

nel n 2 X nel n 2 X nel n  e o X  e T  ext e o X  e T  ext  e o X  e T  ext  e o ^ nþ1 ^ nþ1 f3 nþ1 þ ^ nþ1 ^ nþ1 w w w f int u f1 i nþ1 þ f2 i nþ1 : ¼

e¼1

e¼1

i¼1 e¼1

ð83Þ

i¼1 e¼1

The assembled matrix of internal forces and vectors of external forces are denoted liked in Hughes [12], i.e. nel   e  ^nþ1 Þ ¼ A f int u ^ nþ1 ; F int ðu e¼1





F ext 2 nþ1

nel

¼ A

h

e¼1

 e f2ext 1 nþ1

þ



F ext 1



 nþ1

nel  e ¼ A f1ext nþ1 ;

 e f2ext 2 nþ1

e¼1

i

;



F ext 3

 nþ1

nel

¼ A

e¼1

h      e i e f3ext 1 nþ1 þ f3ext 2 nþ1 :

ð84Þ

We also use the following notation:

F ext nþ1 ¼

3 X 

F ext i

 nþ1

ð85Þ

:

i¼1

Then, the discrete variational formulation (83) becomes

h i ^ nþ1 ; w ^ nþ1 Þ ¼ w ^ Tnþ1 F int ðu ^nþ1 Þ  F ext Gðu nþ1 ¼ 0;

^ 2 W had 8w

ð86Þ

Consequently, the derived system of non-linear Eqs. (86) can be written as (76). ^ nþ1 can be solved using the Newton–Raphson method, namely The non linear system (76) for the unknown u

h i1   j ^ u ^ jþ1 ^j ^ jnþ1 ; u G nþ1 ¼ unþ1  bj K nþ1

ð87Þ

^ is given by where bj 2 ð0; 1 and the Jacobian of the function G

^nþ1 Þi;k ¼ ½K ðu

   Z 1 Z 1 h  ^ i ðu nel i ^ nþ1 Þ nel @ fiint ðuÞenþ1 @G T e ¼ A ¼ A ðBe ðn; gÞÞ ½E ep nþ1 Be ðn; gÞ J e ðn; gÞ dndg : ^ nþ1 Þk ^ nþ1 Þk e¼1 e¼1 @ ðu @ ðu 1 1

ð88Þ

742

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

^ contains all of the displacements provided during the discretization procedure. In In the iterative relations (87), vector u ^ can be expressed using the displacements at the nodes on the order to reduce the computational cost, displacement vector u ^ jCu . Therefor, relation (87) becomes fixed edge u

2  3 2  3 2  3 ^ jþ1 ^ jnþ1 ^ jnþ1 u u GE u nþ1 6 6 7 C 7 C 7 1 6  5;  u 5 ¼ 4  u 5  bj K 4  4 j j ^ ^ jþ1 ^ G u u u F nþ1 nþ1 nþ1 K

" K¼

ð89Þ

K

ðK nþ1 ÞE

ðK nþ1 ÞEF

ðK nþ1 ÞTEF

ðK nþ1 ÞF

# ð90Þ

:

On applying the matrix K, expression (89) becomes

"

ðK nþ1 ÞE ðK nþ1 ÞTEF

2  22   3 2  33 3 # ^jþ1 ^j ^j u G u nþ1 ðK nþ1 ÞEF 66 u 6 E nþ1 7 6 nþ1 Cu 77 Cu 7  5:  5  4  55 ¼ bj 4  44  ðK nþ1 ÞF ^j ^jþ1 ^j GF u u u nþ1

nþ1

K

ð91Þ

nþ1

K

The Newton–Raphson algorithm (87) can be simplified if one takes into account the boundary condition for the displace    ^ jþ1 ^ jnþ1 ment u ¼ u ¼ 0. Then, relation (89) can be written as nþ1 Cu

Cu

h    i   ^ jþ1 ^ jnþ1 ^ jnþ1 ðK nþ1 ÞF u  u ¼ bj GF u nþ1

ð92Þ

    ^ jþ1 ^ jnþ1  bj DF K ¼ u u nþ1

ð93Þ

K

K

or K

K





^j where DF K ¼ ðK nþ1 Þ1 F GF unþ1 . The algorithm to solve the non-linear system (76) can now be applied. h

6. Numerical algorithms We are now able to present the Radial Return Algorithms:

n

o

n

o

^ jþ1 p jþ1 jþ1 ¼ RadialReturn ujnþ1 ; a ^ n ; ½ep n ; kn : rjþ1 nþ1 ; anþ1 ; ½e nþ1 ; knþ1

Step 1. Compute the trial elastic stress for prescribed

str;jþ1 nþ1 ¼ E ð2Þ



eð2Þ

jþ1 nþ1



h

epð2Þ

i  n





eð2Þ



jþ1 nþ1

¼ ½Bujnþ1 :

1 rtr;jþ1 str;jþ1 ; nþ1 ¼ P nþ1

;

where P and E ð2Þ are defined by (19) and (24), respectively. qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi qffiffi   ^ tr;jþ1 ¼ str 2 þ ðtrðstr ÞÞ2  2½H ktr þ rY : Step 2. Check the yield condition F nþ1 nþ1 nþ1 nþ1 3 ^ tr;jþ1 6 0, then: If F nþ1 elastic step: tr;jþ1 ^ tr;jþ1 rjþ1 a^ jþ1 nþ1 ¼ rnþ1 ; nþ1 ¼ anþ1 ;

h

epð2Þ

ijþ1 nþ1

¼

h

epð2Þ

itr;jþ1 nþ1

;

else plastic step: proceed to step 3 endif Step 3. Radial Return Algorithm:

strnþ1 ; strnþ1 ¼ strnþ1  a^ trnþ1 ; njþ1 nþ1 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   fstr 2 þ ðtr ðstr ÞÞ2 nþ1 nþ1

½k nþ1

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffi tr 2 kstrnþ1 k2 þ ðtrðstrnþ1 ÞÞ  23½Hknþ1 þ rY  qffiffi ¼ 2 E H þ 3ð1 3 mÞ þ C

jþ1 tr;jþ1 jþ1 tr;jþ1 sjþ1 knþ1 ¼ knþ1 ; ½E ep nþ1 ¼ E ð2Þ nþ1 ¼ snþ1 ;

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

743

Fig. 1. Deformation problem of the trapezoidal panel.

Fig. 2. The distribution of the stress r11 on the trapezoidal plate at t1 ¼ 14:1 s and t2 ¼ 17:3 s, which correspond to the minimum and maximum of the applied force f and the location of the considered material points.

snþ1 ¼ strnþ1  ½k nþ1 PE ð2Þ njþ1 nþ1 ; h

epð2Þ

i nþ1

¼

h

epð2Þ

itr nþ1

a^ nþ1 ¼ a^ trnþ1 þ ½k nþ1 C ½njþ1 rnþ1 ¼ P1 snþ1 ; nþ1 ;

þ ½k nþ1 njþ1 nþ1 ;

jþ1

tr

knþ1 ¼ knþ1 þ ½k jþ1 nþ1 ;

H1 nþ1 and H2nþ1 are defined by (66);

  ½E ep nþ1 ¼ E ð2Þ H1nþ1  H2nþ1 nnþ1  E ð2Þ P nnþ1 þ n?nþ1 E ð2Þ P Ið2Þ

Return. The proposed algorithm is useful for: 1. the integration of expressions that define the vector of internal forces (82)1. This vector is dependent on the solution rnþ1 which has been given by the Radial Return Algorithm;

744

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

Fig. 3. The distribution of the stress rMises on the trapezoidal plate at t1 ¼ 14:1 s and t 2 ¼ 17:3 s, which correspond to the same minimum and maximum of the applied force f given as in Fig. 2.

Fig. 4. The distribution of the stress r11 for the plate with a hole at t 1 ¼ 14:1 s and t2 ¼ 17:3 s, which correspond to the minimum and maximum of the applied force f given as in Fig. 2. and the location of the material points considered.

^ nþ1 , through the Newton–Rhapson 2. the evaluation of the Jacobian K jnþ1 that is necessary to determine the displacements u method, see formula (92). In order to find the Jacobian K jnþ1 , the elasto-plastic moduli ðE ep Þe are requested. These elasto-plastic moduli have been computed by the Radial Return Algorithm; The vectors of external forces (82)2,3,4 can be obtained using the Simpson integration method, independent of the Radial Return Algorithm. General algorithm to solve Problem P: 0

^ 1 ¼ 0; ½ep 0 ¼ 0; ½a0 ¼ 0; ½k0 ¼ 0. 1. Initialize the vectors: n ¼ 0; ½u 2. For n ¼ 0; N; perform a for loop to compute the time increments with the time step Dt; t0 ¼ 0 ; 2.1. Perform a while loop with respect to the index variable j; j ¼ 0; 1; . . . (i.e. apply the Newton–Raphson method). 2.2. For e ¼ 1; nel perform a for loop with respect to the number of elements of the FEM grid.

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

745

Fig. 5. The distribution of the stress rMises for the plate with a hole at t1 ¼ 14:1 s and t2 ¼ 17:3 s, which correspond to the same minimum and maximum of the applied force f given as in Fig. 2.

Fig. 6. The component of the stress r11 versus the strain 11 at the material points given by nodes 210 and 2301 for the trapezoidal plate and nodes 2619 and 1870 for the plate with hole, respectively, obtained using the in-plane stress solution.

746

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751



^ jnþ1 ; ½ep trnþ1 ; atrnþ1 and ktrnþ1 from the global vectors u ^ jnþ1 ; ½ep n ; and kn . Extract the local vectors u Apply the Radial Return Algorithm. Compute the local Jacobian matrix and the local vectors of the internal and external forces. ^ jþ1 Determine u nþ1 via Newton’s method

^ jþ1 ^j ^j u nþ1 ¼ unþ1  bj ½Dunþ1 ;

h  i1   ^ jnþ1 ¼ K u ^jnþ1 ^ jnþ1 : ½ Du G u

  ^ jnþ1 = ^ jnþ1  Compute the error ½ERRjnþ1 ¼ ½Du u  if ½ERRjnþ1 < TOL

^ nþ2 u

^ jþ1 u nþ1 ;

return else j

jþ1

½ep nþ1 ¼ ½ep nþ1 ;

jþ1

½knþ1 ¼ ½knþ1 ;

½rnþ1 ¼ ½rjþ1 nþ1 ;

½anþ1 ¼ ½ajþ1 nþ1

j þ 1 go to 2.2.

7. Numerical simulation In order to exemplify the proposed numerical algorithm for solving the boundary and initial value Problem P, we consider two examples. The elasto-plastic material is considered to be steel DP600, whose parameters are given by [4]

f E ¼ 182; 000 ½MPa; The yield stress

rY ¼ 349:4 ½MPa; m ¼ 0:3; C ¼ 17; 400 ½MPa; H ¼ 50 ½MPa g

ð94Þ

rY ¼ 349:4 ½MPa for DP600, which is used in the automotive industry, is high.

Fig. 7. The component of the stress r22 versus the strain 22 at the material points given by nodes 210 and 2301 for the trapezoidal plate and nodes 2619 and 1870 for the plate with hole, respectively, obtained using the in-plane stress solution.

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

747

  Example 1. Consider a trapezoidal plate domain X ¼ ðx1 ; x2 Þ 2 R2 j0 6 x2 6 200; x42 6 x1 6 100 see Fig. 1. The vertical left edge is fixed, and the bottom and right vertical edges are stress free, i.e. f ¼ 0. The periodic external force f ¼ 50 sin t; t 2 I, is applied on the top horizontal edge, while the internal forces are assumed to vanish b ¼ 0. Thus in X, for the time interval ½0; TÞ, the equilibrium equations and boundary conditions are to be satisfied

8 div r ¼ 0 > > > < rn ¼ 50; sint e1 > rn ¼ 0 > > : u¼0

in X  I; in Cr1  I;   in Cr2 [ Cr3  I; in Cu  I;

ð95Þ

while the initial conditions are vanishing.

Fig. 8. The graphs of the components of the strain e33 and of the plastic strain ep33 , plotted as functions of time at the material points given by nodes 1225 and 1870 for the plate with a hole, using the Simo solution and the in-plane stress solution.

Fig. 9. The component of the stress r11 versus the strain 11 at the material points given by nodes 207 and 1225 for the plate with hole, obtained using the Simo solution and the in-plane stress solution.

748

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

Example 2. Consider the trapezoidal plate domain given in Example 1 contains a circular hole with the radius of 20 mm and center (60 mm, 80 mm), whose boundary is stress free. We also consider the same problem and type of material as those given by Example 1. The four-node elements were employed and the trapezoidal plate domain was divided into 2500 elements for the two Examples. The numerical simulation was realized for the time interval t 2 ½0; 40 s. To emphasize the differences and facilitate the comparison of the material behaviour of the plate with or without a hole, all numerical results have been obtained using similar conditions. The contour plots of the stress component r11 , which develops in the direction of the applied force, are shown on the plate at time t, corresponding to the minimum and maximum values of the applied force, see Figs. 2 and 4. qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi The distribution of the Mises-stress rMises ¼ r211 þ r222 þ 2r212 (a measure of the Cauchy stress intensity) is represented at the same t in Figs. 3 and 5. The behaviour of the material points remains elastic in a large area starting from a zone near the fixed edge of the plate without a hole and spreading to the opposite free edge since rMises is still less than the initial yield stress hold rY .

Fig. 10. The component of the stress r22 versus the strain 22 at the material points given by nodes 207 and 1225 for the plate with a hole, obtained using the Simo solution and the in-plane stress solution.

Fig. 11. The relative difference of the stress component r11 and the difference of the strain component 11 at the material points given by nodes 210 and 2301 on the trapezoidal plate, obtained using the Simo and the in-plane stress solutions and corresponding to the time variation of the applied force f.

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

749

Fig. 12. The relative difference of the stress component r11 and the difference of the strain component 11 at the material points given by nodes 2619 and 1870 for the plate with a hole, which are plotted using the Simo solution and the in-plane stress solution and corresponding to the time variation of the applied force f shown in Fig. 11.

To emphasize the peculiar elasto-plastic behaviour of the plates at specified material points, we fix two pairs of material points, which coincide with the same geometric points in the plates and are identified by appropriate nodes, see Figs. 2 and 4. One pair of nodes are in the vicinity of the fixed edge and the other one is associated with the material points localized at the nodes close to the hole. Both pairs of material points have been chosen in zones of maximum values of rMises , at a fixed time t. The graphs of the stress component r11 versus the strain e11 and r22 versus 22 are presented at the fixed material points (in the vicinity of the fixed edge) of the two plates in Figs. 6 and 7, respectively. The length of the interval of variation of strain decreases with increasing the level of stress reached when passing from one cycle to another. The normal components of the strains, e33 and ep33 attain larger values (of order 6  103 ) at the upper node near the hole (on the left) than at the node close to the fixed edge as can be seen from Fig. 8. The maximum level of these strains is slowly decreasing during the evolution of the applied force at the node near the hole, see Fig. 8. On the contrary, their corresponding maximum values remain nearly constant at the material point near the fixed edge for both plates. We make a comparison between the numerical solutions of the initial and boundary value problem developed herein by considering an in-plane stress state and Simo solution. The graphs of the stress component r11 versus the strain e11 and r22 versus e22 are plotted, at the fixed material points in the vicinity of the hole, in Figs. 9 and 10, respectively. For the plate without a hole, the elastic solution occurs at these corresponding points. Larger values for the stress r11 and strain e11

Fig. 13. The relative difference of the stress component r11 and the difference of the strain component 11 at the material points given by nodes 207 and 1225 for the plate with a hole, which are plotted using the Simo solution and the in-plane stress solution and corresponding to the time variation of the applied force f shown in Fig. 11.

750

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

are predicted by the in-plane solution than the Simo solution at the material points close to the hole, see Fig. 9. In the direction perpendicular to the applied force, larger values for the stress r22 and smaller values for the strain e22 are predicted by the in-plane solution than the Simo solution, see Fig. 10. r11 Þ The relative difference of the stress component r11 , namely ðrðSÞ11 , and the difference between the strain component rY e11 , i.e ðeðSÞ11  e11 Þ, as function of t or the applied periodic force, are illustrated in Figs. 11 and 12. The differences are calculated near the fixed edge for both plates, with the mention that the smaller values occur for the plate without a hole. We note that the material response is completely different for the points located in the vicinity of the loading boundary (to the right) and close to the free boundary (to the left) see Figs. 11 and 12. In Fig. 13, the relative differences are calculated in the vicinity of the hole and attain a maximum value of 0.12, while the strain difference at the lower point remains small and decreases at the upper point. The numerical results presented in this section have been obtained using Mathcad cod [20]. 8. Conclusions In this paper, the problem of the deformation of an elasto-plastic panel, made up from a mixed hardening material and undergoing an in-plane stress state was reformulated such that the compatibility of the strain state with the in-plane stress is satisfied. We proposed a revised Simo algorithm to solve the bi-dimensional elasto-plastic problem. Only when the plane stress state is considered to be associated with a plane strain state, i.e. the strain in the direction perpendicular to the plane is neglected, the algorithmic elasto-plastic tangent moduli provided here can be reduced to the formulae given by Simo [25], namely the appropriate fourth-order tensors are replaced by two scalars. We make a comparison between the numerical solutions of the initial and boundary value problem developed herein by considering an in-plane stress state and Simo solution. All numerical results have been obtained using similar conditions. Larger values for the stresses are generally predicted by the in-plane solution than the Simo solution. We emphasize the differences of the material behaviour of the plate with or without a hole, with the mention that the smaller values occur for the plate without a hole. The smaller differences between the stresses occur for the material points located near the fixed edge for both plates than for the points located in the vicinity of the hole, the stress values are calculated using the Simo and in-plane solutions, respectively. The numerical algorithm presented herein is applicable to plates with a very small thickness, and this is precisely the case in this paper. We remark that the proposed algorithm is not a particular case of the general standard programme used for solving 3D-problems. Moreover it has been shown that the revisited Simo algorithm is easier and more efficient than the general 3D-algorithm. To illustrate how to integrate the numerical algorithm, two numerical applications have been presented. The efficiency of the integration algorithm was shown and the numerical results obtained have been found to be physically consistent. Acknowledgments The authors acknowledges the support from the Ministry of Education Research and Inovation under CNSIS PN2 Programme Idei, PCCE, Contract No. 100/2009. References [1] P.J. Armstrong, C.O. Frederick, A mathematical representation of the multiaxial Bauschinger effect, G.E.G.B. Report RD-B-N, pp. 731, 1996. [2] J. Bathe, Finite Element Procedure, vols. I–III, Prentice Hall, Englewood Cliffs, NJ, 1996. [3] T. Belytschko, K.L. Wing, B. Moran, Nonlinear Finite Elements for Continua and Structures, British Library Cataloguing in Publication Data, Toronto, 2006. [4] G.B. Broggiato, F. Campana, L. Cortese, The Chaboche nonlinear kinematic hardening model: calibration methodology and validation, Meccanica 43 (2) (2008) 115–124. [5] J.L. Chaboche, Constitutive equation for cyclic plasticity and cyclic viscoplasticity, Int. J. Plast. 5 (1989) 247–302. [6] S. Cleja-T ß igoiu, Anisotropic elasto-plastic model for large metal forming deformation processes, Int. J. Form. Process. 10 (1) (2007) 67–89. [7] L.M. Kachanov, Fundamentals of the Theory of Plasticity, Mir Publishers, Moscow, 1974. [8] A.S. Khan, S. Huang, Continuum theory of plasticity, Fundamentals of the Theory of Plasticity, John Wiley & Sons, Inc., 1995. [9] J. Fish, T. Belytschko, A First Course in Finite Elements, John Wiley and Sons Ltd., USA, 2007. [10] W. Han, B.D. Reddy, Mathematical Theory and Numerical Analysis, Springer, New York, 1999. [11] R. Hill, The Mathematical Theory of Plasticity, Oxford University Press, Oxford, UK, 1950. [12] T.J.R. Hughes, The Finite Element Method. Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Englewood, NJ, 1987. [13] T.J.R. Hughes, R.L. Taylor, Unconditionally stable algorithms for quasi-static elasto/viscoplastic finite element analysis, Comput. Struct. 8 (1978) 169– 173. [14] C. Johnson, On plasticity with hardening, J. Math. Anal. Appl. 62 (1978) 325–336. [15] C. Johnson, Numerical Solution of Partial Differential Equation by the Finite Element Method, Cambridge University Press, Cambridge, 2002. [16] R.D. Krieg, S.W. Key, Implementation of a Time Dependent Plasticity Theory into Structural Computer Programs, in: J.A. Stricklin, K.J. Saczalski (Eds.), Constitutive Equations in Visco-plasticity, Computational and Engineering Aspects, ASME, New York, 1976 (AMD-20). [17] G. Duvaut, J.L. Lions, Les Inéquations en Mécanique et en Physique, Dunot, Paris, 1972. [18] J. Lubliner, A maximum-dissipation principle in generalized plasticity, Acta Mech. 52 (1984) 225–237. [19] J.B. Martin, A New Method of Analyzing Stress and Strain in Work-Hardening Plastic Solid, The MIT Press, Cambridge Mars, 1975. [20] B. Maxfield, Essential Mathcad. For Engineering, Science, and Math, second ed., Academic Press (Elsevier), 2009.

S. Cleja-Tß igoiu, N.E. Stoicutßa / Applied Mathematics and Computation 237 (2014) 730–751

751

[21] M. Ortiz, E.P. Popov, Accuracy and stability of integration algorithms for elastoplastic constitutive equations, Int. J. Numer. Methods Eng. 21 (1985) 1561–1576. [22] W. Prager, Recent development in the mathematical theory of plasticity, J. Appl. Phys 20 (1949) 235–241. [23] W. Prager, A new method of analyzing stress and strain in work-hardening plastic solid, J. Appl. Mech 23 (1956) 493–496. [24] J.C. Simo, R.L. Taylor, Return mapping algorithm for plane stress elastoplasticity, Int. J. Numer. Methods Eng. 22 (1986) 649–670. [25] J.C. Simo, T.J.R. Hughes, Computational Inelasticity, Springer, New-York, 1998. [26] M.L. Wilkins, Calculation of elastic plastic flow, in: B. Alder (Ed.), Methods in Computational Physics, vol. 3, Academic Press, New York, 1964, p. 211.