Describing failure in geomaterials using second-order work approach

* Corresponding author.

?

Describing failure in geomaterials using second-order work approach

Fran?ois Nicota,*,Felix Darveba

aInstitut de Recherche en Sciences et Technologies de l'Environnement et de l'Agriculture (IRSTEA),Erosion Torrentielle Neige et Avalanches (ETNA)-Geomechanics Group,Grenoble 38402,France

bUniversite Joseph Fourier (UJF)-Institut National Polytechnique de Grenoble (INPG)-Centre National de la Recherche Scientifique (CNRS),Laboratoire Sols,Solides,Structures et Risques,Grenoble 38041,France

Received 13 March 2015; accepted 30 March 2015 Available online 12 May 2015

* Corresponding author.

Abstract

Geomaterials are known to be non-associated materials.Granular soils therefore exhibit a variety of failure modes,with diffuse or localized kinematical patterns.In fact,the notion of failure itself can be confusing with regard to granular soils,because it is not associated with an obvious phenomenology.In this study,we built a proper framework,using the second-order work theory,to describe some failure modes in geomaterials based on energy conservation.The occurrence of failure is defined by an abrupt increase in kinetic energy.The increase in kinetic energy from an equilibrium state,under incremental loading,is shown to be equal to the difference between the external second-order work,involving the external loading parameters,and the internal second-order work,involving the constitutive properties of the material.When a stress limit state is reached,a certain stress component passes through a maximum value and then may decrease.Under such a condition,if a certain additional external loading is applied,the system fails,sharply increasing the strain rate.The internal stress is no longer able to balance the external stress,leading to a dynamic response of the specimen.As an illustration,the theoretical framework was applied to the well-known undrained triaxial test for loose soils.The influence of the loading control mode was clearly highlighted.It is shown that the plastic limit theory appears to be a particular case of this more general second-order work theory.When the plastic limit condition is met,the internal second-order work is nil.A class of incremental external loadings causes the kinetic energy to increase dramatically,leading to the sudden collapse of the specimen,as observed in laboratory.?2015 Hohai University.Production and hosting by Elsevier B.V.This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Keywords:Failure in geomaterials; Undrained triaxial loading path; Second-order work; Kinetic energy; Plastic limit condition; Control parameter

E-mail address: francois.nicot@irstea.fr (Fran?ois Nicot).

Peer review under responsibility of Hohai University.

http://dx.doi.org/10.1016/j.wse.2015.05.001

1674-2370/?2015 Hohai University.Production and hosting by Elsevier B.V.This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1.Introduction

The notion of failure can be encountered in many fields,irrespective of the scale considered.This notion is essential in material sciences where the failure can be investigated on the specimen scale (the material point).It is also meaningful in civil engineering with regard to preventing or to predicting the occurrence of failure on a large scale.

Ifthedefinitionoffailureseemsmeaningfulinsomecases,at least from a phenomenological point of view,this is not always true,particularly,when considering heterogeneous materials.With regard to a granular assembly on a microscopic scale,failure might be related to the contact opening between initially contacting grains.However,the kinematic investigation of granular materials,along any given loading path,reveals that a large fraction of the contacts open without any visible failure patternobservedonthemacroscopicscale.Thus,theusualview of failure as the breakage of a given material body into two piecescannot beappliedtocomplex,dividedmaterialsmadeup ofanassemblyofelementaryparticlesorsub-systemsdescribed as approximately non-breakable.

For this reason,a mathematical definition of failure has emerged for solid materials.This definition was progressively built upon the plasticity theory,and developed early in the20th century in the field of metallic materials.Failure means that the external stress applied cannot be increased,and that finite strains may develop and progress in a constant stress state.

For rate-independent materials,before failure occurrence along a given loading path,a unique strain state exists under a given external stress applied to the materials.This mathematical definition,applied to the case of a material point,leads to the following equations:

where s is the Cauchy stress tensor,D is the strain rate tensor,and N0is the tangent stiffness operator.The constitutive behavior for this material point is defined by the incremental constitutive relationbetween the strain rate and a suitable objective time derivative of the Cauchy stress expressed with two six-component vectors D and s (Darve,1990; Wan et al.,2011),which requires that

Eq.(2)is the failure criterion,and the equationcorresponds to the associated failure rule.The failure rule gives the strain rate in different directions once failure has occurred.It should be emphasized that the magnitude of the strain rate remains unknown,and only the direction is determined,in accordance with the kernel of the tangent stiffness operator N′.

Basically,the plastic surface splits the stress space into two parts: the inner part and outer part.The inner part (plastic domain)includes the stress states that can be reached by the materials.The stress states located on the plastic surface are therefore referred to as the stress limit states.The classical theory states that failure occurs in a given stress state which is located on the plastic surface.According to this theory,no failure mode is expected to occur in a mechanical state inside the plastic surface.

One famous counter-example is the liquefaction of loose sands along axisymmetric isochoric triaxial paths (Castro,1969; Lade,1992; Darve,1990).During this test,the volume of the specimen is kept constant (isochoric conditions),and a constant axial displacement rate is imposed.This experiment shows that the curve giving the changes in the deviatoric stress q,defined as the difference between the axial stress s1and the lateral stress s3,against the mean effective pressure p0passes through a maximum value,as shown in Fig.1.If the test is strain-controlled by an imposed constant axial strain rate,the test can be pursued beyond the deviatoric stress peak until the collapse of the specimen: both the deviatoric stress q and the mean effective pressure p0decrease and tend toward zero.This is the well-known liquefaction phenomenon.Otherwise,when the deviatoric stress peak is reached,if an infinitesimal axial load is added,i.e.,the strain control is replaced with a stress control,then a sudden failure occurs.Clearly,such experimental evidence evokes the notion of failure.According to the classical theory,the fact that the deviatoric stress peak (point P in Fig.1)remains strictly inside the plastic domain means that no failure is supposed to occur at this point.However,the experiment shows that a proper loading control (namely,a stress control)can cause a sudden failure,leading to the collapse of the specimen.

Fig.1.Typical undrained triaxial behavior of loose sands.

As a consequence,the classic theory is not general enough,and a variety of failure modes may finally occur strictly within the plastic surface.Plastic failure,detected by Eq.(2),is a particular failure mode.Other failure modes may be encountered as well before the plastic limit is reached (such as the failure mode occurring at point P in Fig.1).The detection of these failure modes requires a novel framework,which should lead to a different criterion from that given in Eq.(2).

Before proceeding with a description of the novel framework's construction,we have to propose a clear definition of failure based on phenomenological arguments.The weakness of the classical theory lies in the fact that it is not based on a physical definition of failure,but rather on a mathematical concept of the stress limit state.The approach presented in this paper allows the recovery of the notion of a limit state in a broader way,as a consequence of the theory,but not as a basic definition.

As far as non-viscous materials are concerned,we conceive that failure is related to a transition from a quasistatic regime toward a dynamic regime,giving rise to a sudden acceleration of the material points: the kinetic energy of the system evolves from a nil value to a strictly positive one.Thus,the failure is not a state,but a transition (bifurcation)from a quasistatic regime with a nil value of kinetic energy toward a dynamic regime with a non-zero value of kinetic energy.The failure occurs in a given mechanical state,which will be described as being potentially unstable: an increase in kinetic energy may take place under the loading conditions (Nicot et al.,2009,2012).

In conclusion,the following definition around the notion of failure can be proposed: a material point in a given mechanical (stressestrain)state after a given loading history is described as being mechanically unstable as loading conditions lead to a bifurcation from a quasistatic regime toward a dynamic regime.This transition corresponds to a failure mode of the material.

2.Second-order work theory

2.1.Formulation on a large scale

An attempt at definition of failure was made in the previous section,asitrelatedtoatransition(bifurcation)fromaquasistatic regimetowardadynamicregime.Inwhatfollows,weinvestigate the conditions in which the kinetic energy of a non-viscous material system,in equilibrium at a given time,may increase.

For this purpose,a system made up of a given material,with a volume V0and a surface boundary S0,initially in a configuration C0(reference configuration),is considered.With a loading history,the system is in a strained configuration C,with a volume V and a surface boundary S,in equilibrium under a prescribed external load.Each material point in the volume V0is transformed into a material point in the volume V (Fig.2).All the material points in the volume V0are displaced,along with the deformation of their geometrical properties,including the surface vector,area,and volume.During this transformation,the material is likely to undergo rigid body motion,along with pure strain induced by stretching and spinning deformations.If large amounts of strain take place,both the initial configuration C0and current configuration C cannot be confused.

We introduce the transformation c relating each material point x of the current configuration C to the corresponding material point X of the initial configurationThe continuity of the material ensures that c is bijective.By means of the transformation,any field ~geXT of the initial position X can be transformed into the field g(x)of the current position,with

As c is bijective,the Jacobian determinant J of the tangent lineartransformation~F,withisstrictlypositive.~F isafunctionofthepositionX.Thedisplacementfields~ueXTatthe materialpoint in the initialconfiguration and u(x)at the material point in the current configuration are defined by the relationThus,is the Kronecher symbol.

The kinetic energy of the system of configuration C,in equilibrium at time t,is given by the energy conservation equation expressed in the rate form: where _Ecrepresents the kinetic energy rate of the system.It is convenient to express the integrals in Eq.(3)with respect to the initial configuration by using the transformation c.Recalling thatand using Nanson's formula,which relates the current surface vector ndS to the corresponding surface vector n0dS0in the initial configuration,Eq.(3)is written as follows: which yields,after some transformations (Nicot and Darve,2007),the following:

S0V0 where P denotes the Piola-Kirchoff stress tensor of the first kind,and

The advantage of the formulation given in Eq.(5)is that all integrals are written with respect to a fixed domain.Thus,differentiation of Eq.(5)gives,after simplification (Nicot and Darve,2007; Nicot et al.,2007),the following: where s?P,n0denotes the stress distribution applied to the initial (reference)configuration.Furthermore,for any time increment Dt,the second-order Taylor expansion of the kinetic energy reads

Combining Eqs.(6)and (8),and ignoring the o(Dt)term,finally yieldEq.(9)introduces explicitly the so-called internal secondorder work that is expressed through a semi-Lagrangian formalism as follows (Hill,1958):

Fig.2.Transformation of a material system and surface vectors of initial and current configurations.

The terminology of the internal second-order work is justified by the fact that Eq.(10)introduces the internal stress and strain variables P and~F.The first termon the right-hand side of Eq.(9)involves both stresses and displacements acting on the boundary of volume V0.This secondorder boundary term is therefore an external second-order work,and will be denoted hereafter as

Thus,for the system in an equilibrium configuration at time t,Eq.(9)indicates that the increase in the kinetic energy of the system,over a small time range from t to ttDt,can be obtained by

In particular,when the internal second-order work is nil,any loading condition,such thatensures thatan outburst in kinetic energy is expected.As an illustration,examples will be given in sections 3 and 4.2.2.Formulation on a small scale

The particularization of this theoretical framework to the case of homogeneous material specimens is worth mentioning,because it corresponds to the laboratory specimen scale.Parallelepiped-like specimens subjected to a prescribed force or displacement on each side surface,directing both stress and strain fields,are examined.Investigating this elementary scale can be useful in the interpretation of the derived experimental results.Material specimens are considered homogeneous when both strain and stress fields are macro-homogeneous,in the sense given by Hill (1967): the external forces applied to the boundary of the specimen are derived from the average stress tensor,and the displacements of each point on the boundary are derived from the average strain tensor.When such conditions are met,the homogenous specimen is also referred to as a representative element volume (REV).

Strictly speaking,a homogeneous specimen is not a material point,but actually a structural system with boundary conditions.As an illustration regarding granular materials,we can mention the fluctuating motion of grains that may exist even though the boundary conditions are kept constant.On average,over the whole specimen,the mechanical imbalance of grains is nil,but locally it is not nil.However,thanks to the macro-homogeneity,both strain and stress states are fully characterized by forces and displacements measured on the boundary.

Let a parallelepiped specimen be considered.Each side face i (i?1,2,3)admits a normal vector Nithat coincides with the direction of velocity viof a fixed reference frame.The initial area of face i is denoted by Ai,and the initial length of the corresponding side is denoted by Li,as shown in Fig.3.The subscript 1 refers to the axial direction,whereas the subscripts 2 and 3 refer to the two lateral directions perpendicular to the axial direction (Fig.3).When a static condition is assigned to face i,it is convenient to introduce a resultant external force fiacting on this face,which is set to be normal to the face considered.The uniform external Lagrangian stress vector distributionacting on face i related to fiis also introducedThe displacement of each side face i,along the direction xi,is denoted by Ui.No tangential displacement is assumed to take place.When a kinematic condition is assigned to face i,the resultant external force fi(orthe stress distribution s0i)acting on this side corresponds to the external loading that must be applied to produce the prescribed displacement Ui.

In these conditions,the displacement ~u of any point (x1,x2,x3)is

Fig.3.Parallelepiped specimen and definition of axes.

Then,we have,under homogeneous conditions,Finally,combining Eqs.(11),(14),and (15)leads to

Eq.(16)indicates that the increase in kinetic energy,from an equilibrium state,is basically related to the imbalance between the external force rate and the internal stress rate.The side face i of the specimen is subjected to the force rateon the external side,and to the force ratedirected by the internal stress,on its internal side.When a mechanical state is reached,in the continuum mechanics sense,imbalances may exist locally,but on average,the equilibrium state is reached; the external force rate equals the internal stress rate,i.e.,on face i.In this case,we verify,according to Eq.(16),that the kinetic energy of the system does not evolve.

It is essential to distinguish displacements and forces acting on the boundary of the specimen,with strain and stress acting within the specimen.During a quasistatic evolution of the specimen,along successive equilibrium states,the internal stress tensors within the specimen derived from internal forces applied to the sides are balanced with the external forces.This is sound until the specimen fails: if inertial effects take place,and the external stress is not balanced by the internal stress,a heterogeneous strain field may develop within the specimen.

3.Interpretation of diffuse failure along undrained triaxial loading paths

A homogeneous,parallelepiped specimen is considered throughout this section (Fig.3).The loading applied to the side face i of this specimen is supposed to occur in the normal direction.Each face i is subjected to a displacement Uialong the normal vector Ni,and the resultant force fiis also oriented along Ni.Neither tangential force nor shear strain is directed on the face.The Piola-Kirchoff stress tensor P is therefore diagonal.The same holds for the tensor_~F,which contains only diagonal terms,i.e.,Thus,it is convenient to replace the componentsandof diagonal tensorsandwith components,respectively,such that?,and

In the axisymmetric undrained triaxial test,both the axial displacement rate (_U1>0)and the volume of the specimen are keptconstant.Thecurrentvolumerate_V isgivenbytheintegralUsing the relationtogether with the second Green formula,finally yields the following:

where S0is the surface boundary of the parallelepiped specimen.Axisymmetric conditions with respect to axial direction mean that U2/L2?U3/L3and P2?P3(or f2/A2?f3/A3).Noting that A2?L1L3and A3?L1L2,the relation U2/L2? U3/L3is equivalent to U2A2?U3A3,so that

Thus,the isochoric condition reads as follows:

The experimental curve from the undrained triaxial test,as shown in Fig.4,giving the evolution of the internal deviatoric stress (qP?P1P3)versus the axial displacement,passes through a peak for sufficiently loose specimens,where point P shows the peak qP,and point M is in the softening regime.This is reported in abundant literature (Castro,1969; Lade and Pradel,1990; Lade,1992; Biarez and Hicher,1994; Chu et al.,2003).With a strain loading control (_U1>0 and _V?0),the response of the specimen follows a quasistatic evolution,passing through a succession of equilibrium states: Pi?fi/Ai.The curve in Fig.4 can be given in terms of the external deviatoric force qf?f1/A1f3/A3.

Let us focus on the deviatoric stress peak (point P).At this equilibrium point,the internal stress leads to axial stressbP1and lateral stressbP3on the boundary of the specimen,where the superscript ^ represents the peak value.Now,let us imagine that an additional small amount of deviatoric loading Dqfis applied to the specimen at time t,over a time range Dt: _qf?Dqf=Dt.This loading causes the system to evolve in such a way that _U1>0.Under such a condition,the internal deviatoric stress qPcannot exceed the peak valuePrecisely at the peak,

As _V?0,then 2A3_U3?A1_U1.Eq.(16)becomes

which yields,as _qP?0,the following:

Fig.4.Internal deviatoric stress versus axial displacement obtained from undrained triaxial test.

At time t,the system is at rest,with Ec(t)?0.Thus,Ec(ttDt)is strictly positive.At the peak of qP,under the effect of an additional deviatoric loading Dqf,the specimen's kinetic energy increases from zero to a strictly positive value DqfA1_U1Dt=2 over the time interval [t,ttDt].The failure mechanism of the specimen is initiated at the same time that the internal second-order work vanishes.In fact,the internal deviatoric stress qPwithin the specimen no longer balances the external deviatoric force qf: qfincreases fromto t Dqf,whereas qPfollows a constitutive path and decreases from the peak valuealong the descending branch.The unbalanced stress causes the dynamic response of the specimen,characterized by a strictly positive value of Ec(ttDt)in Eq.(21).This is exactly what is observed experimentally (Castro,1969; Lade and Pradel,1990; Lade,1992; Chu et al.,2003; Darve et al.,2004)and numerically when using a discrete element method (Sibille et al.,2008; Darve et al.,2007).

4.Plastic limit state and second-order work theory

The plastic surface corresponds to a set of limit states within the stress space.Whatever loading path is considered,the stress state cannot overpass the plastic surface.Let us consider a drained triaxial loading path.A constant axial displacement rate _U1is prescribed,along with a constant lateral force (f2?f3).For dense specimens,the curve from drained triaxial test,as shown in Fig.5,giving the evolution of f1(or P1?f1/A1)versus U1,passes through a peak,and then tends toward a plateau.Using this loading mode,the response of the specimen follows a quasistatic evolution.Thus,the external forces applied to each face are balanced by the internal stresses: fi?PiAi.Using this loading mode,the axial response of the specimen can be analyzed either in terms of the external force f1or in terms of the internal stress P1.Before the peak,along the ascending branch,the internal second-order work can be expressed as and is therefore strictly positive.

Fig.5.Axial stress versus axial displacement obtained from drained triaxial test.

Let us focus on the axial stress peak (point P).At this point,an external axial forcebf1is applied to the specimen,while the lateral force (f2?f3)is kept constant.Now,let us imagine that an additionalsmall amount of axial loading Df1is applied to the

specimen at time t,over a time range Dt:_ f1?Df1=Dt.This loading causes the system to evolve,in a way such that _U1>0.Under such a condition,the internal axial stress P1cannot exceed the peak valuebP1?bf1=A1.Precisely at the peak,_P1?0.As a consequence,Eq.(16)yields the following:

Finally,

At time t,the system is at rest,with Ec(t)?0.Thus,Ec(ttDt)is strictly positive.The specimen's kinetic energy increases from zero to a positive value Df1_U1Dt=2 over the time interval [t,ttDt].The internal axial stress P1within the specimen no longer balances the external axial force f1: with the increase of f1,P1follows a constitutive path and decreases frombP1along the descending branch.The unbalanced axial stress is responsible for the dynamic response of the specimen,characterized by a strictly positive value of Ec(ttDt)in Eq.(23).

Moreover,the failure mechanism of the specimen is initiated when the plastic limit state is reached,withthe internal second-order work is therefore nil.

Thus,the plastic limit theory appears to be a special situation of a more general second-order work theory.Failure can occur when the plastic limit condition is met.This situation is properly described by the second-order work theory (Darve et al.,2004; Nicot et al.,2009).However,this theory states that failure modes also exist and can be encountered before the plastic limit is reached.The test along the undrained triaxial loading path illustrates this situation perfectly well.Recalling that

5.Closing remarks

This paper has investigated the issue of failure in rateindependent materials.Basically,the occurrence of failure is defined by an abrupt increase in kinetic energy.Starting from this physical evidence,the approach developed in this paper is based on energy conservation,leading to a basic equation that introduces both external and internal second-order works.The increase in kinetic energy from an equilibrium state,under incremental loading,is shown to be equal to the difference between the external second-order work,involving the external loading parameters,and the internal second-order work,involving the constitutive properties of the material.The purpose of this study was to investigate how the externalsecond-order work can be greater than the internal secondorder work,leading to an increase in kinetic energy.The mechanical reason for this involves a distinction between the internal stress within the material and forces or stresses applied to the boundary of the system.When a stress limit state is reached,a certain stress component passes through a maximum value and then may decrease.This feature corresponds to a vanishing or a negative value of the internal second-order work.At this point,or along the descending branch,if a certain additional external loading is applied,the system fails,sharply increasing the strain rates.The internal stress no longer balances the external stress,leading to a dynamic response of the specimen.

This theoretical framework was applied to the well-known undrained triaxial test for loose soils.Although the sudden collapse observed after the deviatoric stress peak remains unclarified by the classic theory for the plastic limit state,the second-order work theory is able to describe such an event perfectly.In addition,the plastic limit theory appears to be a particular case of this more general second-order work theory.When the plastic limit condition is met,the internal secondorder work is nil.Thus,a class of incremental external loadings exists that causes the kinetic energy to increase dramatically.

Acknowledgements

The authors would like to express their sincere thanks to the French Research Network MeGe (Multiscale and Multiphysics Couplings in Geo-environmental Mechanics GDR CNRS 3176/2340,2008e2015)for having supported this work.

References

Biarez,J.,Hicher,P.Y.,1994.Elementary Mechanics of Soil Behaviour,Saturated Remoulded Soils.A.A.Balkema,Rotterdam.

Castro,G.,1969.Liquefaction of sands.In: Harvard Soil Mechanics Series,81.Harvard University Press,Cambridge.

Chu,J.,Leroueil,S.,Leong,W.K.,2003.Unstable behavior of sand and its implication for slope instability.Can.Geotech.J.40(5),873e885.

Darve,F.,1990.The expression of rheological laws in incremental form and the main classes of constitutive equations.In: Darve,F.,ed.,Geomaterials Constitutive Equations and Modelling.Elsevier,Amsterdam,pp.123e148.

Darve,F.,Servant,G.,Laouafa,F.,Khoa,H.D.V.,2004.Failure in geomaterials,continuous and discrete analyses.Comput.Fail.Mech.Geomater.193(27e29),3057e3085.http://dx.doi.org/10.1016/j.cma.2003.11.011.

Darve,F.,Sibille,L.,Daouadji,A.,Nicot,F.,2007.Bifurcations in granular media: macro- and micro-mechanics approaches.Comptes Rendus Mecanique 335(9e10),496e515.http://dx.doi.org/10.1016/j.crme.2007.08.005.

Hill,R.,1958.A general theory of uniqueness and stability in elastic-plastic solids.J.Mech.Phys.Solids 6(3),236e249.http://dx.doi.org/10.1016/0022-5096(58)90029-2.

Hill,R.,1967.The essential structure of constitutive laws for metal composites and polycrystals.J.Mech.Phys.Solids 15(2),79e95.http://dx.doi.org/10.1016/0022-5096(67)90018-X.

Lade,P.V.,Pradel,D.,1990.Instability and flow of granular materials,I: Experimental observations.J.Eng.Mech.116(11),2532e2550.

Lade,P.V.,1992.Static instability and liquefaction of loose fine sandy slopes.J.Geotech.Eng.118(1),51e71.http://dx.doi.org/10.1061/(ASCE)0733-9410(1992)118:1(51).

Nicot,F.,Darve,F.,2007.A micro-mechanical investigation of bifurcation in granular materials.Int.J.Solids Struct.44(20),6630e6652.http://dx.doi.org/10.1016/j.ijsolstr.2007.03.002.

Nicot,F.,Darve,F.,Khoa,H.D.V.,2007.Bifurcation and second-order work in geomaterials.Int.J.Numer.Anal.Methods Geomech.31(8),1007e1032.http://dx.doi.org/10.1002/nag.573.

Nicot,F.,Sibille,L.,Darve,F.,2009.Bifurcation in granular materials: an attempt at a unified framework.Int.J.Solids Struct.46(22e23),3938e3947.http://dx.doi.org/10.1016/j.ijsolstr.2009.07.008.

Nicot,F.,Sibille,L.,Darve,F.,2012.Failure in rate-independent granular materials as a bifurcation toward a dynamic regime.Int.J.Plasticity 29,136e154.http://dx.doi.org/10.1016/j.ijplas.2011.08.002.

Sibille,L.,Donze,F.,Nicot,F.,Chareyre,B.,Darve,F.,2008.Bifurcation detection and catastrophic failure.Acta Geotecnica 3(1),14e24.

Wan,R.G.,Pinheiro,M.,Guo,P.J.,2011.Elastoplastic modelling of diffuse instability response of geomaterials.Int.J.Numer.Anal.Methods Geomech.35(2),140e160.http://dx.doi.org/10.1002/nag.921.