Plasticity and strain localization around a horizontal wellbore drilled through a porous rock formation

International Journal of Plasticity 78 (2016) 114e144

Contents lists available at ScienceDirect

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

Plasticity and strain localization around a horizontal wellbore drilled through a porous rock formation
Spiezia*, Valentina A. Salomoni, Carmelo E. Majorana Nicolo Department ICEA, University of Padova, via Marzolo 9, 35131 Padova, Italy

a r t i c l e i n f o

a b s t r a c t

Article history: Received 16 July 2015 Received in revised form 17 October 2015 Available online 1 December 2015

Predicting plastic deformation and localization band around a wellbore drilled through a porous rock is a challenging task that could have important implications for the prediction of instability and sand production. This study focuses on analyzing stress field, plastic zones, and localized deformations around a horizontal borehole drilled at a great depth through a highly porous rock formation. Several laboratory studies demonstrate that, depending on the loading path, highly porous rocks are susceptible to different failure mechanisms, but most of these mechanisms are mainly due to shear-induced dilation and shear-enhanced compaction. Plasticity models, in conjunction with bifurcation analysis, represent an extremely useful framework for describing such detailed constitutive responses. This paper presents a new elastoplastic constitutive model characterized by two yield surfaces intersecting smoothly, that is able to capture the different failure modes. The model is validated against experimental data for several different porous rocks, and it is then used to determine the stress and strain distributions around a horizontal wellbore using nonlinear finite element analysis. Particular interest is devoted to predicting the condition for the formation of a localized band of intense deformation, elucidating the factors that either prevent or enhance the band initiation. Results of simulations show the key role played by the elastoplastic constitutive model and the effects of the mud pressure, the in-situ stress condition and geometric imperfections in the development and propagation of plastic zone, as well as in the initiation of localization zone. © 2015 Elsevier Ltd. All rights reserved.

Keywords: B. Elastic-plastic material B. Rock C. Finite elements B. Finite strain Wellbore stability

1. Introduction Wellbore instability in the oil and gas industry applications, in particular in deep perforation, continues to be one of the major problems faced by scientists and engineers, that can dramatically increase production costs (Adachi et al., 2013; Udegbunam et al., 2013; Cao et al., 2014; Parsamehr et al., 2015; Zhang et al., 2015). Drilled holes alter significantly the insitu stress distribution, leading to a new stress concentration in the rock formation. Depending on environmental factors, this new stress distribution can lead to compressive failures of the rock, known as stress-induced breakouts, and/or tensile failure, known as drilling-induced tensile fractures. Even though a certain amount of failure is tolerated around the wall, excessive breakouts could lead to catastrophic instabilities, resulting in well loss (Zoback, 2010). To prevent instabilities during perforation, wellbores are temporarily supported by the drilling mud pressure. If instability occurs, the value of the mud pressure needs to be sufficiently high to prevent compressional failure, but it should also be lower than a critical value that would cause tensile failure and unintentional hydraulic fracturing. * Corresponding author. E-mail address: [email protected] (N. Spiezia). http://dx.doi.org/10.1016/j.ijplas.2015.10.013 0749-6419/© 2015 Elsevier Ltd. All rights reserved.

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

115

Preventing wellbore instability is even more important in horizontal and inclined wells, where drilling and mud circulation operations are usually more challenging. Even if a horizontal wellbore is sufficiently stable, it is important to prevent breakouts and sand production. Resultant materials are generally more difficult to remove in a horizontal well while drilling and could slow down the rate of penetration. Traditionally, the stability of a well is determined using models based on linear elasticity, where failure is assumed to occur when the stress along the wall of the hole reaches the failure or tensile strength of the rock (Zoback, 2010). Elastoplastic models can provide a more realistic evaluation of stability, with the advantage of being able to delineate the extent of the damaged region (Roshan and Rahman, 2011; Chen et al., 2012, 2014; Hamid et al., 2014; Shojaei et al., 2014). Clearly, the assessment of stability depends on the capability of the constitutive model to capture the different failure modes. Usually, numerical simulations employ elastoplastic models such as the well-known MohreCoulomb or DruckerePrager (Rahimi and Nygaard, 2014; Elyasi and Goshtasbi, 2015; Zhang et al., 2015). These models, although simple to use and to calibrate, are characterized solely by a shear yield surface and thus they can describe only the dilatant plastic process preceding the brittle failure, the so called ‘shear-induced dilation’ mechanism. On the other hand, these models are not able to describe the compactant plastic mechanism that leads to a ductile failure, usually called ‘shear-enhanced compaction’, that strongly influences the collapse, especially in high porosity rocks. Several experiments conducted on high porosity sandstones in the last two decades show that compaction failure can take place under certain stress conditions, in contrast with the more common dilatant failure (Baud et al., 2004, 2009; Xie and Shao, 2012, 2014; Zhou et al., 2013). Compactant failure typically occurs in porous rocks under relatively high confining pressure, with a failure mode conventionally described as homogeneous cataclastic flow. Experimental studies on boreholes drilled in a cubical specimen of rock with different porosity show that, if the porosity is high (22e25%), the failure is associated with the compactant mechanism. In certain situation the specimen develops a long and thin fracture, originating in the region of highest compressive stress concentrations (Haimson, 2001; Haimson and Kovachic, 2003; Buscarnera et al., 2014). Since the compactant plastic mechanism can have an important role in the analysis of wellbore, an appropriate constitutive elastoplastic model is necessary. This model must capture both the dilatant and compactant behavior, and the transition between these two failure criteria (Wong and Baud, 2012; Baud et al., 2015). The scope of this study is twofold. The first is to derive a new elastoplastic constitutive model for high porosity rocks. The outstanding aspect is that the model has been derived in order to describe the main mechanical characteristics, but taking also into account the aspects related to parameter identification and numerical implementation, in the attempt to get the best compromise between these different instances. The second is to evaluate how the plastic zones develop around a horizontal wellbore and which are the conditions for the inception of localizations band. The article is organized in the following sections. Section two briefly reviews the basic theoretical framework for the solution of a nonlinear boundary-value problem with the finite element method, highlighting the conditions for instability according to the bifurcation theory. Section three describes the main aspects of the elastoplastic behavior of high porous rock, which are then used to develop the equations of the constitutive model. The subsequent section explains how to identify the parameters of the constitutive model and demonstrates the capability to reproduce laboratory tests. The final section discusses the results obtained from the numerical simulations of a horizontal wellbore in terms of stress and strain distributions. In the numerical analysis, particular interest is devoted to predicting the conditions for the formation of a localized band of intense deformation, discussing the factors that either enhance or prevent its initiation. The following rules apply for the notation throughout the paper: bold-face letters denote matrices and vectors; the symbol ‘,’ denotes an inner product of two vectors (e.g. a,b ¼ aibi), or a single contraction of adjacent indices of two tensors (e.g. c,d ¼ cijdjk); the symbol ‘:’ denotes an inner product of two second-order tensors (e.g. c:d ¼ cijdij), or a double contraction of adjacent indices of tensor of rank two and higher (e.g. C : εe ¼ Cijkl εekl ); the symbol ‘5’ denotes a juxtaposition, e.g. (a 5 b)ij ¼ aibj. For any symmetric second-order tensor a and b we have (a 5 b)ijkl ¼ aijbkl; (a 4 b)ijkl ¼ ajlbik; and (a . b)ijkl ¼ ailbjk. 2. Basic theoretical framework This part of the work briefly recalls the theoretical background which allows to solve numerically an elastoplastic boundary-value problem, including the condition that leads to instability, with consequent band initiation. Due to the intrinsic structure characterized by a large number of voids, highly porous rocks may experience large deformations under the loading process. Furthermore, former research in the field demonstrates that hyperelastic formulation should be used for problems where shear banding and localization may be encountered (McAuliffe and Waisman, 2015) and that finite deformation effects enhance strain localization (Borja, 2002). Therefore, the complete problem is formulated in the context of the finite strain continuum mechanics, with an hyperelasto-plastic constitutive model. For the sake of briefness, only the fundamental aspect will be recalled, and the reader is referred to Borja (2013) for a detailed description of the subject. Let f : B 0 /B be the function that maps the body from the initial configuration B 0 to the deformed configuration B , and let F be the deformation gradient for any material point X in B 0 , expressed in coordinate form

F aA ¼

vfa ðXÞ vxa ¼ : vXA vXA

(1)

116

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

The deformation of each point of the body can be described by the symmetric left CauchyeGreen tensor b ¼ F,FT. In addition, the stress status can be described by the symmetric Kirchhoff stress t, related to the Cauchy tensor s via the relationship

t ¼ Js;

(2)

where J ¼ det(F) is the Jacobian. In a usual quasi-static problem, the equilibrium condition of the body is ensured by the variational form of the linear momentum balance equation, written with respect to the initial undeformed configuration as

Z

Z ðgrad h : t  r0 h,GÞdV 

B0

h,tdA ¼ 0;

(3)

vB t0

where r0G is the reference body force vector, t is the nominal traction vector in vB t0 3vB 0 , n is the unit vector on vB t0 and h is the weighting function. Several formulations have been proposed in the literature to describe the elastoplastic behavior in finite strains. The model here presented (see Section 3) adopts a formulation based on the multiplicative decomposition of the deformation gradient, which gives rise to the so-called multiplicative plasticity theory (Lee, 1969; Simo, 1992; Borja et al., 2003). Unlike the infinitesimal theory where the deformation tensor is decomposed additively, the corner stone of the finite deformation multiplicative plasticity is the multiplicative decomposition of the deformation gradient

F ¼ F e $F p ;

(4)

where Fe and Fp are, respectively, the elastic and plastic part of the deformation gradient. Ignoring non-mechanical power and kinetic energy production, balance of energy and the use of the second law of thermodynamics lead to the following reduced dissipation inequality

 t2

vJ e $b vbe

 :dþ2

vJ e $b : vbe



1  L v be $be1 2

  0;

(5)

_ 1 is the where be ¼ Fe,FeT is the elastic left CauchyeGreen tensor, d ¼ symm(l) is the rate of deformation tensor, l ¼ F,F overall spatial velocity gradient and J is the stored energy function. The first term of Eq. (5) yields the constitutive relationship

t¼2

vJ e $b ; vbe

(6)

while the second term, after imposing that F ¼ 0 and that the postulate of maximum dissipation holds, yields

1 vG e $b  L v be ¼ g_ 2 vt

(7)

where g_ is the plastic multiplier, F ¼ 0 is the yield function and G ¼ 0 is the plastic potential. The term g_ and F satisfy the _ ¼ 0). By assuming isotropy (Simo, 1992), the elastic left CauchyeGreen tensor be KuhneTucker conditions (g_  0, F  0, gF and the Kirchhoff stress tensor t can be spectrally decomposed as

be ¼

3   X 2 leA nðAÞ 5nðAÞ ; A¼1

t¼

3 X

bA nðAÞ 5nðAÞ ;

(8)

A¼1

where ðleA Þ2 and bA are respectively the eigenvalues of the tensors be and t. With the assumption of frame indifference, the free energy function depends solely on the material point X and the logarithmic elastic principal stretches, defined as

  εeA ¼ ln leA ;

A ¼ 1; 2; 3:

(9)

Therefore, JðX; be Þ ¼ jðX; εeA Þ and the principal Kirchhoff stress b can be derived by means of Eq. (6) along each direction

    vj X; εeA bA X; εeA ¼ ; vεeA

A ¼ 1; 2; 3:

(10)

The rate form of the constitutive equation must be integrated so as to quantify the local stresses and deformations at discrete time stations. The aim is to determine the statically admissible value of be, given the relative deformation gradient f ¼ vx/vxn at a specific time step. The relevant rate equation to integrate is the time derivative of be, which is written in the predictor-corrector format as follows

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

117

e b_ ¼ l$be þ be $lT  ðL v be Þ : |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl}

(11)

corrector

predictor

The algorithm consists of two steps. The first step is the incremental counterpart of the predictor step, computed by freezing the plastic flow

beTr ¼ f $ben $f T ;

(12)

ben

where is the elastic left CauchyeGreen tensor and the relative deformation gradient with respect to the deformed configuration at time tn. The second step is the incremental plastic corrector in presence of plastic flow, which is obtained by integrating the Lie derivative in Eq. (7). In spectral form and in the space of the principal elastic logarithmic stretches, the discretized flow rule reads

εeA ¼ εeTr A  Dg

vG : vbA

(13)

eTr e The term εeTr A is the Trial value of εA , given by the spectral decomposition of b , and Dg is the discretized plastic multiplier _ Due to the nonlinear form of the balance law, arising from both the material and geometric assumptions, Eq. (3) needs to be g. linearized. The most important term arising from the linearization is the expression of the tangent stiffness matrix, which reads

Z K¼

gradh : a : graddudV

(14)

B0

where

a ¼ a  t.1

(15)

and (t . 1)ijkl ¼ tildjk. It is important to underline that, due to the assumption of finite deformation, the tangent stiffness depends also from the current stress status, expressed through the tensor t. The algorithmic tangent tensor a, required for the numerical solution of the linear momentum balance equation, is obtained using the spectral decomposition previously presented and, in agreement with Ogden (1984), Borja et al. (2003), reads

a¼

3 X 3 X

aAB mðAAÞ 5mðBBÞ þ

A¼1 B¼1

3 X X A¼1 BsA

bB  bA eTr lB

eTr  lA

! 

eTr

eTr

lB mðABÞ 5mðABÞ þ lA mðABÞ 5mðBAÞ



(16)

eTr

where aAB ¼ vbA =vεeB , 2εeB ¼ lnðlB Þ and m(AB) ¼ n(A) 5 nB. The tangent operator a plays a central role not only for the solution of the boundary-value problem, but also to detect the onset of instability in the material response. In fact, deformation bands are generally regarded as a result of material instability, involving the bifurcation theory from continuum mechanics (Borja, 2002; Borja and Aydin, 2004; Borja, 2004; Nicot et al., 2010). This theory identifies a stress configuration at which the solution could be non-unique. The localization condition can be written with respect to the current configuration as

b

h0 a$v ¼ 0; h2

Tr

b

Tr

aij ¼ nk aikjl nl :

(17)

where h0 and h are the band thickness respectively in the initial and in the current configuration, a takes the form of Eq. (15), v is the velocity jump and a is the Eulerian acoustic tensor. Non trivial solutions are then possible when det(a) ¼ 0 for some critical band orientation n in the current configuration. Defining the localization function as b

F ¼ inf ðdet aÞ;

(18)

b

n

the inception of a deformation band can be inferred from the initial vanishing of F . 3. Constitutive model Highly porous rocks are characterized by a complex elastoplastic behavior, which is worth being investigated accurately, in order to derive a reliable constitutive model. Elastoplastic models may be derived following a phenomenological or a

118

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

micromechanical approach. The phenomenological approach is based on the macroscopic behavior of the material observed on the basis of available experimental data, and do not require detailed information about the microstructure and its evolution at appropriate scales. On the other hand, micro-mechanical models can properly consider some important additional features of material behavior, for example, impact of porosity, inherent damage-friction coupling, unilateral effects and related continuity at the openingeclosing transition of cracks, volumetric dilatation induced by frictional sliding, hysteretic loop related to local hardening kinetics (Zhu et al., 2010; Shen et al., 2012; Lin et al., 2012; Nicot et al., 2014; Shen and Shao, 2015), but they require a more sophisticated formulation. The paper aims to develop an efficient constitutive model for the study of wellbores, where the macroscopic mechanical behavior plays a predominant role. Hence in this research a phenomenological approach is used. First, this section briefly summarizes the mechanical response revealed by experimental studies conducted on porous rocks. Then, the laboratory observations will form the basis to develop the equations of the constitutive model. Finally, this section describes the numerical implementation of the proposed model. The model is meant to take into account all the fundamental phenomenological aspects, and takes particular care of choosing the most efficient solutions from a numerical point of view. 3.1. Behavior of high porosity rocks In order to develop and evaluate an appropriate constitutive model for highly porous rocks, it is important to understand the main aspects of the complex macro/microscopic behavior under different loading conditions. Under hydrostatic compression, several studies conducted on sandstones of varying porosity (Baud et al., 2006; Zhang et al., 2013; Xie and Shao, 2014) show a mechanical response similar to that reported in Fig. 1. Typically, the response of a porous sandstone is nonlinear, even in the elastic regime, especially for low pressure values. Initial compaction is associated with the rearrangement of cement and loosely bonded grains. Then it is associated with grain crushing and the curve tends to be approximately linear. As the pore space is tightened by elastic deformation, the sample becomes progressively stiffer, as evidenced by a decrease in compressibility. As the hydrostatic loading increases, the sample reaches a point where it suddenly becomes more compliant, showing a dramatic increase in compaction. After this point, corresponding to the grain crushing pressure C*, the rock experiences inelastic deformations, due to grain crashing and pore collapse. Finally, after a considerable amount of porosity has been crushed out, the rock begins to harden. This inflection point occurs at a wide range of effective pressures, depending mainly on the porosity and on the average grain size. If a sample is loaded beyond this inflection point and then unloaded, the permanent compaction is significant, confirming the occurrence of inelastic deformations. Under axisymmetric compression (Wong et al., 1990, 1997) high porosity sandstones have two mainly different responses. When the confining pressure is relatively low, shear induced dilation is observed after volume compaction has occurred. The specimen ultimately fails via shear localization. This behavior is common for most kinds of geomaterials, and suggests the use of a shear yield surface to model them, like the well-known MohreCoulomb or DruckerePrager yield surface. However, when the confining pressure is relatively high, the response is similar to that of the hydrostatic loading previously described. The grain crushing begins at a lower mean stress and a larger porosity reduction occurs for a small increase in mean stress. This response, usually called ‘shear enhanced compaction’ (Wong et al., 1990), reports a larger reduction in grain crushing pressure for smaller confining pressures. To capture this behavior, the shear yield surface needs to be delimited by a so-called cap surface, to describe the compactive plastic response. Increasing the compressive pressure, samples show an hardening

Fig. 1. Hydrostatic pressure versus porosity reduction for high porosity rock. C* indicates the yielding pressure.

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

119

behavior, which is more significant at high values of the confining pressure, thus suggesting that the cap yield surface may expand. In conclusion, a loading path which begins at a lower confining pressure leads to the shear yield surface, while a loading path beginning at a higher confining pressure intersects the cap yield surface. These two different situations are shown schematically in Fig. 2. At intermediate confining pressures, a transitional regime exists where both yield surfaces (and consequently both failure mechanisms) are active (Wong and Baud, 2012; Weng and Ling, 2013). 3.2. Equations of the constitutive model 3.2.1. Hyperelastic response As reported briefly in the previous section, porous rocks have a non linear response, even if elastic, especially for low level of the mean normal stress. The sample becomes stiffer as far as the hydrostatic pressure increases, suggesting a dependence of the bulk modulus from the applied pressure. The model employs an hyperelastic law which was introduced originally by Houlsby (1985), and then extensively used by several other authors (Borja and Tamagnini, 1998; Callari et al., 1998; Houlsby et al., 2005; Yamakawa et al., 2010; Borja et al., 2013) to simulate the non-linear elastic behavior of mainly clay and sand, and here extended for the simulation of porous rocks. The essential feature of this elastic law is a variable elastic bulk and shear modulus, which both depend on the volumetric deformations. The stressestrain relation is given in terms of invariants and can be derived from the assumption of a free energy function. Consider the volumetric and deviatoric elastic invariants defined as

εev ¼ εe ,d;

εes ¼

rffiffiffi 2 e jje jj; 3

1 ee ¼ εe  εev d 3

(19)

where εe ¼ ðεe1 ; εe2 ; εe3 Þ is the vector of the elastic logarithmic principal stretches defined in Eq. (9) and d ¼ (1,1,1). Next, consider a class of two-invariants free energy functions of the form

    b εev ; εes ¼ j e εev þ 3me εe2 j ; 2 s

(20)

  e εe ¼ P0 kexpðUÞ; j v

(21)

where

  U ¼  εev  εev0 k:

e ev Þ represents the stored energy function for isotropic loading. The parameters involved are the elastic The term jðε compressibility modulus k and the elastic volumetric strain εev0 at a mean normal effective stress of P0. The term

ae me ¼ m0 þ j k

(22)

represents the elastic shear modulus and is computed as the sum of a constant term m0 and a variable term, that depends on the elastic volumetric strain through the constant coefficient a. If a ¼ 0, then the elasticity model is defined by a variable elastic bulk modulus and a constant elastic shear modulus. Consider now the vector of the principal Kirchhoff stress b ¼ (b1,b2,b3) and define P and Q, respectively the mean normal and deviatoric Kirchhoff stress invariants, given by

Fig. 2. Two-yield surfaces model for high porosity rocks. The loading path A (low confining pressure) intersects the shear yield surface, associated with dilatant behavior. The loading path B (high confining pressure) intersects the cap yield surface, associated with compactant behavior. The axes represent the mean normal stress p and the norm of the second invariant of the deviatoric stress q.

120

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

P¼

1 b$d; 3

Q¼

rffiffiffi 3 jjsjj; 2

s ¼ b  Pd:

(23)

The relation between stress and strain invariants derives directly from the derivation of the stored energy function with respect to the deformation, and reads

P¼

b vj 3a e 2 ; ¼ P expðUÞ 1 þ ε 0 vεev 2k s

(24a)

Q¼

b vj ¼ 3ðm0  aP0 expUÞεes : vεes

(24b)

3.2.2. Yield surface Depending on the stress condition, porous rocks show two different plastic responses, mainly characterized by shear dilation or compaction. To describe this behavior there are two approaches (Issen, 2002): (i) a single smooth yield surface or (ii) a multi yield surface, where two (or more) functions intersect. As regard single surface models, a pioneering contribution was given by Carroll (1971), who proposed a type of critical state plasticity model using a parabolic yield surface; Grueschow and Rudnicki (2005) suggested to use an elliptical yield surface, with varying axis, to describe compaction. These models give good results with respect to the compaction side, but are usually coarse in reference to the dilatant side, due to the symmetry of the surface. Non-symmetric yield functions (e.g. Jefferies, 1997; Nova et al., 2003; Aubertin and Li, 2004; Marinelli and Buscarnera, 2015; Nguyen et al., 2014; Le et al., 2015), require a more complex formulation, involving more parameters, sometimes arduous to be defined from laboratory test. On the other hand, the plastic surface can be described by using two independent continuously differentiable yield functions (Issen and Challa, 2005). Usually this kind of models are based on the combination of a MohreCoulomb or DruckerePrager yield surface in the zone of shear deformation and a cap surface in the zone of compaction. One of the first and most popular model was proposed by Di Maggio and Sandler (1971), and it consists of an elliptical cap that intersects at the point of horizontal tangency the shear-failure surface. The lack of a smooth transition between the two surfaces, if on one hand it simplifies the formulation, on the other hand gives less accurate results (Fossum and Fredrich, 2000). In fact, the requirement that the cap surface intersects the shear surface at the point of horizontal tangency prevents pre-failure dilatant deformation, contradicting experimental observations. Furthermore, from a numerical point of view, two downsides arise: first, the apex produces an accumulation point of the Return Mapping algorithm, that can give inaccurate results. Second, at the intersection point the flow direction is not well defined, and must me computed as combination of vectors. To circumvent these shortcomings, the two surfaces may be connected by an additional transition function (ABAQUS, 2013; Shin and Kim, 2015), but more parameters need to be introduced, however. This study proposes a yield surface which is characterized by two yield functions that intersect smoothly, without the introduction of a transition function, thus resulting in a continuous surface. In doing so, a combined model based on a linear and an elliptical yield surface is adopted. The linear part of the yield surface F 1 corresponds to the shear failure zone and it is based on the DruckerePrager formulation (Drucker and Prager, 1952)

F 1 ðP; Q Þ ¼ Q  mP  c0 ¼ 0;

(25)

where m is the slope of the linear yield surface and c0 is the intersection of the surface with the vertical axis in a (P,Q) plot. For the compaction failure zone, an elliptical yield surface F 2 is adopted, resembling the Di Maggio and Sandler (1971) failure envelope

F 2 ðP; Q ; Pi Þ ¼ B2 ðP  Pi Þ2 þ A2 Q 2  A2 B2 ¼ 0;

(26)

where A and B are respectively the minor and major semiaxis of the ellipse and Pi is the centroid. The intersection between the two functions is defined as the point in which F 1 is tangent to F 2 , ensuring that the two functions produce a unique surface without any angular point, as shown in Fig. 3. The point (P*, Q*) where the two functions are tangent is given by

mA ffi: P  ¼ Pi  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 B A2 þ m2

(27)

If P > P*, the dilatant plastic mechanism is activated, otherwise if P < P* the compactant plastic mechanism is activated. 3.2.3. Hardening law Experimental evidence shows that, after a significant amount of porosity has been crushed, the rock begins to harden, becoming almost incompressible. In fact, considering the stressestrain relation in Fig. 1 for a hydrostatic compression test, the

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

121

Fig. 3. The two-yield surfaces model. The elliptical yield surface is tangent to the linear surface in the point defined by (P*, Q*).

curve tends to assume an asymptotic vertical trend. In the proposed model, the elliptical surface can contract and expand according to the following criteria, as shown in Fig. 4: 1. The width of the ellipse remain constant, i.e. A ¼ A0 where A0 is the initial horizontal semi-axis of the ellipse (Carroll, 1971). 2. The elliptical surface can slide along the horizontal axis, remaining always tangent to the linear surface, hence

B¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m2 P 2i  A2 m2 þ 2mc0 Pi þ c20 :

(28)

This condition guarantees that the yield surface is always continuous. 3. The variation of Pi, i.e. the variation of the centroid of the ellipse, depends on the accumulated plastic volumetric strain according to

 Pi ¼ Pi0

ε ε  εpv

r ;

(29)

Fig. 4. The two-yield surfaces model. The elliptical surface expands remaining always tangent to the linear surface. The complete yield surface is subdivided in three regions, depending on the yield and the potential functions.

122

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

where Pi0 is the initial value, ε* is the volumetric deformation at ultimate compaction and r is the exponent that controls the rate of volumetric hardening (Stefanov et al., 2011). This proposed hardening law is in agreement with the experimental observations: as soon as the plastic volumetric deformation achieves a maximum value (corresponding to the collapse of the porous skeleton) the further plastic deformation can develop only in shear mode. Even though some hardening is observed in porous rock associated with the dilatant surface, this is significantly smaller if compared with the hardening along the compactant side (Wong et al., 1990). Therefore, the model does not take into account hardening along the linear side, and the linear yield surface is fix throughout the deformation process. 3.2.4. Flow rule The model used in this study is characterized for the shear plastic mechanism by a potential function of the form

G ¼ Q  mP  c ¼ 0:

(30)

The non-associative law is obtained by adopting, as the flow potential, a DruckerePrager yield function with the frictional angle replaced by a smaller dilatancy angle, hence assuming that m < m, to avoid the often excessive dilatancy predicted by the associative rule (Borja, 2013). For the cap side of the plastic surface, the flow rule is associative, using the elliptical function as plastic potential. The smooth continuity between the two potential functions is guaranteed defining a transition zone. Let's introduce the e Þ in which the plastic potential G and the yield function F have the same normal, and hence the same flow vector. e Q point ðP; 2 Imposing the tangent condition between the elliptical function and the potential function we obtain

mA e ¼ P  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi: P i  2 B A2 þ m2

(31)

e > P > P , and in which the plastic The transition zone is defined as the portion of the elliptical yield surface bounded by P potential is defined by G. This formulation ensures that along the complete yield surface there are no discontinuities with respect to the flow vector. Summarizing, the complete plastic surface is subdivided in three regions, depending on the yield and potential functions, as reported in Table 1. The first region is characterized by a linear yield function and a non-associative flow rule. The second region, which guarantees the smooth transition, is characterized by the elliptic yield function and the non-associative flow rule. The third region is characterized by the elliptic function and the associative flow rule. 3.3. Numerical implementation The model permits a fully implicit numerical integration via a classical return mapping scheme performed in the strain invariant space (Borja and Tamagnini, 1998), leading to a system of non linear equations with three unknowns. As usual, it is assumed that the updated displacements are given, which implies that the elastic trial strain εeTr and the solution at time step tn are known. The task is to compute the updated stress state b and the discrete plastic multiplier Dg for a given displacement increment. According to Borja and Tamagnini (1998), consider the following local residual equations generated by the applied strain increment Dε

8 e 9 < εv  εeTr v þ DgvP G = r ¼ rðxÞ ¼ εes  εeTr ; s þ DgvQ G ; : F

8 e 9 < εv = x¼ εe : : s ; Dg

(32)

The aim is to dissipate the residual vector r by finding the solution vector x* using a Newton's method

xkþ1 ¼ xk þ dxk ;

Ak dxk ¼ rk ;

Ak ¼

vrk ; k)k þ 1; vxk

(33)

where k is the iteration counter. Depending on the loading path, different yield functions and plastic potential will be involved. In particular, three different Return Mapping schemes have to be implemented, based on Eq. (32), according to the Table 1 The three plastic regions. Description

Linear non-associative Elliptical non-associative Elliptical associative

Mean stress

Yield function

Plastic potential

P

F

G

P > P* e > P > P P e P
F1 F2 F2

G G F2

Abbr.

Lin. NA Ell. NA Ell. A

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

123

three regions in Table 1. A closed form expression for the consistent tangent operator can be written for the three different set of nonlinear equations. The details can be found in Appendix A. If the elliptical surface is activated, then at each iteration the value of Pi needs to be updated, according to Eq. (29). The p current value of the volumetric plastic strain εv can be easily computed by assuming the additive decomposition of the logarithmic strains p εv ¼ εev þ εpv ¼ εeTr v þ εv;n ;

(34)

where εpv;n is the plastic volumetric deformation at time step tn. The consistent tangent operator, required in Eq. (16), is defined in principal direction as

aep ¼

vb : vεeTr

(35)

The details of computation can be found in Appendix B. The summary of the procedure to select the appropriate Return Mapping algorithm is shown schematically in the flowchart of Fig. 5. The complete procedure to compute the update state of the stress and the consistent tangent operator is summarized in Table 2. As shown in Fig. 4, the linear yield surface and the plastic potential functions create corners on the dilatant side of the hydrostatic axis. In this case, i.e. when the return mapping on the linear surface gives as result a negative value of Q, a standard procedure for apex formulation need to be implemented. The reader can refer to de Souza Neto et al. (2009) for further details. 4. Validation of the constitutive model This section of the paper addresses two issues: the procedure to identify the parameters of the constitutive model and the comparison of the model against experimental data, thus showing the capability and the accuracy to reproduce different laboratory tests. The complete set of parameters describing the constitutive law is listed in Table 3, subdivided in hyperelastic response parameters, compactive plastic response parameters and dilatant plastic response parameters. As far as the first group of parameters (I), k, P0 and εev0 identify the hydrostatic hyperelastic response, since they correlate the elastic volumetric deformation εev with the pressure applied to the specimen according to Eq. (24a), when εes vanishes. Therefore, they can be assessed considering the elastic part of the curve obtained by an hydrostatic laboratory test on the rock. The value of k can be obtained recalling that Ke ¼ P/k, where K e is the bulk modulus of the rock and can be inferred as the tangent of the hydrostatic curve in a ðP; εev Þ diagram. As first hypothesis, it is convenient to assume εev0 ¼ 0, and then compute P0 through the value of the initial bulk modulus K0e ¼ P0 =k, which it can be inferred as the tangent in εev ¼ 0 of the hydrostatic curve. The parameters m0 and a, adopted to compute the elastic shear modulus me, control the deviatoric deformation with respect to the applied stress. In particular, the parameter a describes the pressure-dependence of the elastic shear modulus me. Since the mean normal

Fig. 5. Procedure for the selection of the correct Return-Mapping algorithm.

124

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144 Table 2 Numerical procedure to compute the updated stress state and the consistent tangent operator at the element level. Implicit elastic predictor/return-mapping algorithm Elastic deformation predictor: be Tr ¼ f nþ1 ,ben ,f Tnþ1 . P Spectral decomposition: be Tr ¼ 3A¼1 ðleA Tr Þ2 nTr ðAÞ 5nTr ðAÞ . e Tr Principal elastic logarithmic strains: εTr A ¼ lnðlA Þ. Deformation invariants εev Tr and εes Tr according to Eq. (19). Elastic stress predictor PTr and QTr according to Eq. (24). Check if yielding: F 2 ðP Tr ; Q Tr ; Pi;n Þ  0? No: elastic step, go to (9.). Yes: continue. 7. Check if yielding F 1 ðP Tr ; Q Tr Þ  0? (a) Yes: Linear yield surface. Solve Lin. NA, compute εev , εes , Dg and set Pi ¼ Pi,n.  Check correctness of algorithm: P  P*? No: elliptical yield surface, go to (8.). Yes: continue.  Check cone apex: Q  0? No: Apex algorithm, stop. Yes: correct, go to (10.). (b) No: check plastic point: PTr  P*? No: elastic step, go to (9.). Yes: continue. 8. Elliptical yield surface. e  For every iteration k check: P  P? No: Solve Ell. NA. Yes: Solve Ell. A.  Update Pi according to Eqs. (34) and (29).  Compute εev , εes , Dg, go to (10.). 9. Elastic step. Set ðεev ; εes Þ ¼ ðεev Tr ; εes Tr Þ, Pi ¼ Pi,n and (P, Q) ¼ (PTr, QTr). 10. Principal elastic logarithmic strain εeA and Kirchhoff stress bA. 11. Consistent tangent operator aAB. 12. Updated left elastic CauchyeGreen tensor: P be ¼ 3A¼1 ðexpðεeA ÞÞ2 nTr ðAÞ 5nTr ðAÞ .

1. 2. 3. 4. 5. 6.

Table 3 Parameters for the constitutive model. Group

Description

Parameters

I II III

Hyperelastic response Plastic response (compactive side) Plastic response (dilatant side)

k, P0, εev0 , m0, a Pi0 , A, ε*, r m, c0, m

stress does not vary significantly in common applications, the shear modulus can be assumed constant, independent from the pressure variation. The value of m0 can be estimated from an averaged bulk modulus through the relation m0 ¼ 3(1  2n)Ke/ 2(1 þ n), where n is the Poisson's ratio of the rock. As regards the second group of parameters (II), the compactive plastic response is characterized by Pi0 and A, which describe the initial position and the size of the elliptical yield surface. They both can be evaluated by plotting the initial yield stress in a (P, Q) diagram, for different loading paths, and by interpolating the experimental data with an ellipse. Fig. 6 represents experimental data on the initial yield stress (corresponding to pore collapse) and the evolution of the yield stress as function of the plastic volumetric strain for different sandstones. The data are interpolated by different ellipse, with constant minor semiaxis A and variable major semiaxis B, as assumed for the proposed constitutive model. The pore collapse pressure C  ¼ Pi0 þ A, corresponding to the inflection point in the hydrostatic curve, can provide further information to define Pi0 and A. The parameters ε* and r describe the hardening law of the compactant side, and they can be calibrated considering the plastic part of the curve in the hydrostatic test, and minimizing the differences between the computed data and the experimental data. This can be done rigorously, by assuming a range of values for the two parameters and by running a set of simulations, and finally selecting the values of the parameters that minimize a target function (for example the mean difference between the experimental/computed values). As regards as the third group of parameters (III), m and c0 determine the plastic dilatant behavior. Again, they can be assessed by interpolating the experimental data of the dilatant plastic deformation on a (P, Q) diagram, for different loading paths. Alternatively, if the cohesion c and the angle of internal friction f are known, the parameters m and c0 can be estimated imposing that the linear yield surface matches the MohreCoulomb yield criterion, by forcing both criteria to predict the same collapse loads under plane strain conditions. In this case (Chen and Mizuno, 1990) the constants m and c0 read

pffiffiffi 3tanf m ¼  3 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi; 9 þ 12tan2 f

c0 ¼

pffiffiffi 3 3 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c: 9 þ 12tan2 f

(36)

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

125

Fig. 6. Evolution of the elliptical surface as function of the initial yield stress C* and the plastic volumetric strain for (a) Bentheim, (b) Berea and (c) Adamswiller sandstones. Experimental data are taken from Wong and Baud (2012).

Fig. 7. Loading path with confining pressure equal to (a) 70 MPa and (b) 250 MPa.

126

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

The parameter m, that controls the non-associative flow rule for the dilatant side of the yield surface, is usually assumed to be lower than m. This fact reflects the experimental evidence, in the sense that the dilatancy angle is lower than the frictional angle. In the last two decades, several experiments have been conducted on porous sandstones, with particular interest on the evaluation of the compaction mechanism and the transition from a dilatant to a compactant behavior (Wong et al., 1997, 2001; Vajdova and Wong, 2003; Baud et al., 2006). The standard tests consist on hydrostatic and triaxial experiments, where different axial and confining pressures are applied, monitoring the porosity change and the axial deformation. The experimental results for hydrostatic and triaxial tests have been compared with the numerical results obtained with the developed constitutive model. The procedures to reproduce numerically the experimental test with the proposed constitutive model are explained in detail in Appendix C. Let us consider first a hydrostatic compression experiment, where the porosity reduction is measured as a function of the mean stress. Fig. 8 represents the results obtained from the model compared with the experimental data of a hydrostatic compaction test for several sandstones. The value of the parameters adopted for the different rocks are listed in Table 4. It can be observed that the numerical results match quite satisfactorily the experimental data. In particular the model is able to capture the non linear response of the elastic phase, the hydrostatic yield pressure (identified by the knee in the curve) and the hardening behavior during the plastic compaction. Especially, this model turns to be particularly accurate for the rocks where the elastic or plastic response is highly non linear. The constitutive model can describe accurately also the unloading deformation process. Let us now consider a triaxial compression experiment, where the porosity reduction is measured as a function of the mean stress, for different values of confining pressure Pc. Fig. 7 represents two examples of loading paths, starting from a different confining pressure. The figure shows that for low values of the confining pressure the loading path intersects the dilatant yield surface, instead for high values it intersects the cap surface, determining a progressive expansion according to the hardening law. The transition region is reached when the loading path intersects the compactant surface, which, under continued loading is pushed outward until the other surface is also intersected, thus activating both yield surfaces. The smooth transition between the two surfaces guarantees that also the transition between the two mechanisms is smooth, without any discontinuity between the two processes. The model therefore reflects the main concepts of the brittleeductile transition theory. Fig. 9 shows the results obtained, comparing the experimental data and the prediction obtained by the model for Berea, Bentheim and Darley Dale sandstones. The model can capture well the plastic mechanisms, in particular with respect to the value of the stress yielding and the hardening behavior of the rock. Among the three porous rocks, the model fits better the experimental data for Berea and Bentheim. These two rocks are characterized by a sharp transition between the elastic and the elastoplatic behavior, with a dramatic decrease of porosity after the yielding pressure has been reached. The results are less accurate for the Darley Dale sandstone, for which the transition from the elastic to the plastic response is smoother. The accuracy of the experimental data is crucial to identify correctly the parameters of the model. A high deviation of the experimental data may compromise the correct identification of the parameters of the model. Comparing the values among the different rocks listed in Table 4, most of the parameters are characterized by a small range of values; therefore these data can be taken as reference also for other type of porous rocks, and slightly mitigate an eventual error in the experimental data. The parameters with the largest variation range are Pi0 and A, which change significantly depending on the rock type. Therefore, the identification of these two parameters is particularly important. 5. Numerical simulations The equilibrium equation and the constitutive model proposed have been implemented in a non linear finite element code, to compute the stress status, the elastoplastic deformations and the condition for band initiation around a horizontal wellbore. First of all, the non linear finite element code has been validated taking as reference a laboratory test conducted on a cylindrical sample of Bentheim sandstone, with a borehole drilled along the axis. The details of the experimental test can be found in the work by Dresen et al. (2010). The specimen was equipped with two pairs of orthogonally oriented strain-gages and four single strain-gages oriented in a circumferential direction, glued onto the cylindrical surface of the rock, as shown in Fig. 10. As for as the parameters for the constitutive model, the data in Table 4 for Bentheim sandstone have been adopted. Fig. 11 shows the comparison between the experimental and numerical results, in terms of circumferential (radial) strain. It can be observed that there is a good agreement between the measured and the computed deformations. This study investigates a horizontal wellbore drilled in Campos Basin field (Coelho et al., 2005), a reservoir located 290 Km offshore Brazil. The water depth is about 1900e2400 m and the reservoir is located 5000 m under the sea bed, below a salt layer that may reach 2000 m. The formation is characterized by high porosity (20e30%) and low permeability (1e10 mD). An experimental investigation was conducted to define the mechanical properties of the rock formation, testing samples with classical laboratory experiments. The parameters emerging from the investigations are reported in Table 5. As discussed in the previous section, the availability of laboratory data allows for the identification of the model's parameters. Fig. 12 reproduces the hydrostatic compression test conducted on a sample from the drilled rock formation. As explained before, these experimental data have been used to identify k, P0, εev0 , Pi0 , A, ε* and r. The shear modulus m0 has been assumed constant, and its value has been estimated from the bulk modulus, assuming n ¼ 0.15. The parameters to define the linear yield function,

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

127

Fig. 8. Comparison between experimental data and model simulation of a hydrostatic test for (a) Boise, (b) Berea, (c) Bentheim, (d) St. Peter and (e) Adamswiller sandstones. p is the Cauchy mean stress and Dn is the porosity reduction. Experimental data are taken from Wong and Baud (2012).

i.e. m and c0, have been assumed so as to approximate the MohreCoulomb yield surface through Eq. (36). The cohesion c of the rock and the angle of internal friction f are listed in Table 5. Since this kind of rocks usually do not exhibit a significant dilation, the dilatancy angle has been assumed equal to 5 , and m was estimated again using the first of Eq. (36). Table 6 summarizes the parameters used for the constitutive model in the numerical simulations. As regards the in-situ geostatic stresses measured in the reservoir production region, the effective vertical stress is equal to sV ¼ 32.1 MPa. The effective principal horizontal stresses are equal in both principal directions, sH ¼ sh ¼ 9.0 MPa. The openhole wellbore radius is R ¼ 4.2500 ¼ 107.95 mm. The wellbore axis is parallel to the direction of the principal horizontal far-field stress sh. Since a wellbore is geometrically a 1D structure, it is reasonable to assume that deformations are constrained in the plane perpendicular to the axis. Hence, plane strain conditions are used, so the principal stresses acting on the plane of the borehole section are the vertical and the horizontal stresses. The effects of the gravitational load are already accounted for by the in-situ effective stress field, so the density of the material is ignored. The domain size considered for the numerical

128

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

Table 4 Parameters for the constitutive model. Parameter

Sandstone Boise

Berea

Bentheim

St. Peter

Adamswiller

k (MPa) P0 (MPa) εev0 m0 (MPa)

0.015 10 0 500 0 40 30 0.22 1.5 1.0 100 0.7

0.015 12 0 500 0 190 180 0.1 0.8 0.6 100 0.4

0.010 15 0 500 0 210 200 0.12 0.4 0.6 120 0.4

0.030 10 0 500 0 200 120 0.18 1.2 1.0 145 0.8

0.015 10 0 500 0 100 100 0.15 1.5 0.7 70 0.6

a Pi0 (MPa) A (MPa) ε* r m c0 (MPa) m

simulation is ten times the radius of the bore. Since the problem is doubly symmetric, the discretized domain represents only one quarter of the complete geometry. The finite element geometry, with the boundary conditions, is provided in Fig. 13. The mesh consists of quadrilateral serendipity elements with 8 nodes and 9 Gauss point. The number of elements is equal to 1536 (512 elements are inside the hole), with 32 elements along the wall of the bore, as shown in Fig. 14.

Fig. 9. Mean stress versus volumetric strain for triaxial compression experiments on (a) Berea, (b) Bentheim and (c) Darley Dale sandstones. The confining pressures are indicated by numbers (in MPa) next to each curve. Experimental data are taken from Baud et al. (2006).

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

129

Fig. 10. Section of the cylindrical sample of Bentheim sandstone with a borehole. Geometric details of the samples can be found in Dresen et al. (2010).

Fig. 11. Comparison of experimental data and numerical simulation of the evolution of the circumferential (radial) strain [%] in a cylindrical specimen with a hole. The confining pressure is equal to (a) 92 MPa, (b) 95 MPa and (c) 98 MPa. Experimental data are taken from Dresen et al. (2010).

Table 5 Available material parameters (Coelho et al., 2005). Young's Modulus

Poisson ratio

Cohesion

Friction angle

E (MPa)

n

c (MPa)

f ( )

1200

0.15

8.5

42

Fig. 12. Hydrostatic compression test on Campos Basin field (Coelho et al., 2005).

130

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

Table 6 Parameters data for high porosity rock of Campos Basin field. Parameter

k

P0 (MPa)

εev0

m0

a

Pi0 (MPa)

Value

0.01

5.0

0.01

600

0

25

Parameter

A (MPa)

ε*

r

m

c0 (MPa)

m

Value

18

0.1

1.5

1.08

10

0.15

Fig. 13. Plane strain domain of a quarter of the borehole. (a) In-situ configuration: only far field stress applied without the hole. (b) Drilling configuration: far field stress and mud pressure, with the hole.

The analysis consists of two stages: first the vertical and horizontal stresses are applied in undisturbed conditions, i.e. as the rock formation is still intact, in order to simulate the in-situ stress before the drilling process. Subsequently, the elements corresponding to the borehole are removed from the domain, to simulate the process of drilling, and the mud pressure is applied. The numerical simulations investigate two aspects: (i) the stresses and plastic deformations, depending on the value of the mud pressure and (ii) the conditions for the inception of localized bands. 5.1. Stresses and plastic deformations The first analysis considers different values of the pressure DP ¼ Pm  Pp applied on the internal wall of the cavity. Pm is the pressure of the mud used to support the wellbore and Pp is the pore pressure surrounding the hole. Except for special cases,

Fig. 14. Discretization of the domain with quadrilateral serendipity elements with 8 nodes.

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

131

the minimum mud pressure corresponds to the pore pressure in the rock formation so that the well does not flow while drilling. The maximum mud pressure should not exceed the total minimum far field stress, in order to avoid unintentional hydraulic fracturing. Hence, we assumed a range of pressure spanning from DP ¼ 0 to DP ¼ 8 MPa. The aim of the analysis is to assess the stress configuration around the bore, the deformations and the plastic mechanism. The most interesting results are shown in Fig. 15, which represents the results obtained in terms of radial and circumferential stress, for DP ¼ 0,4,6,8 MPa. At increasing mud pressure we observe a decrease of the radial effective stress (Fig. 15, left column), which leads to a higher compressive stress on the rock along the wall of the wellbore. On the other hand, the distribution of the effective circumferential stress (Fig. 15, right column) does not change significantly depending on the

Fig. 15. Radial (left) and circumferential (right) effective stress (MPa) for different value of pressure (a) DP ¼ 0, (b) DP ¼ 4, (c) DP ¼ 6 and (d) DP ¼ 8.

132

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

applied mud pressure. The stress field around the hole is always compressive (negative), so that the rock formation is never subjected to a tensile stress that may cause a tensile fracture. The stress distribution has a strong implication in the extension of plastic zones around the wellbore, as shown in Fig. 16 that plots the results in terms of volumetric and deviatoric plastic strains. The plastic zone is always located in correspondence of the maximum effective stress, i.e. along the direction of the minimum far field stress. At increasing mud pressure we observe that the zone involved by the plastic deformations becomes smaller, with consequent increase of stability of the wellbore. It is interesting to observe that, increasing the value of the mud pressure, not only the extension of the plastic zone changes, but also the plastic mechanism involved. This can be observed by focusing on the volumetric plastic deformation

Fig. 16. Volumetric (left) and deviatoric (right) plastic strain for different values of pressure (a) DP ¼ 0, (b) DP ¼ 4, (c) DP ¼ 6 and (d) DP ¼ 8.

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

133

Fig. 17. Volumetric (left) and deviatoric (right) plastic strain considering (a) a compactant-dilatant elasto-plastic model and (b) only a dilatant elasto-plastic model with DP ¼ 8.

(Fig. 16, left column), recalling that positive volumetric plastic deformations are associated with dilatant mechanism, while negative volumetric plastic deformations are associated with compactant mechanism. For low values of the mud pressure (DP ¼ 0 MPa (a)) we observe a wide plastic zone along the bore wall, in the direction of the minimum far field stress. In this area the plastic volumetric strain along the wall is positive, therefore the stress is such that the plastic mechanism is dilatant. At the same time, it can be observed a negative plastic volumetric strain zone, associated with the compactant mechanism, situated more inside with respect to the wall of the bore. The higher the values of the mud pressure, the narrower is the area characterized by the dilatant plastic mechanism, while the area with a negative plastic volumetric deformation increases. For high values of the mud pressure (DP ¼ 8 MPa (d)) the area associated with the dilatant mechanism disappears, and the plastic zone is only due to the compactant mechanism. This is in complete agreement with the mechanical behavior of porous rocks observed experimentally, where, under an increasing confining pressure on the specimen, the plastic mechanism changes from dilatant to compactant. To underline this aspect and quantify the differences, we performed a numerical simulations with a pure dilatant plastic constitutive model. The model has the same features of the constitutive law proposed in this work, beside the fact that it does not have a compactant surface which limits the hydrostatic pressure. Therefore, the constitutive model can be considered as a classical DruckerePrager constitutive law. Fig. 17 shows the results in terms of volumetric and deviatoric deformations, with an applied pressure equal to DP ¼ 8 MPa. As it can be clearly observed, in particular focusing on the volumetric deformations contour in the left side, the results differ significantly in terms of extension of the plastic zones and in terms of the plastic mechanism involved. This is due to the fact that, with the pure dilatant plastic model, there are no limits for the effective mean normal stress, which can increase arbitrarily. Hence, this model underestimates the extension of the plastic zone around the wellbore. 5.2. Conditions for band localization In this second section, we discuss the conditions for the initiation of a band of intense deformation around the wellbore, by analyzing the localization function. The localization function F is calculated according to Eq. (18), by computing the minimum value of the determinant of the Eulerian acoustic tensor defined in Eq. (17). Since the plane of the shear band is well defined, the search for the stationary point of F can be made by a sweep over half a unit circle. We defined localization when one or more Gauss points have detected the first negative incursion of the localization function. In this study we focused on the influence of (i) the in-situ stress condition and (ii) the presence of a small imperfection around the wall of the wellbore.

134

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

Fig. 18. Localization function for different value of the vertical in-situ stress (a) sV ¼ 32.1, (b) sV ¼ 33.9, (c) sV ¼ 37.1 MPa with balanced drilling (DP ¼ 0 MPa).

Fig. 19. Localization function for the vertical in-situ stress sV ¼ 33.9 (inception of the localization) (a) 672, (b) 1536 and (c) 1632 elements with balanced drilling (DP ¼ 0 MPa).

5.2.1. In-situ stress condition This set of analysis investigates the in-situ stress condition leading to the inception of a localized band. To address this issue, we incremented the vertical in-situ stress with a constant horizontal stress. Fig. 18 represents the results in terms of localization function, obtained with a mud pressure value equal to DP ¼ 0 MPa. When the in-situ vertical stress is equal to sV ¼ 33.9 MPa (case (b)), we observe the initiation of a deformation band, arising from the bore wall. For an effective vertical stress equal to sV ¼ 37.1 MPa (case (c)) the presence of the band is evident. Remark 1. It is important to underline that the simulations are valid only up to the bifurcation point and slightly beyond it, otherwise the numerical results will exhibit mesh dependence, and finite element enhancement techniques are needed (Andrade and Borja, 2006; Sanavia et al., 2006; Lazari et al., 2015). In fact, the finite element solutions presented above do suffer from a mild pathological mesh dependence once the onset of localization has been detected. However, it is not the aim of this research to describe the post-bifurcation behavior but rather the definition of the conditions leading to the initiation of a band. It has been demonstrated that the inception of a shear band, in the case of drained localization, is not mesh-dependent (see for example Andrade and Borja (2007)). To strengthen this aspect, numerical simulations with different mesh refinements were performed. Fig. 19 shows the results in terms of localization function, obtained with the same boundary conditions but with different discretization of the domain around the hole. In all the cases the inception of the localization is obtained for the same value of the vertical in-situ stress sV ¼ 33.9, and the origin of the shear band is very similar.

Fig. 20. Localization function for different values of the vertical in-situ stress (a) sV ¼ 32.1, (b) sV ¼ 35.2, (c) sV ¼ 37.7 MPa with balanced drilling (DP ¼ 2 MPa).

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

135

Fig. 21. Discretized domain with weak elements, to reproduce an imperfection along the wall of the wellbore.

As the three meshes predict a vary similar onset of localization, this proves that the mesh does not influence the inception of bands. We obtained similar results by decreasing the value of the horizontal in-situ stress, with a constant effective vertical stress. The numerical simulations suggest that, for the rock formation under study, a ratio between the vertical and horizontal stress of 4.1 is the threshold to develop strain localization. We also investigated the influence of a variation of the mud pressure on the band initiation. Fig. 20 shows the most relevant results obtained by increasing the vertical effective stress with a constant mud pressure DP ¼ 2 MPa. The initiation of the first localization band is detected when the vertical stress is increased to the value of sV ¼ 35.2 MPa. Comparing this result with the condition for band initiation when no overpressure is applied at the wall, we conclude that an increase of the mud pressure slightly prevents the formation of the bands. 5.2.2. Presence of a small imperfection So far we considered a perfectly circular wellbore, thus assuming a regular geometry in the finite element analysis. It is very unlikely that a wellbore has a perfect circular shape, in particular when drilled in a formation of high porosity rocks. During drilling operations breakouts may occur, usually along the direction of the minimum far field stress, but these breakouts are tolerated as long as the stability of the well is not a concern. This section shows how the presence of an imperfection in the direction of minimum in-situ stress can influence the inception of bands around the hole. At this scope, we modified the discretized domain as shown in Fig. 21, by inserting some weak elements with a ten times lower bulk modulus. This is equivalent to assume that, while drilling, not only the elements corresponding to the hole are removed, but also some other elements (those with a significantly lower bulk modulus). In this fashion, we simulated the partially removed elements around the wall of the hole. The vertical in-situ stress was incremented progressively, while the horizontal in-situ stress and the mud overpressure were maintained constant, respectively sh ¼ 9.0 MPa and DP ¼ 0 MPa. Fig. 22 shows that, in the presence of an imperfection, the formation localizes at a significant lower value of the effective vertical stress. When sV ¼ 25.2 MPa (case (a)) the first bifurcation points appear, leading to a shear band of intense deformation. Increasing the vertical stress we observe that the band propagated even further. The origin of the band is in correspondence of the imperfection around the wall, and it propagates along the vertical direction. To address the influence of the mud pressure

Fig. 22. Localization function for different values of the vertical in-situ stress (a) sV ¼ 25.2 (band initiation) and (b) sV ¼ 27.8 with applied pressure DP ¼ 0 MPa assuming an imperfection around the wall.

136

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

Fig. 23. Localization function for different value of the vertical in-situ stress (a) sV ¼ 25.2 (band initiation) and (b) sV ¼ 28.4 with applied pressure DP ¼ 2 MPa assuming an imperfection around the wall.

on band localization with an imperfection along the wall, we applied a radial pressure and we increased progressively the vertical effective stress. Fig. 23 shows the most relevant results for DP ¼ 2 MPa applied to the borehole wall. In this case we observe that the mud pressure does not influence significantly the band initiation. In fact, the effective vertical stress associated with band initiation is again sV ¼ 25.2 MPa. The diagram in Fig. 24 summarizes all the results obtained as concerns the band initiation. The graph shows for each value of the mud pressure the corresponding value of the vertical stress leading to a band initiation. We observe that, with an imperfection, whatever the value of the mud pressure is, the vertical stress for the initiation of a localized band is always significantly lower. Hence, the presence of imperfections can play an important role in the occurrence of the localization condition around a borehole, with relevant consequences in the well stability. 6. Conclusions This paper presented a study about the plasticity and the conditions for the initiation of strain localization around a horizontal wellbore, drilled through a high porosity rock formation, in a region characterized by high in-situ stresses. Due to the particular mechanical behavior of porous rocks, an elastoplastic model has been developed and implemented. In particular, the model has been proposed and validated to capture the dilatant and compactant plastic mechanisms. Boundaryvalue problems simulating a horizontal wellbore have been solved using a nonlinear finite element method in the assumption of finite strains. The results of the study show that under certain conditions a wellbore can experience both compactant and dilatant plastic mechanisms. Therefore, an appropriate constitutive model is fundamental to capture correctly the plastic behavior of the rock formations. As demonstrated by the simulations, the in-situ stress condition and the mud pressure have a strong impact as regards both the plastic zones and the band localization. Furthermore, the results show that the presence of an imperfection along the wall drastically influences the failure of the formation.

Fig. 24. Vertical stress sV corresponding to band initiation for different values of the mud pressure DP.

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

137

The analyses presented in this paper allow a better understanding of the elastoplastic behavior experienced by the rock around a borehole. These results can form the basis to practical considerations to limit the sand production and to prevent the instability of the well. Due to the extreme impact that this topic has in real applications, more research is encouraged. In this work we focused the condition for band initiation but also the research on the post-bifurcation behavior is of paramount importance, and it is still an open issue. Furthermore, the study of plasticity and localization around a non circular holes with an accurate constitutive model should be further investigated. Acknowledgment Financial support for nine months was provided to the first author by the Fondazione Gini, while this research was being conducted at Stanford University. The authors wish to thank Prof. R.I. Borja of Stanford University for the much appreciated motivation and support to the research. Appendix A. Derivation of the Jacobian matrix in Eq. (33) The 3  3 Jacobian matrix Ak in Eq. (33) for return mapping iteration is obtained by differentiating the component of r with respect to the variables x defined in Eq. (32). For simplicity, the superscript k denoting the value at k-th iteration is abbreviated in the following derivation. Depending on the return mapping algorithm, we will compute three different matrices, i.e. ANA Lin , A e ANA , A . The Hessian matrix of j , D ¼ VV j reads Ell Ell

De11 ¼ 



P0 3a e 2 expðUÞ 1 þ εs 2k k

(A.1a)

De22 ¼ 3ðm0  aP0 expUÞεes

(A.1b)

3P0 aεes expðUÞ k

De12 ¼ De21 ¼

(A.1c)

For the linear non-associative return mapping, the matrix ANA Lin takes the form

2

ANA Lin

1 0 ¼4 De11 vP F 1 þ De21 vQ F 1

0 0 De12 vP F 1 þ De22 vQ F 1

3 vP G vQ G 5; 0

(A.2)

where vP G ¼ m, vQ G ¼ 1, vP F 1 ¼ m and vQ F 1 ¼ 1. For the elliptical non-associative return mapping, the matrix ANA Ell takes the form

2

ANA Ell

1 0 ¼4 De11 vP F 2 þ De21 vQ F 2 þ Kp vPi F 2

0 0 De12 vP F 2 þ De22 vQ F 2

3 vP G vQ G 5; 0

(A.3)

2 2 2 where vP F 2 ¼ 2B2 ðP  Pi Þ, vQ F 2 ¼ 2A2 Q and vPi F 2 ¼ vB vPi ððP  Pi Þ  A Þ  2b ðP  Pi Þ with plastic hardening modulus and reads 2

vB2 vPi

¼ 2mðmPi þ c0 Þ. Kp is the

  Pi0 rε ε r1 ; x x2

Kp ¼ vεev Pi ¼

(A.4)

with x ¼ ðε  εev Tr þ εev  εpv;n Þ. To compute the Jacobian matrix for the elliptical associative return mapping, let us introduce the Hessian matrix of F 2 ðP; Q Þ with Pi fixed

H ¼ VVF 2 ¼

H11 H21

H12 H22

"

¼

v2PP F 2 v2QP F 2

v2PQ F 2 v2QQ F 2

#

2 ¼2 B 0

0 ; A2

(A.5)

and define the matrix G ¼ HDe. For the elliptical associative return mapping, the matrix AAEll take the form

2 AAEll

6 ¼6 4

  1 þ Dg G11 þ Kp v2PPi F 2   Dg G21 þ Kp v2QPi F 2 De11 vP F 2 

where v2PPi F 2 ¼ 2

vB2 vPi

þ De21 vQ F 2

ðP  Pi Þ  B2



DgG12

vP F 2

7 7 vQ F 2 5 ;

1 þ DG22 De12 vP F 2

þ De22 vQ F 2

and v2QPi F 2 ¼ 0.

3

þ Kp vPi F 2

0

(A.6)

138

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

Appendix B. Derivation of the consistent tangent operator aAB in Eq. (35) The consistent tangent operator is obtained taking the derivative of the principal Kirchhoff stress with respect to the trial logarithmic principal strain,

aep ¼

vb : vεeTr

(B.1)

In case of elasticity, εeTr ≡ εe and the tangential elastic tensor ae can be computed taking advantage of the invariant decomposition. The derivation of the stress invariants with respect of the strain invariants reads

vP ¼ De11 d þ vεe

rffiffiffi 2 e b; D n 3 12

vQ ¼ De21 d þ vεe

rffiffiffi 2 e b; D n 3 22

(B.2a)

(B.2b)

b (see Appendix A). b ¼ ee =jjee jj and De is the symmetric Hessian matrix of j where n The consistent tangent operator is computed using Eq. (B.2) and reads (Borja and Andrade, 2006)

 De11

aAB ¼

rffiffiffi   2Q 2Q 2 e 2 bA n b þ De21 n b A dB þ b B Þ þ De22 n bA n bB:  e dA dB þ ðd  n D d n 9εs 3 12 A B 3εes AB 3

(B.3)

In case of elastoplasticity, consider the following strain derivative

a

ep

vb vP :¼ e Tr ¼ d5 e Tr þ vε vε

rffiffiffi rffiffiffi b 2 vQ 2 vn b5 n Q 5 e Tr ; þ 3 3 vεe Tr vε

(B.4)

where

     b vn vðee =ee Þ v ee Tr ee Tr 1 1 b 5n b : ¼ ¼ ¼ e Tr I  d5d  n e Tr e Tr e Tr 3 vε vε e vε

(B.5)

Substituting (B.5) in (B.4), and using the elements of the matrix De to enforce the chain rule, we have

a

ep

rffiffiffi      e e  2 vεe vεe 2Q 1 e vεv e vεs b 5 De21 v þ De22 s þ b 5n b : n d5d  n I  ¼ d5 D11 e Tr þ D12 eTr þ 3 3 3εeTr vε vε vεeTr vεeTr s

(B.6)

The task is then reduced to determine the strain derivatives of the invariants εev and εes , which are obtained from the discretized flow rule (13) expressed in term of invariants

  vεev v vG eTr ; ε ¼  Dg vP vεeTr vεeTr v

  vεes v vG eTr : ε ¼  Dg vQ vεeTr vεeTr s

(B.7) p

Eq. (B.7) can be written 2  2 operator Dp, with components Dij , which maps the pffiffiffiffiffiffiffiffi in a more efficient way introducing the b onto the derivatives with respect to εeTr of the elastic strain invariants εev and εes basis vectors d and 2=3 n

vεev ¼ Dp11 d þ vεe Tr

rffiffiffi 2 p b; D n 3 12

vεes ¼ Dp21 d þ vεe Tr

rffiffiffi 2 p b D n 3 22

(B.8)

Finally substituting (B.8) in (B.6) and defining

Dep ¼ De Dp ;

(B.9)

we obtain the desired consistent tangent operator ep

aAB ¼



ep

D11 

rffiffiffi   2Q 2 ep 2Q 2 bA n b B þ Dep b A dB þ b B Þ þ Dep b n b : d d þ ðdAB  n n n D12 dA n B A 21 eTr 3 3 22 A B 9εs 3εeTr s

(B.10)

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

139

The close form of a (where the superscript ep has been omitted for the sake of lightness) depends on the specific return mapping algorithm, therefore three different operators need to be derived. In order to reduce the derivatives to their lowest order, we expand (B.7). The expansion of the first of Eq. (B.7) reads

b11

  vεev vεes vG vDg þ b ¼ c d  ; 12 1 vP vεeTr vεeTr vεeTr

(B.11)

while the expansion of the second of Eq. (B.7) reads

vεe vεe b21 e vTr þ b22 e sTr ¼ c2 d þ vε vε

rffiffiffi   2 vG vDg b n 3 vQ vεe Tr

(B.12)

The coefficients bij and ci will be computed in the remaining part of the section, leading to a close form for the operator Dp. The strain-gradient of Dg is obtained from the overall consistency condition

vF ¼ 0: vεeTr

(B.13)

Let us start from the linear non associative return mapping algorithm. In this case

b11 ¼ vεev r1 ¼ 1;

b12 ¼ vεes r1 ¼ 0;

(B.14)

b21 ¼ vεev r2 ¼ 0;

b22 ¼ vεes r2 ¼ 1;

(B.15)

c1 ¼ 1;

c2 ¼ 0:

(B.16)

Substituting in (B.11), (B.12) we obtain

vG vDg vεev ¼d ; vP vεeTr vεeTr

vεes ¼ vεeTr

rffiffiffi v  vDg 2 b G n : vQ vεeTr 3

(B.17)

In this case, Eq. (B.13) reads

vF 1 ðP; Q Þ vF 1 vP vF 1 vQ ¼ þ ¼ 0; vP vεeTr vQ vεeTr vεeTr

(B.18)

which can be written as

vF 1 vεev vεes ¼ d1 eTr þ d2 eTr ¼ 0; eTr vε vε vε

(B.19)

where

d1 ¼ De11 vP F 1 þ De21 vQ F 1 ;

d2 ¼ De12 vP F 1 þ De22 vQ F 1 :

(B.20)

Substituting (B.17) into (B.19) and solving for vDg/vεeTr we obtain

vDg ¼ a1 d þ a2 vεeTr

rffiffiffi 2 b; n 3

(B.21)

where a1 ¼ d1/e, a2 ¼ d2/e and e ¼ d1 vP G þ d2 vQ G. Finally, the last step involves back substitution. Inserting (B.21) into (B.17) and rearranging the terms we have p

DLinNA ¼



1  a1 vP G a2 vP G : a1 vQ G 1  a2 vQ G

(B.22)

For elliptical yield surface with non-associative flow rule, the procedure to compute DpEllNA is very similar to what derived so far. In this case the expression (B.17) remains unchanged, while (B.13) reads

vF 2 ðP; Q ; Pi Þ vF 2 vP vF 2 vQ vF 2 vPi ¼ þ þ ¼ 0; vP vεeTr vQ vεeTr vPi vεeTr vεeTr which can be written as

(B.23)

140

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

vF 2 vεe vεe vF ¼ d1 e vTr þ d2 e sTr þ KpTr 2 ¼ 0; vPi vεe Tr vε vε

(B.24)

where

d1 ¼ De11 vP F 2 þ De21 vQ F 2 þ Kp vP F 2 ;

d2 ¼ De12 vP F 2 þ De22 vQ F 2

(B.25)

and KpTr ¼ vεev Tr Pi ¼ Kp . Substituting (B.17) into (B.24) and solving for vDg/vε we obtain the same expression in (B.21) where now a1 ¼ ðd1 þ KpTr vPi F 2 Þ=e, a2 ¼ d2/e and e as before derived. Again, inserting (B.21) into (B.17) and rearranging the terms we have eTr

p



DEllNA ¼

1  a1 vP G a2 vP G : a1 vQ G 1  a2 vQ G

(B.26) p

Finally, let us compute the operator DEllA for elliptical yield surface with associative flow rule. Also in this case, the derivation follows the same steps as already done for the other algorithms. In this case we have

  b11 ¼ 1 þ Dg G11 þ Kp v2PPi F 2 ;   b21 ¼ Dg G21 þ Kp v2QPi F 2 ;

b12 ¼ DgG12 ;

(B.27)

b22 ¼ 1 þ DgG22 ;

(B.28)

  c1 ¼ 1  Dg G11 þ KpTr v2PPi F 2 ;

  c2 ¼ Dg G11 þ KpTr v2QPi F 2 :

(B.29)

The strain-gradient of Dg is obtained from the overall consistency condition (B.24), which again can be written as (B.24) with d1 and d2 as in (B.20). Solving (B.24) for Dg/vεeTr we obtain the same expression as (B.21), where

i. h a1 ¼ d1 ðb22 c1  b12 c2 Þ þ d2 ðb11 c2  b21 c1 Þ þ e bKpTr vPi F 2 e; pffiffiffiffiffiffiffiffi 2=3ðd2 b11  d1 b12 Þ=e;     e ¼ d1 b22 vP F 2  b12 vQ F 2 þ d2 b11 vQ F 2  b21 vP F 2 ; a2 ¼

(B.30) (B.31) (B.32)

with e b ¼ b11 b22  b21 b12 . Inserting (B.21) into (B.17) and rearranging the terms, we obtain

b

p   D11 ¼ b22 ðc1  a1 vP F 2 Þ  b12 c2  a1 vQ F 2 ;

(B.33a)

  pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi p D12 ¼ b12  1 þ 3=2a2 vQ F 2  3=2b22 a2 vQ F 2 ;

(B.33b)

b b

p   D21 ¼ b11 c2  a1 vQ F 2  b21 ðc1  a1 vP F 2 Þ;

(B.33c)

  pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi p D22 ¼ b11 1  3=2a2 vQ F 2 þ 3=2b21 a2 vP F 2 :

(B.33d)

b

In conclusion

b

DpEllA ¼

Dp : e b

(B.34)

Appendix C. Numerical simulation of hydrostatic compaction and triaxial tests This section describes the procedure to simulate numerically a hydrostatic compaction and triaxial test, and how to obtain numerical results comparable with experimental results. The advantage of this procedure is to avoid the solution of a complete boundary-value to reproduce the experiments.

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

141

Appendix C.1. Hydrostatic compaction test Let us assume the deformation process as governed by a scalar function c(t) > 0, c(0) ¼ 0 which maps the point from the undeformed configuration X to the current configuration x

xA ¼ ð1  cðtÞÞXA ;

A ¼ 1; 2; 3:

(C.1)

The deformation gradient F at time tn can be evaluated as

2 1  cn vxn 4 ¼ Fn ¼ 0 vX 0

0 1  cn 0

3 0 0 5: 1  cn

(C.2)

Hence, the relative deformation gradient f is

f ¼

vxnþ1 ¼ F nþ1 F 1 n : vxn

(C.3)

The role of c(t) can be regarded as a sequence of increments Dc. Then (C.2) and (C.3) take the form

2

3 1  ðn þ 1ÞDc 0 0 5; F nþ1 ¼ 4 0 1  ðn þ 1ÞDc 0 0 0 1  ðn þ 1ÞDc 3 2 1 þ ðn þ 1ÞDc 0 0 7 6 1 þ nDc 7 6 7 6 1 þ ðn þ 1ÞDc 7 6 f ¼6 7; 0 0 7 6 1 þ nDc 7 6 4 1 þ ðn þ 1ÞDc 5 0 0 1 þ nDc

(C.4)

(C.5)

where n indicates the number of the increment. The complete loading and unloading path can be simulated increasing and decreasing the function c(t). For each increment n, it is possible to define the deformation gradient F and compute the Jacobian J ¼ det(F). The total volumetric strain at each deformation increment is related to the Jacobian J according to

εv ¼

3 X

εA ¼ lnðJÞ;

(C.6)

A¼1

The porosity reduction Dn, defined as the difference between the initial porosity of the rock n0 and the porosity n at a certain deformation increment, can be then obtained from the total volumetric deformation. Assuming that the initial undeformed volume V0 ¼ 1 for the sake of simplicity, the deformed volume due to the hydrostatic compression process is equal to V ¼ 1  Dn. Hence the relation between εv and Dn follows from

εv ¼ lnðJÞ ¼ lnðV=V0 Þ ¼ lnð1  DnÞ:

(C.7)

The numerical procedure consists of computing through the Return Mapping the Kirchhoff stress for a given deformation increment, associated to a porosity reduction with Eq. (C.7). The Kirchhoff stress is then converted into Cauchy stress, to compare the computed quantity with laboratory tests. Appendix C.2. Triaxial test The procedure to reproduce numerically a triaxial test is stress driven. The value of P and Q ¼ 3(P  Pc) are increased by discrete steps, with P > Pc and Pc fixed, and for each increment the total volumetric strain εv is computed. This value is then converted in terms of porosity, according to Eq. (C.7). For any given point (P,Q), lying inside the elastic domain, the actual value of εev is computed according to Eq. (24). Since the stress point belongs to the elastic domain, εev corresponds to the total volumetric deformation εv. If the loading path intersects the yield surface in the linear portion, then the surface cannot expand anymore and the limit load, associated with a shear failure, is identified. Otherwise, if the loading path intersects the yield cap, the plastic surface expands, according to the hardening law. In this second case, for each load increment the value of Pi is computed, updating the elliptical plastic surface such that F 2 ¼ 0. Subsequently, the volumetric plastic deformation εpv corresponding to the stress status is computed, which, summed to the elastic counterpart εev , provides the total volumetric deformation εv. This procedure is summed up in the Table C.7.

142

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

Table C.7 Numerical procedure to simulate a triaxial test. Algorithm for triaxial test 1. Initial stress state and hardening parameter: (Pn þ 1, Qn þ 1) and Pin . 2. Check stress status: F 1 ðP nþ1 ; Q nþ1 Þ and F 2 ðP nþ1 ; Q nþ1 ; Pin Þ  0?  If No: elastic status. (a) Solve inversely Eq. (24) and compute εev and εes . (b) Update variables: εv ¼ εev and Pinþ1 ¼ Pin .  If Yes: plastic status. - Check loading path: P  P*? - If No: elliptical surface. (a) Solve inversely F 2 ðP nþ1 ; Q nþ1 ; Pinþ1 Þ ¼ 0 and compute Pinþ1 . (b) Solve inversely Eq. (29) and compute εpv . (c) Update variables: εv ¼ εev þ εpv . - If Yes: linear surface. No expansion allowed.

Appendix D. Notations Tables D.8 and D.9 summarize the notation used in the paper.

Table D.8 Summary of the nomenclature used through the paper. Nomenclature 1 a aep A A b be(Tr) B B c0 c C* F Fe Fp f F ð1;2Þ F G, G h0, h J k K Ke Lv m  m n n P Q P0 Pi P* e P r r s t v x X

Identity tensor Eulerian acoustic tensor Consistent tangent operator in principal direction Minor semiaxis of the ellipse Consistent tangent operator for the Return Mapping algorithm Left CauchyeGreen Elastic left CauchyeGreen (Trial) Major semiaxis of the ellipse Body configuration Intersection of the yield surface F 1 Intersection of the plastic potential function G Pore collapse pressure Deformation gradient Elastic deformation gradient Plastic deformation gradient Relative deformation gradient Yield function (1: Dilatant, 2: Compactant) Localization function Plastic potential function Band thickness Jacobian Elastic compressibility modulus Tangent stiffness matrix Elastic bulk modulus Lie derivative Slope of the yield surface F 1 Slope of the plastic potential function G Eigenvector of tensor be and t Band orientation Mean normal effective stress Deviatoric effective stress Mean normal effective stress at εev ¼ εev0 Centroid of the ellipse Tangent point of F 1 with F 2 Tangent point of F 2 with G Coefficient for the volumetric hardening Local residual equation of the Return Mapping algorithm Deviatoric stress tensor Nominal traction vector Velocity jump Spatial point Material point

b

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

143

Table D.9 Summary of the greek symbols used through the paper. Greek symbols

a a bA d Dg

eðTrÞ

εA eðTrÞ εv eðTrÞ εs p εv εps εev0 ε*

g h

eðTrÞ

lA

me m0 n f r0G s t J j

b j e j

Coefficient for elastic shear modulus Algorithmic tangent tensor Eigenvalues of the tensor t Ones vector Discretized plastic multiplier Logarithmic elastic principal stretches (Trial) Elastic volumetric deformation (Trial) Elastic deviatoric deformation (Trial) Plastic volumetric deformation Plastic deviatoric deformation Elastic volumetric strain at P ¼ P0 Volumetric deformation at ultimate compaction Plastic multiplier Weighting function Square root of the eigenvalues of the tensor be (Trial) Elastic shear modulus Constant elastic shear modulus Poisson ratio Motion of the body B Reference body force Cauchy stress tensor Kirchhoff stress tensor Free energy defined by the tensor be Free energy defined in logarithmic elastic principal stretches Free energy in the invariant space Free energy for isotropic loading

References ABAQUS, 2013. Theory Manual 6.13. Dassault Systmes Simulia Corp. Adachi, J., Nagy, Z.R., Sayers, C.M., Smith, M., Becker, D.F., 2013. Drilling adjacent to salt bodies: definition of mud weight window and pore pressure using numerical models and fast well planning tool. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. http://dx.doi.org/10. 2118/159739-MS. Andrade, J.E., Borja, R.I., 2006. Capturing strain localization in dense sands with random density. Int. J. Numer. Meth. Eng. 67, 1531e1564. Andrade, J.E., Borja, R.I., 2007. Modeling deformation banding in dense and loose fluid-saturated sands. Finite Elem. Anal. Des. 43, 361e383. Aubertin, M., Li, L., 2004. A porosity-dependent inelastic criterion for engineering materials. Int. J. Plast. 20, 2179e2208. Baud, P., Klein, E., Wong, T.-f, 2004. Compaction localization in porous sandstones: spatial evolution of damage and acoustic emission activity. J. Struct. Geol. 26, 603e624. Baud, P., Vajdova, V., Wong, T.-f, 2006. Shear-enhanced compaction and strain localization: inelastic deformation and constitutive modeling of four porous sandstones. J. Geophys. Res. 111, B12401. Baud, P., Vinciguerra, S., David, S., Cavallo, A., Walker, A., Reuschle, T., 2009. Compaction and failure in high porosity carbonates: mechanical data and microstructural observations. Pure Appl. Geophys. 166, 869e898. Baud, P., Zhu, W., Vinciguerra, S., Wong, T.-f, 2015. Micromechanics of dilatancy and compaction in basalt. EGU General Assem. 2015. Borja, R.I., 2002. Bifurcation of elastoplastic solids to shear band mode at finite strain. Comput. Methods. Appl. Mech. Eng. 191, 5287e5314. Borja, R.I., 2004. Computational modeling of deformation bands in granular media. ii. numerical simulations. Comput. Methods. Appl. Mech. Eng. 193, 2699e2718. Borja, R.I., 2013. Plasticity. Springer, Berlin-Heidelberg. Borja, R.I., Andrade, J.E., 2006. Critical state plasticity. Part vi: Meso-scale finite element simulation of strain localization in discrete granular materials. Comput. Methods Appl. Mech. Engrg. 195, 5115e5140. Borja, R.I., Aydin, A., 2004. Computational modeling of deformation bands in granular media. i. geological and mathematical framework. Comput. Methods. Appl. Mech. Eng. 193, 2667e2698. Borja, R.I., Tamagnini, C., 1998. Cam-clay plasticity part iii: extension of the infinitesimal model to include finite strains. Comput. Methods. Appl. Mech. Eng. 155, 73e95. Borja, R.I., Sama, K.M., Sanz, P.F., 2003. On the numerical integration of three-invariant elastoplastic constitutive models. Comput. Methods Appl. Mech. Eng. 192, 1227e1258. Borja, R.I., Song, X., Rechenmacher, A.L., Abedi, S., Wu, W., 2013. Shear band in sand with spatially varying density. J. Mech. Phys. Solids 61, 219e234. Buscarnera, G., Das, A., Laverack, R.T., 2014, December. Plasticity-based interpretation of compaction banding in a high-porosity carbonate rock. In: AGU Fall Meeting Abstracts, vol. 1, p. 4321. Callari, C., Auricchio, F., Sacco, E., 1998. A finite-strain cam-clay model in the framework of multiplicative elasto-plasticity. Int. J. Plast. 14, 1155e1187. Cao, Y., Deng, J., Yu, B., Tan, Q., Ma, C., 2014. Analysis of sandstone creep and wellbore instability prevention. J. Nat. Gas Sci. Eng. 19, 237e243. http://dx.doi. org/10.1016/j.jngse.2014.05.013. Carroll, M.M., 1971. A critical state plasticity theory for porous reservoir rock. Recent Adv. Mech. Struct. Contin. 97, 935e950. Chen, W.-F., Mizuno, E., 1990. Nonlinear Analysis in Soil Mechanics. Theory and Implementation. Elsevier. Chen, S.L., Abousleiman, Y.N., Muraleetharan, 2012. Closed-form elastoplastic solution for wellbore problem in strain hardening/softening rock formations. Int. J. Geomech. 12, 494e507. Chen, S., Abousleiman, Y., Abass, H., 2014. An analytical elasto-plastic analysis for stability of axisymmetric wellbore. In: Materials Technology; Petroleum Technology, vol. 5. ASME. http://dx.doi.org/10.1115/OMAE2014-24619. V005T11A026. Coelho, L.C., Soares, A.C., Ebecken, N.F.F., Alves, J.L.D., Lau, L., 2005. The impact of constitutive modeling of porous rocks on 2-d wellbore stability analysis. J. Petrol. Sci. Eng. 46, 81e100. Di Maggio, F.L., Sandler, I.S., 1971. Material model for granular soils. J. Eng. Mech. Div. 97, 935e950. Dresen, G., Stanchits, S., Rybacki, E., 2010. Borehole breakout evolution through acoustic emission location analysis. Int. J. Rock Mech. Min. Sci. 47, 426e435.

144

N. Spiezia et al. / International Journal of Plasticity 78 (2016) 114e144

Drucker, D.C., Prager, W., 1952. Soil mechanics and plasticity analysis of limit design. Quart. J. Appl. Math. 10, 157e162. Elyasi, A., Goshtasbi, K., 2015. Using different rock failure criteria in wellbore stability analysis. Geomech. Energy Environ. 2, 15e21. http://dx.doi.org/10. 1016/j.gete.2015.04.001. Fossum, A.F., Fredrich, J.T., 2000. Cap plasticity models and compactive and dilatant pre-failure deformation. In: Proceedings of the 4th North American Rock Mechanics Symposium, A.A. Balkema,Rotterdam, pp. 1169e1176. Grueschow, E., Rudnicki, J.W., 2005. Elliptical yield cap constitutive modeling for high porosity sandstone. Int. J. Solids Struct. 42, 4574e4587. Haimson, B., 2001. Fracture-like borehole breakouts in high-porosity sandstone: are they caused by compaction bands? Phys. Chem. Earth. Part A Solid Earth Geod. 26, 15e20. Haimson, B., Kovachic, J., 2003. Borehole instability in high-porosity berea sandstone and factors affecting dimensions and shape of fracture-like breakouts. Eng. Geol. 69, 219e231. Hamid, O.H., Rahim, Z., Ba-Wazir, O., 2014. Wellbore stability analysis using acoustic radial profiles and an elastoplastic model. In: Abu Dhabi International Petroleum Exhibition and Conference. Society of Petroleum Engineers. http://dx.doi.org/10.2118/171838-MS. Houlsby, G.T., 1985. The use of a variable shear modulus in elastic-plastic for clays. Comput. Geotech. 1, 3e13. Houlsby, G.T., Amorosi, A., Rojas, E., 2005. Elastic moduli of of soils dependent on pressure: a hyperelastic formulation. Geotechnique 55, 383e392. Issen, K.A., 2002. The influence of constitutive models on localization conditions for porous rock. Eng. Fract. Mech. 69, 1891e1906. Issen, K.A., Challa, V., 2005. Strain localization conditions in porous rock using a two-yield surface constitutive model. Geol. Soc. Lond. 245, 433e452. Jefferies, M.G., 1997. Nor-sand: a simple critical state model for sand. Geotechnique 47, 255e272. Lazari, M., Sanavia, L., Schrefler, B., 2015. Local and non-local elasto-viscoplasticity in strain localization analysis of multiphase geomaterials. Int. J. Numer. Anal. Methods Geomech. 39, 1570e1592. http://dx.doi.org/10.1002/nag.2408. http://doi.wiley.com/10.1002/nag, 2408. Le, T.M., Fatahi, B., Khabbaz, H., 2015. Numerical optimisation to obtain elastic viscoplastic model parameters for soft clay. Int. J. Plast. 65, 1e21. Lee, E.H., 1969. Elastic-plastic deformations at finite strains. J. Appl. Mech. 36, 1e6. Lin, J., Xie, S., Shao, J.-F., Kondo, D., 2012. A micromechanical modeling of ductile behavior of a porous chalk: formulation, identification, and validation. Int. J. Numer. Anal. Methods Geomech. 36, 1245e1263. http://dx.doi.org/10.1002/nag.1048. http://doi.wiley.com/10.1002/nag, 1048. Marinelli, F., Buscarnera, G., 2015. Parameter calibration for high-porosity sandstones deformed in the compaction banding regime. Int. J. Rock Mech. Min. Sci. 78, 240e252. http://dx.doi.org/10.1016/j.ijrmms.2015.05.004. McAuliffe, C., Waisman, H., 2015. On the importance of nonlinear elastic effects in shear band modeling. Int. J. Plast. 71, 10e31. Nguyen, L.D., Fatahi, B., Khabbaz, H., 2014. A constitutive model for cemented clays capturing cementation degradation. Int. J. Plast. 56, 1e18. Nicot, F., Challamel, N., Lerbet, J., Prunier, F., Darve, F., 2010. Bifurcation and generalized mixed loading conditions in geomaterials. Int. J. Numer. Anal. Methods Geomech. 35 http://dx.doi.org/10.1002/nag.959 n/aen/a. Nicot, F., Sibille, L., Hicher, P.-Y., 2014. Micromacro analysis of granular material behavior along proportional strain paths. Contin. Mech. Thermodyn. 27, 173e193. http://dx.doi.org/10.1007/s00161-014-0347-8. Nova, R., Castellanza, R., Tamagnini, C., 2003. A constitutive model for bonded geomaterials subject to mechanical and/or chemical degradation. Int. J. Numer. Anal. Methods Geomech. 27, 705e732. http://dx.doi.org/10.1002/nag.294. Ogden, R.W., 1984. Nonlinear Elastic Deformations. Ellis Horwood. Parsamehr, H., Mohammadi, S.D., Moarefvand, P., 2015. Numerical modeling of wellbore stability in layered rock masses. Arabian J. Geosci. http://dx.doi.org/ 10.1007/s12517-015-1962-9. Rahimi, R., Nygaard, R., 2014. What Difference Does Selection of Failure Criteria Make in Wellbore Stability Analysis?. Roshan, H., Rahman, S.S., 2011. Analysis of pore pressure and stress distribution around a wellbore drilled in chemically active elastoplastic formations. Rock Mech. Rock Eng. 44, 541e552. http://dx.doi.org/10.1007/s00603-011-0141-x. Sanavia, L., Pesavento, F., Schrefler, B. a, 2006. Finite element analysis of non-isothermal multiphase geomaterials with application to strain localization simulation. Comput. Mech. 37, 331e348. http://dx.doi.org/10.1007/s00466-005-0673-6. Shen, W., Shao, J., 2015. A micromechanical model of inherently anisotropic rocks. Comput. Geotech. 65, 73e79. http://dx.doi.org/10.1016/j.compgeo.2014. 11.016. Shen, W.Q., Shao, J.F., Kondo, D., Gatmiri, B., 2012. A micro-macro model for clayey rocks with a plastic compressible porous matrix. Int. J. Plast. 36, 64e85. http://dx.doi.org/10.1016/j.ijplas.2012.03.006. Shin, H., Kim, J.-B., 2015. Physical interpretations for cap parameters of the modified drucker-prager cap model in relation to the deviator stress curve of a particulate compact in conventional triaxial testing. Powder Technol. 280, 94e102. Shojaei, A., Taleghani, A.D., Li, G., 2014. A continuum damage failure model for hydraulic fracturing of porous rocks. Int. J. Plast. 59, 199e212. Simo, J.C., 1992. Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory. Comput. Meth. Appl. Mech. Eng. 23, 61e112. de Souza Neto, E.A., Peric, D., Owen, D.R.J., 2009. Computational Methods for Plasticity: Theory and Applications. John Wiley and Sons. Stefanov, Y.P., Chertov, M.A., Aidagulov, G.R., Myasnikov, A.V., 2011. Dynamics of inelastic deformation of porous rock and formulation of localized compaction zones studied by numerical modeling. J. Mech. Phys. Solids 59, 2323e2340. Udegbunam, J.E., Aadnø y, B.S., Fjelde, K. K. r, 2013. Uncertainty evaluation of wellbore stability model predictions. In: SPE/IADC Middle East Drilling Technology Conference and Exhibition, vol. 124, p. 15. http://dx.doi.org/10.2118/166788-MS. Vajdova, V., Wong, T.-f, 2003. Incremental propagation of discrete compaction bands: acoustic emission and microstructural observation on circumferentially notched samples of bentheim sandstone. Geophys. Res. Lett. 30, 17e75. Weng, M.-C., Ling, H.I., 2013. Modeling the behavior of sandstone based on generalized plasticity concept. Int. J. Numer. Anal. Methods Geomech. 37, 2154e2169. Wong, T.-f., Baud, P., 2012. The brittle-ductile transition in porous rock: a review. J. Struct. Geol. 44, 25e53. Wong, T.-f., Szeto, H., Zhang, J., 1990. Effect of loading path and porosity on the failure mode of porous rocks. Appl. Mech. Rev. 45, 281e293. Wong, T.-f., David, C., Zhu, W., 1997. The transition from brittle faulting to cataclastic flow in porous sandstones mechanical deformation. J. Geophys. Res. 102, 3009e3025. Wong, T.-f., Baud, P., Klein, E., 2001. Localized failure modes in a compactant porous rock. Geophys. Res. Lett. 28, 2521e2524. Xie, S., Shao, J., 2012. Experimental investigation and poroplastic modelling of saturated porous geomaterials. Int. J. Plast. 39, 27e45. Xie, S., Shao, J., 2014. An experimental study and constitutive modeling of saturated porous rocks. Rock Mech. Rock Eng. 223e234. http://dx.doi.org/10.1007/ s00603-014-0561-5. Yamakawa, Y., Hashiguchi, K., Ikeda, K., 2010. Implicit stress-update algorithm for isotropic cam-clay model based on the subloading surface concept at finite strains. Int. J. Plast. 26, 634e658. Zhang, K., Zhou, H., Shao, J., 2013. An experimental investigation and an elastoplastic constitutive model for a porous rock. Rock Mech. Rock Eng. 46, 1499e1511. http://dx.doi.org/10.1007/s00603-012-0364-5. Zhang, W., Gao, J., Lan, K., Liu, X., Feng, G., Ma, Q., 2015. Analysis of borehole collapse and fracture initiation positions and drilling trajectory optimization. J. Petrol. Sci. Eng. 129, 29e39. http://dx.doi.org/10.1016/j.petrol.2014.08.021. Zhou, H., Bian, H.B., Jia, Y., Shao, J.F., 2013. Elastoplastic damage modeling the mechanical behavior of rock-like materials considering confining pressure dependency. Mech. Res. Commun. 53, 1e8. http://dx.doi.org/10.1016/j.mechrescom.2013.07.008. Zhu, Q., Shao, J., Mainguy, M., 2010. A micromechanics-based elastoplastic damage model for granular materials at low confining pressure. Int. J. Plast. 26, 586e602. http://dx.doi.org/10.1016/j.ijplas.2009.09.006. Zoback, M.D., 2010. Reservoir Geomechanics. Cambridge University Press.