A nonlinear finite element plasticity formulation without matrix inversions

Finite Elements in Analysis and Design 112 (2016) 11–25

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

A nonlinear finite element plasticity formulation without matrix inversions Fabio De Angelis a,n,1, Robert L. Taylor b a b

Department of Structures for Engineering and Architecture, University of Naples Federico II, Naples, Italy Department of Civil and Environmental Engineering, University of California, Berkeley, USA

art ic l e i nf o

a b s t r a c t

Article history: Received 26 July 2015 Received in revised form 3 December 2015 Accepted 8 December 2015

In the present paper a robust computational procedure is proposed for the finite element solution of elasto-plastic boundary value problems in the structural analysis of ductile materials. An implicit numerical scheme is detailed based on a return mapping algorithm. In the present computational scheme the local constitutive equations lead to the solution of a nonlinear scalar equation which ultimately reduces to a single variable algebraic (polynomial) equation in the plastic rate parameter. The proposed integration scheme allows to find the analytical solution of the algebraic equation in closed form. Consequently, in the present algorithmic procedure no iterative method is required to solve the local constitutive problem. The algorithmic scheme presented herein is also characterized by an alternative formulation of the algorithmic counterpart of the plastic consistency condition. Furthermore, in the algorithmic strategy a suitable procedure is presented for the consistent linearization of elasto-plastic boundary value problems. In fact an alternative procedure is proposed herein for the consistent derivation of the tangent operator associated to the algorithmic scheme so that in the present approach the consistent tangent operator is determined without the necessity of computing matrix inversions. The presented computational approach combined with the proposed consistent tangent operator ensures a quadratic rate of asymptotic convergence when an iterative method is pursued for the global solution procedure of elasto-plastic boundary value problems. Numerical applications and illustrative examples are finally reported to show the robustness and effectiveness of the proposed computational procedure in the finite element analysis of elasto-plastic boundary value problems. & 2015 Elsevier B.V. All rights reserved.

Keywords: Finite element method Elastoplasticity Computational procedures

1. Introduction Numerical integration procedures for the analysis of elastoplastic problems appear to have been initially proposed in the form of return mapping algorithms based on an elastic predictor and a plastic corrector scheme introduced by Wilkins [38] and Maenchen and Sack [23]. Extensions to include isotropic and kinematic hardening were presented by Krieg and Key [20], and Krieg and Krieg [21]. Reflections on the accuracy and stability of this class of numerical algorithms were reported e.g. by Krieg and Krieg [21], Ortiz and Popov [26] and Simo [30]. The opportunity of considering a tangent operator consistent with the numerical integration scheme was discussed by Nagtegaal [24] and Simo and n

Corresponding author. E-mail addresses: [email protected] (F. De Angelis), [email protected] (R.L. Taylor). 1 Visiting Scholar at the Department of Civil and Environmental Engineering, University of California, Berkeley, USA. http://dx.doi.org/10.1016/j.finel.2015.12.007 0168-874X/& 2015 Elsevier B.V. All rights reserved.

Taylor [35] with the aim to preserving the quadratic rate of asymptotic convergence typical of Newton's iteration schemes. Among the discussions and proposals for numerical integration schemes we refer to the works presented by, among others, Hughes and Taylor [17], Simo and Govindjee [31], Simo and Govindjee [32], Peric [28], Papadopoulos and Taylor [27], Simo and Hughes [33] and DeAngelis [10]. Numerical integration schemes for applications to non-smooth yield criteria were illustrated by Simo et al. [34] and a tangent operator suitable for general yield criteria was proposed by Alfano et al. [2]. For a comprehensive account see, e.g., Zienkiewicz et al. [40], Wriggers [39], Kojic and Bathe [19] and de Souza Neto et al. [36]. In recent years more elaborate constitutive theories have been proposed to capture the behavior of elastoplastic materials subject to complex loading conditions for the need of reproducing the observed features in the loading–unloading–reloading processes. With reference to the more elaborate constitutive theories standard return mapping algorithms require suitable modifications in order to account for nonlinear hardening effects. Moreover repeated

12

F. De Angelis, R.L. Taylor / Finite Elements in Analysis and Design 112 (2016) 11–25

loading conditions in elastoplasticity boundary value problems often require extensive calculations which necessitate robust and efficient numerical integration schemes for a proper and accurate simulation of the structural behavior. Among the algorithmic procedures suitable for complex and elaborate loading conditions we refer to the works presented by, among others, Doghri [12], Lubliner et al. [22], Auricchio and Taylor [5], Chaboche and Cailletaud [8], Hopperstad and Remseth [16], Nukala [25], Artioli et al. [4], DeAngelis and Taylor [11]. In the present paper a robust algorithmic procedure is proposed for the finite element analysis of small strains elasto-plastic boundary value problems of ductile metals. A return mapping algorithm is adopted to enforce plastic consistency and a nonlinear kinematic hardening mechanism is used to represent the components of the kinematic hardening variables. The proposed algorithmic scheme produces a robust and efficient numerical procedure for the computational simulation of elasto-plastic boundary value problems under elaborate loadings. The present algorithmic scheme has the advantage of reducing the local constitutive equations to a single nonlinear scalar equation. With reference to previously presented proposals in the present approach such nonlinear scalar equation is expressed in a straightforward form since it reduces furtherly to a simple representation of the nonlinear scalar equation namely a single variable algebraic (polynomial) equation in the plastic rate parameter. The present algorithmic approach is therefore simplified with respect to other previously proposed approaches in which the local constitutive equations are also reduced to a nonlinear scalar equation although not to a single variable algebraic (polynomial) equation. In the present algorithmic scheme it is taken advantage of the particularly simple form of the equation since the analytical solutions of the single variable algebraic equation can be found in exact closed form. Consequently, in the present approach no iterative method is required for the solution of the local constitutive problem. In the proposed numerical integration scheme the algorithmic plastic consistency condition is expressed by an alternative form. In addition an alternative procedure is presented for the consistent linearization of elasto-plastic boundary value problems. In fact, a notable feature of the present algorithmic scheme is that the consistent tangent operator is obtained by the following a procedure so that matrix inversions are not required within the present approach. The algorithmic procedure together with the present algorithmic plastic consistency condition and the proposed consistent tangent operator provides a quadratic rate of asymptotic convergence when used with a Newton iterative method for the global solution of elasto-plastic boundary value problems. The effectiveness of the present computational approach is also emphasized for the solution of the local constitutive relations of elasto-plastic boundary value problems. In fact, in the present approach the nonlinear scalar equation is solved in a closed form and thus no numerical methods are required to accelerate or improve the rate of convergence for the solution of the nonlinear scalar equation representing the local constitutive problem. The process of global convergence naturally results in a quadratic rate of asymptotic convergence typical of a Newton iterative method for the global iterative solution procedure of elasto-plastic boundary value problems. Finally, the robustness and the efficiency of the present computational procedure are shown by providing specific numerical examples. The effectiveness of the proposed algorithmic scheme is illustrated for elasto-plastic boundary value problems subject to complex and elaborate monotonic and cyclic loadings. In the numerical examples the robustness of the computational procedure is emphasized by resulting in a stable and accurate numerical solution procedure for elasto-plastic boundary value problems subject to different types of loading conditions.

2. The elastoplastic continuum We consider a continuum body B which occupies the reference configuration Ω
Rn , with 1 r n r 3. Let T
R þ be the time interval of interest, V the space of displacements, D the strain space and S the dual stress space. The compatible strain tensor is denoted by ε ¼ ∇s u : Ω  T -D, where u : Ω  T -V is the displacement and ∇s represents the symmetric part of the gradient operator. The stress tensor σ : Ω  T -S is expressed by

σ ¼ s þp 1;

ð1Þ

where s ¼ dev σ ¼ σ p 1 is the stress deviator, p 1 is the spherical def

part, p ¼ 13 trðσ Þ is the mean stress and 1 is the rank two identity tensor. The relative stress Σ is also introduced as def

Σ ¼ s  α;

ð2Þ

where α is the deviatoric back stress. Likewise, the strain tensor is expressed by

ε ¼ e þ 13 θ 1;

ð3Þ

θ 1 is the strain deviator and θ ¼ trðεÞ where e ¼ dev ε ¼ ε represents the volume change. The deviatoric part of the total strain is additively decomposed into an elastic and a plastic part def

def

 13

e ¼ ee þ ep ;

ð4Þ

e

p

where e is the elastic deviatoric strain and e is the plastic deviatoric strain. The pressure and the volume change are related by the elastic constitutive law p ¼ K θ, where K represents the bulk modulus. The stress deviator and the elastic deviatoric strain are related by the constitutive law s ¼ 2G ee ¼ 2G½e  ep ;

ð5Þ

where G represents the shear modulus. The present procedure is proposed within the context of a J2 von Mises yield function form. Accordingly, we assume the class of J2 material models with a limit equation expressed by the yield criterion qffiffiffi ð6Þ f ðσ ; α; κ Þ ¼ J dev σ  α J  κ ðχ iso Þ ¼ J Σ J  23 ðσ yo þ χ iso Þ r0; qffiffiffi qffiffiffi in which κ ðχ iso Þ ¼ 23 σ y ¼ 23 ðσ yo þ χ iso Þ is the current radius of the yield surface in the deviatoric plane and σyo is the uniaxial yield stress of the virgin material. The conjugate variable related to isotropic hardening is specified as χ iso ¼ R, where R is the increment of the yield stress with respect to the uniaxial yield stress of the virgin material. For linear isotropic hardening it is χ iso ¼ H iso e pq , where Hiso represents ffiffiffi p def R t 2 _p the isotropic hardening modulus and e ¼ 0 3 J e J dt represents def

the equivalent plastic strain. For nonlinear isotropic hardening it is often assumed p

χ iso ¼ Hiso ðe p Þm or R ¼ R1 ð1  e  b e Þ;

ð7Þ

where m, R1 and b are material parameters. The back stress rate is expressed by Prager's law which gives for a linear kinematic hardening behavior [29]

α_ ¼ 23 Hkin e_ p ;

ð8Þ

where Hkin represents the kinematic hardening modulus. In the literature a nonlinear kinematic hardening behavior is often assumed by adopting the model proposed by Armstrong and Frederick [3], see also Frederick and Armstrong [13],

α_ ¼ 23 Hkin e_ p H nl e_ p α;

ð9Þ

where H nl is a non-dimensional material dependent parameter which is null in case of linear kinematic hardening. A better approximation for nonlinear kinematic hardening behavior results in adding several components of the back stress,

F. De Angelis, R.L. Taylor / Finite Elements in Analysis and Design 112 (2016) 11–25

Given Eq. (17), the updated value of the deviatoric stress is expressed by

with different recall constants, see e.g. Chaboche [7],

α¼

M X

αi ;

i¼1

α_ i ¼ 23 Hkin; i e_ p H nl; i e_ p αi :

ð10Þ

For an associative elastoplastic behavior at small strains the evolution flow law is expressed by p e_ ¼ γ_

∂f ∂f ¼ γ_ ¼ γ_ n; ∂σ ∂Σ

ð11Þ

where γ_ is the plastic rate multiplier. The second rank tensor n is defined by def

n¼

s ¼ 2G½e  epn  2Gλn;

ð19Þ

and the updated back stress in discrete form is expressed by qffiffiffi α  αn ¼ 23 Hkin λn  Hnl 23 λα; ð20Þ which can be also written as

α ¼ T λ αn þ 23 H kin T λ λn;

ð21Þ

where we have set

Σ ; JΣJ

ð12Þ

and, accordingly, the equivalent plastic strain rate is expressed by qffiffiffi e_ p ¼ 23 γ_ : ð13Þ To complete the representation of the elastoplastic model the loading–unloading conditions are introduced according to the Kuhn–Tucker optimality conditions

γ_ Z 0;

13

f ðσ ; α; κ Þ r 0;

γ_ f ðσ ; α; κ Þ ¼ 0;

ð14Þ

see e.g. Crisfield [9], Simo and Hughes [33] and Zienkiewicz et al. [40]. The proposed formulation may be investigated to include more hardening functions, whereas the treatment of more general yield functions is outside the scope of the present analysis.

3. Return mapping algorithmic formulation We consider a return mapping discrete formulation based on an elastic prediction and a plastic correction algorithmic scheme. According to a time discrete formulation the time interval of interest ½0; T
R þ is subdivided into a finite number of time steps t n -t n þ 1 . The generic variable at time tn is assumed to be known as evaluated at the previous time step and the value of the variable at time t n þ 1 is to be evaluated according to the discrete form of the evolutive equations. The generic variable u evaluated at time tn will be denoted with the subscript n, so that un ¼ uðt n Þ, whereas to simplify notation the same variable with no subscript will denote the updated value of the variable at time t n þ 1 , so that u ¼ uðt n þ 1 Þ. We also consider a strain driven approach, accordingly the solution at time tn is known and given by the ordered set ðσ n ; εn ; sn ; en ; epn ; e pn ; αn Þ:

ð15Þ

The value of the increment of the displacement field Δu ¼ u un is assigned at time tn, which corresponds to assigning the value of the total strain ε ¼ εn þ ∇s ðΔuÞ at time t n þ 1 . Consequently, we need to compute the solution at time t n þ 1 represented by the ordered set def

ðσ ; ε; s; e; ep ; e p ; αÞ

ð16Þ

according to the discrete forms of the evolutive elastoplasticity laws. In the framework of a backward Euler integration scheme, the evolutive equation for the deviatoric plastic strain rate (11) in discrete form is expressed by ep ¼ epn þ λ n;

ð17Þ

and the equivalent plastic strain rate (13) in discrete form is expressed by qffiffiffi e p ¼ e pn þ 23 λ; ð18Þ where we have set λ ¼ Δγ n þ 1 ¼ γ_ n þ 1 Δt. def

Tλ ¼ 

1þ

1 qffiffiffi 2 3

1  ¼ λ; R H nl λ

Rλ ¼ ð1 þ

qffiffiffi 2 3

H nl λÞ:

ð22Þ

In the above expressions and in the sequel the superscript λ is used to point out a dependence of the variable on the increment of the plastic rate multiplier in the step. It will be noted in the sequel that this feature is typical of a nonlinear kinematic hardening behavior. Accordingly, the updated value of the relative stress is expressed by

Σ ¼ s  α ¼ 2G½e  epn   T λ αn  U λ n; where   h i H U λ ¼ 2G þ 23 kin λ ¼ 2G þ 23 Hkin T λ λ: λ R

ð23Þ

ð24Þ

In the framework of a predictor–corrector scheme, for each time step ½t n ; t n þ 1 , in the first part of the algorithm (elastic predictor) a trial state is computed by assuming a purely elastic behavior in the step. If the trial elastic stress state satisfies the yield condition (6), then it is accepted as the solution at t n þ 1 , otherwise in the second part of the algorithm (plastic corrector) a plastic correction is computed to correct the trial elastic state according to the discrete forms of the evolutive elastoplasticity laws. Consequently, in the generic time step t n -t n þ 1 , the trial elastic state is computed by assuming that the material behavior is purely elastic and the plastic components of the strains are unchanged and set equal to those at time tn (elastic prediction): 8 TR > λ ¼0 > > > > p; TR > ¼ epn e > > > < e p; TR ¼ e p n ð25Þ p TR > > s ¼ 2G½e  en  > > > > > αTR ¼ αn > > : TR Σ ¼ sTR  αTR ¼ sTR  αn : The yielding condition (6) is then checked at time t n þ 1 as a function of the trial values. If the elastic trial state lies within the boundary of the elastic domain then it is an admissible state. In such case the elastic trial state is assumed as the solution at time t n þ 1 and the second part of the algorithm is skipped. Conversely, if the elastic trial state is not admissible then a plastic correction is performed. The increment λ of the plastic rate parameter is computed by enforcing the satisfaction of the nonlinear limit equations at time t n þ 1 as specified in the next section. The normal n to the yield surface is also computed as specified in the next section. Consequently, the updated variables are computed by projecting the trial elastic state onto the boundary of the elastic domain. The variables are therefore updated in terms of the trial elastic state, the increment λ of the plastic rate parameter in the step and the normal n to the yield surface by fulfilling the discrete

14

F. De Angelis, R.L. Taylor / Finite Elements in Analysis and Design 112 (2016) 11–25

forms of the evolutive elastoplasticity laws (plastic correction): 8 ep ¼ ep; TR þ λn > > qffiffiffi > > > p p; TR > þ 23 λ > >e ¼e < ð26Þ s ¼ sTR  2Gλn > > > > α ¼ T λ αTR þ 23 H kin T λ λn > > > > : Σ ¼ ΣTR þ ð1 T λ ÞαTR  U λ n:

the normal to the yield surface nðλÞ and its differential based on the expression (32). In the sequel we will show that this approach leads to a simplified numerical scheme and a simplified expression of the consistent tangent operator.

Consequently, in the generic time step t n -t n þ 1 which requires a plastic correction, after computing the normal n to the yield surface and the increment λ of the plastic rate parameter, the variables are updated at time t n þ 1 according to Eqs. (26).

The limit condition is provided by the condition (6) and enforced at time t n þ 1 by qffiffiffi J ΣðλÞ J  23 σ y ðλÞ ¼ 0: ð33Þ

3.1. Determination of the normal n to the yield surface The updated value of the relative stress is expressed by

Σ ¼ ΣTR þ ð1  T λ ÞαTR  U λ n ¼ sTR  T λ αTR  U λ n:

ð27Þ

In Eq. (27) the term Σ ¼ s  α is the trial value of the relative stress for the case of plasticity with a linear kinematic hardening TR rule. In particular we observe that Σ is independent of λ. In the present context characterized by nonlinear kinematic hardening rules we find convenient to set a trial-like value of the relative stress for plasticity with nonlinear kinematic hardening as TR

TR

TR

ΣλTR ¼ ΣTR þð1 T λ ÞαTR ¼ sTR T λ αTR :

ð28Þ

This is at variance with the case of linear kinematic hardening since now, for nonlinear kinematic hardening, the trial-like value λ of the relative stress ΣTR depends upon λ. Accordingly, the updated value of the relative stress (27) at time t n þ 1 is expressed by λ

λ

Σ ¼ ΣTR  U n:

ð29Þ

The normal n to the yield surface is determined by ∂f ∂f ∂Σ ∂f ∂ðdev σ  αÞ ∂ J Σ J Σ Σ I ¼ ; ¼ ¼ ¼ I ¼ n¼ ∂σ ∂Σ ∂σ ∂Σ ∂σ ∂Σ dev J Σ J dev J Σ J

ð31Þ

that is the normal nðλÞ to the yield surface at time t n þ 1 is evaluated based on the updated value of the relative stress ΣðλÞ at time t n þ 1 . From the above equation we observe that ΣðλÞ is collinear with λ nðλÞ. Accordingly, from Eq. (29) we obtain that ΣTR ðλÞ must also be λ collinear with nðλÞ, and consequently we can write ΣTR ðλÞ ¼ λ J ΣTR ðλÞ J nðλÞ. Therefore the updated normal to the yield surface can also be alternatively expressed by nðλÞ ¼

ΣλTR ðλÞ ; λ J ΣTR ðλÞ J

λ

J ΣðλÞ J ¼ J ΣTR ðλÞ J U λ :

ð34Þ

By considering Eq. (34), the limit condition (33) yields qffiffiffi λ J ΣTR ðλÞ J  U λ  23 σ y ¼ 0:

ð35Þ

λ

The term J ΣTR ðλÞ J is expressed as λ

J ΣTR ðλÞ J ¼ ½Sss  2Ssα T λ þ Sαα ðT λ Þ2 1=2 ;

ð36Þ

where we have set Sss ¼ sTR  sTR ¼ J sTR J 2 Ssα ¼ sTR  αTR ð37Þ

Accordingly, Eq. (35) is specialized to ð30Þ

where  1Þ. It is noted that in plasticity with nonlinear kinematic hardening rules the normal to the yield surface n also depends upon λ. This differs with the standard case of plasticity with linear kinematic hardening rules, in which the tensor n is independent from λ. Accordingly Eq. (30) is expressed by

ΣðλÞ ; J ΣðλÞ J

In this section we show that within the present algorithmic approach the above nonlinear scalar equation may be expressed as an algebraic (polynomial) equation in λ. This equation is to be solved for the increment λ of the plastic rate parameter in the step t n -t n þ 1 . In the previous section it has been shown that both ΣðλÞ and ΣλTR ðλÞ are collinear with nðλÞ. Accordingly, Eq. (29) yields the following relation:

Sαα ¼ αTR  αTR ¼ J αTR J 2 :

Idev ¼ I  13 ð1

nðλÞ ¼

3.2. Nonlinear scalar equation for the determination of the increment λ of the plastic rate parameter

ð32Þ

that is the normal nðλÞ to the yield surface at time t n þ 1 can also be determined based on the trial-like (not updated) value of the λ relative stress ΣTR ðλÞ. In the numerical scheme recently proposed by DeAngelis and Taylor [11], the determination of the updated normal to the yield surface was based on the expression (31). At variance with such approach in the present numerical scheme we will evaluate both

Sss ðRλ Þ2  2Ssα Rλ þ Sαα ¼ 4G2 ðRλ Þ2 λ þ 83 GH kin Rλ λ þ 49 H 2kin λ qffiffiffi qffiffiffi þ4G 23 σ y ðRλ Þ2 λ þ 43 H kin 23 σ y Rλ λ þ 23 σ 2y ðRλ Þ2 : 2

2

2

The yield stress is expressed by  qffiffiffi  qffiffiffi σ y ¼ σ yo þ Hiso e pn þ 1 ¼ σ yo þ Hiso e pn þ 23 λ ¼ σ y;n þ Hiso 23 λ;

ð38Þ

ð39Þ

and, by recalling the expression (22)2 for Rλ , Eq. (38) gives the nonlinear scalar equation gðλÞ ¼ C 1 λ þ C 2 λ þC 3 λ þ C 4 λ þ C 5 ¼ 0; 4

3

2

ð40Þ

where the coefficients of the equation are provided by h i 2 2 2 8 C 1 ¼ 83 G2 H 2nl þ 16 9 GH iso H nl þ 27 H iso H nl ;  qffiffiffi qffiffiffi qffiffiffi qffiffiffi 2 C 2 ¼ 8 23 G2 H nl þ 83 23 GH kin H nl þ 83 23 GH 2nl σ y;n þ 16 3 3 GH iso H nl  qffiffiffi qffiffiffi qffiffiffi þ 89 23 H kin H iso H nl þ 89 23 H 2iso H nl þ 89 23 H iso H 2nl σ y;n ; h 8 8 C 3 ¼ 4G2 þ 83 GH kin þ 49 H 2kin þ 16 3 GH nl σ y;n þ 3 GH iso þ 9 H kin H iso

i 4 2 2 2 2 þ 89 H kin H nl σ y;n þ 49 H 2iso þ 16 9 H iso H nl σ y;n þ 9 H nl σ y;n  3 H nl Sss ;

 qffiffiffi qffiffiffi qffiffiffi qffiffiffi C 4 ¼ 4 23 Gσ y;n þ 43 23 H kin σ y;n þ 43 23 H iso σ y;n þ 43 23 H nl σ 2y;n  qffiffiffi 2 23 H nl ðSss  Ssα Þ ;

F. De Angelis, R.L. Taylor / Finite Elements in Analysis and Design 112 (2016) 11–25

C5 ¼

h

2 3

i

σ 2y;n  ðSss 2Ssα þ Sαα Þ :

ð41Þ

The nonlinear scalar equation (40) represents an algebraic (polynomial) equation which is to be solved for the computation of λ in the generic time step t n -t n þ 1 . Accordingly, one of the advantages of the present algorithmic approach is that the local solution of the constitutive equations reduces to only one nonlinear scalar equation, with the form presented herein. It is to be noted that in the literature other proposals have been presented for algorithmic schemes which also reduce the local solution of the constitutive equations to one nonlinear scalar equation, see e.g. Hartmann et al. [15] and Kobayashi and Ohno [18]. At variance with such proposals, in the present approach the constitutive equations reduce to a straightforward and simple form of nonlinear scalar equation which ultimately boils down to a single variable algebraic (polynomial) equation. Furthermore it is noted that for complex constitutive equations and highly nonlinear scalar equations the solution performed via a Newton iterative method sometimes experiences failures and does not always converge due to the occurring high gradients and the resulting round-off errors. In the sequel it is shown that the herein proposed approach leads naturally to an algorithmic scheme which ensures efficiency and robustness to the computational procedure. In particular it is emphasized that with the present algorithmic scheme no numerical procedures are required to accelerate or improve the rate of convergence for the solution of the nonlinear scalar equation, such as for instance Aitken's process adopted in Kobayashi and Ohno [18]. Moreover, it is emphasized that with the present algorithmic approach the research for the analytical solutions of the single variable algebraic (polynomial) equation (40) can be pursued in a closed form, i.e. with no recourse to iterative methods. In fact, the quartic equation (40) can be solved in a closed form for the real positive roots, see e.g. among others Abramowitz [1], Beyer [6], Hacke [14]. The increment of the plastic rate parameter λ is therefore evaluated, in the step t n -t n þ 1 , as the smallest positive real root of the quartic equation (40). The proposed approach leads naturally to a robust and effective computational solution procedure for elastoplastic structural problems subject to complex loading conditions.

4. Algorithmic plastic consistency condition For the derivation of the algorithmic plastic consistency condition we operate by assuming the updated normal nðλÞ to the yield surface as expressed by Eq. (32) and evaluated as a function of the trial-like value of the relative stress. Likewise, in the limit condition the relative stress at time t n þ 1 is expressed as a function of the trial-like value of the relative stress. In the enforcement of the plastic consistency condition the differentiation of the terms stems from such positions. We note that this approach is at variance with the approach presented by DeAngelis and Taylor [11] in which the normal to the yield surface at time t n þ 1 was based on the expression (31) and evaluated as a function of the updated relative stress at time t n þ 1 and the limit condition and its differentiation were expressed as a function of the updated relative stress. In the sequel we will show that the approach presented herein provides an alternative formulation of the algorithmic plastic consistency condition and an alternative and simplified algorithmic scheme. Given the relation (34), at time t n þ 1 the limit condition (33) is expressed by qffiffiffi λ f ¼ J ΣTR J  U λ  23 σ y ¼ 0: ð42Þ The enforcement of the consistency condition df ¼ 0 provides qffiffiffi qffiffiffi  h i λ 2 2 σ þ H ¼ 0: ð43Þ df ¼ d J ΣTR J  dU λ d y;n iso 3 3λ

15

The first differential of the plastic consistency condition (43) is expressed by h i   ΣλTR λ λ λ λ λ d J ΣTR J ¼ dΣλ J ΣTR J dΣTR ¼ dΣTR ¼ n dΣTR ; λ TR J ΣTR J

ð44Þ

λ

and the term dΣTR is given by h i h i λ dΣTR ¼ d sTR  T λ αTR ¼ dsTR  d T λ αTR h i h i ¼ dsTR  dðT λ ÞαTR þT λ dαTR ¼ dsTR  dλ ðT λ Þdλ  αTR ¼ dsTR þ H λ αTR dλ; in which, respectively, dα dλ T λ is expressed by 2 3 6 dλ T ¼ dλ 4 λ

1þ

1 qffiffiffi 2 3

ð45Þ TR

¼ 0 since α

TR

¼ αn , and the derivative

qffiffiffi  23 H nl 7 qffiffiffi  ¼  Hλ ; 5¼ 2 1 þ 2 23 H nl λ þ 23 H 2nl λ H nl λ

where we have set qffiffiffi qffiffiffi 2 2 qffiffiffi  2 H nl 3 3 H nl qffiffiffi Hλ ¼   ¼  2 ¼ 23 H nl T λ : 2 1 þ 2 23 H nl λ þ 23 H 2nl λ Rλ

ð46Þ

ð47Þ

Given the expression (25)4, it results dsTR ¼ 2G de, and the relation (45) specializes to λ

dΣTR ¼ dsTR þH λ αTR dλ ¼ 2G de þ H λ αTR dλ:

ð48Þ

Accordingly, the first contribution (44) of the plastic consistency condition (43) is expressed by   h i

λ λ d J ΣTR J ¼ n dΣTR ¼ n  2G de þ H λ αTR dλ ¼ 2G n  de  2

qffiffiffi þ H λ ðαn  nÞ dλ ¼ 2G n  de þ 23 H nl T λ ðαn  nÞ dλ: ð49Þ The second contribution of the plastic consistency condition (43) is given by h   i h   dU λ ¼ dλ 2Gλ þ 23 H kin T λ λ dλ ¼ 2G þ dλ 23 H kin T λ  λ  i h i þ 23 H kin T λ dλ ¼ 2G  23 H kin H λ λ þ 23 H kin T λ dλ  qffiffiffi  2  ¼ 2G þ 23 H kin T λ  23 H kin 23 H nl T λ λ dλ: ð50Þ The third contribution of the plastic consistency condition (43) is provided by qffiffiffi qffiffiffi 

2 2 σ þ H ð51Þ ¼ dλ 23 H iso λ dλ ¼ 23 H iso dλ: d y;n iso 3 3λ Accordingly, the plastic consistency condition (43) specializes to  2 h

qffiffiffi df ¼ 2G n  de þ 23 H nl T λ ðαn  nÞ dλ  2G þ 23 H kin T λ qffiffiffi  2   23 H kin 23 H nl T λ λ dλ  23 H iso dλ ¼ 0; ð52Þ and it results  qffiffiffi  2

2G n  de ¼ 2G þ 23 H iso þ 23 H kin T λ  23 H nl T λ ðn  αn Þ qffiffiffi  2   23 H kin 23 H nl T λ λ dλ:

ð53Þ

The above equation can be represented by

dλ ¼ A d n  de ;

ð54Þ

16

F. De Angelis, R.L. Taylor / Finite Elements in Analysis and Design 112 (2016) 11–25

Fig. 3. Tension strip. End stress versus average strain for a loading program in tension–compression with increasing both levels of loading (loading program given in Fig. 2, time increment Δt ¼ 0:05 and t max ¼ 13). The end stress is the traction at the top and bottom boundaries of the section. The average strain is a measure of the strain evaluated as the displacement of the boundary upper edge at vertical axis of symmetry over the height. Fig. 1. Tension strip. Finite element mesh.

Fig. 4. Contour plot of the Mises stress in the strip for t max ¼ 13. Loading program in tension–compression with increasing both levels of loading (loading program given in Fig. 2).

Fig. 2. Cyclic loading program with increasing levels of loading.

where we have set 2G Ad ¼  qffiffiffi qffiffiffi  2  2  λ 2 2 2 2G þ 3 Hiso þ 3 H kin T  3 H nl T λ ðn  αn Þ  23 Hkin 23 H nl T λ λ

ð55Þ We note that the algorithmic plastic consistency condition (54) and the associated discrete algorithmic coefficient (55) differ from the algorithmic plastic consistency condition presented by DeAngelis and Taylor [11] due to the different algorithmic scheme adopted herein. In particular the algorithmic plastic consistency condition and the associated discrete algorithmic coefficient presented by DeAngelis and Taylor [11] depend on α, that is they necessitate the evaluation of the updated value of the back stress at time t n þ 1 . Conversely, in the present algorithmic scheme the expressions of the algorithmic plastic

consistency condition and the associated discrete algorithmic coefficient (55) do not require the evaluation of the updated value of the back stress, since they depend on αn that is the value of the back stress at time tn as evaluated at the previous converged time step. This device simplifies notably the algorithmic scheme and it ensures robustness to the process of convergence and to the overall effectiveness of the computational procedure.

5. Consistent tangent operator In the present section we derive the expression of the consistent tangent operator associated to the proposed algorithmic scheme. One of the advantages of the present algorithmic scheme is that the consistent tangent operator is readily obtained in a simplified form with respect to other algorithmic schemes, see e.g. Auricchio and Taylor [5] and DeAngelis and Taylor [11], in which it is required the inversion of a matrix. Remarkably, in the present algorithmic approach no computation of matrix inversions is requested for the determination of the

F. De Angelis, R.L. Taylor / Finite Elements in Analysis and Design 112 (2016) 11–25

17

Table 1 Tension strip. Relative residual norm values and iterations. Cyclic loading program in tension–compression with increasing both levels of loading (loading program given in Fig. 2, time increment Δt ¼ 0:05 and t max ¼ 13). Relative residual norm values Load step 1 6 7 8 9 10 20 60 100 140 180 220 260

Time

Iteration

1

2

3

4

5

6

5.00E  02 3.00E  01 3.50E  01 4.00E  01 4.50E  01 5.00E  01 1.00E þ00 3.00E þ00 5.00E þ00 7.00E þ00 9.00E þ00 1.10E þ01 1.30E þ01

el el el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl

1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00

4.18E  14 6.85E  14 2.35E 02 6.70E 02 7.37E  02 8.03E  02 6.27E  01 6.34E  01 6.39E  01 6.43E  01 6.47E  01 6.51E  01 6.54E  01

4.38E  03 5.37E  03 7.34E  03 3.68E  03 1.24E  01 1.12E  01 1.13E  01 1.24E  01 1.31E  01 1.35E  01 1.48E  01

3.87E  06 1.29E  05 9.77E 06 4.46E  06 8.07E 03 8.47E  03 4.55E  03 6.44E  03 6.31E  03 8.48E  03 8.71E  03

1.71E  11 1.30E  10 6.11E  11 1.20E  11 1.21E  05 2.23E  05 1.78E  05 2.76E  05 4.04E  05 1.90E  04 7.56E  05

1.11E  10 2.71E  10 2.97E  10 7.05E  10 1.57E  09 3.85E  09 5.97E 09

Table 2 Tension strip. Relative energy norm values and iterations. Cyclic loading program in tension–compression with increasing both levels of loading (loading program given in Fig. 2, time increment Δt ¼ 0:05 and t max ¼ 13). Relative energy norm values Load step 1 6 7 8 9 10 20 60 100 140 180 220 260

Time

Iteration

1

2

3

4

5

6

5.00E  02 3.00E  01 3.50E  01 4.00E  01 4.50E  01 5.00E  01 1.00E þ 00 3.00E þ 00 5.00E þ 00 7.00E þ 00 9.00E þ 00 1.10E þ 01 1.30E þ 01

el el el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl

1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00

9.11E  29 1.64E  28 1.19E  04 8.32E  04 1.80E  03 2.99E  03 4.34E þ00 5.01E þ 00 5.53E þ00 6.06E þ 00 6.55E þ00 7.01E þ 00 7.43E þ00

1.05E  06 6.77E 06 4.38E  06 3.10E  06 1.06E  01 1.36E  01 1.82E  01 2.41E  01 3.09E  01 3.88E  01 4.74E  01

1.74E  12 3.42E  11 9.84E  12 5.26E  12 1.09E  04 1.64E  04 2.94E  04 5.26E  04 8.85E  04 1.42E  03 2.14E  03

3.23E  23 2.38E  21 3.63E  22 3.25E  23 3.57E  10 6.09E  10 1.80E  09 5.25E  09 1.39E  08 3.57E  08 7.16E  08

2.06E  20 3.50E  20 2.21E  19 1.44E  18 8.27E  18 4.07E 17 1.64E  16

Fig. 5. Cyclic loading program with decreasing levels of loading.

Fig. 6. Tension strip. End stress versus average strain for a loading program in tension–compression with decreasing levels of loading (loading program given in Fig. 5, time increment Δt ¼ 0:05 and t max ¼ 13). The end stress is the traction at the top and bottom boundaries of the section. The average strain is a measure of the strain evaluated as the displacement of the boundary upper edge at vertical axis of symmetry over the height.

18

F. De Angelis, R.L. Taylor / Finite Elements in Analysis and Design 112 (2016) 11–25

Fig. 7. Contour plot of Mises stress in the strip for t max ¼ 13. Loading program in tension–compression with decreasing levels of loading (loading program given in Fig. 5). Fig. 9. Tension strip. End stress versus average strain for a loading program in tension–compression with same levels of loading (loading program given in Fig. 8, time increment Δt ¼ 0:05 and t max ¼ 13). The end stress is the traction at the top and bottom boundaries of the section. The average strain is a measure of the strain evaluated as the displacement of the boundary upper edge at vertical axis of symmetry over the height.

Fig. 10. Contour plot of Mises stress in the strip for t max ¼ 13. Loading program in tension–compression with same levels of loading (loading program given in Fig. 8).

Fig. 8. Cyclic loading program with same levels of loading.

consistent tangent operator and thus the determination of the consistent tangent operator is notably simplified. In fact in evaluating the differential of the normal to the yield surface nðλÞ it is customary to differentiate the expression (31), that is the normal to the yield surface at time t n þ 1 is evaluated as a function of the updated value of the relative stress ΣðλÞ at time t n þ 1 . Conversely, in the present algorithmic scheme we assume nðλÞ as expressed by Eq. (32), that is the normal to the yield surface at time t n þ 1 is evaluated as a function of the trial-like value of the relative λ stress ΣTR ðλÞ. We will show that this strategy simplifies the algorithmic procedure and for the determination of the consistent tangent operator the necessity of evaluating inverse matrices is avoided. By assuming Eq. (32) the differential of the normal nðλÞ to the yield surface is provided by λ

dn ¼ dΣλ n dΣTR

ð56Þ

TR

where dΣλ n is expressed by TR d Σλ n ¼ TR

1 λ

J ΣTR J 2

λ

λ

I J ΣTR J  ΣTR 

!

ΣλTR 1 ¼ ðI  n  nÞ; λ λ J ΣTR J J ΣTR J ð57Þ

and, given the expression (48), it results λ

dΣTR ¼ 2G de þ H λ αn dλ: Eq. (56) is therefore expressed by   1 dn ¼ ðI  n  nÞ 2G de þ H λ αn dλ λ J ΣTR J 2G 1 ðI  n  nÞ de þ ðI  n  nÞH λ αn dλ: ¼ λ λ J ΣTR J J ΣTR J

ð58Þ

ð59Þ

By differentiating Eq. (19) we get ds ¼ 2G de  2G dλn 2Gλ dn;

ð60Þ

and, by recalling the algorithmic plastic consistency condition (54) and the differential of the normal to the yield surface (59), we can write "

2G 2G ds ¼ 2G de  2GA d n  de n  2Gλ I de  ðn  nÞ de λ λ J ΣTR J J ΣTR J # 1 λ α d λ þ ð I  n  n ÞH : ð61Þ n λ J ΣTR J

F. De Angelis, R.L. Taylor / Finite Elements in Analysis and Design 112 (2016) 11–25

19

Fig. 13. Contour plot of Mises stress in the strip for t max ¼ 13. Loading program in tension–compression with increasing mean value of loading (loading program given in Fig. 11).

from which we deduce ! ! 2Gλ 2Gλ ds ¼ 2G 1  A de þ 2G  d ðn  nÞ de λ λ J ΣTR J J ΣTR J

Fig. 11. Cyclic loading program with increasing mean value of loading.



2Gλ λ

J ΣTR J

H λ A d ðαn  nÞ de þ

2Gλ λ

J ΣTR J

H λ A d ðn  nÞðαn  nÞ de: ð64Þ

By recalling that ðn  nÞðαn  nÞ ¼ ðn  αn Þðn  nÞ, the above equation specializes to ! " ! 2Gλ 2Gλ A I þ 2G  ds ¼ 2G 1  d ðn  nÞ λ λ J ΣTR J J ΣTR J # 2Gλ λ 2Gλ λ ð65Þ  H A d ðα n  n Þ þ H A d ðn  αn Þðn  nÞ de λ λ J ΣTR J J ΣTR J

Fig. 12. Tension strip. End stress versus average strain for a loading program in tension–compression with increasing mean value of loading (loading program given in Fig. 11, time increment Δt ¼ 0:05 and t max ¼ 13). The end stress is the traction at the top and bottom boundaries of the section. The average strain is a measure of the strain evaluated as the displacement of the boundary upper edge at vertical axis of symmetry over the height.

By considering the algorithmic plastic consistency condition (54) it follows " 2G I de ds ¼ 2G de  2GA d ðn  nÞ de  2Gλ λ J ΣTR J # # " 2G 2Gλ λ A ð α  n Þ de ; ð n  n Þ de  ð I  n  n ÞH  n d λ λ J ΣTR J J ΣTR J ð62Þ and, accordingly, ds ¼ 2G 1  

2Gλ

2Gλ λ

λ

J ΣTR J

J ΣTR J

! de þ 2G

2Gλ λ

J ΣTR J

!  A d ðn  nÞ de

ðI n  nÞH λ A d ðαn  nÞ de;

ð63Þ

and, accordingly, ! " ( ! 2Gλ 2Gλ A Iþ 2G  ds ¼ 2G 1  d λ λ J ΣTR J J ΣTR J # 2Gλ λ þ H A d ðn  αn Þ ðn  nÞ: λ J ΣTR J " # ) 2Gλ λ H A d ðαn  nÞ de:  λ J ΣTR J

ð66Þ

By considering the position (47), the above equation is expressed by ! " ( ! 2Gλ 2Gλ I þ 2G  Ad ds ¼ 2G 1  λ λ J ΣTR J J ΣTR J #  2 2Gλ qffiffiffi 2 þA d T λ 3 H nl ðn  αn Þ ðn  nÞ: λ J ΣTR J ) " #  2 2Gλ qffiffiffi 2 α  nÞ de: ð67Þ ð H  Ad T λ n nl 3 λ J ΣTR J Consequently, in terms of the total stresses and strains and combined with the elastic volumetric terms the consistent tangent operator is represented by " !# ( " ! 2Gλ 2Gλ þ 2G A I  Ddiscr ¼ Kð1  1Þ þ 2G 1  dev d λ λ J ΣTR J J ΣTR J #  2 2Gλ qffiffiffi 2 þA d T λ 3 H nl ðn  αn Þ ðn  nÞ: λ J ΣTR J

20

F. De Angelis, R.L. Taylor / Finite Elements in Analysis and Design 112 (2016) 11–25

"

)

#

2 2Gλ qffiffiffi 2  Ad Tλ 3 H nl ðαn  nÞ ; λ 

approaches the consistent tangent operator (69) is determined without computing matrix inversions. We also emphasize that the present form of consistent tangent operator differs significantly from the expression presented by DeAngelis and Taylor [11], in which the normal nðλÞ to the yield surface was determined based on the expression (31), that is nðλÞ was computed as a function of the updated value of the relative stress ΣðλÞ at time t n þ 1 . Conversely, in the expression (69) of the consistent tangent operator the normal nðλÞ to the yield surface is evaluated based on the expression (32), that is the normal nðλÞ to the yield surface is computed based on the trial-like value of the λ relative stress ΣTR ðλÞ and it does not need to be updated until convergence is attained. Furthermore, it is worthwhile to note that the expression of the consistent tangent operator presented by DeAngelis and Taylor [11] depends on the updated value of the back stress α as evaluated at time t n þ 1 . Remarkably, both the formulation (69) of the consistent tangent operator and the associated coefficients (70) are expressed as a function of the back stress αn which is the back stress at time tn, that is as evaluated at the previous converged time step. Accordingly, one need not to update it until convergence is achieved. These characteristical features of the proposed algorithmic scheme simplify significantly the numerical integration problem by providing robustness and effectiveness to the overall computational procedure. In fact, as detailed in the sequel, it is emphasized that the proposed computational approach preserves the quadratic rate

ð68Þ

J ΣTR J

which can also be expressed by  h  i h   Ddiscr ¼ Kð1  1Þ þ 2G 1  C d Idev þ 2G C d A d

þ B d ðn  αn Þ ðn  nÞ:  B d ðαn  nÞ ;

ð69Þ

where we have set Cd ¼

2Gλ λ

J ΣTR J

1 T λ ¼ λ; R

;

 2 qffiffiffi B d ¼ A d T λ C d 23 H nl ;

Rλ ¼ ð1 þ

qffiffiffi 2 3

H nl λÞ;

2G Ad ¼  qffiffiffi qffiffiffi  2  2 : λ 2 2 2 2G þ 3 Hiso þ 3 H kin T  3 H nl T λ ðn  αn Þ  23 H kin 23 H nl T λ λ

ð70Þ Given the formulations of the coefficients (70), the expression (69) represents the consistent tangent operator associated to the proposed algorithmic scheme. Accordingly, an alternative procedure is proposed for the consistent linearization of elasto-plastic boundary value problems. It is emphasized that the adopted algorithmic strategy has resulted in the expression of the consistent tangent operator in a simplified manner. In fact, at variance with other

Table 3 Tension strip. Relative residual norm values and iterations. Cyclic loading program in tension–compression with increasing mean value of loading (loading program given in Fig. 11, time increment Δt ¼ 0:05 and t max ¼ 13). Relative residual norm values Load step 1 6 7 8 9 10 20 60 100 140 180 220 260

Time

Iteration

1

2

3

4

5

6

5.00E  02 3.00E  01 3.50E  01 4.00E  01 4.50E  01 5.00E  01 1.00E þ 00 3.00E þ 00 5.00E þ 00 7.00E þ 00 9.00E þ 00 1.10E þ 01 1.30E þ 01

el el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl

1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00

4.19E  14 9.53E  03 5.46E  02 7.21E  02 8.00E  02 8.19E  02 6.43E  01 5.86E  01 6.49E  01 5.72E 01 6.56E  01 5.43E  01 6.63E  01

3.99E  05 6.37E  03 1.07E 02 1.01E  02 7.49E  03 1.35E 01 8.07E 02 1.40E  01 7.58E  02 1.58E  01 6.77E 02 1.90E  01

8.03E  10 2.36E  06 9.40E  05 2.00E  05 7.69E  06 8.26E  03 4.29E  03 6.48E  03 1.41E  03 7.49E  03 3.51E  03 9.67E 03

5.16E  12 5.95E  09 2.73E  10 3.86E  11 5.92E  05 8.33E 05 4.25E  05 2.69E  07 5.91E  05 1.40E  04 9.08E  05

2.73E 09 4.64E  10 1.62E  09 4.92E  13 3.51E  09 6.70E 10 8.94E  09

Table 4 Tension strip. Relative energy norm values and iterations. Cyclic loading program in tension–compression with increasing mean value of loading (loading program given in Fig. 11, time increment Δt ¼ 0:05 and t max ¼ 13). Relative energy norm values Load step 1 6 7 8 9 10 20 60 100 140 180 220 260

Time

Iteration

1

2

3

4

5

6

5.00E  02 3.00E  01 3.50E  01 4.00E  01 4.50E  01 5.00E  01 1.00E þ00 3.00E þ00 5.00E þ00 7.00E þ00 9.00E þ00 1.10E þ01 1.30E þ01

el el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl

1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00 1.00E þ 00

9.92E  29 1.70E 05 5.96E  04 1.38E  03 2.77E 03 4.49E  03 6.64E þ 00 2.10E þ 00 7.29E þ 00 1.71E þ 00 8.27Eþ00 1.39E þ 00 9.81E þ00

1.80E  10 2.16E  06 3.01E  05 1.20E  05 7.86E  06 3.48E  01 1.47E  02 3.62E  01 8.78E  03 4.97E  01 5.55E  03 7.53E  01

7.57E  20 4.82E  13 1.26E  09 6.26E  11 2.26E  11 1.20E  03 2.40E  06 1.10E  03 4.29E  07 1.98E  03 1.14E  06 4.18E  03

2.37E  24 5.21E  18 7.51E  21 5.10E  22 2.27E  08 4.05E  10 1.69E  08 8.24E 15 4.86E  08 1.18E  09 1.74E  07

1.93E  17 1.35E 20 9.24E 18 6.94E  27 5.97E  17 3.92E  20 5.37E  16

F. De Angelis, R.L. Taylor / Finite Elements in Analysis and Design 112 (2016) 11–25

21

Fig. 16. Contour plot of the vertical displacements of the membrane at time t max ¼ 13 (cyclic loading program given in Fig. 11).

Fig. 14. Cook's membrane problem. Finite element mesh.

Fig. 17. Contour plot of the horizontal displacements of the membrane at time t max ¼ 13 (cyclic loading program given in Fig. 11).

Fig. 15. Cook's membrane problem. Bending and shearing dominated response. Total load versus vertical displacement of the upper corner of the free end of the membrane (cyclic loading program given in Fig. 11, time increment Δt ¼ 0:05 and t max ¼ 13).

of asymptotic convergence for the global solution problem typical of computationally efficient numerical integration schemes. Accordingly, the present numerical procedure does not require the use of specific numerical methods for increasing the rate of convergence of the local problem since the solution of the nonlinear scalar equation is provided in exact closed form, that is with no recourse to iterative methods. In addition, the consistent tangent operator (69) ensures a quadratic rate of asymptotic convergence for the solution of the global problem. Consequently, the proposed algorithmic scheme naturally leads to a robust and effective

Fig. 18. Contour plot of the 11-stress in the membrane at time t max ¼ 13 (cyclic loading program given in Fig. 11).

22

F. De Angelis, R.L. Taylor / Finite Elements in Analysis and Design 112 (2016) 11–25

Fig. 19. Contour plot of the 22-stress in the membrane at time t max ¼ 13 (cyclic loading program given in Fig. 11).

Fig. 21. Contour plot of the Mises stress in the membrane at time t¼ 5 (cyclic loading program given in Fig. 11).

Fig. 20. Contour plot of the 12-stress in the membrane at time t max ¼ 13 (cyclic loading program given in Fig. 11).

Fig. 22. Contour plot of the Mises stress in the membrane at time t¼ 5.5 (cyclic loading program given in Fig. 11).

computational approach for elasto-plastic boundary value problems subject to complex and elaborate loading conditions.

of the norms are the norm values evaluated at the i-th iteration divided by the corresponding value of the norm at the first iteration. 6.1. Boundary value problem: tension strip

6. Numerical examples In the present section we show that the proposed algorithmic scheme provides a robust and effective computational procedure for the solution of the global problem by preserving the quadratic rate of asymptotic convergence typical of computationally efficient numerical integration schemes. Accordingly, the proposed algorithmic approach represents an efficient computational scheme for elasto-plastic boundary value problems subject to complex loading conditions. The numerical scheme has been implemented into a Finite Element Analysis Program (FEAP), see Zienkiewicz et al. [40] and Taylor [37]. The algorithmic scheme presented in this paper may be adopted with any 2 or 3-d element that accepts a strain-driven approach. In our analyses we use 4-node quadrilateral elements with 2  2 Gaussian quadrature and the mixed approach described in Sec. 2.6.2 of Zienkiewicz et al. [40]. The convergence of the finite element solution is measured in terms of the residual norm and the energy norm, where the energy norm is computed by the residual and the incremental displacement vector. For the i-th load step the relative values

We consider the boundary value problem of a tension strip. The example represents a plane strain problem of a square strip with a circular hole, with a ratio of the radius of the circle over half of the side length equal to r/h ¼0.1. The geometry of the problem is illustrated in Fig. 1. The radius of the circular hole is r ¼10 mm. The dimension of the side of the square section is assumed to be L¼2h ¼200 mm. For symmetry reasons only 1/4 of the strip needs to be considered. The dimension of the side of the section containing 1/4 of the strip is therefore h¼ 100 mm. The adopted mesh consists of 992 nodes and 900 elements. The material properties are assumed to be elastic modulus E ¼ 208 000 MPa, Poisson's ratio ν ¼0.3, yield limit σ y o ¼ 170 MPa, kinematic hardening modulus H kin ¼ 41 080 MPa, nonlinear kinematic hardening parameter H nl ¼ 525, isotropic hardening modulus H iso ¼ 2100 MPa. The load is applied by prescribing the distributed traction at the top and bottom boundaries of the section. The imposed traction has been assigned as qo ¼220 MPa. A proportional load coefficient p(t) amplifies the prescribed traction and describes the loading history according to the time law qðtÞ ¼ pðtÞqo .

F. De Angelis, R.L. Taylor / Finite Elements in Analysis and Design 112 (2016) 11–25

In the sequel the plots of the end stress versus the average strain are illustrated. The end stress is the distributed traction at the top and bottom boundaries of the section and the average strain is a measure of the strain evaluated as the displacement of the boundary upper edge at vertical axis of symmetry over the height. The contour plots of the Mises stress in the strip are also pffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffi reported, where the Mises stress is defined as 3J 2 ¼ 32 J s J 2 . The

Fig. 23. Contour plot of the Mises stress in the membrane at time t max ¼ 13 (cyclic loading program given in Fig. 11).

23

contour plots of the Mises stress are plotted for the whole strip by exploiting symmetry properties.

6.1.1. Cyclic loading program in tension–compression with increasing levels of loading We assume a complex loading program in tension–compression with increasing both levels of loading. The loading history adopted in this example is illustrated in Fig. 2. The corresponding curve of the end stress versus the average strain is reported in Fig. 3, and the related contour plot of the Mises stress in the strip at time t max ¼ 13 is illustrated in Fig. 4. For the given cyclic loading program the iterations required for convergence and the relative residual norm values are reported in Table 1. For a complete evaluation of the convergence the relative energy norm values are also reported in Table 2. The tables clearly illustrate the quadratic rate of asymptotic convergence of the algorithmic procedure.

6.1.2. Cyclic loading program in tension–compression with decreasing levels of loading A complex loading program in tension–compression with decreasing both levels of loading is considered by assuming the loading history described in Fig. 5. The end stress versus the average strain is plotted in Fig. 6 and the final contour plot of the Mises stress in the strip for t max ¼ 13 is illustrated in Fig. 7.

Table 5 Cook's membrane problem. Relative residual norm values and iterations. Cyclic loading program with increasing mean value of loading (loading program given in Fig. 11, time increment Δt ¼ 0:05 and t max ¼ 13). Relative residual norm values Load step 1 6 7 8 9 10 20 60 100 140 180 220 260

Time

Iteration

1

2

3

4

5

6

5.00E  02 3.00E  01 3.50E  01 4.00E  01 4.50E  01 5.00E  01 1.00E þ 00 3.00E þ 00 5.00E þ 00 7.00E þ 00 9.00E þ 00 1.10E þ 01 1.30E þ 01

el el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl

1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00

5.53E  12 3.38E  01 8.02E  01 1.78E  00 2.13E þ 00 4.46E þ 00 8.04E þ00 7.36E þ 00 7.92E þ 00 7.05E þ 00 7.86E þ 00 6.71E þ 00 7.84E þ 00

2.88E  04 2.78E  02 9.93E  02 1.10E  01 4.22E  01 4.54Eþ 00 1.66E þ 00 3.39Eþ 00 1.75E þ00 3.52Eþ 00 1.51E þ 00 4.16E þ00

2.71E  10 2.51E  05 3.38E  04 2.74E  04 3.99E  03 4.89E  01 6.01E  02 2.40E  01 5.43E  02 2.24E 01 2.27E 02 2.97E 01

2.37E  11 4.41E  09 2.85E  09 7.72E  05 3.84E  03 4.68E  03 1.29E  03 1.44E  05 3.82E  03 2.89E  06 2.01E  03

1.46E  11 2.54E  07 4.11E  07 3.20E  08 4.58E  11 1.03E  07 8.11E  11 1.39E  07

7

5.42E  11 1.65E  11

1.24E 10

Table 6 Cook's membrane problem. Relative energy norm values and iterations. Cyclic loading program with increasing mean value of loading (loading program given in Fig. 11, time increment Δt ¼ 0:05 and t max ¼ 13). Relative energy norm values Load step 1 6 7 8 9 10 20 60 100 140 180 220 260

Time

Iteration

1

2

3

4

5

6

5.00E  02 3.00E  01 3.50E  01 4.00E  01 4.50E  01 5.00E  01 1.00E þ 00 3.00E þ00 5.00E þ00 7.00E þ 00 9.00E þ00 1.10E þ 01 1.30E þ01

el el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl el–pl

1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00 1.00E þ00

4.44E  25 2.04E  04 1.26E  03 7.56E  03 1.29E  02 3.92E  02 9.95E þ 00 2.34E þ 00 9.74E þ00 1.73E þ 00 1.01E þ 01 1.31E þ00 1.10E þ 01

6.83E  11 1.17E  06 1.77E 05 2.40E  05 2.08E  04 9.11E  01 1.94E  02 7.75E  01 1.17E  02 9.00E  01 7.08E  03 1.18E þ 00

9.10E  23 9.31E  13 2.18E  10 1.86E  10 2.80E  08 5.62E  03 4.95E  06 3.94E  03 3.09E  06 5.55E  03 5.48E  07 9.58E  03

8.00E  25 3.66E  20 2.04E  20 5.97E 12 2.32E  07 3.03E  08 1.12E  07 1.78E  13 2.34E  07 6.92E  15 6.05E  07

1.91E  25 5.78E  16 1.35E 16 1.24E 16 4.26E  25 4.95E  16 1.32E  24 3.01E  15

7

6.75E  25 5.13E  26

3.79E  24

24

F. De Angelis, R.L. Taylor / Finite Elements in Analysis and Design 112 (2016) 11–25

6.1.3. Cyclic loading program in tension–compression with same levels of loading Let us consider a loading program in tension–compression with same levels of loading by assuming the loading history illustrated in Fig. 8. For the assigned loading program the end stress versus the average strain is described in Fig. 9 and the final contour plot of the Mises stress in the strip for t max ¼ 13 is given in Fig. 10. 6.1.4. Cyclic loading program with increasing mean value of loading We now test the algorithmic procedure by assuming a loading program with increasing mean value of loading and consider the loading history described in Fig. 11. For the assigned loading program the corresponding curve of the end stress versus the average strain is illustrated in Fig. 12 and the final contour plot of the Mises stress in the strip for t max ¼ 13 is plotted in Fig. 13. For the assigned loading program the iterations required for convergence and the relative residual norm values are reported in Table 3. The relative energy norm values are also reported in Table 4. Accordingly, the quadratic rate of asymptotic convergence of the adopted algorithmic procedure is clearly shown also for this type of cyclic loading program. 6.2. Boundary value problem: Cook's membrane problem In order to show the performance of the proposed algorithmic procedure in a bending and shearing dominated response, we consider the numerical test of a tapered panel clamped on one end and subjected to an in-plane shearing load on the free end. The geometry of the problem is illustrated in Fig. 14. The height of the clamped end is 44 mm. The height of the free end is 16 mm. The horizontal projection of the membrane is 48 mm. By setting a reference system centered in the lower corner of the clamped end, the coordinates of the lower corner of the free end are (48, 44), the coordinates of the upper corner of the free end are (48, 60), the coordinates of the upper corner of the clamped end are (0, 44). We assumed 40 elements in x1-direction and 40 elements in x2-direction. Accordingly, the adopted mesh consists of 1600 elements and 1681 nodes. The adopted material properties are elastic modulus E¼208 000 MPa, Poisson's ratio ν ¼0.3, yield limit σ y o ¼ 170 MPa, kinematic hardening modulus H kin ¼ 41080 MPa, nonlinear kinematic hardening parameter H nl ¼ 525, isotropic hardening modulus H iso ¼ 2100 MPa. The loading is applied on the free end of the membrane by prescribing a uniformly distributed tangential load initially directed upwards. The imposed distributed shear has been assigned as so ¼100 MPa. The proportional load coefficient p(t) modifies in time the imposed distributed tangential load. Accordingly, the loading history is described by the time law sðtÞ ¼ pðtÞso . In the sequel we report the plot of the applied total load versus the vertical displacement of the upper corner of the free end of the membrane, the contour plots of the displacements and the contour plots of the stresses in the membrane. 6.2.1. Cyclic loading program with increasing mean value of loading The algorithmic procedure is tested by assigning a cyclic loading program with increasing mean value of the loading. In the numerical test we assume the loading history as described in Fig. 11. For the assigned cyclic loading program the plot of the applied total load versus the vertical displacement of the upper corner of the free end of the membrane is reported in Fig. 15. The contour plot of the vertical displacements of the membrane subjected to the distributed tangential loading on the free end is illustrated in Fig. 16 at time t max ¼ 13. The contour plot of the horizontal displacements of the membrane at time t max ¼ 13 is shown in Fig. 17. The contour plot of the 11-stress in the membrane at time t max ¼ 13 is illustrated in Fig. 18. The contour plot of the 22-stress in the membrane at time t max ¼ 13 is reported in Fig. 19. The contour

plot of the 12-stress in the membrane at time t max ¼ 13 is shown in Fig. 20. The contour plot of the Mises stress in the membrane at time t¼5 is illustrated in Fig. 21. The contour plot of the Mises stress at time t¼5.5 is reported in Fig. 22. The final contour plot of the Mises stress in the membrane at time t max ¼ 13 is shown in Fig. 23. For the examined problem the iterations required for convergence and the relative residual norm values are illustrated in Table 5. The relative energy norm values are reported in Table 6 for a complete evaluation of the convergence. The adopted algorithmic procedure is characterized by a quadratic rate of asymptotic convergence. Consequently, the robustness and effectiveness of the presented algorithmic procedure is clearly shown.

7. Conclusions In recent years constitutive models have been proposed in the literature with the aim of reproducing more refined characteristical material properties in the loading processes. Accordingly, standard return mapping algorithms require the suitable modifications in order to account for complex and elaborate loading conditions by preserving robustness and efficiency in the numerical integration of elasto-plastic boundary value problems. As a matter of fact finite element methods may require extensive computations in the structural analysis of nonlinear elasto-plastic boundary value problems, in particular when analyzing structural problems subject to complex and elaborate loading conditions. In the present paper a robust computational procedure has been illustrated which is well suited for applications to elasto-plastic boundary value problems subject to elaborate loading conditions. The proposed algorithmic scheme reduces the local constitutive problem to only one nonlinear scalar equation. In addition, in the present approach such nonlinear scalar equation has been expressed in a particularly simple form since it has been downsized furtherly to a single variable algebraic equation. The straightforwardness of such nonlinear scalar equation has allowed to finding the analytical solutions of the algebraic equation in a closed form. Accordingly, in the present approach no iterative method is required for the solution of the local constitutive problem. A computational scheme has been presented which is characterized by an alternative expression of the algorithmic plastic consistency condition. The algorithmic scheme is furtherly characterized by an alternative procedure for the consistent linearization of the elasto-plastic boundary value problem. In fact, within the present approach the related expression of the consistent tangent operator is determined so that, remarkably, matrix inversions are not required. Accordingly, one of the advantages of the present computational approach is the determination of the consistent tangent operator in a simplified form without the necessity of computing inverse matrices. It has also been noted that in the proposed approach no numerical method is needed to accelerate or improve the rate of convergence for the local solution problem. In fact in the present algorithmic scheme the solution of the nonlinear scalar equation is provided in exact closed form. Thus, no iterative method is required for the solution of the local constitutive problem. The present computational procedure together with the derived consistent tangent operator has shown to provide a quadratic rate of asymptotic convergence when used with a Newton iterative method for the global solution of nonlinear elasto-plastic boundary value problems subject to complex loading conditions. The process of convergence for the global iterative solution procedure of elasto-plastic boundary value problems naturally results in a quadratic rate of asymptotic convergence typical of computationally effective iterative solution schemes.

F. De Angelis, R.L. Taylor / Finite Elements in Analysis and Design 112 (2016) 11–25

Finally, the efficiency of the proposed computational scheme has been shown with regard to specific numerical examples. In the numerical examples the robustness and the effectiveness of the algorithmic procedure have been illustrated for elasto-plastic boundary value problems subject to different types of loading conditions. In perspective possible future research work is the generalization of the proposed algorithmic framework to other constitutive models. In fact a notable feature of the present algorithmic scheme is the capability of being extensible to more refined elasto-plastic constitutive models. The present algorithmic framework can also comprehend different constitutive models characterized by the same type of kinematic hardening rules, such as elastoviscoplasticity and damage constitutive models.

Acknowledgments The present research work was carried out while the first author was a Visiting Scholar at the Department of Civil and Environmental Engineering, University of California Berkeley. The support of a travel fellowship from the School of Sciences and Technology of the University of Naples Federico II is gratefully acknowledged. The authors wish to thank the reviewers of the paper for their comments and valuable suggestions.

References [1] M. Abramowitz, Solutions of quartic equations, in: M. Abramowitz, I. Stegun (Eds.), Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables, Dover, New York, 1972, pp. 17–18 §3.8.3. [2] G. Alfano, F. DeAngelis, L. Rosati, General solution procedures in elasto/viscoplasticity, Comput. Methods Appl. Mech. Eng. 190 (2001) 5123–5147. [3] P. Armstrong, C. Frederick, A Mathematical Representation of the Multiaxial Baushinger Effect, Technical Report, Central Electricity Generating Board, Berkeley, UK, Report RD/B/N731, 1966. [4] E. Artioli, F. Auricchio, L.B. da Veiga, Second-order accurate integration algorithms for the von-Mises plasticity with a nonlinear kinematic hardening mechanism, Comput. Methods Appl. Mech. Eng. 196 (2007) 1827–1846. [5] F. Auricchio, R. Taylor, Two material models for cyclic plasticity: non-linear kinematic hardening and generalized plasticity, Int. J. Plast. 11 (1995) 65–98. [6] W. Beyer, CRC Standard Mathematical Tables, 28th edition, CRC Press, Boca Raton, FL, 1987. [7] J. Chaboche, A review of some plasticity and viscoplasticity constitutive theories, Int. J. Plast. 24 (2008) 1642–1693. [8] J. Chaboche, G. Cailletaud, Integration methods for complex plastic constitutive equations, Comput. Methods . Appl. Mech. Eng. 133 (1996) 125–155. [9] M. Crisfield, Non-linear Finite Elements of Solids and Structures, vols. 1–2, Wiley and Sons, Chichester, 1991. [10] F. DeAngelis, An internal variable variational formulation of viscoplasticity, Comput. Methods Appl. Mech. Eng. 190 (2000) 35–54. [11] F. De Angelis, R.L. Taylor, An efficient return mapping algorithm for elastoplasticity with exact closed form solution of the local constitutive problem, Eng. Comput. 32 (8) (2015) 2259–2291. [12] I. Doghri, Fully implicit integration and consistent tangent modulus in elastoplasticity, Int. J. Numer. Methods Eng. 36 (1993) 3915–3932. [13] C. Frederick, P. Armstrong, A mathematical representation of the multiaxial Baushinger effect, Mater. High Temp. 24 (1) (2007) 1–26.

25

[14] J.J. Hacke, A simple solution of the general quartic, Am. Math. Mon. 48 (5) (1941) 327–328. [15] S. Hartmann, G. Luhrs, P. Haupt, An efficient stress algorithm with applications in viscoplasticity and plasticity, Int. J. Numer. Methods Eng. 40 (1997) 991–1013. [16] O. Hopperstad, S. Remseth, A return mapping algorithm for a class of cyclic plasticity models, Int. J. Numer. Methods Eng. 38 (1995) 549–564. [17] T. Hughes, R. Taylor, Unconditionally stable algorithms for quasi-static elasto/ viscoplastic finite element analysis, Comput. Struct. 8 (1978) 169–173. [18] M. Kobayashi, N. Ohno, Implementation of cyclic plasticity models based on a general form of kinematic hardening, Int. J. Numer. Methods Eng. 53 (2002) 2217–2238. [19] M. Kojic, K. Bathe, Inelastic Analysis of Solids and Structures, Springer, New York, NY, 2010. [20] R. Krieg, S. Key, Implementation of a time dependent plasticity theory into structural computer programs, in: J. Stricklin, K. Saczalski (Eds.), Constitutive Equations in Viscoplasticity: Computational and Engineering Aspects, ASME (AMD-20), New York, 1976. [21] R. Krieg, D. Krieg, Accuracies of numerical solution methods for the elastic– perfectly plastic model, J. Press. Vessel Technol. Trans. ASME 99 (1977) 510–515. [22] J. Lubliner, R. Taylor, F. Auricchio, A new model of generalized plasticity and its numerical implementation, Int. J. Solids Struct. 30 (1993) 3171–3184. [23] G. Maenchen, S. Sack, The tensor code, in: B. Alder (Ed.), Methods in Computational Physics, vol. 3, Academic Press, New York, 1964, pp. 181–210. [24] J. Nagtegaal, On the implementation of inelastic constitutive equations with special reference to large deformation problems, Comput. Methods Appl. Mech. Eng. 33 (1982) 469–484. [25] P. Nukala, A return mapping algorithm for cyclic viscoplastic constitutive models, Comput. Methods Appl. Mech. Eng. 195 (2006) 148–178. [26] M. Ortiz, E. Popov, Accuracy and stability of integration algorithms for elastoplastic constitutive relations, Int. J. Numer. Methods Eng. 21 (1985) 1561–1576. [27] P. Papadopoulos, R. Taylor, On the application of multi-step integration methods to infinitesimal elasto-plasticity, Int. J. Numer. Methods Eng. 37 (1994) 3169–3184. [28] D. Peric, On a class of constitutive equations in viscoplasticity: formulation and computational issues, Int. J. Numer. Methods Eng. 36 (1993) 1365–1393. [29] W. Prager, Recent developments in the mathematical theory of plasticity, J. Appl. Phys. 20 (3) (1949) 235–241. [30] J. Simo, Nonlinear stability of the time-discrete variational problem of evolution in nonlinear heat conduction, plasticity and viscoplasticity, Comput. Methods Appl. Mech. Eng. 88 (1991) 111–131. [31] J. Simo, S. Govindjee, Exact closed-form solution of the return mapping algorithm in plane stress elasto-viscoplasticity, Eng. Comput. 5 (3) (1988) 254–258. [32] J. Simo, S. Govindjee, Non-linear B-stability and symmetry preserving return mapping algorithms for plasticity and visco-plasticity, Int. J. Numer. Methods Eng. 31 (1991) 151–176. [33] J. Simo, T. Hughes, Computational Inelasticity, Springer-Verlag, New York, 1998. [34] J. Simo, J. Kennedy, S. Govindjee, Non-smooth multisurface plasticity and viscoplasticity, loading/unloading conditions and numerical algorithms, Int. J. Numer. Methods Eng. 26 (1988) 2161–2185. [35] J. Simo, R. Taylor, Consistent tangent operators for rate-independent elastoplasticity, Comput. Methods Appl. Mech. Eng. 48 (1985) 101–118. [36] E. de Souza Neto, D. Peric, D. Owen, Computational Methods for Plasticity: Theory and Applications, Wiley-Blackwell, Chichester, 2008. [37] R. Taylor, A Finite Element Analysis Program, Technical Report, University of California at Berkeley, User Manual, Version 8.4 〈http://www.ce.berkeley.edu/ feap〉, 2014. [38] M. Wilkins, Calculation of elastic plastic flow, in: B. Alder (Ed.), Methods in Computational Physics, vol. 3, Academic Press, New York, 1964, pp. 211–263. [39] P. Wriggers, Nonlinear Finite Element Methods, Springer-Verlag, Berlin, 2008. [40] O. Zienkiewicz, R. Taylor, D. Fox, The Finite Element Method for Solid and Structural Mechanics, 7th edition, Elsevier, Oxford, 2013.