- Email: [email protected]
An improved numerical integration algorithm for elastoplastic constitutive equations
Computers and Geotechnics 64 (2015) 1–9
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Research Paper
An improved numerical integration algorithm for elastoplastic constitutive equations Diogo L. Cecílio a,⇑, Philippe R.B. Devloo a,1, Sônia M. Gomes b, Erick R.S. dos Santos c, Nathan Shauer a,1 a
Universidade Estadual de Campinas, FEC, A. Albert Eintein, 951 Cidade Universitária Zeferino Vaz, Campinas, SP CEP 13083-970, Brazil Universidade Estadual de Campinas, IMECC, Campinas, SP, Brazil c Petrobras, Cenpes, Rio de Janeiro, RJ, Brazil b
a r t i c l e
i n f o
Article history: Received 24 April 2014 Received in revised form 3 September 2014 Accepted 18 October 2014
Keywords: Spectral decomposition Plasticity Geomechanics Return-mapping algorithm
a b s t r a c t A simplified methodology is proposed for elastoplastic calculations which holds for associative models. It is based on the representation of the elastoplastic model based on a rotation of the principal stresses and the fact that, in such system of coordinates, the direction that minimizes the square of a distance between a trial stress and the plastic surface has the same direction as the plastic deformation evolution. Such an approach allows for the elastoplastic calculation of complex models to be simpler and more efficient computationally. The proposed methodology is verified by the application to the elastoplastic model of Sandler–DiMaggio. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction The mechanical behavior of granular materials is very complex, involving plasticity by hydrostatic pressure, differences in the resistance for triaxial compression and triaxial stress, and porosity dependency, among other factors. In order to capture these effects, advanced constitutive models are necessary. Elastoplastic constitutive models are based on a yield criterion, a flow rule and a hardening law, giving rise to an initial value problem with restrictions. Due the complexity of realistic constitutive models, stiff non-linear systems are generated and the use of efficient numerical integration methods are required. An approach for the computation of elastoplastic problems is described in the book [1], Chapter 7, where the numerical integration is divided into two main steps: the elastic trial step and the plastic corrector step (or return-mapping algorithm). If trial stress computed in the first step fails to verify the plastically admissible condition, it is projected onto the yield surface by the returnmapping algorithm. In the present paper we propose a simplified implementation of the plastic corrector step procedure, which holds for associative
⇑ Corresponding author. Tel.: +55 19 3521 2396. E-mail addresses: [email protected] (D.L. Cecílio), [email protected] (P.R.B. Devloo), [email protected] (S.M. Gomes), [email protected] (E.R.S. dos Santos), [email protected] (N. Shauer). 1 Tel.: +55 19 3521 2396. http://dx.doi.org/10.1016/j.compgeo.2014.10.013 0266-352X/Ó 2014 Elsevier Ltd. All rights reserved.
models, reducing the number non-linear equations to be solved and improving efficiency. Using an iterative method for solving a system of non-linear elastoplastic equations, the convergence behavior is strongly dependent on the choice of variables to represent the residual vector. In this sense, an elastoplastic model can be better represented in terms of the principal stresses. Furthermore, in plastic calculations numerical instabilities can occur, in which small variations in the representation of floating point numbers may decide the feasibility of pursuing the calculations. Therefore, in addition to modeling the problem in the space of the principal stresses, we propose a distance function that has to be minimized in order to determine the projection of the plastic deformation on the yield surface. This geometric plastic calculation is possible thanks to the representation of the yield surface within a Haigh–Westergaard cylindrical coordinates [2], and the application of a rotation to a state of intermediate pressure. This procedure avoids the work in threedimensional space, with a reduction in the number of unknowns of the elastoplastic problem. Thus, the convergence of the iterative method can be easily obtained, stabilizing the calculation, with a consequent significant savings in computational time. To illustrate the effect of using this new approach, we consider the elastoplastic model described in the article [3], where the closure of the cap and the surface fault model are modified to allow the adjustment of the surface so that it is possible to characterize the material in triaxial compression and triaxial traction. For the results presented here, the implementation uses the NeoPZ library
2
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
[4,5], which is an object-oriented programming environment providing a framework for developing finite element simulations, augmented with special classes required for the integration of elastoplastic constitutive models. 2. Constitutive elastoplastic model The total deformation tensor e can be divided into two parts: e ¼ ee þ ep , an elastic part ee and a plastic part ep . The free energy u is also divided into portions of elastic ue ðe ep Þ and plastic contributions up ðaÞ, in which a is the internal damage variable. The e @@uee , in which q is the speelastic law establishes the tensor r ¼ q cific mass in the reference configuration. The plastic portion is not related to the strain state of the material; instead, it is related to the history of irreversible dissipative processes to which the material was submitted based on three fundamental axioms: an yield criterion, a flow rule, and a hardening law. Yield criterion. Describes the transition between the elastic and plastic domains through a plasticity function U ¼ Uðr; AÞ, where A ¼ q@ up =@ a is the hardening thermodynamic force. The plasticity function assumes non-positive values in an elastic basis and null values in a plastic basis. Flow rule. Assumes the existence of a plastic potential function W ¼ Wðr; AÞ, which specifies how the plastic deformation tensor ep evolves in a plasticity process e_ p ¼ c_ N, in which Nðr; AÞ ¼ @ W=@ r is the flow direction, and cðtÞ is a plastic multiplier. Hardening law. Specifies how the internal damage variable evolves a_ ¼ c_ H, in which Hðr; AÞ ¼ @ W=@A is the hardening modulus. In summary, the elastic–plastic constitutive model is formed by the following initial value problem: given the initial values ep ðt0 Þ and aðt0 Þ and the history of the infinitesimal deformation tensor eðtÞ; t 2 ½t0 ; T, to find the functions that define the plastic deformation tensor ep ðtÞ, the internal damage variable aðtÞ and a plastic multiplier c_ ðtÞ that meet the constitutive elastoplastic equations
e_ p ¼ c_ N a_ ¼ c_ H
ð1Þ
with the restrictions c_ ðtÞ P 0; UðrðtÞ; AðtÞÞ 0, c_ ðtÞUðrðtÞ; AðtÞÞ ¼ 0 in each (pseudo) instant t 2 ½t 0 ; T. 2.1. Algorithm for solving the incremental elastoplastic constitutive problem For the integration of elastoplastic non-linear systems, the use efficient numerical integration methods is required. Using the implicit Euler method at a step of (pseudo) time ½tn ; tnþ1 of a loading cycle, given a deformation state en and the corresponding plastic deformation ep;n and the internal state variable an at t n , for a prescribed incremental strain De, then the plastic deformation ep;nþ1 , the internal variable anþ1 and Dc at the next step are obtained as a solution of the problem that consists of the incremental non-linear system of equations
ee;nþ1 ¼ ee;n þ De DcN nþ1 anþ1 ¼ an þ DcH nþ1 for the unknowns
Dc P 0;
ð2Þ
ee;nþ1 ; anþ1 and Dc, subjected to the restrictions
Uðrnþ1 ; AÞ 0;
DcUðrnþ1 ; AÞ ¼ 0
ð3Þ
As shown in [1], the imposition of restrictions suggests a procedure for solving the problem in two major steps. It begins with a purely elastic predictor process (elastic trial step), with Dc ¼ 0. In this case, trial elastic strain eetrial ¼ ee;n þ De and internal variables
atrial ¼ an are defined. Then rtrial is calculated according to eetrial , and the corresponding Uðrtrial ; AÞ is given. If Uðrtrial ; AÞ 6 0, already a valid solution to the system is reached, and the variables are updated by the trial ones. Otherwise, a plastic corrector step (also known as plastic return-mapping scheme) is performed reformulating the incremental problem searching ee;nþ1 ; anþ1 and Dc satisfying
ee;nþ1 ¼ eetrial DcNðrnþ1 ; AÞ anþ1 ¼ atrial þ DcHðrnþ1 ; AÞ Dc > 0; Uðrnþ1 ; AÞ ¼ 0
ð4Þ ð5Þ ð6Þ
Next, the plastic strain is updated
ep;nþ1 ¼ ep;n þ De Dee For the specific resolution of the initial elastoplastic value problem, five main classes of tools are available at NeoPZ environment: Tensorial: implements tensor in tree dimensions. Elastic response: implements the elastic response of an isotropic material. Yield criterion: implements the plastic function U, the plastic flow vector N, and the hardening modulus H. Thermodynamic hardening force: implements the calculation of the thermodynamic force A. Incremental stress calculation: implements the solution of the elastoplastic initial value problem using the Newton’s method. 3. Plastic return-mapping scheme using rotated principal stresses In the study of plasticity, instead of using the six stress independent components for the geometric representation of a state of stress at a point, a simplified alternative is to calculate the principal stresses r ¼ ½r1 ; r2 ; r3 T as coordinates. This space is called Haigh–Westergaard stress space (HW). Furthermore, the constitutive law may be simplified by the introduction of a new coordinate system of rotated principal variables, similar to the decompositions defined in [6,7]. 3.1. Haigh–Westergaard stress space According to [2], the stress tensor r may be represented in terms of the principal stresses sorted in descending order r1 > r2 > r3 , which are given by the equations
qffiffi
2
3
2 p1ffiffi n þ 3 q cosðbÞ 3 3 r1 7 6 q ffiffi 6 7 7 1 2 2p 7 r¼ 6 4 r2 5 ¼ 6 6 pffiffi3 n þ 3q cos b 3 7 5 4 qffiffi r3 2 p1ffiffi n þ q cos b þ 23p 3 3
2
ð7Þ
in terms of the hidrostatic and deviatoric components, and the Lode angle
I1 n ¼ pffiffiffi ; 3
q¼
pffiffiffiffiffiffiffi 2J 2 ;
b¼
! pffiffiffi 1 3 2 J3 cos1 3 2 J 3=2 2
ð8Þ
which is only valid for the b 2 0; p3 ; I1 ; J 2 and J3 being the first invariant of the stress tensor, the second and third invariants of the deviatoric stress tensor, respectively. The stress strain relation is given by
2
3
2
3
e1 r1 7 7 1 6 e¼6 4 e2 5 ¼ ðDHW Þ 4 r2 5 e3 r2
ð9Þ
3
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
S ðaÞ ¼ fr : U ðr ; AðaÞÞ ¼ 0Þg
where
3 2 K þ 4G K 2G K 2G 3 3 3 7 6 ¼ 4 K 2G K þ 4G K 2G 3 3 3 5 K 2G K 2G K þ 4G 3 3 3
DHW
ð10Þ
K and G being the Bulk Modulus and the Shear Modulus, respectively. Conversely, the representation of the first invariant of the stress tensor I1 , the second invariant of the deviatoric tensor J 2 , and the Lode Angle b are given by the relations
I1 ¼ r1 þ r2 þ r3
ð11Þ
J 2 ¼ 1=3 r21 þ r22 r2 r3 þ r23 r1 ðr2 þ r3 Þ
ð12Þ
2
e;nþ1 ee ¼ DcN ;nþ1 trial e
ð21Þ
where N ;nþ1 ¼ @ U ðr;nþ1 ; Aðanþ1 ÞÞ=@ r . Therefore, for all possible variations dr in the tangent plane at r;nþ1 2 S ðanþ1 Þ,
1 e;nþ1 0 ¼ de ; ee ¼ dr ; DRHW rtrial r;nþ1 trial e
¼ dr ; rtrial r;nþ1 RHW
3
ð13Þ 3.2. Rotated Haigh–Westergaard space
where for vectors a and b in the RHW space, we define the internal product
ha ; b iRHW ¼ a ; D1 RHW b
ð23Þ
Property (22) means that r;nþ1 is the projected point in 2 S ðanþ1 Þ minimizing the distance induced by the RHW internal product of the trial stress. That is, if
Similarly to the procedure used in [6,7], a rotated Haigh–West T ergaard (RHW) representation r ¼ r1 ; r2 ; r3 may be considered, which is given by the relation
( ) 2 2 2 1=2 a1 b1 a2 b2 a3 b3 dða ; b Þ :¼ þ þ 3K 2G 2G qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ ha b ; a b iRHW
2
then
3
2
3
r1 r1 6 7 6 7 4 r2 5 ¼ R4 r2 5 r3 r3
ð14Þ
where
3
p1ffiffi p1ffiffi p1ffiffi 3 3 6 q3ffiffi 7 6 2 7 1 R ¼ 6 3 pffiffi p1ffiffi 7 6 65 4 p1ffiffi 0 p1ffiffi2 2
ð15Þ
In this rotated coordinate system, the invariants expressions become simpler
pffiffiffi I1 ¼ 3r1
r þr 2 2
2 3
2 b ¼ arctan r3 =r2
e
1 2 3
3
ð25Þ
d rtrial ; r;nþ1 ¼ minr 2S ðanþ1 Þ d rtrial ; r
2
r
1 2 2
r
e
ð17Þ
S ¼ fr ðr 1 ; r 2 Þ : ðr 1 ; r2 Þ 2 I g
3K 6 ¼4 0 0
0 2G 0
3 0 7 0 5 2G
Therefore, the square distance of a given rtrial to a point r ðr1 ; r2 Þ 2 S can be expressed as function of two parameters 2
d
3 ð19Þ
rtrial ; r ¼ d rtrial ; r1 ; r2
Given
e ee trial ; atrial , compute the stress rtrial ¼ DRHW etrial and do:
1. Projection: compute r;nþ1 ðr 1 ; r 2 Þ 2 S by minimizing the RHW ;nþ1 distance of rtrial to the admissibility surface S . That is, by solving the nonlinear system of two equations:
were
2
Motivated by the previous considerations, the incremental plastic corrector step (4)–(6) may be performed by a simplified procedure. In the following, it is described for perfectly plastic material and plastic material with hardening, separately.
ð16Þ
ð18Þ
6 7 7 1 6 4 e 5 ¼ ðDRHW Þ 4 r 5
DRHW
ð24Þ
3.4.1. Case 1: Perfectly plastic material In such case, the three dimensional surface S does not depend on a. Let it be parametrized by two parameters
The stress–strain relation is also simplified by
2
3.4. Solving the incremental plastic corrector step
2
J2 ¼
ð22Þ
! 1 ð2r1 r2 r3 Þðr1 2r2 þ r3 Þðr1 þ r2 2r3 Þ b ¼ arccos 3=2 3 2 r2 þ r2 r2 r3 þ r2 r1 ðr2 þ r3 Þ 1
where U ðr ; AÞ denotes the plastic function expressed in terms of the RHW principal stresses r . Using the RHW coordinate system, the incremental equation for the plastic law (4) reads
ð20Þ
3.3. Incremental flow rule for associative models as a closest point projection in the RHW space For associative plastic models, with W ¼ U, the law of plastic flow may be interpreted via a distance minimization of the trial stress to the admissibility surface [8,9]. For a given variable a, let the admissibility surface be defined by
@d
rtrial ; r1 ; r2
¼0 @r 1 @d rtrial ; r 1 ; r 2 ¼0 @r 2
ð26Þ
ð27Þ
2. If needed, Dc can be computed based on the equation
ee;nþ1 ee;nþ1 ¼ DcN ;nþ1 trial
ð28Þ
by taking the ratio of the norm of both tensors. Equation (28) can be used to check the validity of the projected point.
4
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
3.4.2. Case 2: Plastic material with hardening In the presence of a hardening rule, points on the plastic surface depends on three parameters: one parameter a to define the position of the surface, and two parametrization parameters r 1 and r 2 defining a point on it. Precisely, the target stress is of the form r;nþ1 ¼ r;nþ1 ðr1 ; r2 ; anþ1 Þ. Therefore, the return-mapping scheme is described as a nonlinear system of three equations:
@d rtrial ; r1 ; r 2 ; anþ1 ¼0 @r 1 @d rtrial ; r1 ; r 2 ; anþ1 ¼0 @r 2 anþ1 atrial þ DcH ðr1 ; r2 ; anþ1 Þ ¼ 0
ð29Þ ð30Þ pffiffiffiffi Fig. 1. Original Sandler–DiMaggio yield profile in the I1 ; J2 plane.
ð31Þ
where H ðr1 ; r 2 ; aÞ denotes the hardening modulus expressed in terms of the rotated principal stresses. The value of Dc, as given in (28), is a function of the three unknowns. There are some important implications of this simplified procedure: It reduces the resolution of the non-linear system (4)–(6), with seven equations, to a non-linear system with at most three equations, corresponding to the minimization formality for computing r;nþ1 , and a single equation for the variable a. In case the plasticity surface does not depend upon the Lode angle (e.g. Von Mises or classic Sandler–DiMaggio models), the closest point to the admissibility surface will have the same b angle. This means that r;nþ1 can be found by varying only one parameter. For some plasticity surfaces, the closest point can be computed analytically. This is the case for the Mohr Coulomb criterion (projection on a plane) or the Von Mises criterion (projection on a circle). For those cases, the return-mapping scheme is reduced to the resolution of a single nonlinear equation. When the guess of the projected point is sufficiently close to the target point, the distance function is a convex function of the variables which parametrize the surface. Therefore, algorithms oriented to minimizing the distance of a point to a surface are more stable than the Newton’s method for the resolution of a general non-linear system of equations. In the next section, this procedure is detailed for a generalized Sandler–DiMaggio plastic model.
4. Application: generalized Sandler–DiMaggio plastic model
where
pffiffiffiffi pffiffiffiffi F ðI Þ f 1 F 1 I1 ; J2 ; b ¼ J2 CðbÞ !2 pffiffiffiffi pffiffiffiffi I L 2 J 2 CðbÞ 1 þ 1 F 2 I1 ; J 2 ; L; b ¼ RF f ðLÞ F f ðLÞ
U¼
pffiffiffiffi F 1 I1 ; J2 ; b ; I1 > L pffiffiffiffi F 2 I1 ; J 2 ; L; b ; L P I1 P X
ð34Þ
with
F f ðsÞ ¼ A C expðBsÞ
ð35Þ
and
X ¼ L RF f ðLÞ
ð36Þ
Here the factor CðbÞ is given by
CðbÞ ¼
1 1 ð1 þ sinð3bÞÞ þ ð1 sinð3bÞÞ 2 w
ð37Þ
where w 2 ½7=9; 9=7. The original Sandler–DiMaggio model refers to w ¼ 1, such that CðbÞ ¼ 1. As reported in [3], the values of X and the plastic volumetric strain are related by
epv ¼ W½expðDXÞ 1
ð38Þ
The constants A; B; C; D; R and W are material properties, obtained from laboratory tests. Using representation in principal directions, Fig. 2 shows the aspect of the yield surface for w ¼ 7=9. Cross sections of the yield surfaces by the plane p, for w ¼ 1; w ¼ 7=9 and w ¼ 9=7 are plotted in Fig. 3 to illustrate the effect of different values of CðbÞ on their formats. In the original Sandler–DiMaggio formulation [3], the model is assumed to be associative, with W ¼ U. Furthermore, the hardening modulus and the flow direction are related by
The original Sandler–DiMaggio plastic model was proposed in [3]. It was initially conceived for granular soils, and currently it is widely applied in the oil industry to represent the behavior of rocks at depth. Let L be a parameter whose value is generally negative. The plasticity function U is piecewisely defined by a failure function pffiffiffiffi pffiffiffiffi F 1 I1 ; J 2 , and a cap elliptical function F 2 I1 ; J 2 ; L . A typical 2D profile of the corresponding yield surface is plotted in Fig. 1. Note that L specifies the distance of the origin to the center of the ellipsis. Inspired by the articles [10–12], to control the plastic surface format, so that the response of the material under triaxial compression may be different from its response in triaxial tensile, we consider a generalized Sandler–DiMaggio model, introducing an additional dependence on the Lode angle b. Precisely, the plastic function is defined by the expression
(
ð33Þ
ð32Þ Fig. 2. 3D generalized Sandler–DiMaggio yield surface with w ¼ 7=9.
5
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
H ¼ trðNÞ
ð39Þ
Using this property, the equation for the damage variable at the plastic correction step can be reformulated as
a_ ¼ c_ trðNÞ ¼ trðe_ p Þ ¼ e_ pv
ð40Þ
When using the classical Euler backward algorithm, the incremental plastic volumetric change leads to the following expression
DI 1 3K
Depv ¼
ð41Þ
which, through the relations (36) and (38), defines the evolution of Lepv Þ. This derivation allows the elimination of Dc as a variable in the plastic correction step. Fig. 4. Yield surface represented in the RHW stress space.
4.1. RHW parametrization of the yield surface In the RHW space, the generalized Sandler–DiMaggio yield surface is r1 -axially oriented, as illustrated in Fig. 4. It is formed by two pieces. There is one perfectly plastic portion S 1 , where F 1 ¼ 0, and the cap part S 2 ðLÞ, corresponding to F 2 ¼ 0. The stresses on the failure surface S 1 can be represented in terms of the variables r 1 ¼ I1 and r2 ¼ b. That is,
2
3
2
3
I1ffiffi p
3 r1 pffiffi 6 7 6 7 6 cosðbÞ 2F f ðI1 Þ 7 r ¼ 4 25 6 CðbÞ 7 4 5 pffiffi r3 2F ðI Þ sinðbÞ f 1
ð42Þ
return-mapping should be designed to consider the two cases described in the previous section. The RHW square distance from a given rtrial to a point r on the failure surface S 1 has the form
2 d rtrial ; r ¼
r
þ
Similarly, the RHW representation of any point in the cap surface S 2 ðLÞ can be parametrized in terms of two parameters r 1 ¼ h and r 2 ¼ b. Indeed,
3
r 6 7 4r 5 ¼ r 1 2 3
2
rtrial;1 pI1ffiffi3
CðbÞ
2
ð43Þ
4.2. Plastic return-mapping scheme using rotated principal stresses Since the yield surface for the Sandler–DiMaggio model is a combination of a perfect plastic part S 1 and a cap S 2 ðLÞ, the plastic
trial;3
þ
rtrial;2 cosðbÞ
pffiffi 2 2F f ðI1 Þ sinðbÞ CðbÞ
2G
pffiffi 2 2F f ðI1 Þ CðbÞ
2G :¼ d1
rtrial ; I1 ; b
ð44Þ
Therefore, the procedure to find rproj ¼ rproj ðI1 ; bÞ 2 S 1 that minimizes the distance to rtrial requires the solution of the system
2
3
LþR F f ðLÞ cosðhÞ pffiffi 3 7 6 pffiffi 7 6 6 cosðbÞ 2F f ðLÞ sinðhÞ 7 7 6 CðbÞ 5 4 pffiffi 2F ðLÞ sinðhÞ sinðbÞ fCðbÞ
3K
2
6 R1 ¼ 4
@d1 ðrtrial ;I1 ;bÞ
3
0 7 ¼ 5 @d1 ðrtrial ;I1 ;bÞ 0 @I1
ð45Þ
@b
On the other hand, the RHW square distance from rtrial to a point r 2 S 2 ðLÞ is given by pffiffi 2 2 r1;trial LþR F f pðLÞffiffi3cosðhÞ r2;trial cosðbÞ 2F fCðLÞsinðhÞ ðbÞ d2 rtrial ;h; b; L ¼ þ 3K 2G pffiffi 2 2F f ðLÞsinðhÞ r3;trial sinðbÞ CðbÞ þ :¼ d2 rtrial ; h; b; L ð46Þ 2G Furthermore, Eq. (41), rewritten as
I1
rtrial I1 ðr Þ ¼ 3K epv ðLÞ epv ðLtrial Þ
ð47Þ
is used to update simultaneously the position of the yield surface, defined by the variable L. Remembering the RHW representation (16), the invariant I1 ðr Þ, for r 2 S2 ðLÞ, has the expression
I1 ðr Þ ¼ I1 ðh; LÞ ¼ L þ RF f ðLÞ cosðhÞ
ð48Þ
As a conclusion, given Ltrial and rtrial , the procedure to compute L and rproj 2 S 2 ðLÞ that minimizes the RHW distance to rtrial requires the application of Newton’s method to find the values of h; b and L such that
3 2 3 @d2 ðrtrial ;h;b;LÞ 0 @h 7 6 7 6 7 6 R2 ¼ 6 @d2 ðrtrial ;h;b;LÞ 7 ¼ 4 0 5 5 4 @b 0 resðh; LÞ 2
ð49Þ
where Fig. 3. Yield surface cross sections for w ¼ 7=9; w ¼ 1 and w ¼ 9=7.
resðh; LÞ ¼ I1 rtrial I1 ðh; LÞ 3K epv ðLÞ epv ðLtrial Þ ¼ 0
ð50Þ
6
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
Consequently, from an algorithmic point of view, the returnmapping scheme is performed in the following steps: e Given ee trial ; Ltrial , compute the stress rtrial ¼ DRHW etrial . If I1 rtrial < Lproj , then compute L and rproj 2 S 2 ðLÞ by solving the system (49), and accept them as updated variables. Otherwise, compute rproj 2 S 1 by solving the system (45). – If I1
rproj P Lproj , then accept rproj and Lproj as updated
variables. – Else, compute L and rproj on the intersection of S 1 and S2 ðLÞ by solving a system obtained from (49) by eliminating the first equation, and fixing the third one for h ¼ p=2 (i.e., by reducing the system to two equations for the two unknowns b; L), and accept them as updated variables.
Table 1 Material parameters for the McCormic Ranch Sand [3]. Parameter
Value
E (Young’s modulus) m (Poisson’s ratio) A B
100 ksi 0.25 (dimensionless) 0.25 ksi
C D R W
w
1
0.67 ksi 0.18 ksi
1
0.67 ksi 2.5 (dimensionless) 0.066 (dimensionless) 1 (dimensionless)
A systematic way for computing the derivative of the projected stress as a function of the total strain is included in Appendices A–C. 5. Example The example reproduces the test documented in [13]. A uniaxial compressive strain is applied in the axial direction and zero strain in the other directions, followed by uniaxial strain in the opposite direction. Fig. 5 shows the comparative results between the work of [13] and the present study. The material parameters that were used in this example are listed in Table 1. The load path is also illustrated in Fig. 6, where the stress state is represented as a funcpffiffiffiffi tion of I1 and J 2 . It is graphically recognized that during the compression, the return-mapping scheme involves the projection of the stress points on the cap function defined by F 2 ¼ 0, increasing the plastic damage variable. During the unloading path, the stress point is project on the surface F 1 ¼ 0. Initially, on the path A–B, the material loses volume due the void closure and, thus, it yield triggering the cap function. Then on the path B–C, during unloading the material response is elastic until it reaches a critical point and fails in shear, triggering the failure function, that is represented by the path C–D. Fig. 7 represents the evolution of the residual normal when solving the non-linear systems of equations to project rtrial on S 1 and S 2 ðLÞ, showing that the number of iterations needed for convergence is less than 5.
Fig. 6. Load path of the experiment represented as a function of I1 and
pffiffiffiffi J2 .
Fig. 7. Residual norm per iteration for each load step.
6. Conclusion
Fig. 5. Comparison between results of uniaxial loading conduced by [13] and the results obtained in this work.
An efficient formulation/implementation is presented for any associative elastoplastic models that can be written in terms of the invariants of the stress tensor. It is based on the formulation of the plastic return-mapping in terms of the invariants of the stress tensor and of the minimization of a distance function to the admissibility surface. The improvement obtained with the formulation is that the Newton method is applied to a smaller system of equations and that part of these equations amount to the minimization of a convex function When applied to a generalized Sandler–DiMaggio model, the methodology leads to an iterative procedure which converges in less than 5 iterations.
7
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
Appendix A. Systems to be solved and square distance derivatives
@ðI1 ; bÞ @R1 ¼ J1 1 @ rtrial @ rtrial
A.1. d1 ðI1 ; bÞ derivatives
2 2 @ d1 @R1 @I @ r ¼ 4 12 1 @ d1 @ rtrial
2
d
rtrial ; r ¼
rtrial;1 pI1ffiffi3
rtrial;2
pffiffi 2 2F f ðI1 Þ cosðbÞ CðbÞ
@I1 @d1 ðrtrial ;I1 ;bÞ
2 @ r1 @I1
:¼ d1
2G
@d1 ðrtrial ;I1 ;bÞ
@b@ r1
2G
pffiffi 2 2F f ðI1 Þ CðbÞ
rtrial;3 sinðbÞ
þ
6 R1 ¼ 4
þ
3K
2
2
rtrial ; I1 ; b
ðA:1Þ
6 J1 ¼ 4
7 5
ðA:2Þ
A.2. d2
@ 2 d1 ðrtrial ;I1 ;bÞ
@I21
@I1 @b
@ 2 d1 ðrtrial ;I1 ;bÞ
@ 2 d1 ðrtrial ;I1 ;bÞ
@b@I1
@b2
r
d2 ðrtrial ;h;b;LÞ ¼
3 7 5
@ r1 @b @ r2 @b @ r3 @b
5
ðB:5Þ
3 7 7 5
ðB:6Þ
@d2 ðrtrial ;h;b;LÞ
B.2. Surface F 2 ðh; bLÞ
@ rproj @ rproj @ðh; b; LÞ ¼ @ rtrial @ðh; b; LÞ @ rtrial
ðB:7Þ
3K
r3;trial sinðbÞ
2
r
þ
pffiffi 2 2F f ðLÞsinðhÞ 2;trial cosðbÞ CðbÞ
2G
pffiffi 2 2F f ðLÞsinðhÞ CðbÞ
ðA:4Þ
2G 3
@ 2 d2 ðrtrial ;h;b;LÞ
@h2 6 6 @ 2 d r ;h;b;L ð Þ 2 trial J2 ¼ 6 6 @b@h 4
@R2 @R2 drtrial þ dðh; b; LÞ þ Oðrtrial ; h; b; LÞ ¼ 0 @ rtrial @ðh; b; LÞ
@h@b @ 2 d2 ðrtrial ;h;b;LÞ
@resðh;LÞ @h
@b2 @resðh;LÞ @b
@ 2 d2 ðrtrial ;h;b;LÞ
1 @ðh; b; LÞ @R2 @R2 ¼ @ rtrial @ðh; b; LÞ @ rtrial
ðB:10Þ
@ðh; b; LÞ @R2 ¼ J1 2 @ rtrial @ rtrial
ðB:11Þ
@h@L
7 7 7 5
@ 2 d2 ðrtrial ;h;b;LÞ 7 @b@L
@resðh;LÞ : @L
resðh; LÞ ¼ I1 rtrial I1 ðh; LÞ 3K epv ðLÞ epv ðLtrial Þ ¼ 0
@ 2 d2
6 @h@r1 6 2 ¼ 6 @ d2 4 @b@ r1 @ DL @ r1
2 @ r1
ðA:6Þ
@h
6 @ rproj ¼ 6 @ r2 @ðh; b; LÞ 4 @h @ r3 @h
ðA:7Þ
Appendix B. Derivatives of the projected stress with respect to trial stresses B.1. Surface F 1 ðI1 ; bÞ
@ 2 d2 @h@ r2
@ 2 d2 @h@ r3
@ 2 d2 @b@ r2
@ 2 d2 @b@ r3
@ DL @ r2
@ DL @ r3
@ r1 @b
@ r1 @L
@ r2 @b
@ r2 @L
@ r3 @b
@ r3 @L
3 7 7 7 5
ðB:1Þ
ðB:12Þ
3 7 7 5
ðB:13Þ
Appendix C. Tangent operator The tangent operator to be calculated is given by the derivative of the projected stress tensor by the trial deformation tensor, and, by the chain rule, it is shown that:
@ rproj @ rproj rtrial ¼ @e @ rtrial @ e
For surface F 1 ðI1 ; bÞ the chain rule gives
@ rproj @ rproj @ðI1 ; bÞ ¼ @ rtrial @ðI1 ; bÞ @ rtrial
2 @R2 @ rtrial
3
ðB:9Þ
Rearranging,
ðA:5Þ
@ 2 d2 ðrtrial ;h;b;LÞ
ðB:8Þ
and perform a Taylor expansion on it, we have,
@h 7 6 7 6 R2 ¼ 6 @d2 ðrtrial ;h;b;LÞ 7 5 4 @b resðh; LÞ
2
@ 2 d1 @b@ r3
If we take the derivative of the distance below
LþRF f ðLÞcosðhÞ pffiffi 1;trial 3
r
þ
ðA:3Þ
derivatives
2
@ 2 d1 @b@ r2
3
R2 ¼ 0
trial ; h; b; L
@ 2 d1 @I1 @ r3
For surface F 2 ðh; b; LÞ the chain rule gives
@ 2 d1 ðrtrial ;I1 ;bÞ
@ r3 @I1
@ 2 d1 @I1 @ r2
3
@b
2
6 @ r2 @ rproj ¼6 @ðI1 ; bÞ 4 @I1
ðB:4Þ
ðC:1Þ
And by straightforward derivation, the tangent operator can be calculated by the following formula:
If we take the derivative of the residual equation R1 ¼ 0, and perform a Taylor expansion on it, we have,
X
@ rproj @ rproj;i ¼ ðv i v i Þ ðv j v j Þ C @e @ rtrial;j
@R1 @R1 drtrial þ dðI1 ; bÞ þ Oðrtrial ; I1 ; bÞ ¼ 0 @ rtrial @ðI1 ; bÞ
For the flowing calculations it is convenient to deduce a formulation for the tangent operator in terms of a matrix. To achieve this, first we multiply the tangent operator by a deformation de:
ðB:2Þ
Rearranging,
1 @ðI1 ; bÞ @R1 @R1 ¼ @ rtrial @ðI1 ; bÞ @ rtrial
ðB:3Þ
@ rproj de ¼ @e
X
ðC:2Þ
@ rproj;i ðv i v i Þ ðv j v j Þ 2lI4 þ kI2 I2 de @ rtrial;j ðC:3Þ
8
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
And with some mathematical manipulation Eq. (C.4) is obtained:
@ rproj de ¼ @e
X @ rproj;i h @ rtrial;j
i 2l Tr½ðv j v j Þde þ k I2 de ðv i v i Þ ðC:4Þ
In Voigt notation, the stress and deformation tensors are represented by six element vectors as shown below:
rxx 1 C 1 B B rxy C rxz B C B rxz C B C ryz C ¼ A B C B ryy C rzz B C @ ryz A rzz
1
0
1 0 0 B C dexx ¼ @ 0 0 0 A 0
ðC:5Þ
deyy
0 B ¼ @0 0
1
0
1
ðC:6Þ In Eq. (C.4), de is replaced by each of these matrices. The resultant matrix of each of these six operations represents the rate of variation of each of the stress tensor components in relationship with one of the deformations. Each column of the tangent matrix is composed by one of these operations. Next it is shown the detailed calculation of each of these columns:
@ rproj Column 1 : dexx @e
!
X @ rproj;i h ij
@ rtrial;j
2l
v 2j
Column 2 :
@ rproj dexy @e
X @ rproj;i ! 2 2l @ rtrial;j ij
Column 3 :
@ rproj dexz @e
!
Column 4 :
@ rproj deyy @e
Column 5 :
@ rproj deyz @e
@ rproj Column 6 : dezz @e in winch, k ¼ 1; 2; 3.
vj
k
1
C hx A 0
ðC:8Þ
And given the deformation tensor in expressed in principal directions:
1
e1 0 0 C e¼B @ 0 e2 0 A 0 0 e3
ðC:9Þ
1
e1 ðe1 e2 Þhz ðe1 e3 Þhy B C e2 ðe2 e3 Þhx A RT eR ¼ @ ðe1 e2 Þhz ðe1 e3 Þhy ðe2 e3 Þhx e3
ðC:10Þ
Thus for any deformation, in case of terms out of the principal axes, we have:
0 1 0 0 0 1 B C B C dexy ¼ @ 1 0 0 A dexz ¼ @ 0 0 0 A 0 0 0 0 0 1 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 C B C B C 1 0 A deyz ¼ @ 0 0 1 A dezz ¼ @ 0 0 0 A 0 0 0 1 0 0 0 1
0
0 hx
hy
0
Using this notations, the tangent matrix is a second order tensor of size 6 6. To calculated this matrix, first it is defined six tensor unitary directions:
0
B R ¼ I þ @ hz hy
hz
The rotation of this tensor, already elimination the order terms is given by:
0
rxx rxy r r¼B @ xy ryy rxz ryz
0
0
C.1. Tangent matrix in Voigt notation
0
0
1
þk
i
vi vi
h¼
X @ rproj;i 2 2l v j3 v j1 v i v i @ rtrial;j ij i X @ rproj;i h ! 2l v 2j2 þ k v i v i @ r trial;j ij
ðC:13Þ
In the same way as before, it will be calculated the rate of variation of each of the project stress tensor components in relation with a unitary strain direction. Each column of the correction of the tangent matrix is composed by one of these operations. Next it is shown the detailed calculation of each of these columns:
Column 1 : dexx ! Column 2 : dexy
ðC:7Þ Column 6 : dezz
X rproj;i rproj;j
v v
ðv v þ v v Þ
j 1 i1 i j j i ei ej ij X rproj;i rproj;j ! v j1 v i2 þ v j2 v i1 ðv i v j þ v j v i Þ e e i j ij X rproj;i rproj;j ! v j1 v i3 þ v j3 v i1 ðv i v j þ v j v i Þ e e i j ij X rproj;i rproj;j ! v j2 v i2 ðv i v j þ v j v i Þ ei ej ij X rproj;i rproj;j ! v j2 v i3 þ v j3 v i2 ðv i v j þ v j v i Þ e e i j ij X rproj;i rproj;j ! v j3 v i3 ðv i v j þ v j v i Þ e e i j ij
ðC:14Þ
!
is the k component of the eigenvector
rproj , the correction
rproj;i rproj;j v j ðde v i Þ ðv i v j þ v j v i Þ ei ej
Column 5 : deyz
X @ rproj;i 2 2l v j3 v j2 v i v i @ r trial;j ij i X @ rproj;i h ! 2l v 2j3 þ k v i v i @ rtrial;j ij
ðC:12Þ
And as v i and v j are eigenvector of both e and em this position is given by:
Column 4 : deyy
1
ðC:11Þ
ðei ej Þ
Column 3 : dexz
vj vj vi vi 2
v j ðde v i Þ ¼ hðei ej Þ v j ðde v i Þ
for i ¼ 1; 2 , j ¼ 2; 3 and i – j. C.3. Particular case of repeated eigenvalues
vj,
with
C.2. Correction of the eigenvectors variation In the deduction presented, it was not considered the variation of the eigenvectors between rproj and rtrial . This variation can be represented by a rotation, and its inclusion in the tangent matrix is demonstrated here. A rotation tensor can be seen as the identity tensor plus a part relative to the rotation in each of the axes as below:
In case of repeated eigenvalues, the correction becomes a relation 00 in the flowing equation present in the calculation of each of the columns;
rproj;i rproj;j ei ej
ðC:15Þ
The resolution of this problem will be exemplified for case i ¼ 1; j ¼ 2, with rproj;1 ¼ rproj;2 and e1 ¼ e2 . First, it can be proved that:
rtrial;1 rtrial;2 ¼ 2l e1 e2
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
Which is applied in Eq. (C.15).
References
rproj;1 rproj;2 2l rtrial;1 rtrial;2
ðC:16Þ
And so, the flowing limit is calculated with dr ! 0, and imposing a variation in the direction fdr; dr; 0g:
pr
@ r2 rproj;1 þ @@rrproj;1 fdr; dr; 0g rpr 2 þ @ rtrial fdr; dr; 0g trial lim 2l dr!0 ðrtrial;1 þ drÞ ðrtrial;2 drÞ
ðC:17Þ As
rproj;i ¼ rproj;j and rtrial;i ¼ rtrial;j :
lim 2l
@ rproj;1 @ rtrial;1
@r
@r
@r
dr @ rproj;1 dr @rproj;2 dr þ @ rproj;2 dr trial;2
trial;1
trial;2
ðC:18Þ
2dr
dr!0
And cancelling the terms dr, it is defined the final correction for the case i ¼ 1; j ¼ 2, in which rproj;1 ¼ rproj;2 and e1 ¼ e2 .
l
@ rproj;1 @ rproj;1 @ rproj;2 @ rproj;2 þ @ rtrial;1 @ rtrial;2 @ rtrial;1 @ rtrial;2
ðC:19Þ
Finally, the flowing substitution must be made in the calculation of the correction for each of the columns in case of repeated eigenvalues:
9
rproj;i rproj;j @ rproj;i @ rproj;i @ rproj;j @ rproj;j ¼l þ ei ej @ rtrial;i @ rtrial;j @ rtrial;i @ rtrial;j
ðC:20Þ
[1] Souza Neto E, Peric D, Owen DRJ. Computational Methods for Plasticity. John Wiley & Sons Ltd; 2008. [2] Chen WF, Han DJ. Plasticity for Structural Engineers. Springer-Verlag; 1988. [3] DiMaggio FL, Sandler I. Material model for granular soils. Eng Mech Div 1971;97(3):935–50. [4] Devloo PRB. PZ: an object oriented environment for scientific programming. Comput Methods Appl Mech Eng 1997;150(1-4):133–53. http://dx.doi.org/ 10.1016/S0045-7825(97)00097-2. [5] Devloo PRB. Object oriented tools for scientific computing. Eng Comput 2000;46(2):203–14. http://dx.doi.org/10.1016/0045-7949(93)90185-G. [6] Lainé E, Vallée C, Fortuné D. Nonlinear isotropic constitutive laws: choice of the three invariants, convex potentials and constitutive inequalities. Int J Eng Sci 1999;37(15):1927–41. http://dx.doi.org/10.1016/S0020-7225(99)00006-3. [7] Borja RI. Plasticity Modeling & Computation. Springer Science & Business Media; 2012. [8] Armero F, Pérez-Foguet A. On the formulation of closest-point projection algorithms in elastoplasticity – Part I: The variational structure. Int J Numer Methods Eng 2002;53(2):297–329. http://dx.doi.org/10.1002/nme.278. [9] Armero F, Pérez-Foguet A. On the formulation of closest-point projection algorithms in elastoplasticity – Part II: Globally convergent schemes. Int J Numer Methods Eng 2002;53(2):331–74. http://dx.doi.org/10.1002/nme.279. [10] Foster CD, Regueiro RA, Fossum AF, Borja RI. Implicit numerical integration of a three-invariant, isotropic/kinematic hardening cap plasticity model for geomaterials. Comput Methods Appl Mech Eng 2005;194:5109–38. http:// dx.doi.org/10.1016/j.cma.2005.01.001. [11] A.F. Fossum, J.T. Fredrich, Cap plasticity models and compactive and dilatant pre-failure deformation, in: 4th North American Rock Mechanics Symposium, Seatle, United States of America, 2000. [12] A.F. Fossum, R.M. Brannon, The Sandia Geomodel: Theory and Users Guide, SAND Report, Sandia National Laboratories, 2004. [13] Sandler IS, Rubin D. An algorithm and a modular subroutine for the cap model. Int J Numer Anal Methods Geomech 1979;3(2):173–86. http://dx.doi.org/ 10.1002/nag.1610030206.
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Research Paper
An improved numerical integration algorithm for elastoplastic constitutive equations Diogo L. Cecílio a,⇑, Philippe R.B. Devloo a,1, Sônia M. Gomes b, Erick R.S. dos Santos c, Nathan Shauer a,1 a
Universidade Estadual de Campinas, FEC, A. Albert Eintein, 951 Cidade Universitária Zeferino Vaz, Campinas, SP CEP 13083-970, Brazil Universidade Estadual de Campinas, IMECC, Campinas, SP, Brazil c Petrobras, Cenpes, Rio de Janeiro, RJ, Brazil b
a r t i c l e
i n f o
Article history: Received 24 April 2014 Received in revised form 3 September 2014 Accepted 18 October 2014
Keywords: Spectral decomposition Plasticity Geomechanics Return-mapping algorithm
a b s t r a c t A simplified methodology is proposed for elastoplastic calculations which holds for associative models. It is based on the representation of the elastoplastic model based on a rotation of the principal stresses and the fact that, in such system of coordinates, the direction that minimizes the square of a distance between a trial stress and the plastic surface has the same direction as the plastic deformation evolution. Such an approach allows for the elastoplastic calculation of complex models to be simpler and more efficient computationally. The proposed methodology is verified by the application to the elastoplastic model of Sandler–DiMaggio. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction The mechanical behavior of granular materials is very complex, involving plasticity by hydrostatic pressure, differences in the resistance for triaxial compression and triaxial stress, and porosity dependency, among other factors. In order to capture these effects, advanced constitutive models are necessary. Elastoplastic constitutive models are based on a yield criterion, a flow rule and a hardening law, giving rise to an initial value problem with restrictions. Due the complexity of realistic constitutive models, stiff non-linear systems are generated and the use of efficient numerical integration methods are required. An approach for the computation of elastoplastic problems is described in the book [1], Chapter 7, where the numerical integration is divided into two main steps: the elastic trial step and the plastic corrector step (or return-mapping algorithm). If trial stress computed in the first step fails to verify the plastically admissible condition, it is projected onto the yield surface by the returnmapping algorithm. In the present paper we propose a simplified implementation of the plastic corrector step procedure, which holds for associative
⇑ Corresponding author. Tel.: +55 19 3521 2396. E-mail addresses: [email protected] (D.L. Cecílio), [email protected] (P.R.B. Devloo), [email protected] (S.M. Gomes), [email protected] (E.R.S. dos Santos), [email protected] (N. Shauer). 1 Tel.: +55 19 3521 2396. http://dx.doi.org/10.1016/j.compgeo.2014.10.013 0266-352X/Ó 2014 Elsevier Ltd. All rights reserved.
models, reducing the number non-linear equations to be solved and improving efficiency. Using an iterative method for solving a system of non-linear elastoplastic equations, the convergence behavior is strongly dependent on the choice of variables to represent the residual vector. In this sense, an elastoplastic model can be better represented in terms of the principal stresses. Furthermore, in plastic calculations numerical instabilities can occur, in which small variations in the representation of floating point numbers may decide the feasibility of pursuing the calculations. Therefore, in addition to modeling the problem in the space of the principal stresses, we propose a distance function that has to be minimized in order to determine the projection of the plastic deformation on the yield surface. This geometric plastic calculation is possible thanks to the representation of the yield surface within a Haigh–Westergaard cylindrical coordinates [2], and the application of a rotation to a state of intermediate pressure. This procedure avoids the work in threedimensional space, with a reduction in the number of unknowns of the elastoplastic problem. Thus, the convergence of the iterative method can be easily obtained, stabilizing the calculation, with a consequent significant savings in computational time. To illustrate the effect of using this new approach, we consider the elastoplastic model described in the article [3], where the closure of the cap and the surface fault model are modified to allow the adjustment of the surface so that it is possible to characterize the material in triaxial compression and triaxial traction. For the results presented here, the implementation uses the NeoPZ library
2
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
[4,5], which is an object-oriented programming environment providing a framework for developing finite element simulations, augmented with special classes required for the integration of elastoplastic constitutive models. 2. Constitutive elastoplastic model The total deformation tensor e can be divided into two parts: e ¼ ee þ ep , an elastic part ee and a plastic part ep . The free energy u is also divided into portions of elastic ue ðe ep Þ and plastic contributions up ðaÞ, in which a is the internal damage variable. The e @@uee , in which q is the speelastic law establishes the tensor r ¼ q cific mass in the reference configuration. The plastic portion is not related to the strain state of the material; instead, it is related to the history of irreversible dissipative processes to which the material was submitted based on three fundamental axioms: an yield criterion, a flow rule, and a hardening law. Yield criterion. Describes the transition between the elastic and plastic domains through a plasticity function U ¼ Uðr; AÞ, where A ¼ q@ up =@ a is the hardening thermodynamic force. The plasticity function assumes non-positive values in an elastic basis and null values in a plastic basis. Flow rule. Assumes the existence of a plastic potential function W ¼ Wðr; AÞ, which specifies how the plastic deformation tensor ep evolves in a plasticity process e_ p ¼ c_ N, in which Nðr; AÞ ¼ @ W=@ r is the flow direction, and cðtÞ is a plastic multiplier. Hardening law. Specifies how the internal damage variable evolves a_ ¼ c_ H, in which Hðr; AÞ ¼ @ W=@A is the hardening modulus. In summary, the elastic–plastic constitutive model is formed by the following initial value problem: given the initial values ep ðt0 Þ and aðt0 Þ and the history of the infinitesimal deformation tensor eðtÞ; t 2 ½t0 ; T, to find the functions that define the plastic deformation tensor ep ðtÞ, the internal damage variable aðtÞ and a plastic multiplier c_ ðtÞ that meet the constitutive elastoplastic equations
e_ p ¼ c_ N a_ ¼ c_ H
ð1Þ
with the restrictions c_ ðtÞ P 0; UðrðtÞ; AðtÞÞ 0, c_ ðtÞUðrðtÞ; AðtÞÞ ¼ 0 in each (pseudo) instant t 2 ½t 0 ; T. 2.1. Algorithm for solving the incremental elastoplastic constitutive problem For the integration of elastoplastic non-linear systems, the use efficient numerical integration methods is required. Using the implicit Euler method at a step of (pseudo) time ½tn ; tnþ1 of a loading cycle, given a deformation state en and the corresponding plastic deformation ep;n and the internal state variable an at t n , for a prescribed incremental strain De, then the plastic deformation ep;nþ1 , the internal variable anþ1 and Dc at the next step are obtained as a solution of the problem that consists of the incremental non-linear system of equations
ee;nþ1 ¼ ee;n þ De DcN nþ1 anþ1 ¼ an þ DcH nþ1 for the unknowns
Dc P 0;
ð2Þ
ee;nþ1 ; anþ1 and Dc, subjected to the restrictions
Uðrnþ1 ; AÞ 0;
DcUðrnþ1 ; AÞ ¼ 0
ð3Þ
As shown in [1], the imposition of restrictions suggests a procedure for solving the problem in two major steps. It begins with a purely elastic predictor process (elastic trial step), with Dc ¼ 0. In this case, trial elastic strain eetrial ¼ ee;n þ De and internal variables
atrial ¼ an are defined. Then rtrial is calculated according to eetrial , and the corresponding Uðrtrial ; AÞ is given. If Uðrtrial ; AÞ 6 0, already a valid solution to the system is reached, and the variables are updated by the trial ones. Otherwise, a plastic corrector step (also known as plastic return-mapping scheme) is performed reformulating the incremental problem searching ee;nþ1 ; anþ1 and Dc satisfying
ee;nþ1 ¼ eetrial DcNðrnþ1 ; AÞ anþ1 ¼ atrial þ DcHðrnþ1 ; AÞ Dc > 0; Uðrnþ1 ; AÞ ¼ 0
ð4Þ ð5Þ ð6Þ
Next, the plastic strain is updated
ep;nþ1 ¼ ep;n þ De Dee For the specific resolution of the initial elastoplastic value problem, five main classes of tools are available at NeoPZ environment: Tensorial: implements tensor in tree dimensions. Elastic response: implements the elastic response of an isotropic material. Yield criterion: implements the plastic function U, the plastic flow vector N, and the hardening modulus H. Thermodynamic hardening force: implements the calculation of the thermodynamic force A. Incremental stress calculation: implements the solution of the elastoplastic initial value problem using the Newton’s method. 3. Plastic return-mapping scheme using rotated principal stresses In the study of plasticity, instead of using the six stress independent components for the geometric representation of a state of stress at a point, a simplified alternative is to calculate the principal stresses r ¼ ½r1 ; r2 ; r3 T as coordinates. This space is called Haigh–Westergaard stress space (HW). Furthermore, the constitutive law may be simplified by the introduction of a new coordinate system of rotated principal variables, similar to the decompositions defined in [6,7]. 3.1. Haigh–Westergaard stress space According to [2], the stress tensor r may be represented in terms of the principal stresses sorted in descending order r1 > r2 > r3 , which are given by the equations
qffiffi
2
3
2 p1ffiffi n þ 3 q cosðbÞ 3 3 r1 7 6 q ffiffi 6 7 7 1 2 2p 7 r¼ 6 4 r2 5 ¼ 6 6 pffiffi3 n þ 3q cos b 3 7 5 4 qffiffi r3 2 p1ffiffi n þ q cos b þ 23p 3 3
2
ð7Þ
in terms of the hidrostatic and deviatoric components, and the Lode angle
I1 n ¼ pffiffiffi ; 3
q¼
pffiffiffiffiffiffiffi 2J 2 ;
b¼
! pffiffiffi 1 3 2 J3 cos1 3 2 J 3=2 2
ð8Þ
which is only valid for the b 2 0; p3 ; I1 ; J 2 and J3 being the first invariant of the stress tensor, the second and third invariants of the deviatoric stress tensor, respectively. The stress strain relation is given by
2
3
2
3
e1 r1 7 7 1 6 e¼6 4 e2 5 ¼ ðDHW Þ 4 r2 5 e3 r2
ð9Þ
3
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
S ðaÞ ¼ fr : U ðr ; AðaÞÞ ¼ 0Þg
where
3 2 K þ 4G K 2G K 2G 3 3 3 7 6 ¼ 4 K 2G K þ 4G K 2G 3 3 3 5 K 2G K 2G K þ 4G 3 3 3
DHW
ð10Þ
K and G being the Bulk Modulus and the Shear Modulus, respectively. Conversely, the representation of the first invariant of the stress tensor I1 , the second invariant of the deviatoric tensor J 2 , and the Lode Angle b are given by the relations
I1 ¼ r1 þ r2 þ r3
ð11Þ
J 2 ¼ 1=3 r21 þ r22 r2 r3 þ r23 r1 ðr2 þ r3 Þ
ð12Þ
2
e;nþ1 ee ¼ DcN ;nþ1 trial e
ð21Þ
where N ;nþ1 ¼ @ U ðr;nþ1 ; Aðanþ1 ÞÞ=@ r . Therefore, for all possible variations dr in the tangent plane at r;nþ1 2 S ðanþ1 Þ,
1 e;nþ1 0 ¼ de ; ee ¼ dr ; DRHW rtrial r;nþ1 trial e
¼ dr ; rtrial r;nþ1 RHW
3
ð13Þ 3.2. Rotated Haigh–Westergaard space
where for vectors a and b in the RHW space, we define the internal product
ha ; b iRHW ¼ a ; D1 RHW b
ð23Þ
Property (22) means that r;nþ1 is the projected point in 2 S ðanþ1 Þ minimizing the distance induced by the RHW internal product of the trial stress. That is, if
Similarly to the procedure used in [6,7], a rotated Haigh–West T ergaard (RHW) representation r ¼ r1 ; r2 ; r3 may be considered, which is given by the relation
( ) 2 2 2 1=2 a1 b1 a2 b2 a3 b3 dða ; b Þ :¼ þ þ 3K 2G 2G qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ ha b ; a b iRHW
2
then
3
2
3
r1 r1 6 7 6 7 4 r2 5 ¼ R4 r2 5 r3 r3
ð14Þ
where
3
p1ffiffi p1ffiffi p1ffiffi 3 3 6 q3ffiffi 7 6 2 7 1 R ¼ 6 3 pffiffi p1ffiffi 7 6 65 4 p1ffiffi 0 p1ffiffi2 2
ð15Þ
In this rotated coordinate system, the invariants expressions become simpler
pffiffiffi I1 ¼ 3r1
r þr 2 2
2 3
2 b ¼ arctan r3 =r2
e
1 2 3
3
ð25Þ
d rtrial ; r;nþ1 ¼ minr 2S ðanþ1 Þ d rtrial ; r
2
r
1 2 2
r
e
ð17Þ
S ¼ fr ðr 1 ; r 2 Þ : ðr 1 ; r2 Þ 2 I g
3K 6 ¼4 0 0
0 2G 0
3 0 7 0 5 2G
Therefore, the square distance of a given rtrial to a point r ðr1 ; r2 Þ 2 S can be expressed as function of two parameters 2
d
3 ð19Þ
rtrial ; r ¼ d rtrial ; r1 ; r2
Given
e ee trial ; atrial , compute the stress rtrial ¼ DRHW etrial and do:
1. Projection: compute r;nþ1 ðr 1 ; r 2 Þ 2 S by minimizing the RHW ;nþ1 distance of rtrial to the admissibility surface S . That is, by solving the nonlinear system of two equations:
were
2
Motivated by the previous considerations, the incremental plastic corrector step (4)–(6) may be performed by a simplified procedure. In the following, it is described for perfectly plastic material and plastic material with hardening, separately.
ð16Þ
ð18Þ
6 7 7 1 6 4 e 5 ¼ ðDRHW Þ 4 r 5
DRHW
ð24Þ
3.4.1. Case 1: Perfectly plastic material In such case, the three dimensional surface S does not depend on a. Let it be parametrized by two parameters
The stress–strain relation is also simplified by
2
3.4. Solving the incremental plastic corrector step
2
J2 ¼
ð22Þ
! 1 ð2r1 r2 r3 Þðr1 2r2 þ r3 Þðr1 þ r2 2r3 Þ b ¼ arccos 3=2 3 2 r2 þ r2 r2 r3 þ r2 r1 ðr2 þ r3 Þ 1
where U ðr ; AÞ denotes the plastic function expressed in terms of the RHW principal stresses r . Using the RHW coordinate system, the incremental equation for the plastic law (4) reads
ð20Þ
3.3. Incremental flow rule for associative models as a closest point projection in the RHW space For associative plastic models, with W ¼ U, the law of plastic flow may be interpreted via a distance minimization of the trial stress to the admissibility surface [8,9]. For a given variable a, let the admissibility surface be defined by
@d
rtrial ; r1 ; r2
¼0 @r 1 @d rtrial ; r 1 ; r 2 ¼0 @r 2
ð26Þ
ð27Þ
2. If needed, Dc can be computed based on the equation
ee;nþ1 ee;nþ1 ¼ DcN ;nþ1 trial
ð28Þ
by taking the ratio of the norm of both tensors. Equation (28) can be used to check the validity of the projected point.
4
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
3.4.2. Case 2: Plastic material with hardening In the presence of a hardening rule, points on the plastic surface depends on three parameters: one parameter a to define the position of the surface, and two parametrization parameters r 1 and r 2 defining a point on it. Precisely, the target stress is of the form r;nþ1 ¼ r;nþ1 ðr1 ; r2 ; anþ1 Þ. Therefore, the return-mapping scheme is described as a nonlinear system of three equations:
@d rtrial ; r1 ; r 2 ; anþ1 ¼0 @r 1 @d rtrial ; r1 ; r 2 ; anþ1 ¼0 @r 2 anþ1 atrial þ DcH ðr1 ; r2 ; anþ1 Þ ¼ 0
ð29Þ ð30Þ pffiffiffiffi Fig. 1. Original Sandler–DiMaggio yield profile in the I1 ; J2 plane.
ð31Þ
where H ðr1 ; r 2 ; aÞ denotes the hardening modulus expressed in terms of the rotated principal stresses. The value of Dc, as given in (28), is a function of the three unknowns. There are some important implications of this simplified procedure: It reduces the resolution of the non-linear system (4)–(6), with seven equations, to a non-linear system with at most three equations, corresponding to the minimization formality for computing r;nþ1 , and a single equation for the variable a. In case the plasticity surface does not depend upon the Lode angle (e.g. Von Mises or classic Sandler–DiMaggio models), the closest point to the admissibility surface will have the same b angle. This means that r;nþ1 can be found by varying only one parameter. For some plasticity surfaces, the closest point can be computed analytically. This is the case for the Mohr Coulomb criterion (projection on a plane) or the Von Mises criterion (projection on a circle). For those cases, the return-mapping scheme is reduced to the resolution of a single nonlinear equation. When the guess of the projected point is sufficiently close to the target point, the distance function is a convex function of the variables which parametrize the surface. Therefore, algorithms oriented to minimizing the distance of a point to a surface are more stable than the Newton’s method for the resolution of a general non-linear system of equations. In the next section, this procedure is detailed for a generalized Sandler–DiMaggio plastic model.
4. Application: generalized Sandler–DiMaggio plastic model
where
pffiffiffiffi pffiffiffiffi F ðI Þ f 1 F 1 I1 ; J2 ; b ¼ J2 CðbÞ !2 pffiffiffiffi pffiffiffiffi I L 2 J 2 CðbÞ 1 þ 1 F 2 I1 ; J 2 ; L; b ¼ RF f ðLÞ F f ðLÞ
U¼
pffiffiffiffi F 1 I1 ; J2 ; b ; I1 > L pffiffiffiffi F 2 I1 ; J 2 ; L; b ; L P I1 P X
ð34Þ
with
F f ðsÞ ¼ A C expðBsÞ
ð35Þ
and
X ¼ L RF f ðLÞ
ð36Þ
Here the factor CðbÞ is given by
CðbÞ ¼
1 1 ð1 þ sinð3bÞÞ þ ð1 sinð3bÞÞ 2 w
ð37Þ
where w 2 ½7=9; 9=7. The original Sandler–DiMaggio model refers to w ¼ 1, such that CðbÞ ¼ 1. As reported in [3], the values of X and the plastic volumetric strain are related by
epv ¼ W½expðDXÞ 1
ð38Þ
The constants A; B; C; D; R and W are material properties, obtained from laboratory tests. Using representation in principal directions, Fig. 2 shows the aspect of the yield surface for w ¼ 7=9. Cross sections of the yield surfaces by the plane p, for w ¼ 1; w ¼ 7=9 and w ¼ 9=7 are plotted in Fig. 3 to illustrate the effect of different values of CðbÞ on their formats. In the original Sandler–DiMaggio formulation [3], the model is assumed to be associative, with W ¼ U. Furthermore, the hardening modulus and the flow direction are related by
The original Sandler–DiMaggio plastic model was proposed in [3]. It was initially conceived for granular soils, and currently it is widely applied in the oil industry to represent the behavior of rocks at depth. Let L be a parameter whose value is generally negative. The plasticity function U is piecewisely defined by a failure function pffiffiffiffi pffiffiffiffi F 1 I1 ; J 2 , and a cap elliptical function F 2 I1 ; J 2 ; L . A typical 2D profile of the corresponding yield surface is plotted in Fig. 1. Note that L specifies the distance of the origin to the center of the ellipsis. Inspired by the articles [10–12], to control the plastic surface format, so that the response of the material under triaxial compression may be different from its response in triaxial tensile, we consider a generalized Sandler–DiMaggio model, introducing an additional dependence on the Lode angle b. Precisely, the plastic function is defined by the expression
(
ð33Þ
ð32Þ Fig. 2. 3D generalized Sandler–DiMaggio yield surface with w ¼ 7=9.
5
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
H ¼ trðNÞ
ð39Þ
Using this property, the equation for the damage variable at the plastic correction step can be reformulated as
a_ ¼ c_ trðNÞ ¼ trðe_ p Þ ¼ e_ pv
ð40Þ
When using the classical Euler backward algorithm, the incremental plastic volumetric change leads to the following expression
DI 1 3K
Depv ¼
ð41Þ
which, through the relations (36) and (38), defines the evolution of Lepv Þ. This derivation allows the elimination of Dc as a variable in the plastic correction step. Fig. 4. Yield surface represented in the RHW stress space.
4.1. RHW parametrization of the yield surface In the RHW space, the generalized Sandler–DiMaggio yield surface is r1 -axially oriented, as illustrated in Fig. 4. It is formed by two pieces. There is one perfectly plastic portion S 1 , where F 1 ¼ 0, and the cap part S 2 ðLÞ, corresponding to F 2 ¼ 0. The stresses on the failure surface S 1 can be represented in terms of the variables r 1 ¼ I1 and r2 ¼ b. That is,
2
3
2
3
I1ffiffi p
3 r1 pffiffi 6 7 6 7 6 cosðbÞ 2F f ðI1 Þ 7 r ¼ 4 25 6 CðbÞ 7 4 5 pffiffi r3 2F ðI Þ sinðbÞ f 1
ð42Þ
return-mapping should be designed to consider the two cases described in the previous section. The RHW square distance from a given rtrial to a point r on the failure surface S 1 has the form
2 d rtrial ; r ¼
r
þ
Similarly, the RHW representation of any point in the cap surface S 2 ðLÞ can be parametrized in terms of two parameters r 1 ¼ h and r 2 ¼ b. Indeed,
3
r 6 7 4r 5 ¼ r 1 2 3
2
rtrial;1 pI1ffiffi3
CðbÞ
2
ð43Þ
4.2. Plastic return-mapping scheme using rotated principal stresses Since the yield surface for the Sandler–DiMaggio model is a combination of a perfect plastic part S 1 and a cap S 2 ðLÞ, the plastic
trial;3
þ
rtrial;2 cosðbÞ
pffiffi 2 2F f ðI1 Þ sinðbÞ CðbÞ
2G
pffiffi 2 2F f ðI1 Þ CðbÞ
2G :¼ d1
rtrial ; I1 ; b
ð44Þ
Therefore, the procedure to find rproj ¼ rproj ðI1 ; bÞ 2 S 1 that minimizes the distance to rtrial requires the solution of the system
2
3
LþR F f ðLÞ cosðhÞ pffiffi 3 7 6 pffiffi 7 6 6 cosðbÞ 2F f ðLÞ sinðhÞ 7 7 6 CðbÞ 5 4 pffiffi 2F ðLÞ sinðhÞ sinðbÞ fCðbÞ
3K
2
6 R1 ¼ 4
@d1 ðrtrial ;I1 ;bÞ
3
0 7 ¼ 5 @d1 ðrtrial ;I1 ;bÞ 0 @I1
ð45Þ
@b
On the other hand, the RHW square distance from rtrial to a point r 2 S 2 ðLÞ is given by pffiffi 2 2 r1;trial LþR F f pðLÞffiffi3cosðhÞ r2;trial cosðbÞ 2F fCðLÞsinðhÞ ðbÞ d2 rtrial ;h; b; L ¼ þ 3K 2G pffiffi 2 2F f ðLÞsinðhÞ r3;trial sinðbÞ CðbÞ þ :¼ d2 rtrial ; h; b; L ð46Þ 2G Furthermore, Eq. (41), rewritten as
I1
rtrial I1 ðr Þ ¼ 3K epv ðLÞ epv ðLtrial Þ
ð47Þ
is used to update simultaneously the position of the yield surface, defined by the variable L. Remembering the RHW representation (16), the invariant I1 ðr Þ, for r 2 S2 ðLÞ, has the expression
I1 ðr Þ ¼ I1 ðh; LÞ ¼ L þ RF f ðLÞ cosðhÞ
ð48Þ
As a conclusion, given Ltrial and rtrial , the procedure to compute L and rproj 2 S 2 ðLÞ that minimizes the RHW distance to rtrial requires the application of Newton’s method to find the values of h; b and L such that
3 2 3 @d2 ðrtrial ;h;b;LÞ 0 @h 7 6 7 6 7 6 R2 ¼ 6 @d2 ðrtrial ;h;b;LÞ 7 ¼ 4 0 5 5 4 @b 0 resðh; LÞ 2
ð49Þ
where Fig. 3. Yield surface cross sections for w ¼ 7=9; w ¼ 1 and w ¼ 9=7.
resðh; LÞ ¼ I1 rtrial I1 ðh; LÞ 3K epv ðLÞ epv ðLtrial Þ ¼ 0
ð50Þ
6
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
Consequently, from an algorithmic point of view, the returnmapping scheme is performed in the following steps: e Given ee trial ; Ltrial , compute the stress rtrial ¼ DRHW etrial . If I1 rtrial < Lproj , then compute L and rproj 2 S 2 ðLÞ by solving the system (49), and accept them as updated variables. Otherwise, compute rproj 2 S 1 by solving the system (45). – If I1
rproj P Lproj , then accept rproj and Lproj as updated
variables. – Else, compute L and rproj on the intersection of S 1 and S2 ðLÞ by solving a system obtained from (49) by eliminating the first equation, and fixing the third one for h ¼ p=2 (i.e., by reducing the system to two equations for the two unknowns b; L), and accept them as updated variables.
Table 1 Material parameters for the McCormic Ranch Sand [3]. Parameter
Value
E (Young’s modulus) m (Poisson’s ratio) A B
100 ksi 0.25 (dimensionless) 0.25 ksi
C D R W
w
1
0.67 ksi 0.18 ksi
1
0.67 ksi 2.5 (dimensionless) 0.066 (dimensionless) 1 (dimensionless)
A systematic way for computing the derivative of the projected stress as a function of the total strain is included in Appendices A–C. 5. Example The example reproduces the test documented in [13]. A uniaxial compressive strain is applied in the axial direction and zero strain in the other directions, followed by uniaxial strain in the opposite direction. Fig. 5 shows the comparative results between the work of [13] and the present study. The material parameters that were used in this example are listed in Table 1. The load path is also illustrated in Fig. 6, where the stress state is represented as a funcpffiffiffiffi tion of I1 and J 2 . It is graphically recognized that during the compression, the return-mapping scheme involves the projection of the stress points on the cap function defined by F 2 ¼ 0, increasing the plastic damage variable. During the unloading path, the stress point is project on the surface F 1 ¼ 0. Initially, on the path A–B, the material loses volume due the void closure and, thus, it yield triggering the cap function. Then on the path B–C, during unloading the material response is elastic until it reaches a critical point and fails in shear, triggering the failure function, that is represented by the path C–D. Fig. 7 represents the evolution of the residual normal when solving the non-linear systems of equations to project rtrial on S 1 and S 2 ðLÞ, showing that the number of iterations needed for convergence is less than 5.
Fig. 6. Load path of the experiment represented as a function of I1 and
pffiffiffiffi J2 .
Fig. 7. Residual norm per iteration for each load step.
6. Conclusion
Fig. 5. Comparison between results of uniaxial loading conduced by [13] and the results obtained in this work.
An efficient formulation/implementation is presented for any associative elastoplastic models that can be written in terms of the invariants of the stress tensor. It is based on the formulation of the plastic return-mapping in terms of the invariants of the stress tensor and of the minimization of a distance function to the admissibility surface. The improvement obtained with the formulation is that the Newton method is applied to a smaller system of equations and that part of these equations amount to the minimization of a convex function When applied to a generalized Sandler–DiMaggio model, the methodology leads to an iterative procedure which converges in less than 5 iterations.
7
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
Appendix A. Systems to be solved and square distance derivatives
@ðI1 ; bÞ @R1 ¼ J1 1 @ rtrial @ rtrial
A.1. d1 ðI1 ; bÞ derivatives
2 2 @ d1 @R1 @I @ r ¼ 4 12 1 @ d1 @ rtrial
2
d
rtrial ; r ¼
rtrial;1 pI1ffiffi3
rtrial;2
pffiffi 2 2F f ðI1 Þ cosðbÞ CðbÞ
@I1 @d1 ðrtrial ;I1 ;bÞ
2 @ r1 @I1
:¼ d1
2G
@d1 ðrtrial ;I1 ;bÞ
@b@ r1
2G
pffiffi 2 2F f ðI1 Þ CðbÞ
rtrial;3 sinðbÞ
þ
6 R1 ¼ 4
þ
3K
2
2
rtrial ; I1 ; b
ðA:1Þ
6 J1 ¼ 4
7 5
ðA:2Þ
A.2. d2
@ 2 d1 ðrtrial ;I1 ;bÞ
@I21
@I1 @b
@ 2 d1 ðrtrial ;I1 ;bÞ
@ 2 d1 ðrtrial ;I1 ;bÞ
@b@I1
@b2
r
d2 ðrtrial ;h;b;LÞ ¼
3 7 5
@ r1 @b @ r2 @b @ r3 @b
5
ðB:5Þ
3 7 7 5
ðB:6Þ
@d2 ðrtrial ;h;b;LÞ
B.2. Surface F 2 ðh; bLÞ
@ rproj @ rproj @ðh; b; LÞ ¼ @ rtrial @ðh; b; LÞ @ rtrial
ðB:7Þ
3K
r3;trial sinðbÞ
2
r
þ
pffiffi 2 2F f ðLÞsinðhÞ 2;trial cosðbÞ CðbÞ
2G
pffiffi 2 2F f ðLÞsinðhÞ CðbÞ
ðA:4Þ
2G 3
@ 2 d2 ðrtrial ;h;b;LÞ
@h2 6 6 @ 2 d r ;h;b;L ð Þ 2 trial J2 ¼ 6 6 @b@h 4
@R2 @R2 drtrial þ dðh; b; LÞ þ Oðrtrial ; h; b; LÞ ¼ 0 @ rtrial @ðh; b; LÞ
@h@b @ 2 d2 ðrtrial ;h;b;LÞ
@resðh;LÞ @h
@b2 @resðh;LÞ @b
@ 2 d2 ðrtrial ;h;b;LÞ
1 @ðh; b; LÞ @R2 @R2 ¼ @ rtrial @ðh; b; LÞ @ rtrial
ðB:10Þ
@ðh; b; LÞ @R2 ¼ J1 2 @ rtrial @ rtrial
ðB:11Þ
@h@L
7 7 7 5
@ 2 d2 ðrtrial ;h;b;LÞ 7 @b@L
@resðh;LÞ : @L
resðh; LÞ ¼ I1 rtrial I1 ðh; LÞ 3K epv ðLÞ epv ðLtrial Þ ¼ 0
@ 2 d2
6 @h@r1 6 2 ¼ 6 @ d2 4 @b@ r1 @ DL @ r1
2 @ r1
ðA:6Þ
@h
6 @ rproj ¼ 6 @ r2 @ðh; b; LÞ 4 @h @ r3 @h
ðA:7Þ
Appendix B. Derivatives of the projected stress with respect to trial stresses B.1. Surface F 1 ðI1 ; bÞ
@ 2 d2 @h@ r2
@ 2 d2 @h@ r3
@ 2 d2 @b@ r2
@ 2 d2 @b@ r3
@ DL @ r2
@ DL @ r3
@ r1 @b
@ r1 @L
@ r2 @b
@ r2 @L
@ r3 @b
@ r3 @L
3 7 7 7 5
ðB:1Þ
ðB:12Þ
3 7 7 5
ðB:13Þ
Appendix C. Tangent operator The tangent operator to be calculated is given by the derivative of the projected stress tensor by the trial deformation tensor, and, by the chain rule, it is shown that:
@ rproj @ rproj rtrial ¼ @e @ rtrial @ e
For surface F 1 ðI1 ; bÞ the chain rule gives
@ rproj @ rproj @ðI1 ; bÞ ¼ @ rtrial @ðI1 ; bÞ @ rtrial
2 @R2 @ rtrial
3
ðB:9Þ
Rearranging,
ðA:5Þ
@ 2 d2 ðrtrial ;h;b;LÞ
ðB:8Þ
and perform a Taylor expansion on it, we have,
@h 7 6 7 6 R2 ¼ 6 @d2 ðrtrial ;h;b;LÞ 7 5 4 @b resðh; LÞ
2
@ 2 d1 @b@ r3
If we take the derivative of the distance below
LþRF f ðLÞcosðhÞ pffiffi 1;trial 3
r
þ
ðA:3Þ
derivatives
2
@ 2 d1 @b@ r2
3
R2 ¼ 0
trial ; h; b; L
@ 2 d1 @I1 @ r3
For surface F 2 ðh; b; LÞ the chain rule gives
@ 2 d1 ðrtrial ;I1 ;bÞ
@ r3 @I1
@ 2 d1 @I1 @ r2
3
@b
2
6 @ r2 @ rproj ¼6 @ðI1 ; bÞ 4 @I1
ðB:4Þ
ðC:1Þ
And by straightforward derivation, the tangent operator can be calculated by the following formula:
If we take the derivative of the residual equation R1 ¼ 0, and perform a Taylor expansion on it, we have,
X
@ rproj @ rproj;i ¼ ðv i v i Þ ðv j v j Þ C @e @ rtrial;j
@R1 @R1 drtrial þ dðI1 ; bÞ þ Oðrtrial ; I1 ; bÞ ¼ 0 @ rtrial @ðI1 ; bÞ
For the flowing calculations it is convenient to deduce a formulation for the tangent operator in terms of a matrix. To achieve this, first we multiply the tangent operator by a deformation de:
ðB:2Þ
Rearranging,
1 @ðI1 ; bÞ @R1 @R1 ¼ @ rtrial @ðI1 ; bÞ @ rtrial
ðB:3Þ
@ rproj de ¼ @e
X
ðC:2Þ
@ rproj;i ðv i v i Þ ðv j v j Þ 2lI4 þ kI2 I2 de @ rtrial;j ðC:3Þ
8
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
And with some mathematical manipulation Eq. (C.4) is obtained:
@ rproj de ¼ @e
X @ rproj;i h @ rtrial;j
i 2l Tr½ðv j v j Þde þ k I2 de ðv i v i Þ ðC:4Þ
In Voigt notation, the stress and deformation tensors are represented by six element vectors as shown below:
rxx 1 C 1 B B rxy C rxz B C B rxz C B C ryz C ¼ A B C B ryy C rzz B C @ ryz A rzz
1
0
1 0 0 B C dexx ¼ @ 0 0 0 A 0
ðC:5Þ
deyy
0 B ¼ @0 0
1
0
1
ðC:6Þ In Eq. (C.4), de is replaced by each of these matrices. The resultant matrix of each of these six operations represents the rate of variation of each of the stress tensor components in relationship with one of the deformations. Each column of the tangent matrix is composed by one of these operations. Next it is shown the detailed calculation of each of these columns:
@ rproj Column 1 : dexx @e
!
X @ rproj;i h ij
@ rtrial;j
2l
v 2j
Column 2 :
@ rproj dexy @e
X @ rproj;i ! 2 2l @ rtrial;j ij
Column 3 :
@ rproj dexz @e
!
Column 4 :
@ rproj deyy @e
Column 5 :
@ rproj deyz @e
@ rproj Column 6 : dezz @e in winch, k ¼ 1; 2; 3.
vj
k
1
C hx A 0
ðC:8Þ
And given the deformation tensor in expressed in principal directions:
1
e1 0 0 C e¼B @ 0 e2 0 A 0 0 e3
ðC:9Þ
1
e1 ðe1 e2 Þhz ðe1 e3 Þhy B C e2 ðe2 e3 Þhx A RT eR ¼ @ ðe1 e2 Þhz ðe1 e3 Þhy ðe2 e3 Þhx e3
ðC:10Þ
Thus for any deformation, in case of terms out of the principal axes, we have:
0 1 0 0 0 1 B C B C dexy ¼ @ 1 0 0 A dexz ¼ @ 0 0 0 A 0 0 0 0 0 1 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 C B C B C 1 0 A deyz ¼ @ 0 0 1 A dezz ¼ @ 0 0 0 A 0 0 0 1 0 0 0 1
0
0 hx
hy
0
Using this notations, the tangent matrix is a second order tensor of size 6 6. To calculated this matrix, first it is defined six tensor unitary directions:
0
B R ¼ I þ @ hz hy
hz
The rotation of this tensor, already elimination the order terms is given by:
0
rxx rxy r r¼B @ xy ryy rxz ryz
0
0
C.1. Tangent matrix in Voigt notation
0
0
1
þk
i
vi vi
h¼
X @ rproj;i 2 2l v j3 v j1 v i v i @ rtrial;j ij i X @ rproj;i h ! 2l v 2j2 þ k v i v i @ r trial;j ij
ðC:13Þ
In the same way as before, it will be calculated the rate of variation of each of the project stress tensor components in relation with a unitary strain direction. Each column of the correction of the tangent matrix is composed by one of these operations. Next it is shown the detailed calculation of each of these columns:
Column 1 : dexx ! Column 2 : dexy
ðC:7Þ Column 6 : dezz
X rproj;i rproj;j
v v
ðv v þ v v Þ
j 1 i1 i j j i ei ej ij X rproj;i rproj;j ! v j1 v i2 þ v j2 v i1 ðv i v j þ v j v i Þ e e i j ij X rproj;i rproj;j ! v j1 v i3 þ v j3 v i1 ðv i v j þ v j v i Þ e e i j ij X rproj;i rproj;j ! v j2 v i2 ðv i v j þ v j v i Þ ei ej ij X rproj;i rproj;j ! v j2 v i3 þ v j3 v i2 ðv i v j þ v j v i Þ e e i j ij X rproj;i rproj;j ! v j3 v i3 ðv i v j þ v j v i Þ e e i j ij
ðC:14Þ
!
is the k component of the eigenvector
rproj , the correction
rproj;i rproj;j v j ðde v i Þ ðv i v j þ v j v i Þ ei ej
Column 5 : deyz
X @ rproj;i 2 2l v j3 v j2 v i v i @ r trial;j ij i X @ rproj;i h ! 2l v 2j3 þ k v i v i @ rtrial;j ij
ðC:12Þ
And as v i and v j are eigenvector of both e and em this position is given by:
Column 4 : deyy
1
ðC:11Þ
ðei ej Þ
Column 3 : dexz
vj vj vi vi 2
v j ðde v i Þ ¼ hðei ej Þ v j ðde v i Þ
for i ¼ 1; 2 , j ¼ 2; 3 and i – j. C.3. Particular case of repeated eigenvalues
vj,
with
C.2. Correction of the eigenvectors variation In the deduction presented, it was not considered the variation of the eigenvectors between rproj and rtrial . This variation can be represented by a rotation, and its inclusion in the tangent matrix is demonstrated here. A rotation tensor can be seen as the identity tensor plus a part relative to the rotation in each of the axes as below:
In case of repeated eigenvalues, the correction becomes a relation 00 in the flowing equation present in the calculation of each of the columns;
rproj;i rproj;j ei ej
ðC:15Þ
The resolution of this problem will be exemplified for case i ¼ 1; j ¼ 2, with rproj;1 ¼ rproj;2 and e1 ¼ e2 . First, it can be proved that:
rtrial;1 rtrial;2 ¼ 2l e1 e2
D.L. Cecílio et al. / Computers and Geotechnics 64 (2015) 1–9
Which is applied in Eq. (C.15).
References
rproj;1 rproj;2 2l rtrial;1 rtrial;2
ðC:16Þ
And so, the flowing limit is calculated with dr ! 0, and imposing a variation in the direction fdr; dr; 0g:
pr
@ r2 rproj;1 þ @@rrproj;1 fdr; dr; 0g rpr 2 þ @ rtrial fdr; dr; 0g trial lim 2l dr!0 ðrtrial;1 þ drÞ ðrtrial;2 drÞ
ðC:17Þ As
rproj;i ¼ rproj;j and rtrial;i ¼ rtrial;j :
lim 2l
@ rproj;1 @ rtrial;1
@r
@r
@r
dr @ rproj;1 dr @rproj;2 dr þ @ rproj;2 dr trial;2
trial;1
trial;2
ðC:18Þ
2dr
dr!0
And cancelling the terms dr, it is defined the final correction for the case i ¼ 1; j ¼ 2, in which rproj;1 ¼ rproj;2 and e1 ¼ e2 .
l
@ rproj;1 @ rproj;1 @ rproj;2 @ rproj;2 þ @ rtrial;1 @ rtrial;2 @ rtrial;1 @ rtrial;2
ðC:19Þ
Finally, the flowing substitution must be made in the calculation of the correction for each of the columns in case of repeated eigenvalues:
9
rproj;i rproj;j @ rproj;i @ rproj;i @ rproj;j @ rproj;j ¼l þ ei ej @ rtrial;i @ rtrial;j @ rtrial;i @ rtrial;j
ðC:20Þ
[1] Souza Neto E, Peric D, Owen DRJ. Computational Methods for Plasticity. John Wiley & Sons Ltd; 2008. [2] Chen WF, Han DJ. Plasticity for Structural Engineers. Springer-Verlag; 1988. [3] DiMaggio FL, Sandler I. Material model for granular soils. Eng Mech Div 1971;97(3):935–50. [4] Devloo PRB. PZ: an object oriented environment for scientific programming. Comput Methods Appl Mech Eng 1997;150(1-4):133–53. http://dx.doi.org/ 10.1016/S0045-7825(97)00097-2. [5] Devloo PRB. Object oriented tools for scientific computing. Eng Comput 2000;46(2):203–14. http://dx.doi.org/10.1016/0045-7949(93)90185-G. [6] Lainé E, Vallée C, Fortuné D. Nonlinear isotropic constitutive laws: choice of the three invariants, convex potentials and constitutive inequalities. Int J Eng Sci 1999;37(15):1927–41. http://dx.doi.org/10.1016/S0020-7225(99)00006-3. [7] Borja RI. Plasticity Modeling & Computation. Springer Science & Business Media; 2012. [8] Armero F, Pérez-Foguet A. On the formulation of closest-point projection algorithms in elastoplasticity – Part I: The variational structure. Int J Numer Methods Eng 2002;53(2):297–329. http://dx.doi.org/10.1002/nme.278. [9] Armero F, Pérez-Foguet A. On the formulation of closest-point projection algorithms in elastoplasticity – Part II: Globally convergent schemes. Int J Numer Methods Eng 2002;53(2):331–74. http://dx.doi.org/10.1002/nme.279. [10] Foster CD, Regueiro RA, Fossum AF, Borja RI. Implicit numerical integration of a three-invariant, isotropic/kinematic hardening cap plasticity model for geomaterials. Comput Methods Appl Mech Eng 2005;194:5109–38. http:// dx.doi.org/10.1016/j.cma.2005.01.001. [11] A.F. Fossum, J.T. Fredrich, Cap plasticity models and compactive and dilatant pre-failure deformation, in: 4th North American Rock Mechanics Symposium, Seatle, United States of America, 2000. [12] A.F. Fossum, R.M. Brannon, The Sandia Geomodel: Theory and Users Guide, SAND Report, Sandia National Laboratories, 2004. [13] Sandler IS, Rubin D. An algorithm and a modular subroutine for the cap model. Int J Numer Anal Methods Geomech 1979;3(2):173–86. http://dx.doi.org/ 10.1002/nag.1610030206.