Effects of nonuniform initial stress state on apparent fracture toughness

Journal Pre-proofs Effects of nonuniform initial stress state on apparent fracture toughness I.A. Garagash, A.A. Osiptsov PII: DOI: Reference:

S0013-7944(19)31092-6 https://doi.org/10.1016/j.engfracmech.2019.106837 EFM 106837

To appear in:

Engineering Fracture Mechanics

Received Date: Revised Date: Accepted Date:

31 August 2019 19 December 2019 20 December 2019

Please cite this article as: Garagash, I.A., Osiptsov, A.A., Effects of nonuniform initial stress state on apparent fracture toughness, Engineering Fracture Mechanics (2019), doi: https://doi.org/10.1016/j.engfracmech. 2019.106837

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier Ltd.

Eects of nonuniform initial stress state on apparent fracture toughness I.A. Garagasha,b , A.A. Osiptsova a Multiphase b The

Systems Lab, Skolkovo Institute of Science and Technology (Skoltech), Russia Schmidt Institute of Physics of the Earth of the Russian Academy of Sciences, Moscow, Russia

Abstract We analyze the nonuniform stress state in the vicinity of the fracture tip and its impact on apparent fracture toughness. It is shown that prescribing the boundary conditions at the deformed fracture faces allows one to take into account the eect of nonuniform tectonic stresses acting in the plane of the fracture. We conclude that the fracture may take an ellipsoidal form (aspect ratio departs from unity) as a result of interactions with anisotropic distribution of initial stresses and nonuniform growth in longitudinal and vertical directions. An expression for the eective stress intensity coecient (fracture toughness) is obtained as a function of the nonuniform initial stress state. A simplified asymptotic approximation to the solution for the eective apparent fracture toughness is given, which may be further implemented into the near-tip asymptotics used in the form of a tip element in the numerical implementation of a Planar3D approach to hydraulic fracture simulation. Keywords: fracturing, initial stress state, stress intensity coecient, fracture tip, fracture shape 1. Introduction Hydraulic fracturing technology is designed with the help of simulators, which are based on continuum mechanics models of fracture growth [1, 2, 3, 4]. Some of the most advanced simulators are based on Planar3D approach, which involves fracture propagation criterion at each point of the moving boundary. In hydraulic fracture simulations, it is necessary to take into account the deformation of a reservoir under the action of pressure on fracture surface, flow inside the fracture, and the process of fracture growth [5]. The criterion of fracture propagation is a boundary condition of a specific type, which governs the solution of the global problem of hydraulic fracture growth [1]. One of the most wide spread conditions is the criterion of fracture propagation, which is based on the Linear Elastic Fracture Mechanics (LEFM) [6]. According to LEFM, the stress state at the fracture tip is characterized by the stress intensity coecient. A rapid fracture propagation begins when the stress intensity coecient reaches a certain critical value, which was called fracture toughness and proposed as a material property of the medium [7]. The increase in fracture toughness leads to the decrease in fracture height and the increase in fracture length. Vice versa, with the decrease in fracture toughness the fracture becomes wide and short. Fracture toughness is determined experimentally on core samples without taking into account the initial stress state in the reservoir. Hence, a particular open question of interest is to which extent the fracture toughness depends on the initial stress state. The propagation criterion is typically based on a certain asymptotic solution obtained for the vicinity of the fracture tip in a specific regime. With the development of Irwin-Grith’s Email address: [email protected] (A.A. Osiptsov) Preprint submitted to Engineering Fracture Mechanics

fracture mechanics [8, 7], several parameters have been proposed for characterizing the resistance to fracture growth [9], and fracture toughness is the most widely accepted quantity, which is related with the stress intensity factor of an equilibrium tensile crack: K

= P(li ) = const:

where P is the loading parameter and the function  depends on fracture dimensions. Here, we would refer to the work [10], where this eect has been discussed in application to hydraulic fracturing problems. Particular formulas with respect to li were proposed, e.g. [11]. In [12], there is a universal correction formula for KI C (theoretically evaluated for fast cracks) that shows dependence on the crack velocity. This is relevant for an initial stage of growth, when fractures are short. In addition we refer to [13], where it was shown that K1C (energy release rate) should be a function of velocity. The stress intensity factor K uniquely defines the stress field around the fracture tip: Kc gig ()

i j = p

2r

where gi j () are angular functions. Thus, Kc can be regarded as the driving force of the fracture. Globally, fracture toughness is a combination of elastic and surface-formation properties, while at the micro-scale it relates to nonlinear separation process at the tip. Accurate determination of fracture toughness is of critical importance for comprehensive fracture propagation models capturing the entire parametric space, in which toughness-related energy release, viscous dissipation and leak-o compete on multiple length scales [14]. In particular, the solution [14] gives a uniformly applicable smooth asympJanuary 2, 2020

totic solution varying from toughness-dominated to viscositydominated regimes. The state-of-the art solution for the neartip region with a variety of eects is given in [5]. Impact of the anisotropy of fracture toughness on the propagation of Planar 3D hydraulic fracture was studied in [15], where toughness anisotropy was highlighted as a possible hydraulic fracture height containment mechanism as well as the need for its careful characterization beyond measurements in the directions of divider and arrester. This is a particular justification for the present work, where we will evaluate anisotropy of fracture toughness due to nonuniform initial stress. Fracture orientation depends on the direction of maximum and minimum principal stresses in the reservoir. The rate of fracture propagation in horizontal and vertical directions is generally dierent, so fractures frequently have non-circular form and may be assumed as having elliptic form (see Fig. 1, where the fractures are shown as ellipsoids). One may suppose that it is due to dierent conditions at the fracture tip moving in vertical and horizontal directions (zones I and II in Fig. 1, b). In what follows we will assume that the curvature of the fracture boundary in zones I and II can be neglected. As we discussed above, fracture propagation depends on apparent fracture toughness K1c [17]. The value of K1c is determined experimentally on rock samples (cores) with cuts [18]. In [19], it was shown that fracture toughness of the rock has an eect on characteristic size of the grains and tensile strength of the material. Initial stresses, acting in the direction normal to the fracture plane in a granular medium, lead to the increase of the fracture toughness [20]. Fracture propagation is also affected by the viscous shear stresses exerted from the fluid on the fracture faces [21, 13, 22]. Yet the question regarding the stresses acting in the direction parallel to the fracture plane is still open and being actively discussed [23]. According to the classic solution to the problem of a plane cut tension [24] the parameter K1c should not depend on the initial stress state in the direction of the cut. In what follows, we will demonstrate that by taking into account the non-zero fracture opening and by specifying the boundary conditions at the deformed fracture faces one may include into consideration the eects of initial confining stress in the fracture plane on the stress state near the tip, which leads to the decrease in fracture toughness. Note that usually taking into account the actual location of the boundaries does not aect the stress state, with exception for stability problems [25] and deformation of anisotropic media [26]. However, in the vicinity of the fracture tip (which is a singular point) the fracture boundary is subject to finite shear and strain, hence the perturbations of the boundary should be taken into account in solution to the problem.

Overall we consider a problem of a fracture in an unbounded elastic plane with initial stresses. Following the Sneddon solution [27], we will consider an elliptic fracture, with its shape governed by extension stresses and other parameters. Applying known solutions for an elliptic hole to the cavity, we will find the distribution of stresses at its boundary [28] and near the tip [29] under the influence of the initial stress parallel to the fracture plane. In what follows we will calculate the incremental increase in the potential energy of strain [17], which will finally allow us to generalize the Grith’s failure criterion for the fracture under the initial stresses and to obtain the Irwin’s eective stress intensity coecient. We understand that a hydraulic fracture evolves from the toughness-dominated regime (when the fracture is small) into a viscosity-dominated regime (when large). We assume for simplicity that the fracture is contained in the vertical direction by a fracture toughness contrast between the layers and that, assuming the vertical section is growing slowly in height, that section can be assumed to be uniformly pressurized. This assumption has been used historically in the literature on Pseudo3D (P3D) fracture models as one way to model height growth of each P3D cell. We split the problem into two. One problem corresponds to the triaxial compression under the pressure q, whereas the second problem – to the dry fracture loaded at infinity by the tractions p1 q and q p2 (Fig. 2). According to the solution from the linear elasticity theory [24], the boundaries of the fracture under the load q p2 have the following displacements: w a

=

uy (0; x) a

r

= (1 )

q

p2

1

G

x2 ; a2

(1)

where G is Young’s modulus and  – Poisson’s ratio. In this case the cut is transformed into an ellipse (Fig. 3), since the formula (1) corresponds to the equation of an ellipse: x2 a2

+

w2 b2

= 1;

b = (1

)

q

p2 G

a:

(2)

Typically, when simulating the fracture growth within linear elastic fracture mechanics, the boundary conditions are formulated on unperturbed initial boundaries of the fracture, hence expression (1) does not depend on horizontal stresses p1 q. However, as a result of the action of tractions q p2 the fracture opens, and then the action of the loads p1 q can be considered on new deformed boundaries. For doing so, we will consider an elastic hole in elastic medium in the coordinates  and , which is subject to horizontal loads p1 q (Fig. 4). We will now use the known solution to the problem of an elliptic hole in an otherwise unbounded elastic medium [28] and write the solution for circular stresses at the boundary of the ellipse  = 0 the form:

2. Problem Formulation We consider a plane problem of a fracture of the length 2a embedded in an unbounded elastic medium, under the action of horizontal loads p1 and vertical loads p2 (Fig. 2). The fracture contains fluid under the pressure q. Since the fracture is hydraulically induced, we have q > p2 and p1 > p2 .

( )=0

= (p1

q)

sinh(20 ) exp(20 ) cos 2 + 1 : cos 2 cosh(20 )

(3)

At the top of the ellipse  = 0, the stress  is determined by 2

the relation: ( )=0 ;=0

= (p1

q)

sinh(20 ) exp(20 ) + 1 : 1 cosh(20 )

(Fig. 6). For doing so, following [17], we will first calculate the potential energy of an elastic deformation W1 accumulated under the action of the loading q p2 :

(4)

Taking into account that cch(0 ) = a;

csh(0 ) = b;

c2

= a2

b2 ;

W1 (5)

=

(q p2 )2 a2 2G

(1

)

(12)

Then, taking into account that the fracture under the loading (q p2 ) takes an elliptic form (see formula (2)), we will find the contribution W2 into the potential energy of deformation under the loading (p1 q). In the paper [17] there was a derivation presented for the expression of the potential energy of deformation W of an elastic plane with the elliptic hole for a two-axial extension (see Fig. 7 below).

we obtain that the following tensile stress is acting at the fracture tip: (yy ) x=a;y=0 = ( )=0 ;=0 = p1 q: (6) Note that according to Eq. (3), when moving from the fracture tip to the center the stress  decreases and reaches zero at " # 1 1 + 2 sinh(0 ) cosh 0  = arccos ; (7) 2 1 + 2 sinh 0 (sinh 0 + cosh 0 ) reaching the magnitude of (yy ) x=0;y=b

= ( )= ;==2 = 0

(p1

q)

(8)

at  = =2. Thus the action of confining stresses (p1 q) results in the formation of tensile stresses at the fracture tip. This is particularly clear from Fig. 5, a on an example of a fracture of the length a = 10 m for the material with the Young’s modulus G = 3
109 Pa and Poisson’s ratio  = 0:2 for x  a. Vertical stresses yy to the right from the fracture tip at x  a on the axis y = 0 are determined by the relation [29]: (yy ) xa;y=0

"

(x=a)(x2 =a2 1 + 2b2 =a2 ) + (x2 =a2 1 + b2 =a2 )3=2 # (1 + m)2 (m 1)2 (m + 1)(3s2 m) (9) + ; 2(s2 m) 2(s2 m)3

= (p1 +

q) 1

Figure 7: Biaxial extension of an elliptic cavity (hole).

For an ellipse in the conditions of a plane strain W

r

s=

x 2B

+ m=

x 2 1 m; B = (a + b); 2B 2 a b q p2 ; b = (1 ) a: a+b G

(

= (

= p1

4G

h

i

")2 + 4(1 "2 )m + (1 + ")2 (1 + m2 ) n2 ;

(1 ) 2(1

(13)

m= 1

(10)

b a

!

1+

b a

! 1

a 2

b a

; n = (1 + ):

Taking into account that

Distribution of the stress  =(p1 q) at the fracture tip for x  a, calculated according to formula (9), is shown in Fig. 5, b. We note that, according to the distribution of loads shown in Fig. 2, the full stress at the fracture tip takes the form: f ull  )=0

 p2

where



f ull yy ) x=a;y=0

=

b a and

= (1 ) "=

(11)

q

p2

(14)

G

p1 q q p2

(15)

In accordance with (13) for a narrow ellipse b=a obtain the expression:

Thus the stress yy is acting at the fracture tip, which is accelerating fracture propagation. From which it follows that the fracture has an elliptic (non-circular) form, as the initial tectonic stresses are dierent near the fracture boundary in the vicinity of zones I and II (Fig. 1). Consider the eect of the initial stress state on the stress intensity coecient. First of all, consider the generalization of energy criterion (Grith) to the case of a hydraulic fracture

W2

=

(p1 q)2 (q p2 )a2 4G2

(1

 )2 :



1 we

(16)

Fracture propagation begins when the incremental increase in the surface energy, which grows due to the creation of new rock surface (as a result of failure during crack propagation). The 3

increase in surface energy is balanced by the potential strain energy. The balance condition for the creation of a new segment a + da of a fracture has a form:

@ (W + W2 + ) = 0; @a 1

account the nonuniform initial stress state of the medium. Here P=

(17)

4 G 3 (1 )

!2

(p1

p2 )

A=2



p + a(p1

p2 )2 ; 

G

(p1

p2 ) ;

1  #3 2 G Q=2 (p1 p2 ) 3 1   2p  G a (p1 p2 ) (p1 p2 )2 3 1  2G 2 K : 1  1c "

where

 = 4a

(18)

is the fracture surface energy, – surface energy density. Incorporating the expressions for W and  into (17), we obtain failure criterion in the form: "

p2 )2 1 + (1

(q

)

(p1 q)2 2G(q p2 )

#

=

4G

: a(1 )

In order to evaluate the eect of the dierence (p1 p2 ) on the strength, we will consider the dierence between the Irwin’s coecient and the eective fracture toughness coecient:

(19)

Note that equation (19) coincides with the Grith criterion, if the second term in parentheses is neglected. We transform the relation (19) as follows: s

1 + (1

K1

p

) a

[(p1

p

p2 ) K1 = a]2 2GK1

r

=


K1c =

p = (q p2 ) a is the stress intensity coecient. where

(21)

The following magnitude is the material constant, which is characterizing apparent fracture toughness: r

1 + (1

K1

p

 ) a

[(p1

p

p2 ) K1 = a]2 2GK1

= K1c

"

2G (1 )

#

(p1

p

p2 ) K12 + a(p1 p2 )2 K1

2G 2 K =0 1  1c (23)

where r

=2

2

0

66 1 B P B Q cos 6664 arccos BBB@ 3 3 2

3 P

!3=2 13 C 77 C C 77 C C A75

A 3

(27)



sandstone and siltstone 0:4  K1c



limestone 0:4  K1c



mudstone (argillite) 0:3  K1c

 1:65 MPa
m1=2

 0:95 MPa
m1=2  1:2 MPa
m1=2

As shown above, the fracture toughness coecient (25) depends on the dierence between the maximum and the minimum initial stresses p1 p2 . A conclusion can be drawn from (25) that the eective fracef f ture toughness K1c depends on the fracture half-length a and decreases with fracture growth. The qualitative nature of varief f ation of the ratio K1c =K1c is shown in Fig. 9. Curves 1,2,3 are plotted for  = 0:2 and K1c = 5
105 Pa
m1=2 (a) and K1c = 1
106 Pa
m1=2 (b). G = 6
109 Pa (curve 1), G = 1:2
1010 Pa (2), G = 1:8
1010 Pa (3), respectively. The dierence of initial stresses is related to the initial stress state of a heavy elastic half-space by the formula

Solving this equation, we will represent the strength condition as ef f K1 = K1c ; (24)

ef f K1c

:

3. Discussion (22)

From this relation we drive the cubic equation for K1 : K13 +2

K1c

In Fig. 8, we plot the variation of the fracture coecient increment
K1c (in percents) for K1c = 5
105 Pa
m1=2 (a) and K1C = 1
106 Pa
m1=2 (b). Dierent curves in Fig. 8 are plotted for a = 25m,  = 0:2, G = 6
109 Pa (curve 1), G = 1:2
1010 Pa (curve 2), and G = 1:8
1010 Pa (curve 3).

4G 1 

We will now derive the Irwin’s fracture toughness coecient. Rewrite Eq. (20) as s

ef f K1c

K1c

Assume that the maximum initial stress p1 is equal to the overburden layer weight (density  = 2500kg=m3 and layer height H = 3000m), and the stress dierence varies in the range 1
107 Pa  (p1 p2 )  5
107 Pa. In calculations, we will assume that typical values of the parameter K1c for oil reservoirs are in the following ranges [18]:

4G ; 1  (20)

K1

(26)

(p1

(25)

p2 ) = p1

1 2 1 

at a depth H = 3000m. The initial fracture toughness is specified as K1c = 5
105 Pa
m1=2 (Fig. 9, a) and K1c = 1
106 Pa
m1=2

is the eective fracture toughness coecient, which takes into 4

(Fig. 9, b). It is clear that with fracture growth the apparent fracture toughness rapidly decreases, but then the rate of decrease slows down. The higher the initial fracture toughness, the less pronounced is this eect. The increase in the shear modulus has a similar eect. If we talk about the gravitating fluid-saturated half-space, then the stress dierence will be decreased by the magnitude of the pore pressure and will take the form p1

p2

= (p1

p)

when moving from layer to layer. At the same time, the maximum vertical stress p1 (overburden weight) varies with depth linearly, as expected. Thus, looking back at formula (25) and Fig. 8 one may assert that fracture toughness of the rock depends significantly on the distribution of minimum confining stress, and, hence, varies noticeably with depth, which will have a profound impact on the law of fracture propagation in the toughness-dominated regime that is typical for the earlier stages of hydraulic fracture propagation. For example, when simulated in the Planar3D approach, the condition of the fracture propagation K1 = K1c is enforced in the form of a tip element [1, 34, 14], where the solution near the tip is given by the square-root asymptotics and then matched (actually, glued) to the global solution for the finite fracture in terms of the width. Hence, variation of K1c with depth will affect the rate of fracture propagation in horizontal and vertical directions thereby influencing the fracture shape (aspect ratio), e.g., as predicted within the Planar3D formulation [35, 36]. In practical applications, fracture height growth control is one of the major issues [37], so a more advanced model for the neartip asymptotics governing anysotropic fracture propagation in vertical and horizontal direction due nonuniform initial stress state will add value to the development of fracturing simulators used in the field for operation design purposes.

1 2 1 

ef f . In this case, the ratio K1c =K1c will be smaller (blue curves in Fig. 10) than for the dry reservoir (red curves). In this case the ratio is smaller and it manifests itself with the increase in the initial apparent fracture toughness (compare Fig. 10 a and b). We note that in the real rock the eect of the initial stress state should be considered with account for the increase in the initial value of K1c , which is due to the growth of the plastic dissipation zone in the vicinity of the fracture tip with the increase in the fracture half-length [10, 30]. Using uniform pressure inside the fracture is consistent with toughness dominated hydraulic fracture growth [31]. The change to the eective toughness predicted by this work by including the additional stress from far-field parallel stresses will impact significantly the initial stages of the fracture propagation in the toughness-dominated regime until the transition to the viscosity-dominated regime. One may argue that since nearly all field-scale hydraulic fracturing treatments are, however, in the fluid viscosity dominated propagation regime and in that regime the fracture growth is insensitive to fracture toughness, the present correction to the apparent fracture toughness may not have a noticeable impact on the growth of real finite fractures in oilfield services applications. This is a reasonable viewpoint, but we would like to refer to the chain of works on the impact of morphology of the initial (small) fractures, which begin propagation from perforation tunnels and grow in the toughness-dominated regime, interact with pre-existing natural fractures (if the target rock is shales, then with bedding planes; in other words, rock fabric [32]). The impact of this morphology (spatial shape) on the subsequent fracture growth in the viscosity-dominated regime is profound. We refer to the issue of near-wellbore tortuosity [33], which needs to be modelled, apparently, within the toughness-dominated propagation regime, which justifies the need for accurate consideration of the apparent fracture toughness presented in this work. In order to illustrate the impact of the parameter (p1 p2 ) on hydraulic fracture propagation, we will refer to the case shown in [1], which was related to the reservoir with two oil-bearing layers of sandstone: the upper layer (Pay Zone 1) at the depth 492–509 m, and the lower layer (Pay Zone 2) at the depth 628– 655 m (Fig. 11). Between these two layers there is an alternating array of shale and sandstone layers with varying geomechanics and transport properties. Downhole measurements provided the vertical profile of minimum confining stress p2 (Fig. 11). As it is clear from the Figure, the stresses p2 change abruptly with depth

Figure 11: Distribution of minimum in situ confining stress  vs. depth (Fig. 9 from [38]).

The problem of elasticity theory for the equilibrium of a fracture under the action of pressure is nonlinear in the vicinity of the fracture tip. Nevertheless, in the fracture theory the linear solution is used successfully [27] and the singularity in the fracture tip is characterized by the stress intensity coecient. We note that homogeneous initial stresses are acting in the fracture plane (as in the case of geomechanics problems), whereas these stresses do not enter into the solution of the classical problem. Since the full solution of the nonlinear fracture propagation problem appears to be too complex, the Sneddon solution could be considered as a first approximation, and further improvements (next-order terms) could be constructed based on 5

this first-order approximation. As a basis for the second-order approximation, in our view, one should take the solution from the first-order approximation – the transformation of a cut into an ellipse. For determination of the associated variations of the stress state one may use the classic solution for an elliptic hole in an elastic plane [39]. As a result of this approach, one can determine the contribution of initial stresses into the potential energy of deformations and evaluate their impact on the stress intensity coecient. We note that in this approach, as in the first approximation (the Sneddon solution [27]), finite deformations and rotations near the fracture tip are not taken into account. If one would aim at taking into account finite deformations, then it would take us beyond linear elasticity theory in the first approximation and require considering plastic flow of the medium, since elastic finite deformations are not typically realized in nature. Hence, analytical solution cannot be found in this case and one needs to use numerical methods.

unity [43]. It should be pointed out that such a fracture when pressurized uniformly will have higher stress intensity along the short axis, at its top and bottom sides [44]. Thus, fracture growth resulting in an elliptical profile must be justified by invoking other factors such as horizontal layering and associated stress contrasts, modulus contrasts, and bedding planes and their eect on fracture height growth [43]. Numerical simulation further demonstrates that the contrast in geomechanics parameters (layering, bedding planes) and stress intensity coecient in vertical and horizontal directions results in heightcontainment [45] and, hence, non-radial (e.g., elliptic) shape, which justifies the assumption for a non-radial elliptic shape of a fracture made in the problem formulation section of the present paper. ef f 3.3. Simplification of formula for K1c for engineering applications Here we would like to revisit the derivation of formula (25) for the eective apparent fracture toughness and simplify it for most practical applications. Assume that there is a small parameter (p1 p2 )=G  1, so the cubic equation can be significantly simplified and reduced to a quadratic equation, so the solution can be represented in the form 1 :

3.1. Toughness-dominated vs. viscosity-dominated regimes. The analysis presented in this paper has been performed for a fluid driven crack while it propagates in the toughnessdominated regime. On the other hand, in case of the viscosity dominated regime, this influence maybe even stronger, but it should be a subject of a separate research. Indeed the lag of a fluid with high viscosity from the tip of a propagating fracture diminishes the impact of the stress intensity coecient on the process of rock failure [31]. However, asymptotic solutions can be constructed, which allow one to neglect the lag [40]. In [41], it was demonstrated that in many practical cases the lag is important only at very early stages of the hydraulic fracturing treatment and it gradually vanishes, especially in the cases of high initial stresses at large depths. It should be noted that modeling of hydraulic fracturing on large rock outcrops demonstrates that the stress intensity coecient governs the fracturing process even in the case of injection of a high-viscosity fluid [1]. Contemporary numerical models of fracture propagation allow one to use the stress intensity coecient as an input parameter [3]. Note that even in the viscosity-dominated regime there is an action of initial confining stress in the fracture plane and it affects the propagation of a fracture tip. Hence, the eect of the initial stresses will remain even in the so-called zero-toughness case [42]. Evaluation of this eect goes beyond the scope of the present paper and will be considered in a separate study. As well as the impact on fracture redirection as a result of tangential stresses [21].

ef f K1c

= p

K1c 1+

2

+

; =

1  4GK1c

p

a(p1

p2 )2

(28)

In Fig. 12, we plot the cross-comparison of the numerical results by the full solution (25) and by the asymptotic approximation (28). Comparison clearly shows that for practical ranges of governing parameters the approximate solution matches the full solution with graphical accuracy, and, hence, can be eectively used in applications, e.g., in fracture tip element modules for the Planar3D approach to numerically modeling the propagation of a finite-size fracture. 4. Conclusions In this work, we developed a model for the eective fracture toughness coecient related to the initial nonuniform stress distribution. This study investigates the eect of the initial stress state on the apparent fracture toughness. We consider a vertical plane strain fracture in a strongly anisotropic stress field, such that the vertical stress is much larger than the horizontal stress. First, an elliptical opening profile is computed on the basis of linear elasticity and the assumption of uniform pressure inside the fracture. Then, this width profile is used to compute stress concentration around the elliptical hole caused by the large vertical stress. Correction to elastic energy is computed based on the stresses. This, in turn, is then translated into the value of the apparent toughness that becomes a function of the initial stresses. It is shown that the initial confining stress may significantly decrease fracture toughness of the reservoir. It is most

3.2. Nonuniform fracture growth in vertical and horizontal directions resulting in an elliptic shape It should be explicitly emphasized that even under uniform initial conditions in terms of the initial stress the solutions for stress intensity factors along the leading edge of an elliptical fracture (in profile) are dierent. Frequently when observing natural fractures and hydraulically-driven fractures it is noted that the shape is not radial and aspect ratio departs from

1 This remarkable approximation was proposed by Prof. Gennady Mishuris at the revision stage of the preparation of this paper for publication.

6

pronounced for soft media with the largest contrast between the smallest and the largest principal stresses. Since the stress difference (p1 p2 ) varies along the fracture boundary, the fracture propagation speed is dierent along the fracture boundary and, hence, the fracture shape departs from radial symmetry. A simplified asymptotic approximation to the solution for effective apparent fracture toughness is given, which is applicable for practical ranges of engineering applications (e.g. hydraulic fracturing technology). The expression for the eective fracture toughness obtained in this work may be then incorporated, e.g., into Planar3D models of commercial hydraulic fracturing simulators used for the design of fracturing jobs in the field.

[14] E. Dontsov, A. Peirce, A non-singular integral equation formulation to analyse multiscale behaviour in semi-infinite hydraulic fractures, Journal of Fluid Mechanics 781. [15] H. Zia, B. Lecampion, W. Zhang, Impact of the anisotropy of fracture toughness on the propagation of planar 3d hydraulic fracture, International Journal of Fracture 211 (1-2) (2018) 103–123. [16] N. Bahrami, D. Pena, I. Lusted, Well test, rate transient analysis and reservoir simulation for characterizing multi-fractured unconventional oil and gas reservoirs, Journal of Petroleum Exploration and Production Technology 6 (4) (2016) 675–689. [17] G. Sih, H. Liebowitz, Mathematical theories of brittle fracture, Fracture 2 (1968) 67–190. [18] M. J. Economides, K. G. Nolte, et al., Reservoir stimulation, Vol. 2, Prentice Hall Englewood Clis, NJ, 1989. [19] D. O. Potyondy, P. Cundall, A bonded-particle model for rock, International journal of rock mechanics and mining sciences 41 (8) (2004) 1329– 1364. [20] L. Huang, J. Liu, F. Zhang, E. Dontsov, B. Damjanac, Exploring the influence of rock inherent heterogeneity and grain size on hydraulic fracturing using discrete element modeling, International Journal of Solids and Structures. [21] M. Perkowska, A. Piccolroaz, M. Wrobel, G. Mishuris, Redirection of a crack driven by viscous fluid, International Journal of Engineering Science 121 (2017) 182–193. [22] M. Wrobel, A. Piccolroaz, P. Papanastasiou, G. Mishuris, Redirection of a crack driven by viscous fluid taking into account plastic eects in the process zone, Geomechanics for Energy and the Environment (2019) 100147. [23] M. Wrobel, G. Mishuris, A. Piccolroaz, On the impact of tangential traction on the crack surfaces induced by fluid in hydraulic fracture: Response to the letter of am linkov. int. j. eng. sci.(2018) 127, 217–219, International Journal of Engineering Science 127 (2018) 220–224. [24] I. N. Sneddon, D. S. Berry, The classical theory of elasticity, in: Elasticity and Plasticity/Elastizit¨at und Plastizit¨at, Springer, 1958, pp. 1–126. [25] M. A. Biot, Mechanics of incremental deformations, 1965. [26] I. A. Garagash, On the development of cracks in anisotropic bodies with initial stresses [in russian], Proc. Acad. Sci. Kazakh SSR 3 (1980) 16–21. [27] I. N. Sneddon, The distribution of stress in the neighbourhood of a crack in an elastic solid, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 187 (1009) (1946) 229–260. [28] S. Timoshenko, J. Goodier, Theory of Elasticity, McGraw-Hill, New York, 1970. [29] C. E. Inglis, Stresses in a plate due to the presence of cracks and sharp corners, Trans Inst Naval Archit 55 (1913) 219–241. [30] T. Hashida, K. Sato, H. Takahashi, Significance of crack opening for determining the growth behavior, 18th Workshop on Geothermal Reservoir Engineering. [31] E. Detournay, Propagation regimes of fluid-driven fractures in impermeable rocks, International Journal of Geomechanics 4 (1) (2004) 35–45. [32] S. Stanchits, J. Desroches, J. Burghardt, A. Surdi, N. Whitney, Rock fabric influence on hydraulic fracture propagation, in: 77th EAGE Conference and Exhibition 2015, 2015. [33] A. Bunger, B. Lecampion, Four critical issues for successful hydraulic fracturing applications, Tech. rep., CRC Press (2017). [34] H. Zia, B. Lecampion, Explicit versus implicit front advancing schemes for the simulation of hydraulic fracture growth, International Journal for Numerical and Analytical Methods in Geomechanics 43 (6) (2019) 1300– 1315. [35] E. Dontsov, A. Peirce, A multiscale implicit level set algorithm (ilsa) to model hydraulic fracture propagation incorporating combined viscous, toughness, and leak-o asymptotics, Computer Methods in Applied Mechanics and Engineering 313 (2017) 53–84. [36] H. Zia, B. Lecampion, Pyfrac: A planar 3d hydraulic fracture simulator, arXiv preprint arXiv:1908.10788. [37] H. Gu, E. Siebrits, et al., Eect of formation modulus contrast on hydraulic fracture height containment, in: International Oil & Gas Conference and Exhibition in China, Society of Petroleum Engineers, 2006. [38] J. Adachi, E. Siebrits, A. Peirce, J. Desroches, Computer simulation of hydraulic fractures, International Journal of Rock Mechanics and Mining Sciences 44 (5) (2007) 739–757. [39] N. Muskhelishvili, Some basic problems of the mathematical theory of

Acknowledgements The authors are grateful to Prof. Gennady Mishuris for fruitful discussions on the results of this work and for proposing a simplified approximate expression for the eective apparent toughness based on the full solution. Startup funds of Skolkovo Institute of Science and Technology are gratefully acknowledged by A. Osiptsov. This work received financial support from the Ministry of Education and Science of the Russian Federation (project No. 14.581.21.0027, unique identifier RFMEFI58117X0027; project No. AAAAA17-117051110248-3). References [1] J. Adachi, E. Siebrits, A. Peirce, J. Desroches, Computer simulation of hydraulic fractures, International Journal of Rock Mechanics and Mining Sciences 44 (5) (2007) 739–757. [2] E. Detournay, Mechanics of hydraulic fractures, Annual Review of Fluid Mechanics 48 (2016) 311–339. [3] B. Lecampion, A. Bunger, X. Zhang, Numerical methods for hydraulic fracture propagation: a review of recent trends, Journal of natural gas science and engineering 49 (2018) 66–83. [4] A. A. Osiptsov, Fluid mechanics of hydraulic fracturing: a review, Journal of petroleum science and engineering 156 (2017) 513–535. [5] D. I. Garagash, E. Detournay, J. I. Adachi, Multiscale tip asymptotics in hydraulic fracture with leak-o, Journal of Fluid Mechanics 669 (2011) 260–297. [6] J. R. Rice, Mathematical analysis in the mechanics of fracture, Fracture: an advanced treatise 2 (1968) 191–311. [7] G. R. Irwin, Analysis of stresses and strains near the end of a crack transversing a plate, Trans. ASME, Ser. E, J. Appl. Mech. 24 (1957) 361– 364. [8] A. A. Grith, Vi. the phenomena of rupture and flow in solids, Philosophical transactions of the royal society of london. Series A, containing papers of a mathematical or physical character 221 (582-593) (1921) 163–198. [9] B. Lawn, D. Marshall, Hardness, toughness, and brittleness: an indentation analysis, Journal of the American ceramic society 62 (7-8) (1979) 347–350. [10] D. Liu, B. Lecampion, D. I. Garagash, Propagation of a fluid-driven fracture with fracture length dependent apparent toughness, Engineering Fracture Mechanics (2019) 106616. [11] D. Tromans, Crack propagation in brittle materials: relevance to minerals comminution, Int. J. Res. Rev. Appl. Sci 13 (2012) 406–427. [12] L. B. Freund, Dynamic fracture mechanics, Cambridge university press, 1998. [13] M. Wrobel, G. Mishuris, A. Piccolroaz, Energy release rate in hydraulic fracture: Can we neglect an impact of the hydraulically induced shear stress?, International Journal of Engineering Science 111 (2017) 28–51.

7

elasticity, Noordho, Groningen 17404. [40] E. Detournay, D. Garagash, The near-tip region of a fluid-driven fracture propagating in a permeable elastic solid, Journal of Fluid Mechanics 494 (2003) 1–32. [41] D. I. Garagash, Propagation of a plane-strain hydraulic fracture with a fluid lag: Early-time solution, International journal of solids and structures 43 (18-19) (2006) 5811–5835. [42] J. Desroches, E. Detournay, B. Lenoach, P. Papanastasiou, J. R. A. Pearson, M. Thiercelin, A. Cheng, The crack tip region in hydraulic fracturing, Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 447 (1929) (1994) 39–48. [43] H. Van Eekelen, et al., Hydraulic fracture geometry: fracture containment in layered formations, Society of Petroleum Engineers Journal 22 (03) (1982) 341–349. [44] H. Nisitani, Y. Murakami, Stress intensity factors of an elliptical crack or a semi-elliptical crack subject to tension, International Journal of Fracture 10 (3) (1974) 353–368. [45] D. A. Chuprakov, R. Prioul, et al., Hydraulic fracture height containment by weak horizontal interfaces, in: SPE hydraulic fracturing technology conference, Society of Petroleum Engineers, 2015.

8

Figure 1: Orientation of hydraulic fractures in the initial stress field. Longitudinal and transverse fractures (a) (sketch from [16]); an individual fracture as an ellipse with half-axes (b).

Figure 2: Distribution of loads exerted on the fracture.

Figure 3: Fracture shape for  = 0:2 and (q

9

p2 )=G

= 3
10

3.

Figure 4: Elliptic hole formed by the loads q

p2 .

Figure 5: Distribution of stress  at the fracture tip for x  a (a) and for x  a (b).

10

Figure 6: Hydraulic fracture in elastic medium.

Figure 8: variation of the fracture coecient increment
K1c for a (curve 1), G = 1:2
1010 Pa (2), G = 1:8
1010 Pa (3).

ef f

= 25m,  = 0:2 and K1c = 5
105 Pa
m1=2 (a) and K1c = 1
106 Pa
m1=2 (b). G = 6
109 Pa

Figure 9: Variation of the ratio K1c =K1c with increasing fracture length for  (curve 1), G = 1:2
1010 Pa (2), G = 1:8
1010 Pa (3).

= 0:2 and K1c = 5
105 Pa
m1=2 (a) and K1c = 1
106 Pa
m1=2 (b). G = 6
109 Pa

11

Figure 10: The decrease in the stress intensity coecient K1c =K1c for K1c = 5
105 Pa
m0:5 (a) and K1c correspond to fluid-saturated half-space with hydrostatic pressure distribution, red curves – for dry rock. ef f

= 5
16 Pa
m0:5 (b) for G = 1:2
1010 Pa. Blue curves

Figure 12: Comparison of formula (25) shown by red line with the simplified asymptotic formula (28) shown by blue dots suitable for engineering applications, obtained for: K1c = 5
105 Pa
m1=2 , G = 1:2
1010 Pa,  = 0:5 (a,b) and (p1 p2 ) = 5:625
107 Pa (b).

12

   

Analytical model is developed for fracture toughness dependent on non-uniform initial stress state; A simplified approximate formula is obtained for apparent fracture toughness as a function of the difference between vertical and horizontal initial loads; Parametric study is conducted, showing also the dependence of fracture toughness on the length of a growing fracture; Formula may be further used in near-tip asymptotics within the tip element in the numerical implementation of a Planar3D approach to fracture simulation;

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: