Critical softening in Cam-Clay plasticity: Adaptive viscous regularization, dilated time and numerical integration across stress–strain jump discontinuities

Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133

Contents lists available at SciVerse ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Critical softening in Cam-Clay plasticity: Adaptive viscous regularization, dilated time and numerical integration across stress–strain jump discontinuities R. Conti a,⇑, C. Tamagnini b, A. DeSimone a a b

SISSA-International School for Advanced Studies, Via Bonomea 265, 34136 Trieste, Italy Università degli Studi di Perugia, Dipartimento di Ingegneria Civile e Ambientale, Via G. Duranti 93, Perugia, Italy

a r t i c l e

i n f o

Article history: Received 3 August 2012 Received in revised form 12 December 2012 Accepted 7 February 2013 Available online 17 February 2013 Keywords: Cam-Clay Elastoplasticity Viscoplasticity Instability Subcritical softening

a b s t r a c t Within the framework of continuum mechanics, the mechanical behaviour of geomaterials is often described through rate-independent elastoplasticity. In this field, the Cam-Clay models are considered as the paradigmatic example of hardening plasticity models exhibiting pressure dependence and dilation-related hardening/softening. Depending on the amount of softening exhibited by the material, the equations governing the elastoplastic evolution problem may become ill-posed, leading to either no solutions or two solution branches (critical and sub-critical softening). Recently, a method was proposed to handle subcritical softening in Cam-Clay plasticity through an adaptive viscoplastic regularization for the equations of the rate-independent evolution problem. In this work, an algorithm for the numerical integration of the Cam-Clay model with adaptive viscoplastic regularization is presented, allowing the numerical treatment of stress–strain jumps in the constitutive response of the material. The algorithm belongs to the class of implicit return mapping schemes, slightly rearranged to take into account the rate-dependent nature of inelastic deformations. Applications of the algorithm to standard axisymmetric compression tests are discussed. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction Many granular materials, such as rocks, concrete, dense sands and stiff clays, exhibit strain softening when subjected to intense shear deformation. This behaviour, which usually manifests itself as a macroscopic loss of strength after a peak load level has been reached, may lead to both spatial (strain localization) and time (critical softening) discontinuities in the response of the material. The first phenomenon occurs when deformations in solids localize into narrow bands [40], while critical softening is related to the loss of test controllability and the onset of snap-back processes during displacement controlled tests [43,17]. In this condition, jumps in the stress–strain behaviour of the material occur and the response of the material evolves faster than the applied loading conditions. Within the framework of continuum mechanics, the mechanical behaviour of geomaterials is often described through rate-independent elastoplasticity. In this field, the Cam-Clay family of models [32,34,31] are considered as the paradigmatic example of hardening plasticity models exhibiting pressure dependence and dilation-related hardening/softening. As far as rate-independent ⇑ Corresponding author. E-mail address: [email protected] (R. Conti). 0045-7825/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cma.2013.02.002

elastoplasticity is concerned, both localization and critical softening phenomena have been recognized to emerge as local instabilities in the constitutive equations, leading to the ill-posedness of the evolution problem and a non-uniqueness in the incremental response of the material. Strain localization has been widely analyzed in the literature both from a mathematical [30,3] and a numerical [21,24] point of view. The onset of localization is characterized by the loss of hyperbolicity of the dynamical equations of motion, resulting in the fact that the wave speeds vanish or become imaginary. As a consequence, numerical solutions of localization problems carried out by adopting rate-independent plasticity models suffer from pathological mesh dependence and length-scales effects. Many strategies have been proposed in the literature to model strainsoftening behaviour even beyond the onset of localization. Among these regularization techniques, viscoplasticity has been recognized to provide a valuable framework for the analysis of strain localization in solids [21,24,42,29]. One main obstacle in applying viscoplastic regularization techniques to materials exhibiting strong softening response is the possible occurrence of time discontinuities, which emerge in critical softening conditions. As pointed out by many authors [23,19,7], the inception of critical softening in a strain controlled process corresponds to the vanishing of the determinant of the elastoplastic compliance

119

R. Conti et al. / Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133

matrix or, equivalently, to the vanishing of the so-called plastic modulus, K p . In this condition, for a given strain rate e_ , the corresponding stress rate r_ goes to infinity and uniqueness of the solution for the rate-independent elastoplastic evolution problem is lost. It is worth noting that, in the literature, the plastic modulus is always assumed to be strictly positive a priori. This assumption, which ensures the positiveness of the consistency parameter during a plastic process and the well-posedness of the rate-independent evolution problem, clearly poses a restriction on the allowable amount of material softening that the constitutive model can handle. Only recently, the possibility to guarantee the well-posedness of the evolution problem even beyond the onset of critical conditions has been investigated, with reference to Cam-Clay plasticity [9,12,11]. It was demonstrated that, by adopting a Duvaut–Lions viscoplastic regularization [15] for the equations of the rate-independent evolution problem, a unique solution exists in the whole space of the admissible stress states [10]. The heuristic idea is that, in the limit as the viscous regularization parameter, s, goes to zero, two different regimes may occur during plastic loading, depending on the sign of the plastic modulus. For all initial conditions corresponding to which the plastic modulus remains positive, and then the inviscid equations are well-posed, the viscous solution is continuous and tends to the solution of the rate-independent problem as s ! 0. In this regime, which we will refer to as slow dynamics, the response of the material takes place at the same time scale at which the loading conditions evolve. On the other hand, if the plastic modulus vanishes or becomes negative, the limit solution for the viscoplastic problem is discontinuous. To study the evolution of both the stress state and the internal variable along the jump, i.e. to follow the instantaneous response of the material to the applied deformation, a dilated time s :¼ 1s t may be used to rescale the viscoplastic evolution equations and to compute the state of the material at the end of the jump. We will refer to this rescaled evolution problem as the fast dynamics regime. In this work we focus on the numerical treatment of time instabilities related to critical and subcritical softening in Cam-Clay plasticity models. More specifically, we present an algorithm for the numerical integration of the evolution equations of the regularized viscoplastic Modified Cam-Clay model [31]. The proposed algorithm belongs to the class of implicit return mapping schemes, slightly rearranged to take into account the viscous nature of plastic deformations [37]. Two different strategies have been used to integrate the viscous equations in the slow dynamics regime and their rescaled version along the jumps. As far as the slow dynamics is concerned, the well-conditioned structure of the adopted algorithm allows to recover exactly the return-mapping scheme for the rateindependent problem, in the limit as s ! 0. On the other hand, by exploiting the mathematical properties of the governing equations, the same structure is essentially preserved also for the integration in the fast dynamics regime. In this case, however, time discontinuities (of stress, plastic deformations and internal variable) are generated in the limit as s ! 0. This discussion highlights one of the main differences of our approach with respect to viscoplastic regularization schemes already discussed in the literature [21,42]. While in these cases the (regularizing) viscous perturbation is always present, in our case it is invoked only when critical softening conditions occur, and it is switched off otherwise. For this reason, our regularization strategy can be described as adaptive. A comprehensive description of return mapping strategies for the integration of both rate-independent and rate-dependent elastoplastic evolution problems can be found in the literature (see e.g. [36,13]). A key ingredient in return mapping schemes is the use of an operator split strategy for the original evolution problem, leading to the so-called closest-point projection approximation [36]. Briefly, the final stress and internal variables are computed from

an initial trial elastic state by solving the plastic evolution equations whenever the plastic constraint is activated. By adopting an implicit approximation of the governing equations, the evolution problem is then reduced to a system of non-linear algebraic equations, which is solved customarily with a Newton–Raphson iterative procedure. As far as the rate-independent problem is concerned, the unknowns to be computed during the integration are the stress (or elastic strains), the internal variables and the discrete consistency parameter, which can no longer be assumed positive a priori. Just as for the continuum case, well-posedness of standard numerical algorithms is guaranteed in this case only as long as the discrete plastic multiplier is positive. In the rate-dependent problem, the stress state can lie outside the yield surface during a plastic process and the consistency condition is no longer necessary. However, by exploiting the variational structure of the rate-dependent Cam-Clay constitutive equations for the slow dynamics regime, it is possible to show that an analogous Lagrange multiplier emerges naturally from the corresponding return mapping scheme, see Section 5. Moreover, the well-conditioned structure of the algorithm allows us to show that this Lagrange multiplier tends to the rate-independent consistency parameter in the limit as s ! 0. As for the continuum case, jumps in the viscous solution can be detected in the numerical algorithm by the sign of the Lagrange multiplier: as long as the discrete Lagrange multiplier is positive, the return mapping scheme for the slow dynamics remains wellposed and the viscous solution evolves without discontinuities; on the other hand, if the solution converges to a negative value of the Lagrange multiplier during a time step, then a jump occurs in the stress–strain response of the material and the equations of the fast dynamics have to be integrated to follow the evolution of the stress state along the discontinuity; at the end of the jump, the discrete Lagrange multiplier becomes again positive and the solution returns to evolve in slow dynamics conditions. The paper is organized as follows. After the introduction of the main notation conventions (Section 2), in Section 3 we recall the existence and uniqueness conditions for the solution of the rateindependent elastoplastic evolution problem. In Section 4 the Modified Cam Clay model is briefly reviewed and properties of the solution for its viscoplastic regularization are outlined. The return mapping algorithm we propose for the numerical integration of the regularized viscoplastic model is presented in Section 5. Results from numerical tests are given in Section 6, where the ability of the algorithm to handle instabilities due to strain softening is assessed by means of a series of single-element tests. 2. Notation The usual sign convention of soil mechanics (compression positive) is adopted throughout. Following standard notation, bold-face letters denote vectors and second-order tensors, while blackboardbold symbols denote fourth-order tensors. Accordingly, I and I are the second-order and the fourth-order identity tensor respectively. For any two vectors v ; w the scalar product is defined as: v  w :¼ v i wi , and the dyadic product as: ½v  wij :¼ v i wj . Accordingly, for any two second-order tensors X; Y, X  Y :¼ X ij Y ij and ½X  Yijkl :¼ X ij Y kl . The Euclidean norm of a second order tensor X pffiffiffiffiffiffiffiffiffiffi ffi is defined as kXk :¼ X  X. For any vector field v ; rv denotes the spatial gradient of v . In the representation of stress and strain states, the following invariant quantities will be used in the paper:

p :¼

1 trðrÞ; 3

q :¼

rffiffiffi 3 ksk; 2

S :¼ sinð3hÞ :¼

pffiffiffi trðs3 Þ 6 ½trðs2 Þ3=2

ð1Þ

120

R. Conti et al. / Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133

and

ev :¼ e  I; eq

rffiffiffi 2 kek :¼ 3

which forces the state of the material to remain on the yield surface whenever a plastic loading occurs (f_ ¼ 0 if c_ > 0).

ð2Þ

In Eq. (1) and (2) above,

1 s :¼ r  trðrÞI; 3

1 e :¼ e  trðeÞI 3

ð3Þ

are the deviatoric part of the stress and the strain tensor respectively; s2 and s3 are the square and the cube of the deviatoric stress tensor, whose components are given by ðs2 Þij :¼ sik skj and ðs3 Þij :¼ sik skl slj , and h is the Lode angle. 3. Existence and uniqueness of the solution of the rateindependent elastoplastic evolution problem Within the framework of linearized kinematics (small strain assumption), let u 2 R3 be the displacement vector and e :¼ 12 ðru þ r> uÞ the total strain tensor. Moreover, let r 2 S3 be the stress tensor, where S3 is the space of symmetric tensors in R3 , and let z 2 Rnint be a vector of (nint ) stress-like internal variables, which take into account the effects of previous loading history. The basic assumptions of the rate-independent elastoplastic evolution problem can be summarized as follows [36]: (i) Strain rates are decomposed additively into an elastic, reversible part, e_ e , and a plastic, irreversible part, e_ p , as

e_ ¼ e_ e þ e_ p

ð4Þ

(ii) The rate form of the constitutive equation is

r_ ¼ Ce_ e

ð5Þ

@g ðr; zÞ @r

ð6Þ

where g : S3  Rnint # R is the plastic potential, and c_ is a non-negative scalar parameter, customarily defined as consistency parameter or plastic multiplier. (v) The evolution of the internal variables is defined by the hardening law:

z_ ¼ c_ hðr; zÞ 3

nint

where h : S  R # R is a suitable hardening function of the state of the material. During a plastic process, the yield surface may either expand (hardening behaviour), or shrink (softening behaviour), or remain unaltered (perfectly plastic behaviour at critical state), depending on the sign of the hardening modulus (see e.g. [36]):

@f  z_ @z

(vi) The consistency parameter is subjected to the Kuhn–Tucker complementarity conditions

c_ P 0; f 6 0; c_ f ¼ 0

ð8Þ

according to which plastic deformations may occur only for stress states on the yield surface, and to the consistency condition

c_ f_ ¼ 0

1 @f  Ce_ Kp @r

ð10Þ

where:

K p :¼ H þ

@f @g C ¼ H  Hc @r @r

ð11Þ

is the plastic modulus of the material and Hc is the so-called critical hardening modulus. For an incremental plastic loading starting from an initial stress state on the yield surface ð@@fr  Ce_ > 0Þ, Eq. (10), together with the plastic constraints (8)–(9), implies that the solution of the rate-independent elastoplastic evolution problem exists and is unique as long as H  Hc > 0. As a matter of fact, for admissible stress states ðr; zÞ which do not obey this condition, i.e. for critical ðH  Hc ¼ 0Þ and sub-criticalðH  Hc < 0Þ softening conditions, the evolution equations have either no solutions or two solutions, one elastic and one elastoplastic [23]. As recognized by many authors [23,19,7], the condition H ¼ Hc , which from a mathematical point of view corresponds to the vanishing of the determinant of the material compliance tensor [19], is related to the loss of loading path controllability and the onset of snap-back processes, during which discontinuities (jumps) in the stress–strain response of the material occur.

4.1. The Modified Cam-Clay model Within the framework of Critical State soil mechanics, the Modified Cam-Clay model (MCC) was proposed to describe, from a phenomenological point of view, the observed behaviour of fine-grained soils [31]. The model is characterized by a small number of material constants, easily measured in conventional laboratory tests. Despite its simplicity, the MCC model is able to capture, at least qualitatively, the main features of the mechanical behaviour of clays, i.e. pressure sensitivity and hardening/softening behaviour associated with plastic volumetric compaction/dilation. 4.1.1. Yield function The yield function is defined by

ð7Þ

nint

Hðr; zÞ :¼ 

c_ ¼

4. Cam-Clay plasticity

where C is the fourth-order elastic stiffness tensor. (iii) Admissible stress states ðr; zÞ are constrained to lie in the convex set Er :¼ fðr; zÞjf ðr; zÞ 6 0g, where f : S3  Rnint # R is the yield function. We define Er ð zÞ :¼ frjf ðr;  zÞ 6 0g as the section of Er for z ¼  z. (iv) The evolution of plastic strains is defined by the flow rule:

e_ p ¼ c_

During a strain-driven process, c_ is related to the current strain rate by (see e.g. [36]):

ð9Þ

f ðp; q; zÞ ¼

q2 M2

þ pðp  2zÞ

ð12Þ

where M is a material constant defining the slope of the critical state line in the p : q plane, the internal variable z ¼ pc =2 ðnint ¼ 1Þ is the radius along the hydrostatic axis of the ellipsoid representing the yield surface (f ¼ 0), and pc is the preconsolidation pressure (see Fig. 1). Eq. (12), which controls the shape of the yield function in the triaxial compression range, can be generalized to three-dimensional stress state by replacing the constant M with a suitable function MðSÞ of the third invariant [20]. In this paper, MðSÞ is defined as in [39], starting from the work of van Eekelen [16]:

MðSÞ ¼ Mc c1 ð1 þ c2 SÞn

ð13Þ

where:

c1 :¼

in 1h 1=n ; n 1 þ ðcM Þ 2

c2 :¼

1  ðcM Þ1=n 1 þ ðcM Þ1=n

ð14Þ

R. Conti et al. / Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133

121

where G is the constant shear modulus and K is the pressure dependent bulk modulus, defined as:

(a)

( K :¼

^ p=j

ðp P pr Þ

^ ðp < pr Þ p r =j

ð17Þ

As discussed in [6,39], for an isotropic loading process starting from an arbitrary initial isotropic state ðeev 0 ; p0 P pr Þ, the expression adopted for the elastic stiffness matrix leads to the empirical compression law proposed in [18,8]. 4.1.4. Evolution equations for rate-independent Cam-Clay model The constitutive equations and the evolution laws for the internal variables result in the system

8 e_ ¼ e_ e þ e_ p > > > > < r_ ¼ C½e_  e_ p  > e_ p ¼ c_ @@fr > > > : z_ ¼ qc ztrðe_ p Þ

(b)

ð18Þ

which, together with the Kuhn–Tucker and the consistency conditions, defines completely the constrained evolution problem for the rate-independent MCC model. 4.2. Viscoplastic regularization of Modified Cam-Clay plasticity

Fig. 1. Modified Cam Clay model: (a) yield surface and (b) response in isotropic compression.

In Eqs. (14) n is a material constant and cM ¼ Me =M c is the ratio between the values of MðSÞ in triaxial extension ðh ¼ p=6Þ and compression ðh ¼ p=6Þ respectively. 4.1.2. Flow rule and hardening law An associative flow rule is assumed for the plastic strain rate (g ¼ f ), while a non-associative hardening law is assumed for the internal variable, in the form of

z_ ¼ qc z trðe_ p Þ

ð15Þ

^ Þ. The material constants ^ ^ can be derived where qc ¼ 1=ð^ kj k and j from the interpretation of isotropic compression tests, as the interpolation of the bilogarithmic relations between the specific volume v and the mean effective pressure p for virgin compression and unloading–reloading [18,8] (see Fig. 1). From Eqs. (12) and (15), it follows that the elastic domain expands or shrinks during a loading process depending on the sign of z_ . Moreover, as z > 0, hardening and softening behaviour of the model is determined only by the sign of trðe_ p Þ. By comparing Eq. (15) with Eq. (7), the hardening function takes the form:

  @f z_ ¼ c_ hðr; zÞ h :¼ qc z tr @r 4.1.3. Elastic response The isotropic fourth-order elastic stiffness tensor is defined as:

  1 C :¼ KI  I þ 2G I  I  I 3

ð16Þ

To handle instabilities in softening regime, a viscous approximation of Duvaut–Lions type [15] was proposed in [12,11] for the evolution problem (18). In this context, the stress state is no longer constrained to lie on the yield surface during a plastic process, and the resulting unconstrained system of ODEs for the rate-dependent evolution problem takes the form

8 e_ ¼ e_ e þ e_ p > > > < r_ ¼ C½e_  e_ p  > e_ p ¼ 1s A½r  pA ðrÞ > > : z_ ¼ qc z trðe_ p Þ

ð19Þ

where pA ðrÞ is the minimal distance projection of r onto Er ðzÞ in the metric induced by A; s is the regularization parameter and ð1=sÞA is a fluidity-like tensor. We require A to be an isotropic and invertible fourth-order tensor. Remark 4.1. The third equation in (19), defining the rate at which inelastic deformations evolve during a loading process, represents a generalization of the viscoplastic law proposed by Duvaut and Lions [15] for perfect plasticity and successively by Simo et al. [37] for the case of hardening plasticity. Indeed, the formulations proposed by Duvaut, Lions and Simo et al. correspond to the particular choices A ¼ I and A ¼ C1 , respectively. Remark 4.2. Even in the case A ¼ C1 , the regularized evolution Eqs. (19) are slightly different from the viscoplastic format proposed by Simo et al. and adopted in [33,29,42,21]. Indeed, in the present formulation, the hardening law is not defined in terms of the quasi-static solution of the corresponding rate-independent problem. As discussed in the following, this will results in a different algorithm for the numerical solution of the viscoplastic evolution problem. During a generic loading process, by solving the evolution equations (19) in the limit as s ! 0, the viscous dynamics presents three possible regimes [12]:

122

R. Conti et al. / Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133

(a) Elastic regime. This situation occurs in a time interval ½t 1 ; t 2  when the plastic deformation, and hence the internal variable, do not evolve, while the stress is completely determined through the elastic constitutive relation rðtÞ ¼ rðt1 Þþ CðeðtÞ  eðt1 ÞÞ for every t 2 ½t1 ; t2 . A necessary condition for this behaviour to occur is clearly ðrðtÞ; zðt1 ÞÞ 2 Er for every t 2 ½t 1 ; t 2 . (b) Slow dynamics. For a plastic loading ð@f =@ r  Ce_ > 0Þ starting at time t1 from an initial state ðr1 ; z1 Þ :¼ ðrðt 1 Þ; zðt1 ÞÞ on the yield surface, the stress evolves smoothly and Eqs. (18) for the rate-independent problem are recovered as a limit case. This regime occurs when the plastic modulus is positive:

K p ðr1 ; z1 Þ > 0

K p :¼

@f @f @f C  h @r @ r @z

The evolution problem can be studied adopting the standard time scale t and both hardening and softening behaviour may occur. (c) Fast dynamics. For a plastic loading occurring in softening regime, if K p ðr1 ; z1 Þ 6 0, the solution of system (19) is discontinuous and the stress–strain response of the material exhibits a jump. By introducing a rescaled (dilated) time s :¼ t=s, the evolution of both the stress state and the internal variable along the jump can be recovered from the solution of the problem:

8 < dr ¼ CAðpA ðrðsÞÞ  rðsÞÞ ds : dz ¼ q ztr½AðrðsÞ  p ðrðsÞÞÞ A c ds

ð20Þ

The asymptotic values at s ¼ 1 of the solution of Eqs. (20) define the solution of Eqs. (18) before and after the jump time, that is:

lim ðrðsÞ; zðsÞÞ ¼ ðrðt1 Þ; zðt 1 ÞÞ

s!1

ð21Þ

Eqs. (20), called equations of the fast dynamics, are formally obtained by rescaling time in (19) and neglecting all terms of order s. It is shown in [12] that, along the solution of (20), the internal variable is strictly decreasing (softening regime). At the end of the jump, the stress state lies on the þ yield surface, that is f ðrþ 1 ; z1 Þ ¼ 0; moreover, the plastic þ modulus is positive, K p ðrþ ; z Þ 1 1 > 0, and the solution of (19) evolves following either the elastic regime or the slow dynamics. In the latter case, the material exhibits softening behaviour even after the end of the jumps (see Fig. 2).

The following heuristic argument is helpful to clarify the rationale behind Eqs. (20). Using the chain rule:

d d ¼s ds dt

()

d 1 d ¼ dt s ds

ð22Þ

we can rewrite the equations in (19) as:

1 dr 1 de 1 dep ¼ C  C s ds s ds s ds 1 dep 1 ¼ A½r  pA ðrÞ s ds s  p 1 dz 1 de ¼ qc ztr s ds s ds

ð23Þ ð24Þ ð25Þ

Since

C

de ¼ sCe_ ds

we have, from Eqs. (23) and (24):

dr dep ¼ sCe_  C ¼ sCe_  CA½r  pA ðrÞ ds ds

ð26Þ

Now, taking the limit as s ! 0 in Eqs. (25) and (26) and observing that the term sCe_ ! 0 because e_ (i.e., the prescribed total strain rate) is bounded, we recover Eqs. (20). The chain rule (22) shows that we may obtain nontrivial (fast) dynamics in the dilated time s only at jump times, where the time derivative d=dt becomes infinite. Eq. (26)1 and its limit as s ! 0 shows that stress jumps occur at frozen total strain and are entirely due to jumps in plastic deformation. Finally, in the regular regime of slow dynamics, where:

dep dep ¼s ! 0 as ds dt

s!0

Eq. (24) shows that, in the limit as s ! 0, the evolution takes place on the yield surface, where r ¼ pA ðrÞ. As discussed in [35], the transition between elastic and inelastic regimes can be detected through an initial trial elastic step, just as in the rate-independent evolution problem. On the other hand, the transition between slow and fast dynamics, during a plastic process, depends on the sign of the plastic modulus. From Eq. (10), it follows that we can use the sign of the plastic multiplier, which

Fig. 2. Schematic representation of the stress–strain behaviour in Cam-Clay plasticity before, during and after the jump.

123

R. Conti et al. / Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133

is no longer assumed to be positive a priori, to identify the onset of instabilities in the constitutive response of the material. 5. Numerical integration of rate-dependent MCC model Let ½0; T  R be the time interval of interest, and et the strain tensor history, which is assumed to be given in a standard straindriven problem. We assume that the state of the material ðrn ; zn Þ is known at time t n 2 ½0; T. Let Denþ1 ¼ enþ1  en be the incremental strain at time tnþ1 2 ½0; T. The problem to be addressed is to update the state variables ðrnþ1 ; znþ1 Þ through the integration of the viscoplastic constitutive equations, defining either the slow dynamics – Eqs. (19) – or the fast dynamics – Eqs. (20) – evolution problem, with the initial conditions ðr; zÞjt¼tn ¼ ðrn ; zn Þ. While the integration of the slow dynamics equations is quite straightforward using a standard return mapping scheme, two problems of major concern have to be addressed in the numerical algorithm: (i) the choice of a robust scheme for the integration of the fast dynamics equations, and (ii) the definition of the conditions according to which the transition between the two regimes takes place. Accordingly, the integration of the equations for slow and fast dynamics will be treated separately in the following. As pointed out by many authors (see e.g. [26,2]), the classical Newton–Raphson method, applied to the plastic corrector stage in a return mapping scheme, may not converge for initial trial states far from the final solution, if the constitutive equations exhibit strong non-linearities. For the Cam-Clay model this is indeed the case in regions of the stress space where high softening is expected and, particularly, where subcritical softening occurs and the transition between slow and fast dynamics is detected. To handle such non-linearities, a simple substepping procedure based on time step halving [26] is activated whenever the plastic corrector requires more than a prescribed number of iterations to converge.

rtrnþ1 ¼ rn þ Cn Denþ1

ð27Þ

ztr nþ1 ¼ zn

tr tr If fnþ1 :¼ f ðrtr nþ1 ; znþ1 Þ 6 0, the process is declared elastic and the trial state represents the actual final state of the material (see Box 1).

5.1.2. Plastic corrector tr If fnþ1 > 0, the convexity of the yield surface implies that the trial state lies outside the yield locus, and hence inelastic deformations occur. The plastic corrector problem is solved by integrating Eqs. (19) with an implicit backward Euler scheme, at fixed total strain (e_ ¼ 0), taking the trial state as the new initial condition:

Depnþ1 ¼

Dt

s

 nþ1 Þ Aðrnþ1  r

r

e

r

q

Box 1. Numerical algorithm for the elastic predictor scheme in slow dynamics.

ð30Þ

 nþ1 :¼ pA ðrnþ1 Þ is the minimal distance projection of rnþ1 where r onto the updated elastic domain. Substituting (28) into (29) and (30) we get:

rnþ1 ¼ rtrnþ1  znþ1 ¼ ztr nþ1 þ

Dt

s

Dt

 nþ1 Þ Cn Aðrnþ1  r

ð31Þ

qc znþ1 tr½Aðrnþ1  r nþ1 Þ

s

r nþ1 ¼ arg



1 krnþ1  rk2A f ðr;znþ1 Þ60 2 min

 ð33Þ

which, from standard convex analysis arguments [22], implies that

@f @ r ðr nþ1 ;znþ1 Þ

 nþ1 ; znþ1 Þ ¼ 0 Dcnþ1 f ðr

ð34Þ ð35Þ

where Dcnþ1 P 0 is the Lagrange multiplier for the unconstrained minimization problem associated with (33).  nþ1 can be By adopting a standard return mapping scheme, r computed from the trial stress as

r nþ1 ¼ rtrnþ1  Dcnþ1 X

@f @ r ðr nþ1 ;znþ1 Þ

ð36Þ

where X is a fourth-order tensor defined so that Eq. (34) holds. Using (31) and (36) we get

X¼

Dt

s

FA1 ;

with F ¼

s Dt

I þ Cn A

ð37Þ

Note that, as both A and Cn are invertible and s=Dt P 0, we have that F is an invertible fourth-order tensor. Finally, substituting (37) into (36) and (34) into (32), and using (31), we obtain

f (1) Evaluate elastic predictor: rtrnþ1 ¼ rn þ Cn Denþ1 ztrnþ1 ¼ zn (2) Check yield condition: tr tr fnþ1 :¼ f ðrtr nþ1 ; znþ1 Þ tr if fnþ1 < TOLf then Set ðÞnþ1 ¼ ðÞtr nþ1 , CONV=’YES’ and EXIT else GO TO Box 2 end if

ð32Þ

 nþ1 is characterized as the solution of Remark 5.1. By definition, r the minimization problem

5.1. Slow dynamics

5.1.1. Elastic predictor To handle nonlinearities arising from pressure sensitivity of the bulk modulus, a hybrid algorithm is adopted, in which the elastic stiffness tensor is considered as constant and evaluated at the beginning of each step [41,4,1]. Under this assumption, the trial state is given by:

ð29Þ

e

r nþ1 ¼ rnþ1  Dcnþ1 A1 As in classical rate-independent elasto–plasticity [36], a twostep product formula algorithm is constructed by exploiting the additive structure of Eqs. (19). First, the elastic predictor problem is solved by freezing the plastic flow during the time step, and a trial elastic state is obtained. Then, if the plastic constraint is violated, the trial state is taken as the initial condition for the plastic corrector problem.

ð28Þ

p tr nþ1 ¼ nþ1  Cn D nþ1   tr znþ1 ¼ znþ1 þ c znþ1 tr D pnþ1

(38) (39) (40) (41)

where we have defined

! @f @r

:¼ nþ1

@f @ r

 nþ1 ;znþ1 Þ ðr

124

R. Conti et al. / Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133

for simplicity of notation, and

B :¼ B

1 þ s=Dt Dcnþ1 s=Dt

Dc nþ1 :¼

for convenience. In Eq. (39), hnþ1 is the discrete form of the hardening function for the rate-dependent case, given by:

!

hnþ1

@f ¼ qc znþ1 tr @r

@f @r

D :¼ Dr

Bnþ1 ¼

)

A¼1

)

Dtr nþ1

3 X ðAÞ;tr ¼ ðDA Þtr nþ1 mnþ1

)

Tnþ1 ¼

3 X

Remark 5.2. As in [35], system (38)–(41), which defines the numerical algorithm for the plastic corrector problem, is wellconditioned for any value of the viscous parameter s 2 ½0; 1Þ. In particular, one can set s ¼ 0 to recover exactly the rate-independent limit. Remark 5.3. From Eqs. (34) and (38) it follows that, if A ¼ C1 n , r nþ1 is the projection of both the trial stress, rtrnþ1 , and the current stress, rnþ1 , onto the updated elastic domain, in the metric induced by C1 n . However, as already pointed out in Section 4.2, the proposed algorithm does not coincide with the viscoplastic formula nþ1 is not the solution of the ratetion proposed in [37], as r independent evolution problem. Remark 5.4. The Lagrange multiplier in Eqs. (38) and (39), which was formally introduced to define the projection of the current stress on the elastic domain, corresponds, in the limit as s ! 0, to the consistency parameter for the rate-independent case. It is worth noting that, due to the assumption of isotropic behaviour, the system of Eqs. (38)–(41) is redundant [39,38,5]. As a matter of fact, the number of independent stress components  nþ1 Þ can be reduced to six, thus simplifying considerably ðrnþ1 ; r the numerical solution of the plastic corrector problem. One possible strategy to perform such a reduction is considering as independent stress invariants the three principal stresses [5,39].

where BA ; DA and T A are the eigenvalues of B, D and T, respectively. From the isotropy assumption, we get

ðBA Þnþ1

! @f ¼ b1 tr @r

þ 2b2 nþ1

3 X

rA mðAÞ ; mðAÞ :¼ nðAÞ  nðAÞ

! nþ1

ð47Þ

Using the spectral representation, the plastic corrector problem may be formulated in principal stress space, and Eq. (38) transforms into 3 3 X X ðAÞ;tr  ðAÞ  A Þnþ1 m ðr ðrA Þtr nþ1 mnþ1 nþ1 ¼ A¼1

A¼1



3 Dc nþ1 X  ðAÞ ðBA Þnþ1 m nþ1 1 þ s=Dt A¼1

ð48Þ

From Eq. (48), it follows that: ðAÞ;tr  ðAÞ m nþ1 ¼ mnþ1

 A Þnþ1 ¼ ðrA Þtr ðr nþ1 

Dc nþ1 ðBA Þnþ1 1 þ s=Dt

ð49Þ

for A ¼ 1; 2; 3. Note that, since the trial stress is known, so are its principal directions. By applying the same arguments to Eq. (41), the only unknown quantities to be determined are the principal  A Þnþ1 , the internal variable znþ1 and the Lastresses ðrA Þnþ1 and ðr grange multiplier Dc nþ1 . Therefore, the plastic corrector problem in principal stress space can be recast as: (50) (51) (52)

ð43Þ

(53)

A¼1

where rA and nðAÞ are the principal stress and principal stress direction, respectively. If material isotropy is assumed, the yield function depends on r through its principal values, that is f ðr; zÞ ¼ f ðr1 ; r2 ; r3 ; zÞ, and its gradient can be computed as 3 X @f @f ¼ mðAÞ @ r A¼1 @ rA

ð44Þ

Let us introduce for convenience the following fourth-order isotropic tensors 1

B :¼ FA

1

D :¼ F

@f @ rA

tr tr ðDA Þtr nþ1 ¼ d1 trðrnþ1 Þ þ 2d2 ðrA Þnþ1  A Þnþ1  nþ1 Þ þ 2s2 ðr ðT A Þnþ1 ¼ s1 trðr

5.1.3. Plastic corrector in principal stress space We recall the spectral representation

r¼

ðAÞ  nþ1 ðT A Þnþ1 m

A¼1

ð42Þ

Eqs. (38)–(41) represent an uncoupled system of 8 þ 6 non-linear  nþ1 ; znþ1 ; Dc nþ1 and algebraic equations in the 8 þ 6 unknowns r rnþ1 , which can be solved iteratively by Newton’s method.

ð46Þ

A¼1

T :¼ Tr nþ1

3 X ðAÞ  nþ1 ðBA Þnþ1 m

where we have defined

r^ :¼ ðr1 ; r2 ; r3 Þ> r^ :¼ ðr 1 ; r 2 ; r 3 Þ> r^ tr :¼ ðrtr1 ; rtr2 ; rtr3 Þ> ^ :¼ ðB1 ; B2 ; B3 Þ> B ^ :¼ ðD1 ; D2 ; D3 Þ> D T^ :¼ ðT 1 ; T 2 ; T 3 Þ>

¼ b1 I  I þ 2b2 I

¼ d1 I  I þ 2d2 I

ð45Þ

T :¼ F1 Cn A ¼ s1 I  I þ 2s2 I where bi ; di and si can be easily computed from Cn and A, by exploiting some properties of fourth-order isotropic tensors (see Appendix A). By introducing the spectral decomposition of  ð@f =@ rÞnþ1 , rtr nþ1 and rnþ1 , we can define

5.1.4. Iterative solution of plastic corrector equations By exploiting the uncoupled nature of system (50)–(53), a twostep procedure can be used for the plastic corrector problem (see Box 2). Step 1. Solve iteratively the system of five non-linear algebraic equations (50)–(52), which can be rewritten in residual form as:

125

R. Conti et al. / Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133

8 9 8 r 9 Dc ^ nþ1 > ^ nþ1 þ r ^ tr > B  1þsnþ1 r r > > > > nþ1 = D t nþ1 < = < = Dc nþ1 tr rðxÞ :¼ r znþ1 :¼ znþ1 þ znþ1 þ 1þs=Dt hnþ1 > > > > > > : ; : ; fnþ1 fnþ1

ð54Þ

^  nþ1 ; znþ1 ; Dc nþ1 g. The convergence of the algorithm is where x :¼ fr detected at each iteration through the expressions:

8 r ;ðkÞ > kr k < TOLr :¼ TOL  krtr > nþ1 k < nþ1 z;ðkÞ tr krnþ1 k < TOLz :¼ TOL  kznþ1 k > > : ðkÞ 2 < TOLf :¼ TOL  kztr jfnþ1 j nþ1 k

ð55Þ Box: 2. Numerical algorithm for the plastic corrector scheme in slow dynamics.

where TOL is a prescribed tolerance value. Step 2. Compute the rate-dependent viscous solution by applying the explicit Eq. (53). By exploiting the property of Dc nþ1 (see Remark 5.4), we can use the sign of the Lagrange multiplier to determine if critical softening conditions are attained during the plastic process at time tnþ1 . Just as in the continuum problem (see Section 3), indeed, the sign of the plastic modulus, which determines the transition between slow and fast dynamics, is related to the sign of the discrete consistency parameter. Accordingly, at the end of each time step (or substep, if the substepping procedure is activated), the sign of the updated consistency parameter is checked. Depending on the sign of Dc nþ1 and on the state of the material at time t ¼ t n , three conditions may occur: (i) If Dc nþ1 > 0, the process is declared to occur under slow dynamics conditions (no critical softening). In this case the viscous solution represents the actual final state of the material. (ii) If Dc nþ1 < 0 and the state of the material at time t n is plastic, critical softening occurs. The viscous solution is then rejected and the equations of the fast dynamics are solved to find the þ new initial conditions ðrþ for the incremental n ; zn Þ ¼ ðr; zÞjt¼t þ n evolution problem at time tnþ1 (see Section 5.2). (iii) If Dc nþ1 < 0, and the state of the material at time t n is inside the elastic domain, the actual time step must be subdivided in order to determine the portion of the stress increment that lies within the yield surface. This additional procedure, which is necessary to avoid the fast dynamics being started from an initial stress state inside the yield locus, is performed with a simple bisection algorithm (See Box 3). Remark 5.5. As pointed out in [12], for every s > 0 the continuum viscoplastic evolution problem is well posed and a unique solution exists for system (19). As far as the numerical integration is concerned, this would imply the discrete Lagrange multiplier being always positive in the plastic corrector scheme, i.e. Dc nþ1 > 0 for every s > 0. Nonetheless, as will be shown in Section 6, for very small values of the viscosity parameter (s 1) the converged solution in the return mapping scheme may correspond to a negative value of Dc nþ1 , this implying a discontinuity in the viscoplastic solution and a transition between slow and fast dynamics. 5.2. Fast Dynamics Let ð1; þ1Þ be the time interval for the dilated time s. We as sume that the state of the material ðr n ; zn Þ is known at s ¼ 1, i.e. at time t ¼ t  when a discontinuity occurs in the solution of system n (19). The problem to be addressed is to compute the state variables þ þ ðrþ n ; zn Þ at s ¼ þ1, that is at the end of the jump (t ¼ t n ), through the integration of (20). As the internal variable is strictly decreasing, dz=ds is invertible. Therefore, using z as independent variable, from Eqs. (20) we get:

dr CAðpA ðrðzÞÞ  rðzÞÞ ¼ dz qc z tr½AðrðzÞ  pA ðrðzÞÞÞ

þ Note that we can use the condition f ðrþ n ; zn Þ ¼ 0 (see Section 4.2) to detect the end point of the fast dynamics along the integration of (56). To integrate Eq. (56) let us define a decreasing sequence zs > zsþ1 > . . ., where zsþ1 ¼ zs þ Dzsþ1 and Dzsþ1 < 0 is imposed, such that the first and the last term in the series correspond to þ z n (known) and zn (unknown) respectively. At each integration step, the problem to be addressed is then to update the stress state rsþ1 by integrating Eq. (56), with the initial condition rjz¼zs ¼ rs .

ð56Þ

(1) Initialize: ðkÞ ðkÞ ðkÞ ^ ^ tr  nþ1 ¼r znþ1 ¼ ztr Dcnþ1 ¼0 k ¼ 0; r nþ1 ; nþ1 ; (2) Check for iteration number: if k > ITMAX then Set CONV=‘NO’ and EXIT end if (3) Check for convergence: ðkÞ if rnþ1 < TOL then ðkÞ if Dcnþ1 P 0 then ðkÞ xnþ1 ¼ xnþ1 ^ nþ1 þ T ^ nþ1 r^ nþ1 ¼ s=DtD Set CONV=‘YES’ and EXIT ðkÞ else if Dcnþ1 < 0 and fn < TOLf then GO TO Box 3 else GO TO Box 4 end if else GO TO (4) end if (4) Evaluate increment at local iteration k "  #1 ðkÞ @r ðkÞ ðkÞ dxnþ1 ¼  rnþ1 @x nþ1 (5) Update state variables and Lagrange multiplier ðkþ1Þ

ðkÞ

ðkÞ

xnþ1 ¼ xnþ1 þ dxnþ1 (6) Set k

k þ 1 and GO TO (2)

Box: 3. Numerical algorithm for the bisection step. (1) Use bisection method to find a such that: a Dtanþ1 ¼ aDt nþ1 Dt 1 nþ1 ¼ ð1  aÞDt nþ1 1a a _ Deanþ1 ¼ e_ nþ1 Dtanþ1 De1 nþ1 ¼ enþ1 Dt nþ1 and ranþ1 ¼ rn þ Cn Deanþ1 zanþ1 ¼ zn jf ðranþ1 ; zanþ1 Þj 6 TOLf (2) Update variables a Set rn ¼ ranþ1 ; zn ¼ zanþ1 , Denþ1 ¼ De1 nþ1 ðiÞ 1a ðiÞ Set Dt ¼ Dtnþ1 ; t ¼ t n þ aDtnþ1 GO TO Box 1

By using an implicit backward Euler scheme and freezing the elastic parameters at the beginning of the fast dynamics, we get

rsþ1 ¼ rs 

Dzsþ1  sþ1 Þ Cn Aðrsþ1  r hsþ1

ð57Þ

126

R. Conti et al. / Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133

As for the integration of the plastic corrector problem in the slow dynamics regime, we adopt a return mapping scheme to com sþ1 from the initial stress, rs , as: pute r

q σs=σn-

σs+1

f(σn-,zn-)=0

r sþ1 ¼ rs  Dcsþ1 X

σs+1

X ¼

zs=zn

p

-

Δzs+1

ð61Þ

where X is a fourth-order tensor to be defined in order that Eq. (59) must hold. Using (57) and (61) we get

f(σs+1,zs+1)=0 zs+1

@f @ r ðr sþ1 ;zsþ1 Þ



Dzsþ1 1 FA ; hsþ1

F ¼

hsþ1 I þ Cn A Dzsþ1

ð62Þ

where, again, invertibility of F is guaranteed by the fact that hsþ1 =Dzsþ1 > 0. Substituting (62) into (61), and using (57), we obtain

(a)

(63) (64)

q

(65)

σs

σs+1

where, again, we have defined

f(σs,zs)=0

σs+1

! @f @r

f(σs+1,zs+1)=0

:¼ sþ1

@f @ r

 sþ1 ;zsþ1 Þ ðr

for simplicity of notation, and

zs+1 zs Δzs+1

p

(b)



Dc sþ1 :¼

1 þ hsþ1 =Dzsþ1 Dcsþ1 hsþ1 =Dzsþ1

for convenience. Using Eqs. (58) and (59), we can rewrite hsþ1 as:

!



hsþ1 ¼ Dc sþ1 qc zsþ1 tr

q

 Dzsþ1

ð66Þ

sþ1

It follows that Eqs. (63)–(65) represent an uncoupled system of 7 + 6 nonlinear algebraic equations in the 7 + 6 unknowns r sþ1 ; Dc sþ1 and rsþ1 . By exploiting the isotropy assumption for the material behaviour, we can reduce the number of independent stress component  sþ1 Þ to six, thus simplifying considerably the numerical ðrsþ1 ; r solution of the fast dynamics.

σs σn+=σs+1=σs+1 f(σn+,zn+)=0 zn+=zs+1 zs Δzs+1

p f(σs,zs)=0

(c)

Fig. 3. Schematic representation of the return mapping algorithm for fast dynamics: (a) at the inception, (b) during and (c) at the end of the jump. hsþ1

where may be viewed as the discrete counterpart of the hardening function for the fast dynamics, defined as

 sþ1 Þ hsþ1 :¼ qc zsþ1 tr½Aðrsþ1  r

ð58Þ

 sþ1 :¼ pA ðrsþ1 Þ is the minimal distance projection of rsþ1 onto and r  sþ1 we have (see Rethe updated elastic domain. By definition of r mark 5.4),

r sþ1

@f @r

@f ¼ rsþ1  Dcsþ1 A1 @r

ð59Þ

5.2.1. Fast dynamics in principal stress space Let us introduce for convenience the following tensor notation:

 ; zÞI  I þ 2/ 2 I F ¼ / 1 ðr and

 ; zÞI  I þ 2b 2 I B :¼ F A1 ¼ b 1 ðr  ; zÞI  I þ 2d 2 I D :¼ F 1 ¼ d 1 ðr

1

T :¼ F

 1 ð ; zÞI

Cn A ¼ s r

where Dcsþ1 P 0 is the Lagrange multiplier.

 I þ 2s

ð68Þ 2I

By introducing the spectral decomposition of ð@f =@ rÞsþ1 ; rs and r sþ1 , we can define

B :¼ B

@f @r

D :¼ D r

)

Bsþ1 ¼

3 X  ðAÞ ðB A Þsþ1 m sþ1 A¼1

)

Dsþ1

3 X ¼ ðD A Þsþ1 msðAÞ

 sþ1 ;zsþ1 Þ ðr

 sþ1 ; zsþ1 Þ ¼ 0 Dcsþ1 f ðr

ð67Þ

A¼1

ð60Þ

T :¼ T r

)

Tsþ1 ¼

3 X  ðAÞ ðT A Þsþ1 m sþ1 A¼1

ð69Þ

R. Conti et al. / Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133

where B A ; D A and T A are the eigenvalues of B ; D and T , respectively. From the isotropy assumption, we get

ðB A Þsþ1 ¼ b 1;sþ1 tr

! @f @r

þ 2b 2 sþ1

@f @ rA

 1;sþ1 trð sþ1 Þ

¼s

r

Box: 4. Numerical algorithm for the fast dynamics. sþ1

ð70Þ

 2 ð A Þsþ1

þ 2s

r

 sþ1 ; zsþ1 Þ, d 1;sþ1 :¼ d 1 ðr  sþ1 ; zsþ1 Þ where we have defined b 1;sþ1 :¼ b 1 ðr  sþ1 ; zsþ1 Þ. and s 1;sþ1 :¼ s 1 ðr Using this spectral representation and by applying the same arguments adopted in Section 5.1.3, Eqs. (63)–(65) can be rewritten in principal stress space as (71) (72) (73)

where the only unknown quantities to be determined are the prin A Þsþ1 , with A ¼ 1; 2; 3, and the Lagrange cipal stresses ðrA Þsþ1 and ðr multiplier Dc sþ1 . 5.2.2. Iterative solution of fast dynamics equations The following procedure can be used for the solution of the fast dynamics equations (see Box 4 and Fig. 3). Step 1. Update the internal variable zsþ1 ¼ zs þ Dzsþ1 , where we can define Dzsþ1 :¼ fz n . The integration parameter f defines the ’’rate’’ at which the yield surface shrinks during the fast dynamics. In this work we have used f ¼ ½104 ; 102 . Step 2. Update the viscous stresses by solving system (70)–(72) in a two-step procedure. First, solve iteratively Eqs. (70) and (71) , which can be rewritten in residual form as:

rðxÞ :¼

r  rsþ1 fsþ1

( :¼

initial conditions and/or model parameters critical softening instabilities may be detected during a loading path.

!

ðD A Þsþ1 ¼ d 1;sþ1 trðrs Þ þ 2d 2 ðrA Þs ðT A Þsþ1

127

Dc sþ1 sþ1 =Dzsþ1

^ sþ1 þ r ^ s  1þh r

^ sþ1 B

) ð74Þ

(1) Initialization:  ^ s ; zs Þ ¼ ðr ^ ðr n ; zn Þ (2) Update internal variable: Dzsþ1 ¼ fzn zsþ1 ¼ zs þ Dzsþ1 (3) Update stresses (Newton–Raphson algorithm): (i) Initialize: ðkÞ ^ ^  ðkÞ Dcsþ1 ¼0 k ¼ 0; r sþ1 ¼ rs ; (ii) Check for iteration number: if k > ITMAX then Set CONV=‘NO’ and GO TO (4) end if (iii) Check for convergence: ðkÞ if rsþ1 < TOLr then ðkÞ xsþ1 ¼ xsþ1 ^ sþ1 ^ sþ1 þ T ^ sþ1 ¼ h =Dzsþ1 D r sþ1

GO TO (4) else GO TO (iv) end if (iv) Evaluate increment at local iteration k ðkÞ ðkÞ ðkÞ dxsþ1 ¼ ðJsþ1 Þ1 rsþ1 (v) Update state variables and Lagrange multiplier ðkþ1Þ ðkÞ ðkÞ xsþ1 ¼ xsþ1 þ dxsþ1 (vi) Set k k þ 1 and GO TO (ii) (4) Check yield condition ^ sþ1 ; zsþ1 Þ fsþ1 :¼ f ðr if jfsþ1 j < TOLf then þ ^þ ^ ðr n ; zn Þ ¼ ðrsþ1 ; zsþ1 Þ and GO TO Box 1 else if fsþ1 P TOLf then28 Set ðÞs ¼ ðÞsþ1 and GO TO (2) else if fsþ1 < TOLf or CONV=‘NO’ then Set f ¼ rf (0 < r < 1) and GO TO (2) end if

fsþ1

^  sþ1 ; Dc sþ1 g. Then, compute the rate-depenwhere x :¼ fr dent viscous stress by applying the explicit Eq. (73). ^ sþ1 ; zsþ1 Þ ¼ 0, the Step 3. Check the consistency condition. If f ðr stress has returned to the yield surface: the viscous solution represents the state of the material at the end of the jump and we can proceed along the solution of the equations of the slow dynamics with the new ^ sþ1 ; zsþ1 Þ > 0, the viscous initial conditions. If f ðr stress is outside the yield surface and the viscous solution is still in fast dynamics regime. Condition ^ sþ1 ; zsþ1 Þ < 0, which implies that the Lagrange multif ðr plier Dc sþ1 is negative, can occur if the value adopted for Dzsþ1 is too large. In this case the viscous solution obtained is rejected, Dzsþ1 is reduced and the procedure resumes from Step 1. 6. Numerical examples In this section we assess the ability of the proposed algorithm to handle critical softening in Cam Clay plasticity. More specifically, we present the results of a number of single-element tests reporting the convergence properties of the algorithm. Moreover, we study the sensitivity of the numerical results to some integration and physical parameters. Finally, we evaluate under which set of

All the analyses were performed prescribing an error tolerance TOL ¼ 1:0  1010 and a maximum number of iterations of 10 (slow dynamics) and 20 (fast dynamics) in the Newton procedure.

6.1. Convergence and accuracy analyses To evaluate the convergence and accuracy properties of the proposed algorithm, we present in this section a number of isochoric axisymmetric compression tests for which a reference solution was obtained by Dal Maso and DeSimone [9] by integrating the evolution equations using the code Mathematica v.8 [44]. Dal Maso and DeSimone considered a simplified version of the MCC model, in which the (isotropic) elastic stiffness has been assumed constant and characterized by a zero Poisson ratio, so that C ¼ 2GI, and the hardening function is given by:

  @f _ c p tr z_ ¼ qc p trðe_ p Þ ¼ cq @r

ð75Þ

This last equation is slightly different from the hardening law (15) as z_ is assumed proportional to p instead of z. Note that, under this last assumption, equations (50)–(53) for the slow dynamics and (70)–(72) for the fast dynamics represent a coupled system of equations which have to be solved simultaneously.

128

R. Conti et al. / Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133

Table 1 Single-element testing program Run #

s

z0

f

Step #

1

0

1:0  10

6:54  13:48  22:50  33:71

10

2

0

1:0  103

6:54  13:48  22:50  33:71

20

3

0

1:0  103

6:54  13:48  22:50  33:71

100

4

0

1:0  103

6:54  13:48  22:50  33:71

10000

5

3

1:0  103

22:50

10000

3

22:50

10000

1:0  102

1:0  103

22:50

10000

1

1:0  103

22:50

10000

1:0  1010

6

1:0  10

7 8

1:0  10

4

1:0  10

9

0

1:0  104

22:50

10000

10

0

1:0  102

22:50

10000

11

0

5:0  102

22:50

10000

12

0

1:0  101

22:50

10000

The simulations have been performed adopting the following transformation of the stress and strain quantities:

r :¼

1 z r z :¼ pa pa

e :¼

2G e pa

e p :¼

2G p e pa

where pa is the atmospheric pressure. In terms of the above non– dimensional variables, the evolution equations for the state variables and the yield function reduce to:

r_ ¼ e_  e_ p z_ ¼ q c p trðe_ p Þ f ðr ; z Þ ¼ p ðp  2z Þ  ðq =MÞ2 ¼ 0

ð76Þ ð77Þ ð78Þ

where q ¼ qc pa =ð2GÞ. Following Dal Maso and DeSimone [9], the values M ¼ 1 and q s ¼ 1 have been adopted in the simulations. The complete testing programme for this first group of analyses is detailed in Table 1. The initial (isotropic) stress state is defined by r 11 ¼ r 22 ¼ r 33 ¼ 2=3, while a total axial strain increment e 11 ¼ 10 was applied with a rate of e_ 11 ¼ 0:1 in 10, 20, 100 and 10,000 steps. The initial values for the internal variables are z 0 ¼ 6:54; 13:48; 22:50; 33:71. Fig. 4 shows the results from Tests #1 to #4 in the q : p ; q : e q and p : e q planes. Due to the assumed initial conditions, the material response is always characterized by an initial elastic state followed by softening behaviour, until critical conditions are reached. However, while for z 0 ¼ 6:54  13:48 the whole plastic process evolves in the slow dynamics regime, for z 0 ¼ 22:50  33:71 fast dynamics is detected during the plastic loading, resulting in a jump c

Fig. 4. Isochoric axisymmetric compression tests at different values of the internal variable: comparison with literature results.[9]

Fig. 5. Isochoric axisymmetric compression test results: influence of the viscous parameter.

R. Conti et al. / Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133

129

Fig. 6. Isochoric axisymmetric compression test results: influence of the integration parameter f.

solution obtained for z 0 ¼ 22:50 does not show any instability in the plastic regime. Fig. 5 shows the results from Tests #5 to #8, together with Test #4, which is taken as reference. The solution obtained with s ¼ 1:0  1010 is practically identical with the inviscid one and fast dynamics regime is detected during plastic loading. On the other hand, for larger values of the viscous parameter, the regularization introduced inhibits the inception instabilities. Note that the quite different trend exhibited by Tests #6 in the q : p plane is only apparent and reflects the triggering of substepping procedure during the integration process. Fig. 6 shows the results from Tests #4 and #9 to #12, in terms of the stress–strain response of the model during the fast dynamics regime. Fig. 7 shows, for the same tests, the evolution of the La grange multiplier (Dcsþ1 ; Dc sþ1 ) and the yield function (fsþ1 ) along the jump, as a function of the normalized number of integration steps performed during the integration of the fast dynamics. At the end of the fast dynamics we have Dcsþ1 ¼ 0 and fsþ1 ¼ 0, which reflect the fact that the stress state has returned on the yield surface. It is apparent, from Fig. 7, that resolving stress jumps may require many local time steps. However, there is simply no alternative within the available integration algorithms, since they loose meaning in the presence of sub-critical softening. 6.2. Isochoric axisymmetric compression tests at different preconsolidation pressures

Fig. 7. Isochoric axisymmetric compression test results: evolution of the Lagrange multiplier and the yield condition in fast dynamics.

in the stress–strain response of the material. It is worth noting that, while the 100 and 10,000 step solutions are practically identical to each other and to the reference solution, the 10 step

In this example, the element was compressed in isochoric axisymmetric conditions, starting from an isotropic stress state at p0 ¼ 200 kPa, and applying a total axial strain increment e11 ¼ 0:3 with a strain rate e_ 11 ¼ 0:1s1 in 10, 20, 100 and 10 000 steps. The initial values of the preconsolidation pressure are pc0 ¼ 300; 600; 900; 1200 kPa. We assume as material parameters: ^ ¼ 0:013, ^ M ¼ 1:3; G ¼ 1000 kPa; j k ¼ 0:032; pr ¼ 10 kPa. Moreover, we assume s ¼ 0; f ¼ 1:0  103 and A ¼ I. Fig. 8 shows, for the initial conditions considered in this example, the yield surface and the locus at which critical softening is attained in the q : p plane. For pc0 ¼ 300 kPa, corresponding to which the condition K p ðq; p; zÞ ¼ 0 is always attained for stress states inside the yield surface, no instabilities can be detected for any loading condition applied. For pc0 ¼ 600 kPa, critical softening is

130

R. Conti et al. / Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133

Fig. 8. Yield surface and limit locus for critical softening for different values of the initial reconsolidation pressure.

Fig. 9. Isochoric axisymmetric compression test results.

expected along the isochoric compression, as soon as the stress state reach the yield surface. Finally, for pc0 ¼ 900 kPa and 1200 kPa a jump in the stress–strain response of the model can occur only after the plastic loading has been evolved initially in slow dynamics regime. Fig. 9 shows the numerical results in in the q : p; q : eq and p : eq planes. Due to the assumed initial conditions, the material response is characterized by an initial elastic behaviour followed by hardening for pc0 ¼ 300 kPa and softening in the other tests. Critical conditions are always attained and the reduction of the deviatoric stress along the jump is more pronounced as the ratio pc0 =p0 increases. 6.3. Isochoric axisymmetric compression tests at different initial mean pressures In this example, the element was compressed in isochoric axisymmetric conditions, starting from an isotropic stress state at three different values of the mean pressure (p0 ¼ 20; 60; 100 kPa), and applying a total axial strain increment e11 ¼ 0:3 with a

strain rate e_ 11 ¼ 0:1s1 in 10 000 steps. The initial values of the preconsolidation pressure is pc0 ¼ 1200 kPa. We assume as ^ ¼ 0:045; ^ material parameters: M ¼ 1:3, G ¼ 1000 kPa; j k¼ 0:083, pr ¼ 10 kPa. Moreover, we assume f ¼ 5:0  104 and s ¼ 0; 1:0 104 ; 1:0  103 and we use the formulation proposed by both Duvaut–Lions (A ¼ I) and Simo (A ¼ C1 ) for the flow rule in the regularized model. Fig. 10 shows the results in the q : p and q : eq planes, obtained with the Duvaut–Lions viscoplastic flow rule. The model exhibits softening after the yield surface has been reached; moreover, critical conditions are attained for tests starting from p0 ¼ 60 kPa and p0 ¼ 100 kPa. It is worth noting that, for the adopted values of the viscous parameter, no practical differences can be observed between the inviscid solutions and the rate-dependent ones. A completely different behaviour is observed for A ¼ C1 (see Fig. 11), where the regularization parameter affects significantly the plastic regime, by inhibiting the inception of instabilities for the two cases with p0 ¼ 60 kPa and p0 ¼ 100 kPa. On the contrary, when the viscous parameter is set to zero the choice of A does not affect the

R. Conti et al. / Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133

131

Fig. 10. Isochoric axisymmetric compression test results for A ¼ I.

Fig. 11. Isochoric axisymmetric compression test results for A ¼ C1 .

model prediction along the fast dynamics, as can be observed by comparing Figs. 10 and 11, where the rate-independent solutions are always equal. 7. Conclusions Following the theoretical works by Dal Maso and coworkers [9,12,11], in which a viscous regularization technique was proposed to handle critical and subcritical softening in Cam-Clay plasticity, a numerical algorithm for the integration of the regularized viscoplastic equations has been developed for dealing with the occurrence of time discontinuities in computational plasticity. As for the continuum case, the proposed strategy guarantees the existence and uniqueness of the solution for the discrete evolution problem in the whole space of the admissible stress states. In the continuum problem, the mechanical response of the material during a strain driven process can evolve either smoothly (slow dynamics) or by exhibiting time discontinuities (fast dynamics) if critical or subcritical softening is attained. The transition between the two regimes occurs whenever the plastic modulus becomes null or negative, which is not uncommon under low mean

pressure and in strongly dilating materials. Moreover, strong softening leading eventually to critical conditions may occur in bonded granular materials such as cemented soils or weak rocks, see e.g. [25] and references therein. We have shown that in the discrete evolution problem, the inception of such discontinuities during a generic time integration step can be related to the sign of a discrete Lagrange multiplier, which correspond to the consistency parameter in the limit as s ! 0. Whenever the Lagrange multiplier becomes negative, the equations of the fast dynamics are integrated to find the updated state of the material at the end of the jump. This work has focused solely on the constitutive response of the material at the element level, without taking into account the solution of boundary value problems. The MCC model has been selected as the prototype of a family of isotropic hardening plasticity models exhibiting pressure dependence and dilationrelated hardening/softening. However, the numerical strategies developed and the results obtained are by no means limited to MCC alone, and could be extended to more recent and versatile isotropic plasticity models for granular materials, see e.g. [25,14,27].

132

R. Conti et al. / Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133

As viscous regularization techniques have been applied successfully to the analysis of strain localization in solids [28], a future extension of this work will consider the application of the proposed numerical techniques to the solution of boundary value problems. This will enable us to investigate the interaction between time and spatial discontinuities in determining the structural response of the material body during arbitrary loading processes. Acknowledgements The research leading to these results has received specific funding under the ’’Young SISSA Scientists’ Research Projects’’ scheme 2011–2012, promoted by the International School for Advanced Studies (SISSA), Trieste, Italy. Appendix A Let A :¼ a1 I  I þ 2a2 I and B :¼ b1 I  I þ 2b2 I be two fourth-order isotropic tensors. Then it can be easily proved that:

C :¼ A þ B ¼ c1 I  I þ 2c2 I;

c1 ¼ a1 þ b1 ; c2 ¼ a2 þ b2

ðA:1Þ

and

C :¼ AB ¼ c1 I  I þ 2c2 I;

c1 ¼ 3a1 b1 þ 2ða1 b2 þ a2 b1 Þ; c2 ¼ 2a2 b2

ðA:2Þ

and

C :¼ A1 ¼ c1 I  I þ 2c2 I;

a1

c1 ¼ 

2a2 ð2a2 þ 3a1 Þ 1 c2 ¼ 4a2

;

ðA:3Þ

Let E be a second-order tensor. Its spectral representation is:

E¼

3 X EA mðAÞ ;

mðAÞ :¼ nðAÞ  nðAÞ

ðA:4Þ

A¼1

with EA and nðAÞ being the eigenvalues and eigenvectors of E, respectively. Then

S :¼ AE ¼

3 X SA mðAÞ

ðA:5Þ

A¼1

with

SA ¼ a1 trðEÞ þ 2a2 ðEA Þ:

ðA:6Þ

References [1] H. Alawaji, K. Runesson, S. Sture, K. Axelsson, Implicit integration in soil plasticity under mixed control for drained and undrained response, Int. J. Numer. Anal. Methods Geomech. 13 (1992) 737–756. [2] F. Armero, A. Pérez-Foguet, On the formulation of closest-point projection algorithms in elastoplasticity – Part I: The variational structure, Int. J. Numer. Methods Engrg. 53 (2002) 297–329. [3] D. Bigoni, T. Hueckel, Uniqueness and localization-i. associative and nonassociative elasto–plasticity, Int. J. Solids Struct. 28 (1991) 197–213. [4] R. Borja, S. Lee, Cam-Clay plasticity, Part I: Implicit integration of elasto–plastic constitutive relations, Comput. Methods Appl. Mech. Engrg 78 (1990) 49–72. [5] R. Borja, K. Sama, P. Sanz, On the numerical integration of three-invariant elastoplastic constitutive models, Comput. Methods Appl. Mech. Engrg 192 (2003) 1227–1258. [6] R. Borja, C. Tamagnini, Cam-Clay plasticity, Part III: Extension of the infinitesimal model to include finite strains, Comput. Methods Appl. Mech. Engrg 155 (1998) 73–95.

[7] G. Buscarnera, G. Dattola, C. di Prisco, Controllability, uniqueness and existence of the incremental response: A mathematical criterion for elastoplastic constitutive laws, Int. J. Solids Struct. 48 (2011) 1867–1878. [8] R. Butterfield, A natural compression law for soils, Géotechnique 29 (1979) 469–480. [9] G. Dal Maso, A. DeSimone, Quasistatic evolution for cam-clay plasticity: Examples of spatially homogeneous solutions, Math. Models Methods Appl. Sci. 19 (2009) 1643–1711. [10] G. Dal Maso, A. DeSimone, M. Mora, M. Morini, A vanishing viscosity approach to quasistatic evolution in plasticity with softening, Arch. Ration. Mech. Anal. 189 (2008) 469–544. [11] G. Dal Maso, A. DeSimone, F. Solombrino, Quasistatic evolution for cam-clay plasticity: A weak formulation via viscoplastic regularization and time rescaling, Calc. Var. 40 (2011) 125–181. [12] G. Dal Maso, F. Solombrino, Quasistatic evolution for Cam-Clay plasticity: The spatially homogeneous case, Netw. Heter. Media 5 (2010) 97–132. [13] E. de Souza Neto, D. Peric´ , D. Owen, Computational Methods for Plasticity: Theory and Applications, Wiley & Sons, New-York, 2008. [14] A. DeSimone, C. Tamagnini, Stress-dilatancy based modelling of granular materials and extensions to soils with crushable grains, Int. J. Numer. Anal. Methods Geomech. 29 (2005) 73–101. [15] G. Duvaut, J. Lions, Les Inéquations en Méchanique et En Physique, Dunod, Paris, 1972. [16] H. van Eekelen, Isotropic yield surfaces in three dimensions for use in soil mechanics, Int. J. Numer. Anal. Meth. Geomech. 4 (1980) 89–101. [17] M. Hardy, J. Hudson, C. Fairhurst, The failure of rock beams, Part II – experimental studies, Int. J. Rock Mech. Min. Sci. 10 (1973) 69–82. [18] K. Hashiguchi, M. Ueno, Elasto–plastic constitutive laws of granular materials, in: Booktitle Constitutive Equations in Soils, 9th ICSMFE, Spec. Sess. 9, Tokyo, pp. 73–82. [19] S. Imposimato, R. Nova, An investigation on the uniqueness of the incremental response of elastoplastic models for virgin sand, Mech. Cohes. Frict. Mater. 3 (1998) 65–87. [20] R. Lagioia, A. Puzrin, D. Potts, A new versatile expression for yield and plastic potential surfaces, Comput. Geotech. 19 (1996) 171–191. [21] B. Loret, H. Prevost, Dynamic strain localization in elasto–(visco–)plastic solids, Part 1. General formulation and one-dimensional examples, Comput. Methods Appl. Mech. Engrg 83 (1990) 247–273. [22] D. Luenberger, Linear and Nonlinear Programming, Addison-Wesley, Reading, MA, 1984. [23] G. Maier, T. Hueckel, Non-associated and coupled flow-rules of elastoplasticity for rock-like materials, Int. J. Rock Mech. Min. Sci. 16 (1979) 77–92. [24] A. Needleman, Material rate dependence and mesh sensitivity in localization problems, Comput. Methods Appl. Mech. Engrg 67 (1988) 69–85. [25] R. Nova, C. Tamagnini, R. Castellanza, A constitutive model for bonded geomaterials subjected to mechanical and/or chemical degradation, Int. J. Num. Anal. Meth. Geomech. 27 (2003) 705–732. [26] A. Pérez-Foguet, A. Rodriguez-Ferran, A. Huerta, Consistent tangent matrices for substepping schemes, Comput. Methods Appl. Mech. Engrg 190 (2001) 4627–4647. [27] A. Piccolroaz, D. Bigoni, Yield criteria for quasibrittle and frictionals materials: A generalization to surface with corners, Int. J. Solids Struct. 46 (2009) 3587– 3596. [28] H. Prevost, B. Loret, Dynamic strain localization in elasto–(visco–)plastic solids, Part 2. Plane strain examples, Comput. Methods Appl. Mech. Engrg 83 (1990) 275–294. [29] R. Regueiro, C. Foster, Bifurcation analysis for a rate-sensitive, non-associative, three-invariant, isotropic/kinematic hardening cap plasticity model for geomaterials: Part I. Small strain, Int. J. Num. Anal. Meth. Geomech. 35 (2011) 201–225. [30] J. Rice, The localization of plastic deformation, in: Koiter (Ed.), Theoretical and Applied Mechanics, North-Holland Publishing Co., Amsterdam, p. 207–220. [31] K. Roscoe, J. Burland, On the generalized behaviour of wet clay, Engrg. Plast. 48 (1968) 535–609. [32] K. Roscoe, A. Schofield, C. Wroth, On the yielding of soils, Géotechnique 8 (1958) 22–53. [33] K. Runesson, M. Ristinmaa, L. Mähler, A comparison of viscoplasticity formats and algorithms, Mech. Cohes. Frict. Mater. 4 (1999) 75–98. [34] A. Schofield, C. Wroth, Critical State Soil Mechanics, McGraw-Hill, London, 1968. [35] J. Simo, Numerical analysis and simulation of plasticity, in: P. Ciarlet, J. Lions (Eds.), Handbook of Numerical Analysis,Elsevier Science, North Holland, The Netherlands, pp. 183–499. [36] J. Simo, T. Hughes, Computational Inelasticity, Springer, New-York, 1997. [37] J. Simo, J. Kennedy, S. Govindjee, Non-smooth multisurface plasticity and viscoplasticity. Loading/unloading conditions and numerical algorithms, Int. J. Numer. Methods Engrg. 26 (1988) 2161–2185. [38] J. Simo, G. Meschke, A new class of algorithms for classical plasticity extended to finite strains. Application to geomaterials, Comput. Mech. 11 (1993) 253– 278. [39] C. Tamagnini, R. Castellanza, R. Nova, A generalized backward euler algorithm for the numerical integration of an isotropic hardening elastoplastic model for mechanical and chemical degradation of bonded geomaterials, Int. J. Numer. Anal. Meth. Geomech. 26 (2002) 963–1004. [40] I. Vardoulakis, J. Sulem, Bifurcation Analysis in Geomaterials, Blackie Academic & Professional, Glasgow, 1995.

R. Conti et al. / Comput. Methods Appl. Mech. Engrg. 258 (2013) 118–133 [41] J. Vaunat, J. Cante, A. Ledesma, A. Gens, A stress point algorithm for an elastoplastic model in unsaturated soils, Int. J. Plast. 16 (2000) 121–141. [42] W. Wang, L. Sluys, R. de Borst, Viscoplasticity for instabilities due to strain softening and strain-rate softening, Int. J. Numer. Anal. Meth. Geomech. 40 (1997) 3839–3864.

133

[43] W. Wawersik, C. Fairhurst, A study of brittle rock fracture in laboratory compression experiments, Int. J. Rock Mech. Min. Sci. 7 (1970) 561–575. [44] S. Wolfram, The Mathematica Book, Cambridge University Press, New York, 1999.