Exact integration of the von Mises elastoplasticity model with combined linear isotropic-kinematic hardening

International Journal of Plasticity 25 (2009) 1083–1106

Contents lists available at ScienceDirect

International Journal of Plasticity journal homepage: www.elsevier.com/locate/ijplas

Exact integration of the von Mises elastoplasticity model with combined linear isotropic-kinematic hardening Attila Kossa *, László Szabó }egyetem rkp. 5, H-1111 Budapest, Hungary Department of Applied Mechanics, Budapest University of Technology and Economics, Mu

a r t i c l e

i n f o

Article history: Received 10 February 2008 Received in final revised form 4 June 2008 Available online 26 August 2008

Keywords: Exact integration Combined hardening von Mises elastoplasticity Consistent tangent

a b s t r a c t The main purpose of this work is to present two semi-analytical solutions for the von Mises elastoplasticity model governed by combined linear isotropic-kinematic hardening. The first solution (SOLe) corresponds to strain-driven problems with constant strain rate assumption, whereas the second one (SOLr) is proposed for stress-driven problems using constant stress rate assumption. The formulas are derived within the small strain theory Besides the new analytical solutions, a new discretized integration scheme (AMe) based on the time-continuous SOLe is also presented and the corresponding algorithmically consistent tangent tensor is provided. A main advantage of the discretized stress updating algorithm is its accuracy; it renders the exact solution if constant strain rate is assumed during the strain increment, which is a commonly adopted assumption in the standard finite element calculations. The improved accuracy of the new method (AMe) compared with the well-known radial return method (RRM) is demonstrated by evaluating two simple examples characterized by generic nonlinear strain paths. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction Accurate prediction of the nonlinear behavior of real materials is a topic of considerable interest. There are material models for which analytical solutions have been derived, but in most cases of elastoplastic models, the integration of governing constitutive equations can only be performed numerically due to their complex structure. The accuracy of the calculations is strongly dependent on the integration scheme of the constitutive equation. This means that the integration, which is usually performed in a computational environment (such as a finite element code), must be as accurate as * Corresponding author. Fax: +36 1 463 3471. E-mail addresses: [email protected] (A. Kossa), [email protected] (L. Szabó). 0749-6419/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijplas.2008.08.003

1084

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

possible. Therefore, the requirement for analytical solutions of specific material models is obvious. In general, deriving closed-form solutions for elastoplastic constitutive models, including the von Mises model, is quite difficult, because of the highly nonlinear nature of the governing equations. The nonlinear modeling of the hardening rules makes the analytical treatment more complicated or even impossible. For such simple cases as linear hardening models, the analytical solutions are not clarified for all hardening rules, especially not for the combined hardening rule. In this work, we restrict ourselves to linear hardening. Additionally, we provide insight into the analytical studies of nonlinear hardening materials. By far, the most widely used plasticity model is the von Mises elastoplasticity. Depending on the hardening rules, material models with linear hardening can be categorized into four different groups: perfect plasticity (no hardening), linear kinematic hardening, linear isotropic hardening and combined linear hardening (both isotropic and kinematic hardening). Papers that analytically study these models are summarized in the following paragraphs. For the simplest case (perfect plasticity), an analytical solution was presented by Krieg and Krieg (1977). Another analytical solution was also derived by Reuss (1930), where the author solved a system of differential equations of the stress components. Hong and Liu (1997) treated the problem as a two-phase linear system with an on-off switch, and presented an integration method solution. By exact linearization of the stress update procedure, Wei et al. (1996) derived the consistent modulus corresponding to the exact integration formula. For the purely kinematic hardening model, Wang and Chang (1985) proposed an exact formula for the integration of the constitutive equations. Numerical implementation of their method can be found in the work of Szabó and Kovács (1987). The integration scheme presented by Auricchio and Beirão da Veiga (2003) is also an exact solution for purely kinematic hardening. A closed-form solution of purely kinematic hardening and softening is also given in the paper of Yoder and Whirley (1984) in strain space description. To the knowledge of the authors, a fully analytical solution for isotropic hardening cannot be found in the available literature. The solution method presented by Ristinmaa and Tryding (1993) leads to an integral, which cannot be integrated explicitly in general case. Szabó (submitted for publication) overcomes this problem by using the incomplete Beta function in solving the governing equation. This yields a semi-analytical solution for the von Mises elastoplasticity model with linear isotropic hardening. In the paper of Mukherjee and Liu (2003), the authors reformulated the basic governing constitutive equations and discussed different strategies to obtain numerical solutions for these differential equations. For combined hardening, a truncated series solution was presented by Chan (1996). Romashchenko et al. (1999) proposed an analytical solution when the loading is given in the form of multisection polygonal lines in the deviatoric stress space. The solution obtained by Ristinmaa and Tryding (1993) requires numerical integration during the stress update procedure. An approach for this numerical integration can be found in the work of Krieg and Xu (1997). Another exact scheme is presented in the work of Liu (2004), where the author has proposed two numerical schemes to solve the constitutive equation reformulated into an integral formulation. Therefore, the latter schemes also involve numerical solutions in the final derivation. The main objective of this work is to derive a semi-analytical solution (SOLe) for the constitutive equations of the von Mises elastoplasticity model governed by combined linear isotropic-kinematic hardening, using the constant strain rate assumption. The present solution combines the pioneering work of Krieg and Krieg (1977) and the new solution presented by Szabó (submitted for publication). The angle defined between the strain rate tensor and the so-called relative stress deviator is introduced as a variable w(t). If the solution for this angle is known, then the exact expression for the deviatoric stress can be obtained in a linear combination. The governing equations lead to an integral, which defines the variation of the angle w(t). This integral was solved numerically in the work of Ristinmaa and Tryding (1993) and was approximated by Krieg and Xu (1997). However, this integration can be performed, and the solution presented by Szabó (submitted for publication) can be employed here. Finally, we arrive at the formula which implicitly defines the angle w(t) using the incomplete Beta function. An additional semi-analytical solution (SOLr) is also presented for stress-driven problems using constant stress rate assumption. Besides the analytical study of the constitutive equations,

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

1085

the discretized stress update algorithm (AMe) based on the time-continuous SOLe is clarified and the corresponding consistent tangent tensor is also provided. The solution technique proposed by Krieg and Krieg (1977) is also adopted by others in the analysis of other specific material models. For example, Wallin et al. (2008) have proposed an integration scheme for a nonlinear isotropic hardening von Mises model coupled with damage evolution; whereas, Carosio et al. (2000) have employed the technique mentioned above in order to obtain a reference solution for a perfect J2 viscoplasticity model. Loret and Prevost (1986) and Rezaiee-Pajand and Nasirai (2008) have employed the method for Drucker–Prager’s elastoplastic equations. Many papers have been published analyzing numerical integration methods for constitutive equations of elastoplastic solids. Although this paper is mainly concerned with the analytical solution of the governing equations, here we mention some relevant numerical studies. Two of the most widely used schemes are the generalized trapezoidal (GTR) and the generalized midpoint rule (GMR). For the von Mises material model with combined isotropic-kinematic hardening, Ristinmaa (1990) derived the consistent tangent modulus of the GTR and GMR in case of general loading, where the integration path starts from an elastic state and ends in an elastoplastic state. A detailed discussion, and the construction of the consistent tangent modulus of GMR in case of isotropic hardening, is given by Gratacos et al. (1992). Caddemi (1994) has presented an unified treatment of the backward-difference, midpoint, and trapezoidal algorithm for combined hardening. The error involved in backward-difference time integrations of elastoplastic models has been discussed in the paper of Cocchetti and Perego (2003). Auricchio and Beirão da Veiga (2003) proposed a new integration scheme based on the computation of an integration factor for von Mises elastoplasticity model with combined linear hardening. A new exponential based integration algorithm for associative von Mises elastoplasticity model with combined linear isotropic-kinematic hardening has been presented by Artioli et al. (2006). A comprehensive study of four integration methods based on the GMR is given by Artioli et al. (2007). The numerical performance of the GTR is investigated by Yang et al. (2008) by using an advanced soil model. Application of the return map algorithm and the corresponding consistent tangent tensor to nonlinear combined hardening is given by Auricchio and Taylor (1995), where authors have discussed also the generalized plasticity model. Khoei and Jamali (2005) developed a solution method based on the return map algorithm for a multi-surface plasticity model with both isotropic and kinematic hardening. In the work of Wallin et al. (2001), a Runge–Kutta integration scheme is investigated for von Mises materials with isotropic hardening and for von Mises materials with damage evolution coupled to nonlinear mixed hardening. Application of explicit Runge–Kutta methods with error control to general class of elastoplastic models are under consideration in the paper of Hiley and Rouainia (2008). When the problem is considered in finite strain framework, Ponthot (2002) proposed an unified integration algorithm based on the classical radial return method for von Mises materials. Detailed study of the integration of inelastic constitutive models is given in the textbooks of Simo and Hughes (1998), Dunne and Petrinic (2005) and Stein et al. (2004). The paper is organized as follows. The notation for equations will be introduced at the end of this section. Section 2 contains a brief review of the constitutive relations for von Mises elastoplasticity model with combined linear isotropic-kinematic hardening. The new semi-analytical solutions are derived in Section 3, where the solution for problems assuming constant strain rate and constant stress rate will also be presented (SOLe and SOLr). The last part of Section 3 is concerned with the analytical study of nonlinear hardening models. The first part of Section 4 contains the time discretized stress updating algorithm (AMe) corresponding to the SOLe. Various cases depending on the yield condition and the strain increment are discussed. The second part of Section 4 presents the algorithmically consistent tangent tensor. In Section 5, two simple examples are provided in order to illustrate the accuracy of the AMe compared with the well established radial return map algorithm (RRM). Finally, in Appendix A, the detailed derivation of the consistent tangent tensor corresponding to AMe is provided. The following abbreviations are used throughout the paper: SOLe denotes the new semi-analytical solution for strain-driven problems (Section 3.1); SOLr is used for the new semi-analytical solution for stress-driven problems (Section 3.2); AMe stands for the discretized integration scheme (Section 4.1) based on the SOLe; RRM means the radial return method; CT is the computational time and DE denotes differential equation.

1086

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

We employ the following notation. Tensors are denoted by bold-face characters, the order of which is indicated in the text. The tensor product is denoted by , and the following symbolic operations apply: a:b = aijbij and (C:a)ij = Cijklakl with the summation over repeated indices. The prefix tr refers to the pffiffiffiffiffiffiffiffiffiffi trace. The symbol kak ¼ a : a is used to denote a norm of second order tensor a. Second- and fourthorder identity tensors are designated by d and I, respectively. The superposed dot denotes the material time derivative. 2. Framework for small deformation von Mises elastoplasticity model In this section, the constitutive equations of the commonly used von Mises elastoplasticity model with combined linear isotropic-kinematic hardening are briefly reviewed. For a more detailed discussion of constitutive relations see Khan and Huang (1995), Chaboche (2008), Simo and Hughes (1998) or Nemat-Nasser (2004). The rate of change of the strain measure can be decomposed into an elastic and a plastic contribution as

e_ ¼ e_ e þ e_ p :

ð1Þ

In case of linear elasticity, the stress rate tensor is related to the elastic strain rate as

r_ ¼ De : e_ e :

ð2Þ

e

Here, D is the fourth-order elasticity tensor. In the following derivation, linear isotropic elasticity is assumed, so that

De ¼ 2GT þ Kd  d;

ð3Þ

1 d 3

where T ¼ I   d is the fourth-order deviatoric operator tensor, and G and K are the shear and bulk moduli, respectively. The von Mises yield condition for combined isotropic-kinematic hardening is defined by

Fðn; RÞ ¼ knk  RðcÞ 6 0;

ð4Þ

where the function R(c) defines the isotropic hardening behavior in terms of a scalar plastic state variable. The so-called relative stress deviator is defined as

n ¼ s  a;

ð5Þ

1 trrd 3

where s ¼ r  is the deviatoric stress and a represents the back stress describing the shifting of the von Mises yield surface in the deviatoric stress space due to the kinematic hardening. The associated flow rule for the plastic strain rate is given by

e_ p ¼ k_

oF _ n ¼k ; or knk

ð6Þ

where k_ is the plastic multiplier. The linear isotropic and kinematic hardening moduli can be defined as

hiso ¼ Mh;

hkin ¼ ð1  MÞh;

ð7Þ

where h = 2H/3 and H is the constant plastic hardening modulus. The mixed hardening parameter M 2 [0, 1] describes the relation between the isotropic and kinematic hardening parts (for more details see Axelsson and Samuelsson (1979), Chen and Han (1988), Simo and Hughes (1998)). For the case of linear isotropic hardening, the function R(c) takes the form

R ¼ R0 þ chiso ;

ð8Þ

pffiffiffiffiffiffiffiffi where R0 is a material constant related to the initial value of yield stress rY as R0 ¼ 2=3rY . Since c_ ¼ ke_ p k, the plastic flow rule (6) implies that

c¼

Z

t

k_ d~t: 0

ð9Þ

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

1087

The evolution law for the back stress a is defined according to the Prager–Ziegler rule:

a_ ¼ hkin e_ p ¼ hkin k_

n : knk

ð10Þ

The loading/unloading conditions can be expressed in the Kuhn–Tucker form as

k_ P 0;

Fðn; RÞ 6 0;

_ kFðn; RÞ ¼ 0:

ð11Þ

_ using the plastic consistency condition F_ ¼ 0; and the (1)–(9) can be derived: The plastic multiplier k,

2Gn : e_ k_ ¼ : Rð2G þ hÞ

ð12Þ

Finally, combining (2), (6) and (12), the elastoplastic constitutive relations can be written as

r_ ¼ Dep : e_ ;

Dep ¼ De 

4G2 2

R ð2G þ hÞ

n  n;

ð13Þ

where Dep is the so-called elastoplastic or continuum tangent modulus tensor. The constitutive equation defined above can be separated into deviatoric (s) and hydrostatic (p) parts as follows:

_ r_ ¼ s_ þ p; s_ ¼ 2Ge_ 

ð14Þ 4G

2

R2 ð2G þ hÞ

_ ðn : eÞn;

p_ ¼ Ktre_ d;

ð15Þ

where e_ ¼ e_  13 tre_ d is deviatoric strain rate. The rate of the back stress a can be written using (10), (7) and (12):

a_ ¼

2Gð1  MÞh R2 ð2G þ hÞ

_ ðn : eÞn:

ð16Þ

The evolution law for the radius of the yield surface can be expressed by combining (7)–(9) and (12):

R_ ¼

2GMh _ ðn : eÞ: Rð2G þ hÞ

ð17Þ

Eqs. (15)–(17) are the governing equations of von Mises elastoplasticity with combined linear isotropic-kinematic hardening. The expression for the rate of the relative stress deviator can also be expressed by using (15) and (16) in (5):

  2G Mh _ n_ ¼ 2Ge_  2 1  ðn : eÞn: 2G þ h R

ð18Þ

3. Integration of the constitutive equations of the von Mises elastoplasticity model with combined linear isotropic-kinematic hardening 3.1. Time integration of the constitutive equations with constant strain rate assumption (SOLe) The main goal of the time integration procedure is to solve the system of DEs defined by (15)–(17) for a given strain input. In the following analysis, the constant strain rate assumption will be used. We will restrict the analysis to purely plastic loading, i.e., where both initial and final stresses lie on the yield surface. The variable t represents the time measured from the previous equilibrium state denoted by subscript n. Note that constitutive equations defined by (15)–(17) contain the scalar variable _ This implies the introduction of the angle w through the following inner product: ðn : eÞ.

_ cos w ¼ Rkek _ cos w; n : e_  n : e_ ¼ knkkek

ð19Þ

where w is the angle defined between the strain rate e_ and the relative stress deviator n in the deviatoric plane. Therefore, w is restricted to be between 0 6 wn 6 p, where the range 0 6 wn 6 p/2

1088

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

represents plastic loading, while p/2 < wn 6 p yields elastic unloading. Introducing the scalar parameter w allows one to express the final solution in a linear combination with the help of the law of sines. This was first introduced by Krieg and Krieg (1977). Fig. 1 shows a schematic illustration for the angle w in the deviatoric principal stress plane. After some straightforward algebraic manipulation the constitutive relations (15)–(17) can be converted into a system of nonlinear ordinary DEs:

_ sin w 2Gkek w_ ¼  ; R _R ¼ 2Gakek _ cos w;

ð20Þ

where the dimensionless parameter a is expressed by

a¼

Mc ; 1þc

c¼

h : 2G

ð21Þ

Initial values are defined by R(t = 0) = Rn and w(t = 0) = wn. The procedure for solving the system of differential equations given by (20) is the following. Step 1. The parameter R, representing the radius of the yield surface in the deviatoric plane, can be derived from (20) as a function of w:

Z w 1 e 1 ~ d R ¼ a dw; ~ e Rn R wn tan w  a sin wn ; RðwÞ ¼ Rn sin w Z

R

ð22Þ ð23Þ

which was also presented by Ristinmaa and Tryding (1993) and Krieg and Xu (1997), respectively. Step 2. Substituting the time derivative of (23) into (20), one obtains

_ 2Gkek dw ¼ ðsin wÞð1þaÞ : dt Rn ðsin wn Þa

ð24Þ

Integrating Eq. (24) within the time interval [0, t] we have

_ 2Gkek Rn ðsin wn Þa

Z

t

d~t ¼  0

Z

w

~ ð1þaÞ dw ~ ðsin wÞ

wn

Fig. 1. Schematic illustration of the angle w in the deviatoric principal stress plane.

ð25Þ

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

1089

with the constraint

0 6 wðtÞ 6

p

ð26Þ

:

2

Note that constraint (26) implies that plastic loading is considered. Evaluating the integrals of (25) yields the following expression, which implicitly contains the w(t) function:

    _ 4Gkekt 1 a 1 a 2  B cos2 wn ; ;  : a ¼ B cos w; ;  2 2 2 2 Rn ðsin wn Þ

ð27Þ

Here, the incomplete Beta function B(x, m, l) (see Abramowitz and Stegun (1968); Spanier and Oldham (1987)) is defined as

Bðx; m; lÞ ¼

Z

x

sm1 ð1  sÞl1 ds;

0 6 x 6 1:

ð28Þ

0

Further numerical calculations require that (27) be solved for w(t). An efficient method can be found in the paper of Dominici (2003), where the author derived a series expansion for several inverse functions, including the inverse incomplete Beta function. This approach is very accurate and does not require iteration procedure such as the classical Newton scheme. As mentioned earlier, in the paper of Ristinmaa and Tryding (1993), the right-hand side of Eq. (25) was numerically solved, while Krieg and Xu (1997) proposed an approximation for the function w(t). Step 3. In order to compute s(t), the solution for the function n(t) must be determined first. Using _ the law of sines, the relative stress deviator can be expressed as a linear combination of nn and e:

_ nðtÞ ¼ An nn þ Bn e;

ð29Þ

where the constant coefficients An and Bn have the forms:

An ¼

 a1 sin wn ; sin w

Bn ¼

 a Rn sinðwn  wÞ sin wn : _ sin wn kek sin w

ð30Þ

Step 4. Substituting the solution (29) into (15) and using (19), one obtains the following ordinary differential equation for the deviatoric stress s(t):

_ s_ ðtÞ ¼ A_ s nn þ B_ s e;

ð31Þ

where

A_ s ¼ 

_ sin 2w Gkek ; Rn ð1 þ cÞ sin wn

  cos w sin ðwn  wÞ : B_ s ¼ 2G 1  ð1 þ cÞ sin wn

ð32Þ

The integration of (31) cannot be performed without having the w(t) function either in implicit or explicit form. The new implicit expression (27) is adopted here and allows one to explicitly integrate (31), which leads to the following solution for s(t):

_ sðtÞ ¼ sn þ As nn þ Bs e;

ð33Þ

where the constants As and Bs take the following forms:

As ¼ Bs ¼

1 ð1  An Þ; ða  1Þðc þ 1Þ

     a 2Gc Rn cos wn Rn sin wn 1 a 1 a t  B cos2 wn ; ; 1  : As þ B cos2 w; ; 1  _ _ ð1 þ cÞ 2 2 2 2 kek 2ð1 þ cÞkek ð34Þ

These parameters cannot be expressed explicitly without using the new solution (27). Again, under the proposed technique, the solutions (23), (29) and (33) of the constitutive relations (15)–(17) with the constant strain rate assumption, are expressed with the help of the variable w(t), defined by (19). The first two steps in the numerical calculations are to determine wn and w(t), where the former can be computed from (19) using the variables corresponding to the nth state. For the latter, an incomplete

1090

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

beta function must be inverted, for which, an efficient solution technique is proposed by Dominici (2003).It should be noted that the new solution has singularity at wn = 0 (proportional loading). However, the analytical solution for proportional loading can be easily derived using the fact that e_ and n are coaxial. Finally, the hydrostatic part can be integrated directly:

pðtÞ ¼ pn þ tKtre_ d:

ð35Þ

3.2. Time integration of the constitutive equations with constant stress rate assumption (SOLr) In the previous section, strain-driven problems using constant strain rate assumption were considered. In this section, stress-driven problems are discussed and the integration of the constitutive equations is performed using constant stress rate assumption. The solution for this case can also be derived. However, in order to obtain this, the inverse of (15) is required. It has the form:

e_ ¼

1 1 s_ þ 2 ðn : s_ Þn: 2G R h

ð36Þ

Substituting (36) into (18), we have:

ð1  MÞ n_ ¼ s_  ðn : s_ Þn: R2

ð37Þ

Let x be the angle between the deviatoric stress rate and the relative stress. Then, the inner product ðn : s_ Þ can be written as

n : s_ ¼ Rks_ k cos x:

ð38Þ

Next, combining the time derivative of (38), the plastic consistency condition RR_ ¼ n : n_ and Eq. (37) yields the following system of nonlinear ordinary DEs:

ks_ k sin x ; R R_ ¼ ks_ kM cos x:

x_ ¼ 

ð39Þ

The initial values are defined by R(t = 0) = Rn and x(t = 0) = xn. For plastic loading the angle xn is restricted to be between 0 6 xn 6 p/2 and elastic unloading occurs if p/2 < xn 6 p. The system of DEs (39) can be solved in a similar fashion as was shown in the previous section. Step 1. The R(x) function has the form:

RðxÞ ¼ Rn



sin xn sin x

M ð40Þ

:

Step 2. Differentiating (40) with respect to time and combining it with (39), the following nonlinear equation is obtained:

    1 M 1 M 2 2 ¼ B cos x ; x ; ;   B cos ;  : n 2 2 2 2 Rn ðsin xn ÞM 2ks_ kt

ð41Þ

This expression implicitly defines the function x(t). Step 3. Using the law of sines in the deviatoric plane, the function n(t) can be expressed as a linear combination of nn and s_ :

nðtÞ ¼ An nn þ Bn s_ :

ð42Þ

Note that constants in this expression differ from those in the previous section:

An ¼

 M1 sin xn ; sin x

Bn ¼

 M Rn sinðxn  xÞ sin xn : ks_ k sin xn sin x

ð43Þ

_ Step 4. Substituting (42) into (36) one obtains the following for e:

e_ ¼ A_ e nn þ B_ e s_ ;

ð44Þ

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

1091

where

ks_ k sin 2x A_ e ¼ ; GcRn sin xn

  1 sinðxn  xÞ cos x : B_ e ¼ 1þ 2G c sin xn

ð45Þ

Integrating (44) with respect to time yields to the solution:

eðtÞ ¼ en þ Ae nn þ Be s_ ;

ð46Þ

where the parameters Ae and Be are

Ae ¼ Be ¼

1 ðAn  1Þ; 2GcðM  1Þ

     M ð1 þ cÞ Rn cos xn Rn sin xn 1 M 1 M t  B cos2 xn ; ; 1  : Ae  B cos2 x; ; 1  2Gc 2 2 2 2 ks_ k 4Gcks_ k ð47Þ

Solutions derived for stress-driven problems with the constant stress rate assumption have forms which are similar in structure to those presented for the constant strain rate loading case in the previous subsection. The main difference lies in the definition of the angles, defined by (19) and (38), respectively. Here, also, it should be noted that the solution has singularity at xn = 0 (proportional loading). Finally, the solution for the volume strain can be calculated directly by

treðtÞ ¼ tren þ t

1 _ trr: 3K

ð48Þ

3.3. Partial results for nonlinear hardening models Attention is now focused on the study of the von Mises elastoplasticity model governed by nonlinear hardening. The difficulties appearing in the derivation of the analytical solution will be investigated here. The problem is considered in the strain-driven description. In order to derive a possible analytical solution for combined nonlinear hardening models, the isotropic solution should be clarified first. Therefore, we restrict ourselves to purely isotropic nonlinear hardening. As will be seen, the derivation of analytical solution is difficult, if not impossible, even for this case. Consequently, nonlinear combined hardening models are not discussed here. The two most commonly used nonlinear isotropic hardening rules under review are the power law hardening (PLH) and the exponential hardening (EXH):

RðcÞ ¼



R0 þ bcm ;

PLH;

R0 þ R1 ð1  expðmcÞÞ; EXH;

ð49Þ

with the tangents of the hardening functions expressed in terms of R:

8  11=m > R  R0 oR < mb ; PLH; hðRÞ ¼ ¼ b oc > : EXH: mðR0 þ R1  RÞ;

ð50Þ

Here, b, m and R1 are material parameters. In Section 3.1, the solution was derived for combined linear hardening, where the first key step was the reduction of the original constitutive equations defined by (15)–(17) to the system of DEs through the angle w. After that, the next task was to solve (20) for R(t) and w(t), respectively. Therefore, when applying this solution technique to material models with nonlinear hardening rule, first the construction of the system of DEs for R(t) and w(t) is needed. Since purely isotropic hardening is considered, the _ The resulting system of DEs for angle w is defined between the deviatoric stress s and the strain rate e. R(t) and w(t) takes the following form: _ w_ ¼  2GkekR sin w ; _ cos w: R_ ¼ 2Gh kek

ð51Þ

2Gþh

The relation between the radius of the yield surface (R) and the angle w can be derived from (51). For the models considered here, these equations have the form:

1092

where

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

sin w ¼ gðRÞ;

ð52Þ

8 Rn > > PLH; sin wn expðKÞ; > < R gðRÞ ¼  mðR 2GþR Þ > 0 1 Rn RðR  R0  R1 Þ > > ; EXH: sin wn : Rn ðRn  R0  R1 Þ R

ð53Þ

with

K¼

    1   2G R0 m R0 1 1 R0 1 1 B ;1  ; B : ;1 ; mR0 b m m m m R Rn

ð54Þ

Substituting (53) into the second equation of (51) and separating the variables t and R, we have the following expression (assuming constant strain-rate):

Z

R

Rn

e 2G þ hð RÞ e ¼ kekðt _ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d R  t n Þ; e e 2 2Ghð RÞ 1  g ð RÞ

ð55Þ

The integral in (55) is so highly nonlinear that the integration cannot be performed. Therefore, the R(t) and w(t) functions cannot be derived either in explicit, nor in implicit form. This implies that the solution for s(t) cannot be obtained in this analytical manner. However, analytical solutions are derived for the relation between w and R. The w(R) function can be expressed from (52) in closed-form whereas, the R(w) relation only given in implicit form. When purely nonlinear kinematical hardening is considered and the evaluation of the back stress is described, for example, with the Armstrong–Frederick type model (Armstrong and Frederick, 1966; Chaboche, 2008) then the above analytical technique fails, because the constitutive relations cannot be reduced to a (51)-type system of DEs for R(t) and w(t), respectively. The discussion above points out the limitations of the solution technique considered in this paper. However, the new scheme can be applied for nonlinear hardening materials if we approximate the nonlinear characteristics of the hardening rule with piecewise linear segments. 4. Stress updating method and the construction of the consistent tangent modulus 4.1. Stress updating algorithm (AMe) In the previous section, the exact solution of the governing equations was derived in a time-continuous form, meaning that all quantities were expressed in terms of time. In the standard terminology of finite element implementation, the incremental form is required. In this case, the independent variable is commonly the strain increment defined by De ¼ Dt e_ , where Dt = tn+1  tn. For this reason, the solution SOLe presented in Section 3 needs to be converted into discretized form (AMe). Input parameters are: the previous state variables (rn, en, an, Rn), the material constants (G, h, M) and the strain increment (De). Trial elastic stresses are calculated as

atrial nþ1 ¼ an ;

strial nþ1 ¼ sn þ 2GDe;

ntrial nþ1 ¼ nn þ 2GDe:

ð56Þ

During the increment, several different scenarios can occur (see Fig. 2). These cases are discussed in detail in the following paragraphs. The nth state can take place either inside the yield surface or on the yield surface. First, consider the case when nth state is in the elastic domain, i.e., inside the yield surface. Depending on the direction and magnitude of the strain increment, two scenarios are possible. Either the trial state remains in the elastic region (Case A), i.e., the increment is elastic, or goes outside of the yield surface (Case B), i.e., first elastic and then elastoplastic loading will occur. A special case is when the trial state lies on the yield surface. This situation corresponds to Case A. If the nth state lies on the yield surface, then Case B*, Case C or Case D can occur. Case C represents full elastoplastic increment; Case D produces elastic unloading; and Case B* involves first elastic unloading and then purely elastoplastic loading.

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

1093

Fig. 2. Geometrical interpretation of various cases in the stress updating procedure.

Case A: Elastic increment (elastic to elastic). When the relations

Fðnn ; Rn Þ < 0

and

Fðntrial nþ1 ; Rn Þ 6 0

ð57Þ

are valid, the increment is elastic and the new stress state at time tn+1 is calculated by

snþ1 ¼ sn þ 2GDe;

ð58Þ

anþ1 ¼ an ;

ð59Þ

nnþ1 ¼ snþ1  anþ1 ¼ nn þ 2GDe:

ð60Þ

Case B and Case B*: Elastic–elastoplastic transition (elastic prior to elastoplastic). There are two possible cases resulting from this type of increment. The first (Case B) is when the state at time tn is elastic and the trial state is outside of the yield surface, i.e.,

Fðnn ; Rn Þ < 0

and

  F ntrial nþ1 ; Rn > 0:

ð61Þ

The second (Case B*) is when the nth state is on the yield surface, the trial state is outside the yield surface and the angle wn is larger than p/2:

Fðnn ; Rn Þ ¼ 0;

  F ntrial nþ1 ; Rn > 0

and

wn >

p 2

:

ð62Þ

In this latter case, elastic unloading comes prior to elastoplastic loading. For this type of loading the increment is separated into two parts: first elastic loading and then elastoplastic loading. The contact state at which the stress just reaches the yield surface is

sk ¼ sn þ 2GkDe;

ak ¼ an ;

nk ¼ nn þ 2GkDe;

ð63Þ

and the strain increment, which involves elastoplastic loading is calculated by

D~e ¼ ð1  kÞDe:

ð64Þ

From the yield condition nk : nk ¼

k¼

R2n ,

the parameter k can be calculated as

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi nn : De þ ðnn : DeÞ2 þ kDek2 ðR2n  nn : nn Þ 2GkDek2

:

ð65Þ

e is purely elastoplastic, i.e., the solutions derived in From the contact point, the loading given by D~ e the following nonlinear equation Section 3.1 can be employed here. Based on stress state nk and D~ determines the angle wn+1:

    4GkD~ek 1 a 1 a 2 2 ¼ B cos w ; w ; ;   B cos ;  ; nþ1 k a 2 2 2 2 Rn sin wk

ð66Þ

1094

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

where the initial value of w is computed from:

cos wk ¼

nk : De~ n : De :  k Rn kD~ek Rn kDek

ð67Þ

The stress nn+1 can be expressed using the solution (29):

nnþ1 ¼ An nk þ Bn De;

ð68Þ

where the parameters An and Bn are

An ¼



sin wk sin wnþ1

a1

Bn ¼

;

  a Rn sin wk  wnþ1 sin wk : sin wnþ1 kDek sin wk

ð69Þ

Using (33), the solution for the deviatoric stress sn+1 has the form:

snþ1 ¼ sk þ As nk þ Bs De;

ð70Þ

with the following parameters:

As ¼

1 ð1  An Þ; ða  1Þðc þ 1Þ

ð71Þ       a 2Gcð1  kÞ Rn cos wk Rn sin wk 1 a 1 a  B cos2 wk ; ; 1  : Bs ¼  As þ B cos2 wnþ1 ; ; 1  ð1 þ cÞ 2 2 2 2 kDek 2ð1 þ cÞkDek ð72Þ

The radius of the yield surface at the end of the strain increment according to (23) is given by



Rnþ1 ¼ Rn

sin wk sin wnþ1

a

ð73Þ

:

Case C: Elastoplastic increment (elastoplastic to elastoplastic). When the stress state at time tn is on the yield surface, the trial state is outside and 0 6 w 6 p/2, then:

Fðnn ; Rn Þ ¼ 0

and

  F ntrial nþ1 ; Rn > 0:

ð74Þ

In this case, the new stress state also becomes elastoplastic. The stress updating procedure can be simply derived from Case B by substituting k = 0, i.e., variables corresponding to nth state are used instead of contact point values. Case D: Unloading (elastoplastic to elastic). When the stress state at the time tn is elastoplastic and the trial stress takes place inside the yield surface, i.e.,

  Fðnn ; Rn Þ ¼ 0 and F ntrial nþ1 ; Rn < 0;

ð75Þ

then the new stress state at time tn+1 can be calculated according to Case A. Finally, addition of the deviatoric and hydrostatic stresses completes the overall stress updating procedure:

1 rnþ1 ¼ snþ1 þ pnþ1 ¼ snþ1 þ trrnþ1 d ¼ snþ1 þ Ktrenþ1 d: 3

ð76Þ

4.2. Consistent tangent modulus In general, the first step in elastoplastic finite element analysis involves the stress updating procedure after which the global nonlinear system of algebraic equations is solved. The structure of the tangent tensor has a strong influence on the convergence speed of the calculations. For finite element implementations, where usually the Newton–Raphson method is employed, the construction of the consistent tangent modulus tensor is crucial. This yields faster convergence of the numerical scheme, compared to the case when continuum tangent modulus is used. The consistent tangent modulus

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

1095

gives a quadratic rate of asymptotic convergence for the iterative solution (Simo and Hughes, 1998). This, so-called algorithmically consistent tangent modulus, can be obtained by exact linearization of the stress update procedure AMe, which was described in the previous section. The definition for the consistent tangent modulus is based on

ornþ1 : denþ1 ¼ Dcons : denþ1 ; oenþ1 ornþ1 osnþ1 osnþ1 ¼ ¼ þ Kd  d ¼ : T þ Kd  d: oenþ1 oenþ1 o De

drnþ1 ¼

ð77Þ

Dcons

ð78Þ

Here, we only present the consistent tangent modulus for the elastic–elastoplastic transition (Case B, Case B*), because the expression for purely elastoplastic loading (Case C) can be derived by substituting k = 0 and by replacing the contact point variables with those corresponding to the nth state. The derivative of the deviatoric stress with respect to the deviatoric strain increment De at state n+1 is

osnþ1 osk oAs on oBs ¼ þ nk  þ A s k þ De  þ Bs I: oDe oDe oDe o De oDe

ð79Þ

Substituting (79) into (78), the consistent tangent modulus can be expressed as

Dcons ¼ a1 nk  nk þ a2 nk  De þ a3 De  nk þ a4 De  De þ ð2Gkð1 þ As Þ þ Bs ÞT þ Kd  d;

ð80Þ

where scalar parameters a1, a2, a3 and a4 come from the derivatives of As and Bs with respect to De. Since the following parameters were derived by a lengthy, otherwise straightforward algebraic manipulation, here we only present the final results. Detailed calculations are summarized in Appendix A.

  An C3 C1 ;  1 þ c tan wnþ1 tan wk   An C4 C2 a2 ¼  ; 1 þ c tan wnþ1 tan wk

a1 ¼

ð81Þ ð82Þ

a1 C 2 aAs C 22 kDek As C 1 Rn sin wk Rn sin wk ðC 1  An C 3 Þ þ þ þ C 1 Rn sin wk ð1 þ cÞkDek C1 kDek   3 2 aC 2 kDekð2Gc  Bs ð1 þ cÞÞ 2Gc C 1 Rn sin wk þ aC 2 kDek 2Gkð1 þ As Þ þ  k ; 2 ð1 þ cÞRn sin wk cos wk Rn kDek ð1 þ cÞC 2 kDek Rn sin wk

ð83Þ

a2 C 2 aAs C 32 kDek As C 2 Rn sin wk Rn sin wk ðC 2  An C 4 Þ þ 2 þ þ ð1 þ cÞkDek C1 kDek C 1 Rn sin wk   3 2 aC 22 kDekð2Gc  Bs ð1 þ cÞÞ 2Gc C 1 Rn sin wk þ aC 2 kDek 2Gc  Bs ð1 þ cÞ  kþ ; þ ð1 þ cÞC 1 Rn sin wk ð1 þ cÞC 1 kDek2 Rn sin wk ð1 þ cÞkDek2

ð84Þ

a3 ¼

a4 ¼

where further parameters are introduced:

C1 ¼ 

cos wk Rn  2GkkDek

; sin wk cos wk R2n kDek cos wk Rn  2GkkDek C2 ¼ ; sin wk Rn kDek2  An  C3 ¼ 2 C 1 R2n þ 2GðaC 1 kDekRn cos wk ð1  kÞ  k tan wk Þ ; Rnþ1   An R n  C 2 kDekRn þ 2Gð1  kÞ aC 2 kDek2 cos wk  sin wk : C4 ¼ 2 kDekRnþ1

ð85Þ ð86Þ ð87Þ ð88Þ

1096

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

5. Examples The validity, the accuracy and the computational cost of the new method AMe based on SOLe is investigated through two simple numerical test examples, comparing the new results with those obtained using the classical RRM. Note that if a prescribed rectilinear strain path is applied then the new method always gives the exact solution in every linear part. Therefore, a nonlinear strain path is needed in order to analyze the efficiency. All numerical calculations presented herein were performed using the symbolic software MATHEMATICA (Wolfram Research, Inc., 2005). In the first example, the aim is to compute the corresponding stress response in the case when a nonlinear strain path is given. Analytical solutions for general nonlinear strain input do not exist. Consequently, the problem is solved by approximating the strain path with a rectilinear polygon. The second example is performed through a displacement-based finite element analysis. One 3-D hexahedron element is considered under homogenous stress loading defined by combined tension and shear loading. Errors are analyzed at the end of the loading path. This example is a stress-driven problem with a constant stress rate loading. Therefore, the solution SOLr can be employed here, in order to directly obtain the exact solution in one step. 5.1. Example 1: prescribed nonlinear strain input 5.1.1. The problem and the exact solution In order to apply the new method AMe for a problem described with a nonlinear strain path, an approximation of the whole strain path with a multisection polygonal line is needed. Therefore, the proposed technique AMe would yield a solution which can only be considered to be an approximation, where the nonlinear strain input is approximated by superposition of linear segments. The accuracy of the solution depends upon the degree of the approximation. Let the material parameters be the following:

Young modulus : E ¼ 200 ½GPa; Poisson-ratio : m ¼ 0:3; hardening parameter : H ¼ 40 ½GPa;

ð89Þ

mixed hardening parameter : M ¼ 0:5; initial yield stress :

rY ¼ 150 ½MPa:

The strain path is given by the prescribed e11 and e12 components of the strain tensor:

e11 ¼ 0:003ð1  cos tÞ; e12

1 ¼ c12 ¼ 0:003 sin t: 2

ð90Þ ð91Þ

Fig. 3. (a) Prescribed nonlinear strain path in Example 1. (b) Rectilinear approximations of the prescribed nonlinear strain path.

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

1097

Since the constitutive model is rate-independent, the variable t denotes a dimensionless parameter. The interval is considered: 0 6 t 6 3p/2. Fig. 3a displays the whole strain path in the (e11, e12) plane. At the beginning, elastic deformation occurs until the initial yield surface is reached. This point can be calculated from the initial yield condition:

rffiffiffi 2 ksðtÞk ¼ rY ) tel ¼ 0:1884739447: 3

ð92Þ

Elastoplastic loading results in the interval tel < t 6 3p/2. The reference solution in the elastoplastic region was obtained by numerically (fourth-order Runge–Kutta method with 1000 steps) solving the corresponding system of nonlinear DEs described by (15), (17) and (18). These results are collected in Fig. 4. Fig. 4a shows the angle defined by (19) versus the parameter t, Fig. 4b illustrates the variation of the radius of the yield surface, while the non-zero components of the back-stress a and the stress r are displayed in Figs. 4c and d, respectively. 5.1.2. Numerical calculations For calculations using either the new method (AMe) or the RRM, the strain path given by (90) and (91) needs to be divided into increments. Fig. 3b shows two linear approximations of the originally nonlinear strain path. The errors of both methods at a given time are reported as the relative norm of the error between the exact and the computed solution:

kR  Rexact j ½%; Rexact kr  rexact k ¼ 100 ½%; krexact k

EkRk ¼ 100 Ekrk

Ekak ¼ 100

ka  aexact k ½%; kaexact k

ð93Þ ð94Þ

where R, a and r represents the results obtained by applying the specific algorithm, whereas Rexact, aexact and rexact denotes the exact solutions. Errors resulting in the approximations (5 and 20 steps, respectively) are illustrated in Fig. 5. Results obtained using the RRM are indicated with squares; whereas, circles denote the solutions computed with the new method AMe. Filled markers correspond

Fig. 4. Exact solutions in Example 1. (a) Variation of the angle w. (b) Radius of the yield surface. (c) Components of the back stress. (d) Components of the stress.

1098

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

Fig. 5. Numerical results for the AMe and the RRM in Example 1. (a) Comparison of errors in the radius of the yield surface. (b) Comparisons of errors in the back stress. (c) Comparison of errors in the stress tensor.

Table 1 Comparison of the RRM and the AMe in Example 1 Number of segments

Tcomp

5 6 7 8 9 10 15 20 25 30 50

AMe

RRM

1.00 1.19 1.38 1.56 1.75 2.00 2.94 3.94 5.00 6.12 10.56

Error values [%] Ekrk

Ekak

ER

8.27 6.87 5.86 5.10 4.52 4.06 2.67 1.99 1.59 1.32 0.79

11.10 9.22 7.87 6.86 6.07 5.45 3.59 2.68 2.13 1.77 1.06

6.77 5.67 4.91 4.32 3.86 3.49 2.36 1.78 1.43 1.19 0.72

Tcomp

Error values [%] Ekrk

Ekak

ER

3.88 4.57 5.27 6.14 6.89 7.64 11.46 15.45 19.50 23.50 38.83

3.48 2.47 1.84 1.42 1.13 0.92 0.41 0.23 0.15 0.10 0.04

4.68 3.32 2.47 1.90 1.51 1.23 0.55 0.31 0.20 0.14 0.05

3.27 2.31 1.72 1.32 1.05 0.86 0.38 0.22 0.14 0.10 0.04

to the 20, and empty markers to the five step case, respectively. Furthermore, Table 1 shows the required computational time (CT), besides the considered errors, for both methods. Here, Tcomp denotes a dimensionless relative time measure. Tcomp = 1 represents the required CT for the RRM, when the strain path is approximated with five steps. The higher CT of AMe is explained with the more complex structure of the new method. However, in the calculation of the computational cost, other factors should be considered, besides the CT. For example: how many intermediate steps must be stored in order to obtain the required accuracy? In this simple example, using almost the same amount of CPU time, the new method AMe with 8 steps renders roughly the same accuracy as the RRM can yield with only 30 steps.

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

1099

5.2. Example 2: prescribed rectilinear stress loading 5.2.1. The problem and the exact solution The previous example was a strain-driven problem with a prescribed nonlinear strain path. In this second example, the loading is given in the stress space, or more precisely, in the (r = r11, s = r12) plane (see Fig. 6). Since only the r11 and r12 components are controlled, we have a plane stress problem. The whole loading path consists of three different loading sections. Starting from the initial stress-free state (point O), first the initial yield surface is reached (point Y) by uniaxial tension involving only pure elastic deformation. Then, the uniaxial tension is further applied until the point A is reached. This second loading part (Y ? A) represents a proportional elastoplastic loading. Beyond the point A, an additional stress component, namely, s = r12 begins to increase from sB = 0 [MPa] to sB = 300 [MPa]; whereas, r = r11 is kept to be constant (rB = rA = 300 [MPa]). The last loading section (A ? B) is a non-proportional elastoplastic loading. Material parameters are the same as in the first example. The attention is focused on the non-proportional loading part (A ? B). Although, the loading path is linear in the stress space (consequently the stress rate is constant), the components of the corresponding strain are nonlinear; therefore, it cannot be characterized by constant strain rate. The main goal during this analysis is to determine the total strain e, the back-stress a and the radius of the yield surface R at point B. Two solution methods are considered and compared, similar to the previous example. The first one is the RRM, and the second is the new AMe. Note that, if the loadings were characterized by constant strain rate (which is equivalent to a linear loading path in the strain space), then AMe would provide the exact solution using only one loading step. However, in this case we do not have constant strain rate. Therefore, the accuracy of each variable is strongly dependent on the number of loading steps applied in section A ? B. The exact values at point Y are determined according to the elastic Hooke law. Since the loading Y ? A is proportional, the exact values at point A can be calculated simply by using the exact solution corresponding to a proportional loading. At the end of the last loading section (A ? B), the exact solution is determined by the solution SOLr presented in Section 3.1. This can be done, because the loading part A ? B is a linear path in the stress space; therefore, we have constant stress rate. The exact values of e, a and R are summarized in Table 2. Furthermore, the shape and the position of the yield curve in the (r, s) plane at loading points Y, A and B are illustrated in Fig. 6, where the contributions of the exact solution of the total strain components in terms of the applied s stress are also presented. 5.2.2. Numerical calculations In order to obtain information about the performance of the consistent tangent tensor, besides the accuracy of the corresponding stress update algorithm, this example is performed in a standard displacement-based finite element analysis created in MATHEMATICA. One 3-D linear (8-noded) hexahedron element (with 100 [mm] edge dimension) is considered with 2-point Gaussian quadrature.

Fig. 6. Loading path in Example 2, and the corresponding exact solution of the strain components and yield curve.

1100

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

Table 2 Exact values in Example 1 at points Y, A and B, respectively Y

A

B

e11 e22 = e33 e12 = e21

0.75  103 0.225  103 0

5.25  103 2.325  103 0

10.194  103 4.797  103 7.948  103

a11 a22 = a33 a12 = a21

0 0 0

50 [MPa] 25 [MPa] 0

115.924 [MPa] 57.962 [MPa] 79.973 [MPa]

R

122.474 [MPa]

183.712 [MPa]

327.761 [MPa]

Because the loading Y ? A is proportional, the values computed at point A are unaffected by the number of steps applied in section Y ? A. Therefore, even the RRM renders the exact solution using one loading step. As mentioned earlier, the new method cannot be applicable for proportional loading. However, it cannot be deemed as a significant drawback, because for proportional loading, the corresponding exact expression should be used in numerical implementation, in order to save computational costs. The numerical analysis is focused on the non-proportional loading path A ? B, where an iterative solution procedure is involved to obtain the resulting displacement values. Two issues are discussed: (a) The performance of the consistent tangent tensor compared to the continuum tangent tensor, and (b) The accuracies and the computational costs of the RRM and the new AMe. Issue (a) is demonstrated by presenting the required iterations for various cases of load step number (1, 2, 5 and 10), and presenting the energy and residual norms at the end of each iteration in a selected loading step. Whereas, issue (b) is illustrated by comparing the errors in the computed variables (e, a and R) at the end of the loading path (point B), and comparing the computational times of the two methods. The convergence of the iteration in a given load step is measured in terms of the discrete energy norm, which is computed from the incremental nodal displacement vector Ddn+1 and the residual vector Rðdnþ1 Þ as

  ðiÞ ðiþ1Þ ðiÞ DEðdnþ1 Þ :¼ Ddnþ1  R dnþ1 ;

ð95Þ

ðiÞ

where dnþ1 denotes the ith approximation of the nodal displacement corresponding to the actual load ðiÞ step, whereas Ddnþ1 is the incremental nodal displacement at ith iteration (Simo and Hughes, 1998). The termination criteria is defined by ðiÞ

ð1Þ

DEðdnþ1 Þ 6 109 DEðdnþ1 Þ:

ð96Þ

The CT is measured with a dimensionless relative time measure Tcomp. Here, contrary to Example 1, Tcomp = 1 denotes the CT required for the RRM when the loading in A ? B is given in one load step, i.e., path A ? B is not divided into steps. The computed errors are measured in the same manner as in Example 1. Similarly, the error of the strain is defined by

Ekek ¼ 100

ke  eexact k ½%: keexact k

ð97Þ

Table 3 contains the number of iterations for convergence. For illustration purposes, 4 cases are considered: the loading A ? B is divided into 1, 2, 5 and 10 steps, respectively. The obviously better performance of the consistent tangent tensor can be clearly seen in these results. In Table 4, the values of the energy and residual norms at the end of each iteration are displayed, corresponding to the last (5th) step of the case when the loading A ? B is divided into five steps. These results show that, when the iterations are terminated, both the residual and the energy norms have much smaller values if the consistent tangent tensor is used. Table 5 is intended to show the numerical accuracies and the CTs of the two methods, in terms of the number of the applied steps. The errors considered are presented in

1101

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106 Table 3 Iterations for each step using continuum and consistent tangent tensors in Example 2 Step

Number of applied steps 1

1 2 3 4 5 6 7 8 9 10

2

5

10

Dep

Dcons

Dep

Dcons

Dep

Dcons

36

5

22 27

5 4

10 12 15 15 15

5 4 4 4 4

Dep 7 7 8 9 10 10 10 10 10 10

Dcons 4 4 4 4 4 4 3 3 3 3

Four cases are considered: the whole loading path is divided into 1, 2, 5 and 10 steps, respectively.

Table 4 Energy and residual norm values in Example 2 at the end of each iteration of the last (5th) step in the case when the loading path is divided into five steps Iteration

Energy norm [Nmm] Dep

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

5.28 1.01 2.04 4.19 8.67 1.80 3.74 7.78 1.62 3.37 7.00 1.46 3.03 6.30 1.31

Residual norm [N] Dcons

E+3 E+3 E+2 E+1 E+0 E+0 E1 E2 E2 E3 E4 E4 E5 E6 E6

9.14 3.05 7.29 1.83

Dep E+3 E+1 E5 E15

6.19 2.80 1.27 5.80 2.65 1.21 5.50 2.51 1.14 5.22 2.38 1.08 4.96 2.26 1.03

Dcons E+4 E+4 E+4 E+3 E+3 E+3 E+2 E+2 E+2 E+1 E+1 E+1 E+0 E+0 E+0

5.05 1.22 3.89 1.51

E+3 E+1 E5 E9

Figs. 7a–c in log–log plots. As can be seen in Table 5, the AMe requires more CT than the RRM, because of the more complex structure of the integration scheme. However, substantial information can be extracted from the results if we combine the CTs and the resulting errors. For example, if we consider the error Ekek and the corresponding CT of various load step numbers, then we can conclude that the AMe can be faster than the RRM when higher accuracy is required (see Fig. 7d and Table 6). This kind of efficiency of the new method comes from the fact that we need less loading steps in order to achieve the same accuracy. 6. Conclusion The purpose of this paper is twofold. First, it presents semi-analytical solutions for the von Mises elastoplasticity model governed by combined linear isotropic-kinematic hardening. Secondly, it provides the corresponding discretized integration scheme AMe based on the new solution SOLe. The algorithmically consistent tangent tensor is also provided. Depending on the input parameters, both the strain-driven and stress-driven problems are considered and semi-analytical solutions (SOLe

1102

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

Table 5 Comparison of the RRM and the AMe Applied steps

AMe

RRM Tcomp

Riter

Error values [%] Ekak

E k ek 1 2 3 4 5 6 7 8 9 10 15 20 25 30 50 100

1.00 2.49 3.78 5.27 6.87 7.86 9.06 10.85 12.14 13.04 18.90 23.48 29.55 35.53 58.42 115.94

3 7 11 15 19 22 25 29 32 35 49 60 75 90 150 300

13.24 8.11 5.96 4.72 3.92 3.34 2.91 2.58 2.32 2.10 1.44 1.09 0.88 0.74 0.45 0.23

16.17 9.90 7.27 5.77 4.78 4.08 3.56 3.16 2.83 2.57 1.75 1.34 1.08 0.90 0.55 0.28

Tcomp

Riter

ER 1.45 1.07 0.83 0.69 0.59 0.51 0.46 0.41 0.37 0.34 0.24 0.18 0.15 0.13 0.08 0.04

2.69 4.97 7.17 9.06 11.45 13.83 14.92 17.11 19.61 21.40 29.96 37.42 46.08 54.94 90.66 182.12

5 9 13 17 21 25 27 30 33 36 50 61 76 90 150 300

Error values [%] E k ek

Ekak

ER

8.90 3.43 1.86 1.16 0.78 0.56 0.42 0.33 0.26 0.21 0.10 0.05 0.03 0.02 0.01 0.00

10.87 4.19 2.27 1.42 0.96 0.68 0.51 0.40 0.32 0.26 0.12 0.07 0.04 0.03 0.01 0.00

0.94 0.31 0.13 0.08 0.06 0.04 0.03 0.03 0.02 0.02 0.01 0.00 0.00 0.00 0.00 0.00

Riter counts the total number of iterations (sum of step  iterations of step).

Table 6 Comparison of the RRM and the AMe in Example 2, in view of the required accuracy in Ekek Required accuracy in Ekek ½%

10 5 2 1 0.5 0.1

AMe

RRM Required steps

Tcomp

Required steps

Tcomp

2 4 11 22 45 230

2.49 5.27 13.93 25.68 54.04 273.29

1 2 3 5 7 15

2.69 4.97 7.17 11.45 14.92 29.96

and SOLr) are derived for both cases using constant strain and constant stress rate assumptions, respectively. The new discretized integration scheme AMe can be implemented easily into a standard displacement-based finite element code. Since an analytical solution for the inverse incomplete beta function does not exist, the most crucial task during the numerical calculations is to invert an incomplete beta function. An efficient method is proposed in the work of Dominici (2003), where the author derived series expansions for several inverse functions, including the inverse incomplete beta function. In order to obtain information about the validity and the accuracy of the AMe, two simple examples were considered, where the loadings were characterized by non-constant strain rates. The results were compared to those obtained by the well-established and efficient RRM. Resulting from the more complex structure of the new method, it is clearly concluded that the new scheme requires more computational time than the RRM needs. However, if we globally consider the computational time in terms of the required accuracy, then the new method (AMe) can be faster than the RRM. It is explained with their accuracy: the new method can render the same accuracy with less loading steps than the RRM. The presented numerical scheme is not directly applicable for the most widely used nonlinear hardening materials, because it is based on a semi-analytical solution, which cannot, at this time, be derived for nonlinear hardening models. In the case of nonlinear hardening, the governing equations and the resulting system of DEs are discussed, where some new partial results and the difficulties appearing in the derivation of the analytical solution are presented. However, the new scheme can be applied for nonlinear hardening materials if we approximate the nonlinear characteristics of the hardening rule with piecewise linear segments.

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

1103

Fig. 7. Numerical results for the AMe and the RRM in Example 2 in log–log plots. (a) Comparison of errors in the total strain at point B. (b) Comparison of errors in the back stress at point B. (c) Comparison of errors in the stress at point B. (d) The CTs of the two methods considered in terms of the required accuracy in the total strain.

Acknowledgements This research has been supported by the National Development and Research Foundation, Hungary (under Contract: OTKA, K 72572). This support is gratefully acknowledged. Appendix A This appendix presents the detailed derivation of the consistent tangent tensor for the stress update procedure AMe involving elastic–elastoplastic transition (Section 4.1, Case B). Some derivatives required in further calculations are derived first. ok=oDe: Derivatives of the yield condition F(nk, Rn) = 0 with respect to De gives:

2nk :

onk ¼ 0 with oDe

onk ok ¼ 2GkI þ 2GDe  : o De oDe

ð98Þ

By combining the two expressions in (98), we obtain:

ok k ¼ n: oDe ðnk : DeÞ k

ð99Þ

onk/oDe: Using the definition (99) in (98), we have:

onk 2Gk De  nk þ 2GkI: ¼ oDe ðnk : DeÞ

ð100Þ

1104

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

owk/oDe: Recalling definition (67), it follows:

   1 owk o cos wk o cos wk 1 o n k : De ¼ ¼ : sin wk oDe Rn kDek o De owk o De

ð101Þ

After some algebraic manipulations and substitutions, this formula can be written in the following linear combination:

owk ¼ C 1 nk þ C 2 De; o De

ð102Þ

where the scalar parameters C1 and C2 are given by (85) and (86), respectively. own+1/oDe: For simplicity, denote the left and right hand sides of Eq. (66) by z and b1, respectively. Then the derivatives of z and b1 with respect to De are:

oz oz ok oz okDek oz owk ¼ þ þ ; oDe ok oDe okDek oDe owk oDe ob1 ob1 ownþ1 ob1 owk ¼ þ : oDe ownþ1 oDe owk oDe

ð103Þ ð104Þ

Eqs. (103) and (104) are equal; therefore, combining them allows us to express:

ownþ1 ¼ o De



ob1 ownþ1

1 

 oz ok oz okDek oz owk ob1 owk : þ þ  ok oDe okDek oDe owk oDe owk oDe

ð105Þ

After some algebraic manipulations and substitutions, expression (105) reduces to a linear combination of nk and De:

ownþ1 ¼ C 3 nk þ C 4 De; o De

ð106Þ

where the scalar parameters C3 and C4 are given by (87) and (88), respectively. oAs/oDe: The derivative of (71) with respect to De is calculated by

oAs oAs owk oAs ownþ1 ¼ þ : oDe owk oDe ownþ1 oDe

ð107Þ

Using (102), (106) in (107), the derivative of As with respect to De can be expressed as

oAs ¼ a1 nk þ a2 De; o De

ð108Þ

where the scalar parameters a1 and a2 are given in (81) and (82), respectively. oBs/oDe: The derivative of (72) with respect to De is obtained by applying the following chain rule:

oBs oBs ok oBs owk oBs ownþ1 oBs oAs oBs okDek ¼ þ þ þ þ : oDe ok oDe owk oDe ownþ1 oDe oAs oDe okDek oDe

ð109Þ

Finally, we arrive at the following linear form:

oBs ~3 nk þ a4 De; ¼a o De

~3 ¼ a3 þ with a

2Gkð1 þ As Þ ; cos wk Rn kDek

ð110Þ

where the scalar parameters a3 and a4 are given by (83) and (84). The origin of the additional term in ~3 will be clarified later. a After the above derivatives are provided, the derivation of the consistent tangent tensor is presented as follows. Expression (79) can be reformulated using the identity osk/oDe = onk/oDe:

osnþ1 on oAs oBs ¼ ð1 þ As Þ k þ nk  þ De  þ Bs I o De oDe oDe oDe

ð111Þ

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

1105

Substituting (100) in (111) and then multiplying with: T from the right hand side, we obtain

osnþ1 ¼ o De

 nk 

 oAs oBs 2Gkð1 þ As Þ De  nk þ De  :T cos wk Rn kDek o De oDe

þ ðBs þ 2Gkð1 þ As ÞÞT:

ð112Þ

Using (108) and (110) in (112), we arrive at the following expression: osnþ1 oD e

~3 nk þ a4 DeÞÞ : T ¼ ðnk  ða1 nk þ a2 DeÞ þ De  ða 2Gkð1þAs Þ  cos De  nk þ ðBs þ 2Gkð1 þ As ÞÞT: w Rn kDek

ð113Þ

k

Finally, substituting (113) into (78) and collect together the coefficients of the dyad De  nk are the last steps in the construction of the consistent tangent tensor: ornþ1 oenþ1

¼ a1 nk  nk þ a2 nk  De þ a3 De  nk þ a4 De  De þðBs þ 2Gkð1 þ As ÞÞT þ Kd  d:

ð114Þ

References Abramowitz, M., Stegun, I.A., 1968. Handbook of Mathematical Functions. Applied Mathematics Series, vol. 55. Dover Publications, New York. Chapter 5. Armstrong, P.J., Frederick, C.O., 1966. A mathematical representation of the multiaxial Bauschinger effect. Report RD/B/N731, CEGB, Central Electricity Generating Board, Berkeley, UK. Artioli, E., Auricchio, F., Beir ãoda Veiga, L., 2006. A novel ‘optimal’ exponential-based integration algorithm for von-Mises plasticity with linear hardening: theoretical analysis on yield consistency, accuracy, convergence and numerical investigations. Int. J. Numer. Methods Eng. 67, 449–498. Artioli, E., Auricchio, F., Beir ãoda Veiga, L., 2007. Generalized midpoint integration algorithms for J2 plasticity with linear hardening. Int. J. Numer. Methods Eng. 72, 422–463. Auricchio, F., Beirão da Veiga, L., 2003. On a new integration scheme for von-Mises plasticity with linear hardening. Int. J. Numer. Methods Eng. 56, 1375–1396. Auricchio, F., Taylor, R.L., 1995. Two material models for cyclic plasticity: nonlinear kinematic hardening and generalized plasticity. Int. J. Plasticity 11, 65–98. Axelsson, K., Samuelsson, A., 1979. Finite element analysis of elastic–plastic materials displaying mixed hardening. Int. J. Numer. Methods Eng. 14, 211–225. Caddemi, S., 1994. Computational aspects of the integration of the von Mises linear hardening constitutive laws. Int. J. Plasticity 10, 935–956. Carosio, A., Willam, K., Etse, G., 2000. On the consistency of viscoplastic formulations. Int. J. Solids Struct. 37, 7349–7369. Chaboche, J.L., 2008. A review of some plasticity and viscoplasticity constitutive theories. Int. J. Plasticity 1. doi:10.1016/ j.ijplas.2008.03.009. Chan, A.H.C., 1996. Exact stress integration for von Mises elasto-plastic model with constant hardening modulus. Int. J. Numer. Anal. Methods Geomech. 20, 605–613. Chen, W.F., Han, D.J., 1988. Plasticity for Structural Engineers. Springer, New York. Cocchetti, G., Perego, U., 2003. A rigorous bound on error in backward-difference elastoplastic time-integration. Comput. Methods Appl. Mech. Eng. 192, 4909–4927. Dominici, D., 2003. Nested derivatives: a simple method for computing series expansions of inverse functions. Int. J. Math. Math. Sci. 2003, 3699–3715. Dunne, F., Petrinic, N., 2005. Introduction to Computational Plasticity. Oxford University Press, New York. Gratacos, P., Montmitonnet, P., Chenot, J.L., 1992. An integration scheme for Prandtl–Reuss elastoplastic constitutive equations. Int. J. Numer. Methods Eng. 33, 943–961. Hiley, R.A., Rouainia, M., 2008. Explicit Runge–Kutta methods for the integration of rate-type constitutive equations. Comput. Mech. 42, 53–66. Hong, H.K., Liu, C.S., 1997. Prandtl–Reuss elastoplasticity: on–off switch and superposition formulae. Int. J. Solids Struct. 34, 4281–4304. Khan, A.S., Huang, S., 1995. Continuum Theory of Plasticity. John Wiley & Sons, New York. Khoei, A.R., Jamali, N., 2005. On the implementation of a multi-surface kinematic hardening plasticity and its applications. Int. J. Plasticity 21, 1741–1770. Krieg, R.D., Krieg, D.B., 1977. Accuracies of numerical solution methods for the elastic-perfectly plastic model. J. Pressure Vessel Technol. (ASME) 99, 510–515. Krieg, R.D., Xu, S., 1997. Plane stress linear hardening plasticity theory. Finite Elements Anal. Des. 27, 41–67. Liu, C.S., 2004. A consistent numerical scheme for Mises mixed hardening constitutive equations. Int. J. Plasticity 20, 663–704. Loret, B., Prevost, J.H., 1986. Accurate numerical solution for Drucker–Prager elastic–plastic models. Comput. Methods Appl. Mech. Eng. 54, 259–277. Mukherjee, S., Liu, C.S., 2003. Computational isotropic-workhardening rate-independent elastoplasticity. J. Appl. Mech. 70, 644– 648.

1106

A. Kossa, L. Szabó / International Journal of Plasticity 25 (2009) 1083–1106

Nemat-Nasser, S., 2004. Plasticity. A Treatise on Finite Deformation of Heterogeneous Inelastic Materials. Cambridge University Press, Cambridge. Ponthot, J.P., 2002. Unified stress update algorithms for the numerical simulation of large deformation elasto-plastic and elastoviscoplastic processes. Int. J. Plasticity 18, 91–126. Reuss, E., 1930. Berücksichtigung der elastischen Formänderung in der Plastizitätstheorie. Zeits. Angew. Math. Mech. (ZAMM) 10, 266–274. Rezaiee-Pajand, M., Nasirai, C., 2008. On the integration schemes for Drucker–Prager’s elastoplastic models based on exponential maps. Int. J. Numer. Methods Eng. 74, 799–826. Ristinmaa, M., 1990. Tangent Modulus in Finite Element Calculations for Non-Linear Materials. LUTFD2/(TFHT-3033), Div. of Solid Mechanics. Lund University. Ristinmaa, M., Tryding, J., 1993. Exact integration of constitutive equations in elasto-plasticity. Int. J. Numer. Methods Eng. 36, 2525–2544. Romashchenko, V.A., Lepikhin, P.P., Ivashchenko, K.B., 1999. Exact solution of problems of flow theory with isotropic-kinematic hardening. Part 1. Setting the loading trajectory in the space of stress. Strength Mater. 31, 582–591. Simo, J.C., Hughes, T.J.R., 1998. Computational Inelasticity. Springer, Berlin. Spanier, J., Oldham, K.N., 1987. An Atlas of Functions. Springer, Berlin. Stein, E., de Borst, R., Hughes, T.J.R., 2004. Encyclopedia of Computational Mechanics, vol. 2. John Wiley & Sons, New York. Szabó, L., accepted for publication. A semi-analytical integration method for J2 flow theory of plasticity with linear isotropic hardening. Comput. Methods Appl. Mech. Eng. Szabó, L., Kovács, Á., 1987. Numerical implementation of Prager’s kinematic hardening model in exactly integrated form for elastic–plastic analysis. Comput. Struct. 26, 815–822. Wallin, M., Ristinmaa, M., 2001. Accurate stress updating algorithm based on constant strain rate assumption. Comput. Methods Appl. Mech. Eng. 190, 5583–5601. Wallin, M., Ristinmaa, M., 2008. An alternative method for the integration of continuum damage evolution laws. Comput. Mech. 41, 347–359. Wang, X., Chang, L., 1985. Exact integration of constitutive equations of kinematic hardening material and its extended applications. SMIRT-8, Brussels, Proc. Paper L2/3, pp. 65–70. Wei, Z., Peric, D., Owen, D.R.J., 1996. Consistent linearization for the exact stress update of Prandtl–Reuss non-hardening elastoplastic models. Int. J. Numer. Methods Eng. 39, 1219–1235. Wolfram Research, Inc., 2005. Mathematica Edition: Version 5.2. Wolfram Research, Inc., Champaign, IL. Yoder, P.J., Whirley, R.G., 1984. On the numerical implementation of elastoplastic models. J. Appl. Mech. 51, 283–288. Yang, Y., Muraleetharan, K.K., Yu, H.S., 2008. Generalized trapezoidal numerical integration of an advanced soil model. Int. J. Numer. Anal. Methods Geomech. 32, 43–64.